// Copyright (C) 2002, International Business Machines
// Corporation and others.  All Rights Reserved.

#include "CoinPragma.hpp"
#include <stdio.h>
#include <stdarg.h>
#include <stdlib.h>
#include <math.h>
#include "ClpPresolve.hpp"
#include "Idiot.hpp"
#include "CoinTime.hpp"
#include "CoinSort.hpp"
#include "CoinMessageHandler.hpp"
#include "CoinHelperFunctions.hpp"
// Redefine stuff for Clp
#ifndef OSI_IDIOT
#include "ClpMessage.hpp"
#define OsiObjOffset ClpObjOffset
#endif
/**** strategy 4 - drop, exitDrop and djTolerance all relative:
For first two major iterations these are small.  Then:

drop - exit a major iteration if drop over 5*checkFrequency < this is
used as info->drop*(10.0+fabs(last weighted objective))

exitDrop - exit idiot if feasible and drop < this is
used as info->exitDrop*(10.0+fabs(last objective))

djExit - exit a major iteration if largest dj (averaged over 5 checks)
drops below this - used as info->djTolerance*(10.0+fabs(last weighted objective)

djFlag - mostly skip variables with bad dj worse than this => 2*djExit

djTol - only look at variables with dj better than this => 0.01*djExit
****************/

#define IDIOT_FIX_TOLERANCE 1e-6
#define SMALL_IDIOT_FIX_TOLERANCE 1e-10
int 
Idiot::dropping(IdiotResult result,
		    double tolerance,
		    double small,
                    int *nbad)
{
  if (result.infeas<=small) {
    double value=CoinMax(fabs(result.objval),fabs(result.dropThis))+1.0;
    if (result.dropThis>tolerance*value) {
      *nbad=0;
      return 1;
    } else {
      (*nbad)++;
      if (*nbad>4) {
	return 0;
      } else {
	return 1;
      }
    }
  } else {
    *nbad=0;
    return 1;
  }
}
// Deals with whenUsed and slacks
int 
Idiot::cleanIteration(int iteration, int ordinaryStart, int ordinaryEnd,
		      double * colsol, const double * lower, const double * upper,
		      const double * rowLower, const double * rowUpper,
		      const double * cost, const double * element, double fixTolerance,
		      double & objValue, double & infValue)
{
  int n=0;
  if ((strategy_&16384)==0) {
    for (int i=ordinaryStart;i<ordinaryEnd;i++) {
      if (colsol[i]>lower[i]+fixTolerance) {
	if (colsol[i]<upper[i]-fixTolerance) {
	  n++;
	} else {
	  colsol[i]=upper[i];
	}
	whenUsed_[i]=iteration;
      } else {
	colsol[i]=lower[i];
      }
    }
    return n;
  } else {
    printf("entering inf %g, obj %g\n",infValue,objValue);
    int nrows=model_->getNumRows();
    int ncols=model_->getNumCols();
    int * posSlack = whenUsed_+ncols;
    int * negSlack = posSlack+nrows;
    int * nextSlack = negSlack + nrows;
    double * rowsol = (double *) (nextSlack+ncols);
    memset(rowsol,0,nrows*sizeof(double));
#ifdef OSI_IDIOT
    const CoinPackedMatrix * matrix = model_->getMatrixByCol();
#else
    ClpMatrixBase * matrix = model_->clpMatrix();
#endif
    const int * row = matrix->getIndices();
    const CoinBigIndex * columnStart = matrix->getVectorStarts();
    const int * columnLength = matrix->getVectorLengths(); 
    //const double * element = matrix->getElements();
    int i;
    objValue=0.0;
    infValue=0.0;
    for ( i=0;i<ncols;i++) {
      if (nextSlack[i]==-1) {
	// not a slack
	if (colsol[i]>lower[i]+fixTolerance) {
	  if (colsol[i]<upper[i]-fixTolerance) {
	    n++;
	    whenUsed_[i]=iteration;
	  } else {
	    colsol[i]=upper[i];
	  }
	  whenUsed_[i]=iteration;
	} else {
	  colsol[i]=lower[i];
	}
	double value = colsol[i];
	if (value) {
	  objValue += cost[i]*value;
	  CoinBigIndex j;
	  for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
	    int iRow = row[j];
	    rowsol[iRow] += value*element[j];
	  }
	}
      }
    }
    // temp fix for infinite lbs - just limit to -1000
    for (i=0;i<nrows;i++) {
      double rowSave=rowsol[i];
      int iCol;
      iCol =posSlack[i];
      if (iCol>=0) {
	// slide all slack down
	double rowValue=rowsol[i];
	CoinBigIndex j=columnStart[iCol];
	double lowerValue = CoinMax(CoinMin(colsol[iCol],0.0)-1000.0,lower[iCol]);
	rowSave += (colsol[iCol]-lowerValue)*element[j];
	colsol[iCol]=lowerValue;
	while (nextSlack[iCol]>=0) {
	  iCol = nextSlack[iCol];
	  double lowerValue = CoinMax(CoinMin(colsol[iCol],0.0)-1000.0,lower[iCol]);
	  j=columnStart[iCol];
	  rowSave += (colsol[iCol]-lowerValue)*element[j];
	  colsol[iCol]=lowerValue;
	}
	iCol =posSlack[i];
	while (rowValue<rowLower[i]&&iCol>=0) {
	  // want to increase
	  double distance = rowLower[i]-rowValue;
	  double value = element[columnStart[iCol]];
	  double thisCost = cost[iCol];
	  if (distance<=value*(upper[iCol]-colsol[iCol])) {
	    // can get there
	    double movement = distance/value;
	    objValue += movement*thisCost;
	    rowValue = rowLower[i];
	    colsol[iCol] += movement;
	  } else {
	    // can't get there
	    double movement = upper[iCol]-colsol[iCol];
	    objValue += movement*thisCost;
	    rowValue += movement*value;
	    colsol[iCol] = upper[iCol];
	    iCol = nextSlack[iCol];
	  }
	}
	if (iCol>=0) {
	  // may want to carry on - because of cost?
	  while (cost[iCol]<0&&rowValue<rowUpper[i]) {
	    // want to increase
	    double distance = rowUpper[i]-rowValue;
	    double value = element[columnStart[iCol]];
	    double thisCost = cost[iCol];
	    if (distance<=value*(upper[iCol]-colsol[iCol])) {
	      // can get there
	      double movement = distance/value;
	      objValue += movement*thisCost;
	      rowValue = rowUpper[i];
	      colsol[iCol] += movement;
	      iCol=-1;
	    } else {
	      // can't get there
	      double movement = upper[iCol]-colsol[iCol];
	      objValue += movement*thisCost;
	      rowValue += movement*value;
	      colsol[iCol] = upper[iCol];
	      iCol = nextSlack[iCol];
	    }
	  }
	  if (iCol>=0&&colsol[iCol]>lower[iCol]+fixTolerance&&
	      colsol[iCol]<upper[iCol]-fixTolerance) {
	    whenUsed_[i]=iteration;
	    n++;
	  }
	}
	rowsol[i]=rowValue;
      }
      iCol =negSlack[i];
      if (iCol>=0) {
	// slide all slack down
	double rowValue=rowsol[i];
	CoinBigIndex j=columnStart[iCol];
	rowSave += (colsol[iCol]-lower[iCol])*element[j];
	colsol[iCol]=lower[iCol];
	assert (lower[iCol]>-1.0e20);
	while (nextSlack[iCol]>=0) {
	  iCol = nextSlack[iCol];
	  j=columnStart[iCol];
	  rowSave += (colsol[iCol]-lower[iCol])*element[j];
	  colsol[iCol]=lower[iCol];
	}
	iCol =negSlack[i];
	while (rowValue>rowUpper[i]&&iCol>=0) {
	  // want to increase
	  double distance = -(rowUpper[i]-rowValue);
	  double value = -element[columnStart[iCol]];
	  double thisCost = cost[iCol];
	  if (distance<=value*(upper[iCol]-lower[iCol])) {
	    // can get there
	    double movement = distance/value;
	    objValue += movement*thisCost;
	    rowValue = rowUpper[i];
	    colsol[iCol] += movement;
	  } else {
	    // can't get there
	    double movement = upper[iCol]-lower[iCol];
	    objValue += movement*thisCost;
	    rowValue -= movement*value;
	    colsol[iCol] = upper[iCol];
	    iCol = nextSlack[iCol];
	  }
	}
	if (iCol>=0) {
	  // may want to carry on - because of cost?
	  while (cost[iCol]<0&&rowValue>rowLower[i]) {
	    // want to increase
	    double distance = -(rowLower[i]-rowValue);
	    double value = -element[columnStart[iCol]];
	    double thisCost = cost[iCol];
	    if (distance<=value*(upper[iCol]-colsol[iCol])) {
	      // can get there
	      double movement = distance/value;
	      objValue += movement*thisCost;
	      rowValue = rowLower[i];
	      colsol[iCol] += movement;
	      iCol=-1;
	    } else {
	      // can't get there
	      double movement = upper[iCol]-colsol[iCol];
	      objValue += movement*thisCost;
	      rowValue -= movement*value;
	      colsol[iCol] = upper[iCol];
	      iCol = nextSlack[iCol];
	    }
	  }
	  if (iCol>=0&&colsol[iCol]>lower[iCol]+fixTolerance&&
	      colsol[iCol]<upper[iCol]-fixTolerance) {
	    whenUsed_[i]=iteration;
	    n++;
	  }
	}
	rowsol[i]=rowValue;
      }
      infValue += CoinMax(CoinMax(0.0,rowLower[i]-rowsol[i]),rowsol[i]-rowUpper[i]);
      // just change
      rowsol[i] -= rowSave;
    }
    return n;
  }
}

/* returns -1 if none or start of costed slacks or -2 if
   there are costed slacks but it is messy */
 static int countCostedSlacks(OsiSolverInterface * model)
{
#ifdef OSI_IDIOT
  const CoinPackedMatrix * matrix = model->getMatrixByCol();
#else
  ClpMatrixBase * matrix = model->clpMatrix();
#endif
  const int * row = matrix->getIndices();
  const CoinBigIndex * columnStart = matrix->getVectorStarts();
  const int * columnLength = matrix->getVectorLengths(); 
  const double * element = matrix->getElements();
  const double * rowupper = model->getRowUpper();
  int nrows=model->getNumRows();
  int ncols=model->getNumCols();
  int slackStart = ncols-nrows;
  int nSlacks=nrows;
  int i;
  
  if (ncols<=nrows) return -1;
  while (1) {
    for (i=0;i<nrows;i++) {
      int j=i+slackStart;
      CoinBigIndex k=columnStart[j];
      if (columnLength[j]==1) {
	if (row[k]!=i||element[k]!=1.0) {
	  nSlacks=0;
	  break;
	}
      } else {
	nSlacks=0;
	break;
      }
      if (rowupper[i]<=0.0) {
	nSlacks=0;
	break;
      }
    }
    if (nSlacks||!slackStart) break;
    slackStart=0;
  }
  if (!nSlacks) slackStart=-1;
  return slackStart;
}
void
Idiot::crash(int numberPass, CoinMessageHandler * handler,const CoinMessages *messages)
{
  // lightweight options
  int numberColumns = model_->getNumCols();
  const double * objective = model_->getObjCoefficients();
  int nnzero=0;
  double sum=0.0;
  int i;
  for (i=0;i<numberColumns;i++) {
    if (objective[i]) {
      sum += fabs(objective[i]);
      nnzero++;
    }
  }
  sum /= (double) (nnzero+1);
  if (maxIts_==5)
    maxIts_=2;
  if (numberPass<=0)
    // Cast to double to avoid VACPP complaining
    majorIterations_=(int)(2+log10((double)(numberColumns+1)));
  else
    majorIterations_=numberPass;
  // If mu not changed then compute
  if (mu_==1e-4)
    mu_= CoinMax(1.0e-3,sum*1.0e-5);
  if (maxIts2_==100) {
    if (!lightWeight_) {
      maxIts2_=105;
    } else if (lightWeight_==1) {
      mu_ *= 1000.0;
      maxIts2_=23;
    } else if (lightWeight_==2) {
      maxIts2_=11;
    } else {
      maxIts2_=23;
    }
  }
  //printf("setting mu to %g and doing %d passes\n",mu_,majorIterations_);
  solve2(handler,messages);
#ifndef OSI_IDIOT
  double averageInfeas = model_->sumPrimalInfeasibilities()/((double) model_->numberRows());
  if ((averageInfeas<0.01&&(strategy_&512)!=0)||(strategy_&8192)!=0) 
    crossOver(16+1); 
  else
    crossOver(3);
#endif
}
void
Idiot::solve()
{
  CoinMessages dummy;
  solve2(NULL,&dummy);
}
void
Idiot::solve2(CoinMessageHandler * handler,const CoinMessages * messages)
{
  int strategy=0;
  double d2;
  int i,n;
  int allOnes=1;
  int iteration=0;
  int iterationTotal=0;
  int nTry=0; /* number of tries at same weight */
  double fixTolerance=IDIOT_FIX_TOLERANCE;
  int maxBigIts=maxBigIts_;
  int maxIts=maxIts_;
  int logLevel=logLevel_;
  int saveMajorIterations = majorIterations_;
  if (handler) {
    if (handler->logLevel()>0&&handler->logLevel()<3)
      logLevel=1;
    else if (!handler->logLevel())
      logLevel=0;
    else
      logLevel=7;
  }
  double djExit=djTolerance_;
  double djFlag=1.0+100.0*djExit;
  double djTol=0.00001;
  double mu =mu_;
  double drop=drop_;
  int maxIts2=maxIts2_;
  double factor=muFactor_;
  double smallInfeas=smallInfeas_;
  double reasonableInfeas=reasonableInfeas_;
  double stopMu=stopMu_;
  double maxmin,offset;
  double lastWeighted=1.0e50;
  double exitDrop=exitDrop_;
  double fakeSmall=smallInfeas;
  double firstInfeas;
  int badIts=0;
  int slackStart,slackEnd,ordStart,ordEnd;
  int checkIteration=0;
  int lambdaIteration=0;
  int belowReasonable=0; /* set if ever gone below reasonable infeas */
  double bestWeighted=1.0e60;
  double bestFeasible=1.0e60; /* best solution while feasible */
  IdiotResult result,lastResult;
  int saveStrategy=strategy_;
  const int strategies[]={0,2,128};
  double saveLambdaScale=0.0;
  if ((saveStrategy&128)!=0) {
    fixTolerance=SMALL_IDIOT_FIX_TOLERANCE;
  }
#ifdef OSI_IDIOT
  const CoinPackedMatrix * matrix = model_->getMatrixByCol();
#else
  ClpMatrixBase * matrix = model_->clpMatrix();
#endif
  const int * row = matrix->getIndices();
  const CoinBigIndex * columnStart = matrix->getVectorStarts();
  const int * columnLength = matrix->getVectorLengths(); 
  const double * element = matrix->getElements();
  int nrows=model_->getNumRows();
  int ncols=model_->getNumCols();
  double * rowsol, * colsol;
  double * pi, * dj;
#ifndef OSI_IDIOT
  double * cost = model_->objective();
  double * lower = model_->columnLower();
  double * upper = model_->columnUpper();
#else
  double * cost = new double [ncols];
  memcpy( cost, model_->getObjCoefficients(), ncols*sizeof(double));
  const double * lower = model_->getColLower();
  const double * upper = model_->getColUpper();
#endif
  const double *elemXX;
  double * saveSol;
  double * rowupper= new double[nrows]; // not const as modified
  memcpy(rowupper,model_->getRowUpper(),nrows*sizeof(double));
  double * rowlower= new double[nrows]; // not const as modified
  memcpy(rowlower,model_->getRowLower(),nrows*sizeof(double));
  int * whenUsed;
  double * lambda;
  saveSol=new double[ncols];
  lambda=new double [nrows];
  rowsol= new double[nrows];
  colsol= new double [ncols];
  memcpy(colsol,model_->getColSolution(),ncols*sizeof(double));
  pi= new double[nrows];
  dj=new double[ncols];
  delete [] whenUsed_;
  bool oddSlacks=false;
  // See if any costed slacks
  int numberSlacks=0;
  for (i=0;i<ncols;i++) {
    if (columnLength[i]==1)
      numberSlacks++;
  }
  if (!numberSlacks) {
    whenUsed_=new int[ncols];
  } else {
#ifdef COIN_DEVELOP
    printf("%d slacks\n",numberSlacks);
#endif
    oddSlacks=true;
    int extra = (int) (nrows*sizeof(double)/sizeof(int));
    whenUsed_=new int[2*ncols+2*nrows+extra];
    int * posSlack = whenUsed_+ncols;
    int * negSlack = posSlack+nrows;
    int * nextSlack = negSlack + nrows;
    for (i=0;i<nrows;i++) {
      posSlack[i]=-1;
      negSlack[i]=-1;
    }
    for (i=0;i<ncols;i++) 
      nextSlack[i]=-1;
    for (i=0;i<ncols;i++) {
      if (columnLength[i]==1) {
	CoinBigIndex j=columnStart[i];
	int iRow = row[j];
	if (element[j]>0.0) {
	  if (posSlack[iRow]==-1) {
	    posSlack[iRow]=i;
	  } else {
	    int iCol = posSlack[iRow];
	    while (nextSlack[iCol]>=0)
	      iCol = nextSlack[iCol];
	    nextSlack[iCol]=i;
	  }
	} else {
	  if (negSlack[iRow]==-1) {
	    negSlack[iRow]=i;
	  } else {
	    int iCol = negSlack[iRow];
	    while (nextSlack[iCol]>=0)
	      iCol = nextSlack[iCol];
	    nextSlack[iCol]=i;
	  }
	}
      }
    }
    // now sort
    for (i=0;i<nrows;i++) {
      int iCol;
      iCol = posSlack[i];
      if (iCol>=0) {
	CoinBigIndex j=columnStart[iCol];
#ifndef NDEBUG
	int iRow = row[j];
#endif
	assert (element[j]>0.0);
	assert (iRow==i);
	dj[0]= cost[iCol]/element[j];
	whenUsed_[0]=iCol;
	int n=1;
	while (nextSlack[iCol]>=0) {
	  iCol = nextSlack[iCol];
	  CoinBigIndex j=columnStart[iCol];
#ifndef NDEBUG
	  int iRow = row[j];
#endif
	  assert (element[j]>0.0);
	  assert (iRow==i);
	  dj[n]= cost[iCol]/element[j];
	  whenUsed_[n++]=iCol;
	}
	for (j=0;j<n;j++) {
	  int jCol = whenUsed_[j];
	  nextSlack[jCol]=-2;
	}
	CoinSort_2(dj,dj+n,whenUsed_);
	// put back
	iCol = whenUsed_[0];
	posSlack[i]=iCol;
	for (j=1;j<n;j++) {
	  int jCol = whenUsed_[j];
	  nextSlack[iCol]=jCol;
	  iCol=jCol;
	}
      }
      iCol = negSlack[i];
      if (iCol>=0) {
	CoinBigIndex j=columnStart[iCol];
#ifndef NDEBUG
	int iRow = row[j];
#endif
	assert (element[j]<0.0);
	assert (iRow==i);
	dj[0]= -cost[iCol]/element[j];
	whenUsed_[0]=iCol;
	int n=1;
	while (nextSlack[iCol]>=0) {
	  iCol = nextSlack[iCol];
	  CoinBigIndex j=columnStart[iCol];
#ifndef NDEBUG
	  int iRow = row[j];
#endif
	  assert (element[j]<0.0);
	  assert (iRow==i);
	  dj[n]= -cost[iCol]/element[j];
	  whenUsed_[n++]=iCol;
	}
	for (j=0;j<n;j++) {
	  int jCol = whenUsed_[j];
	  nextSlack[jCol]=-2;
	}
	CoinSort_2(dj,dj+n,whenUsed_);
	// put back
	iCol = whenUsed_[0];
	negSlack[i]=iCol;
	for (j=1;j<n;j++) {
	  int jCol = whenUsed_[j];
	  nextSlack[iCol]=jCol;
	  iCol=jCol;
	}
      }
    }
  }
  whenUsed=whenUsed_;
  if (model_->getObjSense()==-1.0) {
    maxmin=-1.0;
  } else {
    maxmin=1.0;
  }
  model_->getDblParam(OsiObjOffset,offset);
  if (!maxIts2) maxIts2=maxIts;
  strategy=strategy_;
  strategy &= 3;
  memset(lambda,0,nrows*sizeof(double));
  slackStart=countCostedSlacks(model_);
  if (slackStart>=0) {
    printf("This model has costed slacks\n");
    slackEnd=slackStart+nrows;
    if (slackStart) {
      ordStart=0;
      ordEnd=slackStart;
    } else {
      ordStart=nrows;
      ordEnd=ncols;
    }
  } else {
    slackEnd=slackStart;
    ordStart=0;
    ordEnd=ncols;
  }
  if (offset&&logLevel>2) {
    printf("** Objective offset is %g\n",offset);
  }
  /* compute reasonable solution cost */
  for (i=0;i<nrows;i++) {
    rowsol[i]=1.0e31;
  }
  for (i=0;i<ncols;i++) {
    CoinBigIndex j;
    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
      if (element[j]!=1.0) {
	allOnes=0;
	break;
      }
    }
  }
  if (allOnes) {
    elemXX=NULL;
  } else {
    elemXX=element;
  }
  // Do scaling if wanted
  bool scaled=false;
#ifndef OSI_IDIOT
  if ((strategy_&32)!=0&&!allOnes) {
    if (model_->scalingFlag()>0)
      scaled = model_->clpMatrix()->scale(model_)==0;
    if (scaled) {
      const double * rowScale = model_->rowScale();
      const double * columnScale = model_->columnScale();
      double * oldLower = lower;
      double * oldUpper = upper;
      double * oldCost = cost;
      lower = new double[ncols];
      upper = new double[ncols];
      cost = new double[ncols];
      memcpy(lower,oldLower,ncols*sizeof(double));
      memcpy(upper,oldUpper,ncols*sizeof(double));
      memcpy(cost,oldCost,ncols*sizeof(double));
      int icol,irow;
      for (icol=0;icol<ncols;icol++) {
	double multiplier = 1.0/columnScale[icol];
	if (lower[icol]>-1.0e50)
	  lower[icol] *= multiplier;
	if (upper[icol]<1.0e50)
	  upper[icol] *= multiplier;
	colsol[icol] *= multiplier;
	cost[icol] *= columnScale[icol];
      }
      memcpy(rowlower,model_->rowLower(),nrows*sizeof(double));
      for (irow=0;irow<nrows;irow++) {
	double multiplier = rowScale[irow];
	if (rowlower[irow]>-1.0e50)
	  rowlower[irow] *= multiplier;
	if (rowupper[irow]<1.0e50)
	  rowupper[irow] *= multiplier;
	rowsol[irow] *= multiplier;
      }
      int length = columnStart[ncols-1]+columnLength[ncols-1];
      double * elemYY = new double[length];
      for (i=0;i<ncols;i++) {
	CoinBigIndex j;
	double scale = columnScale[i];
	for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
	  int irow=row[j];
	  elemYY[j] = element[j]*scale*rowScale[irow];
	}
      }
      elemXX=elemYY;
    }
  }
#endif
  for (i=0;i<ncols;i++) {
    CoinBigIndex j;
    double dd=columnLength[i];
    dd=cost[i]/dd;
    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
      int irow=row[j];
      if (dd<rowsol[irow]) {
	rowsol[irow]=dd;
      }
    }
  }
  d2=0.0;
  for (i=0;i<nrows;i++) {
    d2+=rowsol[i];
  }
  d2*=2.0; /* for luck */
  
  d2=d2/((double) (4*nrows+8000));
  d2*=0.5; /* halve with more flexible method */
  if (d2<5.0) d2=5.0;
  if (djExit==0.0) {
    djExit=d2;
  }
  if ((saveStrategy&4)!=0) {
    /* go to relative tolerances - first small */
    djExit=1.0e-10;
    djFlag=1.0e-5;
    drop=1.0e-10;
  }
  memset(whenUsed,0,ncols*sizeof(int));
  strategy=strategies[strategy];
  if ((saveStrategy&8)!=0) strategy |= 64; /* don't allow large theta */
  memcpy(saveSol,colsol,ncols*sizeof(double));
  
  lastResult=IdiSolve(nrows,ncols,rowsol ,colsol,pi,
		       dj,cost,rowlower,rowupper,
		       lower,upper,elemXX,row,columnStart,columnLength,lambda,
		       0,mu,drop,
		       maxmin,offset,strategy,djTol,djExit,djFlag);
  // update whenUsed_
  n = cleanIteration(iteration, ordStart,ordEnd,
		     colsol,  lower,  upper,
		     rowlower, rowupper,
		     cost, elemXX, fixTolerance,lastResult.objval,lastResult.infeas);
  if ((strategy_&16384)!=0) {
    int * posSlack = whenUsed_+ncols;
    int * negSlack = posSlack+nrows;
    int * nextSlack = negSlack + nrows;
    double * rowsol2 = (double *) (nextSlack+ncols);
    for (i=0;i<nrows;i++)
      rowsol[i] += rowsol2[i];
  }
  if ((logLevel_&1)!=0) {
#ifndef OSI_IDIOT
    if (!handler) {
#endif
      printf("Iteration %d infeasibility %g, objective %g - mu %g, its %d, %d interior\n", 
	     iteration,lastResult.infeas,lastResult.objval,mu,lastResult.iteration,n);
#ifndef OSI_IDIOT
    } else {
      handler->message(CLP_IDIOT_ITERATION,*messages)
	<<iteration<<lastResult.infeas<<lastResult.objval<<mu<<lastResult.iteration<<n
	<<CoinMessageEol;
    }
#endif
  }
  int numberBaseTrys=0; // for first time
  int numberAway=-1;
  iterationTotal = lastResult.iteration;
  firstInfeas=lastResult.infeas;
  if ((strategy_&1024)!=0) reasonableInfeas=0.5*firstInfeas;
  if (lastResult.infeas<reasonableInfeas) lastResult.infeas=reasonableInfeas;
  double keepinfeas=1.0e31;
  double lastInfeas=1.0e31;
  double bestInfeas=1.0e31;
  while ((mu>stopMu&&lastResult.infeas>smallInfeas)||
         (lastResult.infeas<=smallInfeas&&
	 dropping(lastResult,exitDrop,smallInfeas,&badIts))||
	 checkIteration<2||lambdaIteration<lambdaIterations_) {
    if (lastResult.infeas<=exitFeasibility_)
      break; 
    iteration++;
    checkIteration++;
    if (lastResult.infeas<=smallInfeas&&lastResult.objval<bestFeasible) {
      bestFeasible=lastResult.objval;
    }
    if (lastResult.infeas+mu*lastResult.objval<bestWeighted) {
      bestWeighted=lastResult.objval+mu*lastResult.objval;
    }
    if ((saveStrategy&4096)) strategy |=256;
    if ((saveStrategy&4)!=0&&iteration>2) {
      /* go to relative tolerances */
      double weighted=10.0+fabs(lastWeighted);
      djExit=djTolerance_*weighted;
      djFlag=2.0*djExit;
      drop=drop_*weighted;
      djTol=0.01*djExit;
    }
    memcpy(saveSol,colsol,ncols*sizeof(double));
    result=IdiSolve(nrows,ncols,rowsol ,colsol,pi,dj,
		     cost,rowlower,rowupper,
		     lower,upper,elemXX,row,columnStart,columnLength,lambda,
		     maxIts,mu,drop,
		     maxmin,offset,strategy,djTol,djExit,djFlag);
    n = cleanIteration(iteration, ordStart,ordEnd,
		       colsol,  lower,  upper,
		       rowlower, rowupper,
		       cost, elemXX, fixTolerance,result.objval,result.infeas);
    if ((strategy_&16384)!=0) {
      int * posSlack = whenUsed_+ncols;
      int * negSlack = posSlack+nrows;
      int * nextSlack = negSlack + nrows;
      double * rowsol2 = (double *) (nextSlack+ncols);
      for (i=0;i<nrows;i++)
	rowsol[i] += rowsol2[i];
    }
    if ((logLevel_&1)!=0) {
#ifndef OSI_IDIOT
      if (!handler) {
#endif
	printf("Iteration %d infeasibility %g, objective %g - mu %g, its %d, %d interior\n", 
	       iteration,result.infeas,result.objval,mu,result.iteration,n);
#ifndef OSI_IDIOT
      } else {
	handler->message(CLP_IDIOT_ITERATION,*messages)
	  <<iteration<<result.infeas<<result.objval<<mu<<result.iteration<<n
	  <<CoinMessageEol;
      }
#endif
    }
    if (iteration>50&&n==numberAway&&result.infeas<1.0e-4) {
      printf("infeas small %g\n",result.infeas);
      break; // not much happening
    }
    if (lightWeight_==1&&iteration>10&&result.infeas>1.0&&maxIts!=7) {
      if (lastInfeas!=bestInfeas&&CoinMin(result.infeas,lastInfeas)>0.95*bestInfeas)
	majorIterations_ = CoinMin(majorIterations_,iteration); // not getting feasible
    }
    lastInfeas = result.infeas;
    numberAway=n;
    keepinfeas = result.infeas;
    lastWeighted=result.weighted;
    iterationTotal += result.iteration;
    if (iteration==1) {
      if ((strategy_&1024)!=0&&mu<1.0e-10) 
	result.infeas=firstInfeas*0.8;
      if (majorIterations_>=50||dropEnoughFeasibility_<=0.0)
        result.infeas *= 0.8;
      if (result.infeas>firstInfeas*0.9
	  &&result.infeas>reasonableInfeas) {
	iteration--;
	if (majorIterations_<50)
	  mu*=1.0e-1;
	else
	  mu*=0.7;
	bestFeasible=1.0e31;
        bestWeighted=1.0e60;
        numberBaseTrys++;
        if (mu<1.0e-30||(numberBaseTrys>10&&lightWeight_)) {
          // back to all slack basis
          lightWeight_=2;
          break;
        }
	memcpy(colsol,saveSol,ncols*sizeof(double));
      } else {
	maxIts=maxIts2;
	checkIteration=0;
	if ((strategy_&1024)!=0) mu *= 1.0e-1;
      }
    } else {
    }
    bestInfeas=CoinMin(bestInfeas,result.infeas);
    if (iteration) {
      /* this code is in to force it to terminate sometime */
      double changeMu=factor;
      if ((saveStrategy&64)!=0) {
	keepinfeas=0.0; /* switch off ranga's increase */
	fakeSmall=smallInfeas;
      } else {
	fakeSmall=-1.0;
      }
      saveLambdaScale=0.0;
      if (result.infeas>reasonableInfeas||
	  (nTry+1==maxBigIts&&result.infeas>fakeSmall)) {
	if (result.infeas>lastResult.infeas*(1.0-dropEnoughFeasibility_)||
	    nTry+1==maxBigIts||
	    (result.infeas>lastResult.infeas*0.9
	     &&result.weighted>lastResult.weighted
	     -dropEnoughWeighted_*CoinMax(fabs(lastResult.weighted),fabs(result.weighted)))) {
	  mu*=changeMu;
          if ((saveStrategy&32)!=0&&result.infeas<reasonableInfeas&&0) {
	    reasonableInfeas=CoinMax(smallInfeas,reasonableInfeas*sqrt(changeMu));
	    printf("reasonable infeas now %g\n",reasonableInfeas);
	  }
	  result.weighted=1.0e60;
	  nTry=0;
	  bestFeasible=1.0e31;
	  bestWeighted=1.0e60;
	  checkIteration=0;
	  lambdaIteration=0;
#define LAMBDA
#ifdef LAMBDA
	  if ((saveStrategy&2048)==0) {
	    memset(lambda,0,nrows*sizeof(double));
	  }
#else
	  memset(lambda,0,nrows*sizeof(double));
#endif
	} else {
	  nTry++;
	}
      } else if (lambdaIterations_>=0) {
	/* update lambda  */
	double scale=1.0/mu;
	int i,nnz=0;
	saveLambdaScale=scale;
         lambdaIteration++;
         if ((saveStrategy&4)==0) drop = drop_/50.0;
         if (lambdaIteration>4 && 
            (((lambdaIteration%10)==0 && smallInfeas<keepinfeas) ||
             ((lambdaIteration%5)==0 && 1.5*smallInfeas<keepinfeas))) {
           //printf(" Increasing smallInfeas from %f to %f\n",smallInfeas,1.5*smallInfeas);
           smallInfeas *= 1.5;
         }
         if ((saveStrategy&2048)==0) {
	   for (i=0;i<nrows;i++) {
	     if (lambda[i]) nnz++;
	     lambda[i]+= scale*rowsol[i];
	   }
	 } else {
	   nnz=1;
#ifdef LAMBDA
	   for (i=0;i<nrows;i++) {
	     lambda[i]+= scale*rowsol[i];
	   }
#else
	   for (i=0;i<nrows;i++) {
	     lambda[i] = scale*rowsol[i];
	   }
	   for (i=0;i<ncols;i++) {
	     CoinBigIndex j;
	     double value=cost[i]*maxmin;
	     for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
	       int irow=row[j];
	       value+=element[j]*lambda[irow];
	     }
	     cost[i]=value*maxmin;
	   }
	   for (i=0;i<nrows;i++) {
	     offset+=lambda[i]*rowupper[i];
	     lambda[i]=0.0;
	   }
#ifdef DEBUG
           printf("offset %g\n",offset);
#endif
	   model_->setDblParam(OsiObjOffset,offset);
#endif
	 }
        nTry++;
	if (!nnz) {
	  bestFeasible=1.0e32;
	  bestWeighted=1.0e60;
	  checkIteration=0;
	  result.weighted=1.0e31;
	}
#ifdef DEBUG
        double trueCost=0.0;
	for (i=0;i<ncols;i++) {
	  int j;
	  trueCost+=cost[i]*colsol[i];
	}
	printf("True objective %g\n",trueCost-offset);
#endif
      } else {
        nTry++;
      }
      lastResult=result;
      if (result.infeas<reasonableInfeas&&!belowReasonable) {
	belowReasonable=1;
	bestFeasible=1.0e32;
        bestWeighted=1.0e60;
	checkIteration=0;
	result.weighted=1.0e31;
      }
    }
    if (iteration>=majorIterations_) {
      // If not feasible and crash then dive dive dive
      if (mu>1.0e-12&&result.infeas>1.0&&majorIterations_<40) {
	mu=1.0e-30;
	majorIterations_=iteration+1;
	stopMu=0.0;
      } else {
	if (logLevel>2) 
	  printf("Exiting due to number of major iterations\n");
	break;
      }
    }
  }
  majorIterations_ = saveMajorIterations;
#ifndef OSI_IDIOT
  if (scaled) {
    // Scale solution and free arrays
    const double * rowScale = model_->rowScale();
    const double * columnScale = model_->columnScale();
    int icol,irow;
    for (icol=0;icol<ncols;icol++) {
      colsol[icol] *= columnScale[icol];
      saveSol[icol] *= columnScale[icol];
      dj[icol] /= columnScale[icol];
    }
    for (irow=0;irow<nrows;irow++) {
      rowsol[irow] /= rowScale[irow];
      pi[irow] *= rowScale[irow];
    }
    // Don't know why getting Microsoft problems
#if defined (_MSC_VER)
    delete [] ( double *) elemXX;
#else
    delete [] elemXX;
#endif
    model_->setRowScale(NULL);
    model_->setColumnScale(NULL);
    delete [] lower;
    delete [] upper;
    delete [] cost;
    lower = model_->columnLower();
    upper = model_->columnUpper();
    cost = model_->objective();
    //rowlower = model_->rowLower();
  }
#endif
#define TRYTHIS
#ifdef TRYTHIS
  if ((saveStrategy&2048)!=0) {
    double offset;
    model_->getDblParam(OsiObjOffset,offset);
    for (i=0;i<ncols;i++) {
      CoinBigIndex j;
      double djval=cost[i]*maxmin;
      for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
	int irow=row[j];
	djval -= element[j]*lambda[irow];
      }
      cost[i]=djval;
    }
    for (i=0;i<nrows;i++) {
      offset+=lambda[i]*rowupper[i];
    }
    model_->setDblParam(OsiObjOffset,offset);
  }
#endif
  if (saveLambdaScale) {
    /* back off last update */
    for (i=0;i<nrows;i++) {
      lambda[i]-= saveLambdaScale*rowsol[i];
    }
  }
  muAtExit_=mu;
  // For last iteration make as feasible as possible
  if (oddSlacks)
    strategy_ |= 16384;
  // not scaled
  n = cleanIteration(iteration, ordStart,ordEnd,
		     colsol,  lower,  upper,
		     model_->rowLower(), model_->rowUpper(),
		     cost, element, fixTolerance,lastResult.objval,lastResult.infeas);
#if 0
  if ((logLevel&1)==0||(strategy_&16384)!=0) {
    printf(
	    "%d - mu %g, infeasibility %g, objective %g, %d interior\n",
	    iteration,mu,lastResult.infeas,lastResult.objval,n);
  }
#endif
#ifndef OSI_IDIOT
  model_->setSumPrimalInfeasibilities(lastResult.infeas);
#endif
  // Put back more feasible solution
  double saveInfeas[]={0.0,0.0};
  for (int iSol=0;iSol<3;iSol++) {
    const double * solution = iSol ? colsol : saveSol;
    if (iSol==2&&saveInfeas[0]<saveInfeas[1]) {
      // put back best solution
      memcpy(colsol,saveSol,ncols*sizeof(double));
    }
    double large=0.0;
    int i;
    memset(rowsol,0,nrows*sizeof(double));
    for (i=0;i<ncols;i++) {
      CoinBigIndex j;
      double value=solution[i];
      for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
	int irow=row[j];
	rowsol[irow] += element[j]*value;
      }
    }
    for (i=0;i<nrows;i++) {
      if (rowsol[i] > rowupper[i]) {
	double diff=rowsol[i]-rowupper[i];
	if (diff>large) 
	  large=diff;
      } else if (rowsol[i] < rowlower[i]) {
	double diff=rowlower[i]-rowsol[i];
	if (diff>large) 
	  large=diff;
      } 
    }
    if (iSol<2)
      saveInfeas[iSol]=large;
    if (logLevel>2)
      printf("largest infeasibility is %g\n",large);
  }
  /* subtract out lambda */
  for (i=0;i<nrows;i++) {
    pi[i]-=lambda[i];
  }
  for (i=0;i<ncols;i++) {
    CoinBigIndex j;
    double djval=cost[i]*maxmin;
    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
      int irow=row[j];
      djval -= element[j]*pi[irow];
    }
    dj[i]=djval;
  }
  if ((strategy_&1024)!=0) {
    double ratio = ((double) ncols)/((double) nrows);
    printf("col/row ratio %g infeas ratio %g\n",ratio,lastResult.infeas/firstInfeas);
    if (lastResult.infeas>0.01*firstInfeas*ratio) {
      strategy_ &= (~1024);
      printf(" - layer off\n");
    } else {
      printf(" - layer on\n");
    }
  }
  delete [] saveSol;
  delete [] lambda;
  // save solution 
  // duals not much use - but save anyway
#ifndef OSI_IDIOT
  memcpy(model_->primalRowSolution(),rowsol,nrows*sizeof(double));
  memcpy(model_->primalColumnSolution(),colsol,ncols*sizeof(double));
  memcpy(model_->dualRowSolution(),pi,nrows*sizeof(double));
  memcpy(model_->dualColumnSolution(),dj,ncols*sizeof(double));
#else
  model_->setColSolution(colsol);
  model_->setRowPrice(pi);
  delete [] cost;
#endif
  delete [] rowsol;
  delete [] colsol;
  delete [] pi;
  delete [] dj;
  delete [] rowlower;
  delete [] rowupper;
  return ;
}
#ifndef OSI_IDIOT
void
Idiot::crossOver(int mode)
{
  if (lightWeight_==2) {
    // total failure
    model_->allSlackBasis();
    return;
  }
  double fixTolerance=IDIOT_FIX_TOLERANCE;
  double startTime = CoinCpuTime();
  ClpSimplex * saveModel=NULL;
  ClpMatrixBase * matrix = model_->clpMatrix();
  const int * row = matrix->getIndices();
  const CoinBigIndex * columnStart = matrix->getVectorStarts();
  const int * columnLength = matrix->getVectorLengths(); 
  const double * element = matrix->getElements();
  const double * rowupper = model_->getRowUpper();
  int nrows=model_->getNumRows();
  int ncols=model_->getNumCols();
  double * rowsol, * colsol;
  // different for Osi
  double * lower = model_->columnLower();
  double * upper = model_->columnUpper();
  const double * rowlower= model_->getRowLower();
  int * whenUsed=whenUsed_;
  rowsol= model_->primalRowSolution();
  colsol= model_->primalColumnSolution();;
  double * cost=model_->objective();

  int slackEnd,ordStart,ordEnd;
  int slackStart = countCostedSlacks(model_);

  int addAll = mode&7;
  int presolve=0;

  double djTolerance = djTolerance_;
  if (djTolerance>0.0&&djTolerance<1.0)
    djTolerance=1.0;
  int iteration;
  int i, n=0;
  double ratio=1.0;
  double objValue=0.0;
  if ((strategy_&128)!=0) {
    fixTolerance=SMALL_IDIOT_FIX_TOLERANCE;
  }
  if ((mode&16)!=0&&addAll<3) presolve=1;
  double * saveUpper = NULL;
  double * saveLower = NULL;
  double * saveRowUpper = NULL;
  double * saveRowLower = NULL;
  bool allowInfeasible = (strategy_&8192)!=0;
  if (addAll<3) {
    saveUpper = new double [ncols];
    saveLower = new double [ncols];
    memcpy(saveUpper,upper,ncols*sizeof(double));
    memcpy(saveLower,lower,ncols*sizeof(double));
    if (allowInfeasible) {
      saveRowUpper = new double [nrows];
      saveRowLower = new double [nrows];
      memcpy(saveRowUpper,rowupper,nrows*sizeof(double));
      memcpy(saveRowLower,rowlower,nrows*sizeof(double));
      double averageInfeas = model_->sumPrimalInfeasibilities()/((double) model_->numberRows());
      fixTolerance = CoinMax(fixTolerance,1.0e-5*averageInfeas);
    }
  }
  if (slackStart>=0) {
    slackEnd=slackStart+nrows;
    if (slackStart) {
      ordStart=0;
      ordEnd=slackStart;
    } else {
      ordStart=nrows;
      ordEnd=ncols;
    }
  } else {
    slackEnd=slackStart;
    ordStart=0;
    ordEnd=ncols;
  }
  /* get correct rowsol (without known slacks) */
  memset(rowsol,0,nrows*sizeof(double));
  for (i=ordStart;i<ordEnd;i++) {
    CoinBigIndex j;
    double value=colsol[i];
    if (value<lower[i]+fixTolerance) {
      value=lower[i];
      colsol[i]=value;
    }
    for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
      int irow=row[j];
      rowsol[irow]+=value*element[j];
    }
  }
  if (slackStart>=0) {
    for (i=0;i<nrows;i++) {
      if (ratio*rowsol[i]>rowlower[i]&&rowsol[i]>1.0e-8) {
	ratio=rowlower[i]/rowsol[i];
      }
    }
    for (i=0;i<nrows;i++) {
      rowsol[i]*=ratio;
    }
    for (i=ordStart;i<ordEnd;i++) {
      double value=colsol[i]*ratio;
      colsol[i]=value;
      objValue+=value*cost[i];
    }
    for (i=0;i<nrows;i++) {
      double value=rowlower[i]-rowsol[i];
      colsol[i+slackStart]=value;
      objValue+=value*cost[i+slackStart];
    }
    printf("New objective after scaling %g\n",objValue);
  }
#if 0
   maybe put back - but just get feasible ?
  // If not many fixed then just exit
  int numberFixed=0;
  for (i=ordStart;i<ordEnd;i++) {
    if (colsol[i]<lower[i]+fixTolerance) 
      numberFixed++;
    else if (colsol[i]>upper[i]-fixTolerance) 
      numberFixed++;
  }
  if (numberFixed<ncols/2) {
    addAll=3;
    presolve=0;
  }
#endif
  model_->createStatus();
  /* addAll
     0 - chosen,all used, all 
     1 - chosen, all
     2 - all
     3 - do not do anything  - maybe basis
  */
  for (i=ordStart;i<ordEnd;i++) {
    if (addAll<2) {
      if (colsol[i]<lower[i]+fixTolerance) {
	upper[i]=lower[i];
	colsol[i]=lower[i];
      } else if (colsol[i]>upper[i]-fixTolerance) {
	lower[i]=upper[i];
	colsol[i]=upper[i];
      }
    }
    model_->setColumnStatus(i,ClpSimplex::superBasic);
  }
  if ((strategy_&16384)!=0) {
    // put in basis
    int * posSlack = whenUsed_+ncols;
    int * negSlack = posSlack+nrows;
    int * nextSlack = negSlack + nrows;
    for (i=0;i<nrows;i++) {
      int n=0;
      int iCol;
      iCol =posSlack[i];
      if (iCol>=0) {
	if (colsol[iCol]>lower[iCol]+1.0e-8&&
	    colsol[iCol]<upper[iCol]+1.0e-8) {
	  model_->setColumnStatus(iCol,ClpSimplex::basic);
	  n++;
	}
	while (nextSlack[iCol]>=0) {
	  iCol = nextSlack[iCol];
	  if (colsol[iCol]>lower[iCol]+1.0e-8&&
	      colsol[iCol]<upper[iCol]+1.0e-8) {
	    model_->setColumnStatus(iCol,ClpSimplex::basic);
	    n++;
	  }
	}
      }
      iCol =negSlack[i];
      if (iCol>=0) {
	if (colsol[iCol]>lower[iCol]+1.0e-8&&
	    colsol[iCol]<upper[iCol]+1.0e-8) {
	  model_->setColumnStatus(iCol,ClpSimplex::basic);
	  n++;
	}
	while (nextSlack[iCol]>=0) {
	  iCol = nextSlack[iCol];
	  if (colsol[iCol]>lower[iCol]+1.0e-8&&
	      colsol[iCol]<upper[iCol]+1.0e-8) {
	    model_->setColumnStatus(iCol,ClpSimplex::basic);
	    n++;
	  }
	}
      }
      if (n) {
	if (fabs(rowsol[i]-rowlower[i])<fabs(rowsol[i]-rowupper[i]))
	  model_->setRowStatus(i,ClpSimplex::atLowerBound);
	else
	  model_->setRowStatus(i,ClpSimplex::atUpperBound);
	if (n>1)
	  printf("%d basic on row %d!\n",n,i);
      }
    }
  }
  double maxmin;
  if (model_->getObjSense()==-1.0) {
    maxmin=-1.0;
  } else {
    maxmin=1.0;
  }
  if (slackStart>=0) {
    for (i=0;i<nrows;i++) {
      model_->setRowStatus(i,ClpSimplex::superBasic);
    }
    for (i=slackStart;i<slackEnd;i++) {
      model_->setColumnStatus(i,ClpSimplex::basic);
    }
  } else {
    /* still try and put singletons rather than artificials in basis */
    int ninbas=0;
    for (i=0;i<nrows;i++) {
      model_->setRowStatus(i,ClpSimplex::basic);
    }
    for (i=0;i<ncols;i++) {
      if (columnLength[i]==1&&upper[i]>lower[i]+1.0e-5) {
	CoinBigIndex j =columnStart[i];
	double value=element[j];
	int irow=row[j];
	double rlo=rowlower[irow];
	double rup=rowupper[irow];
	double clo=lower[i];
	double cup=upper[i];
	double csol=colsol[i];
	/* adjust towards feasibility */
	double move=0.0;
	if (rowsol[irow]>rup) {
	  move=(rup-rowsol[irow])/value;
	  if (value>0.0) {
	    /* reduce */
	    if (csol+move<clo) move=CoinMin(0.0,clo-csol);
	  } else {
	    /* increase */
	    if (csol+move>cup) move=CoinMax(0.0,cup-csol);
	  }
	} else if (rowsol[irow]<rlo) {
	  move=(rlo-rowsol[irow])/value;
	  if (value>0.0) {
	    /* increase */
	    if (csol+move>cup) move=CoinMax(0.0,cup-csol);
	  } else {
	    /* reduce */
	    if (csol+move<clo) move=CoinMin(0.0,clo-csol);
	  }
	} else {
	  /* move to improve objective */
	  if (cost[i]*maxmin>0.0) {
	    if (value>0.0) {
	      move=(rlo-rowsol[irow])/value;
	      /* reduce */
	      if (csol+move<clo) move=CoinMin(0.0,clo-csol);
	    } else {
	      move=(rup-rowsol[irow])/value;
	      /* increase */
	      if (csol+move>cup) move=CoinMax(0.0,cup-csol);
	    }
	  } else if (cost[i]*maxmin<0.0) {
	    if (value>0.0) {
	      move=(rup-rowsol[irow])/value;
	      /* increase */
	      if (csol+move>cup) move=CoinMax(0.0,cup-csol);
	    } else {
	      move=(rlo-rowsol[irow])/value;
	      /* reduce */
	      if (csol+move<clo) move=CoinMin(0.0,clo-csol);
	    }
	  }
	}
	rowsol[irow] +=move*value;
	colsol[i]+=move;
	/* put in basis if row was artificial */
	if (rup-rlo<1.0e-7&&model_->getRowStatus(irow)==ClpSimplex::basic) {
	  model_->setRowStatus(irow,ClpSimplex::superBasic);
	  model_->setColumnStatus(i,ClpSimplex::basic);
	  ninbas++;
	}
      }
    }
    /*printf("%d in basis\n",ninbas);*/
  }
  bool wantVector=false;
  if (dynamic_cast< ClpPackedMatrix*>(model_->clpMatrix())) {
    // See if original wanted vector
    ClpPackedMatrix * clpMatrixO = dynamic_cast< ClpPackedMatrix*>(model_->clpMatrix());
    wantVector = clpMatrixO->wantsSpecialColumnCopy();
  }
  if (addAll<3) {
    ClpPresolve pinfo;
    if (presolve) {
      if (allowInfeasible) {
	// fix up so will be feasible
	double * rhs = new double[nrows];
	memset(rhs,0,nrows*sizeof(double));
	model_->clpMatrix()->times(1.0,colsol,rhs);
	double * rowupper = model_->rowUpper();
	double * rowlower= model_->rowLower();
	double sum = 0.0;
	for (i=0;i<nrows;i++) {
	  if (rhs[i]>rowupper[i]) {
	    sum += rhs[i]-rowupper[i];
	    rowupper[i]=rhs[i];
	  }
	  if (rhs[i]<rowlower[i]) {
	    sum += rowlower[i]-rhs[i];
	    rowlower[i]=rhs[i];
	  }
	}
	printf("sum of infeasibilities %g\n",sum);
	delete [] rhs;
      }
      saveModel = model_;
      pinfo.setPresolveActions(pinfo.presolveActions()|16384);
      model_ = pinfo.presolvedModel(*model_,1.0e-8,false,5);
    }
    if (model_) {
      if (!wantVector) {
	model_->primal(1);
      } else {
	ClpMatrixBase * matrix = model_->clpMatrix();
	ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
	assert (clpMatrix);
	clpMatrix->makeSpecialColumnCopy();
	model_->primal(1);
	clpMatrix->releaseSpecialColumnCopy();
      }
      if (presolve) {
	pinfo.postsolve(true);
	delete model_;
	model_ = saveModel;
	saveModel=NULL;
      }
    } else {
      // not feasible
      addAll=1;
      presolve=0;
      model_ = saveModel;
      saveModel=NULL;
    }
    if (allowInfeasible) {
      memcpy(model_->rowUpper(),saveRowUpper,nrows*sizeof(double));
      memcpy(model_->rowLower(),saveRowLower,nrows*sizeof(double));
      delete [] saveRowUpper;
      delete [] saveRowLower;
      saveRowUpper = NULL;
      saveRowLower = NULL;
    }
    if (addAll<2) {
      n=0;
      if (!addAll ) {
	/* could do scans to get a good number */
	iteration=1;
	for (i=ordStart;i<ordEnd;i++) {
	  if (whenUsed[i]>=iteration) {
	    if (upper[i]-lower[i]<1.0e-5&&saveUpper[i]-saveLower[i]>1.0e-5) {
	      n++;
	      upper[i]=saveUpper[i];
	      lower[i]=saveLower[i];
	    }
	  }
	}
      } else {
	for (i=ordStart;i<ordEnd;i++) {
	  if (upper[i]-lower[i]<1.0e-5&&saveUpper[i]-saveLower[i]>1.0e-5) {
	    n++;
	    upper[i]=saveUpper[i];
	    lower[i]=saveLower[i];
	  }
	}
	delete [] saveUpper;
	delete [] saveLower;
	saveUpper=NULL;
	saveLower=NULL;
      }
      printf("Time so far %g, %d now added from previous iterations\n",
	     CoinCpuTime()-startTime,n);
      if (addAll)
	presolve=0;
      if (presolve) {
	saveModel = model_;
	model_ = pinfo.presolvedModel(*model_,1.0e-8,false,5);
      } else {
	presolve=0;
      }
      if (!wantVector) {
	model_->primal(1);
      } else {
	ClpMatrixBase * matrix = model_->clpMatrix();
	ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
	assert (clpMatrix);
	clpMatrix->makeSpecialColumnCopy();
	model_->primal(1);
	clpMatrix->releaseSpecialColumnCopy();
      }
      if (presolve) {
	pinfo.postsolve(true);
	delete model_;
	model_ = saveModel;
	saveModel=NULL;
      }
      if (!addAll) {
	n=0;
	for (i=ordStart;i<ordEnd;i++) {
	  if (upper[i]-lower[i]<1.0e-5&&saveUpper[i]-saveLower[i]>1.0e-5) {
	    n++;
	    upper[i]=saveUpper[i];
	    lower[i]=saveLower[i];
	  }
	}
	delete [] saveUpper;
	delete [] saveLower;
	saveUpper=NULL;
	saveLower=NULL;
	printf("Time so far %g, %d now added from previous iterations\n",
	       CoinCpuTime()-startTime,n);
      }
      if (presolve) {
	saveModel = model_;
	model_ = pinfo.presolvedModel(*model_,1.0e-8,false,5);
      } else {
	presolve=0;
      }
      if (!wantVector) {
	model_->primal(1);
      } else {
	ClpMatrixBase * matrix = model_->clpMatrix();
	ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
	assert (clpMatrix);
	clpMatrix->makeSpecialColumnCopy();
	model_->primal(1);
	clpMatrix->releaseSpecialColumnCopy();
      }
      if (presolve) {
	pinfo.postsolve(true);
	delete model_;
	model_ = saveModel;
	saveModel=NULL;
      }
    }
    printf("Total time in crossover %g\n", CoinCpuTime()-startTime);
    delete [] saveUpper;
    delete [] saveLower;
  }
  return ;
}
#endif
/*****************************************************************************/

// Default contructor
Idiot::Idiot()
{
  model_ = NULL;
  maxBigIts_ = 3;
  maxIts_ = 5;
  logLevel_ = 1; 
  logFreq_ = 100;
  maxIts2_ = 100;
  djTolerance_ = 1e-1;
  mu_ = 1e-4;
  drop_ = 5.0;
  exitDrop_=-1.0e20;
  muFactor_ = 0.3333;
  stopMu_ = 1e-12;
  smallInfeas_ = 1e-1;
  reasonableInfeas_ = 1e2;
  muAtExit_ =1.0e31;
  strategy_ =8;
  lambdaIterations_ =0;
  checkFrequency_ =100;
  whenUsed_ = NULL;
  majorIterations_ =30;
  exitFeasibility_ =-1.0;
  dropEnoughFeasibility_ =0.02;
  dropEnoughWeighted_ =0.01;
  // adjust
  double nrows=10000.0;
  int baseIts =(int) sqrt((double)nrows);
  baseIts =baseIts/10;
  baseIts *= 10;
  maxIts2_ =200+baseIts+5;
  maxIts2_=100;
  reasonableInfeas_ =((double) nrows)*0.05;
  lightWeight_=0;
}
// Constructor from model
Idiot::Idiot(OsiSolverInterface &model)
{
  model_ = & model;
  maxBigIts_ = 3;
  maxIts_ = 5;
  logLevel_ = 1; 
  logFreq_ = 100;
  maxIts2_ = 100;
  djTolerance_ = 1e-1;
  mu_ = 1e-4;
  drop_ = 5.0;
  exitDrop_=-1.0e20;
  muFactor_ = 0.3333;
  stopMu_ = 1e-12;
  smallInfeas_ = 1e-1;
  reasonableInfeas_ = 1e2;
  muAtExit_ =1.0e31;
  strategy_ =8;
  lambdaIterations_ =0;
  checkFrequency_ =100;
  whenUsed_ = NULL;
  majorIterations_ =30;
  exitFeasibility_ =-1.0;
  dropEnoughFeasibility_ =0.02;
  dropEnoughWeighted_ =0.01;
  // adjust
  double nrows;
  if (model_)
    nrows=model_->getNumRows();
  else
    nrows=10000.0;
  int baseIts =(int) sqrt((double)nrows);
  baseIts =baseIts/10;
  baseIts *= 10;
  maxIts2_ =200+baseIts+5;
  maxIts2_=100;
  reasonableInfeas_ =((double) nrows)*0.05;
  lightWeight_=0;
}
// Copy constructor. 
Idiot::Idiot(const Idiot &rhs)
{
  model_ = rhs.model_;
  if (model_&&rhs.whenUsed_) {
    int numberColumns = model_->getNumCols();
    whenUsed_ = new int [numberColumns];
    memcpy(whenUsed_,rhs.whenUsed_,numberColumns*sizeof(int));
  } else {
    whenUsed_=NULL;
  }
  djTolerance_ = rhs.djTolerance_;
  mu_ = rhs.mu_;
  drop_ = rhs.drop_;
  muFactor_ = rhs.muFactor_;
  stopMu_ = rhs.stopMu_;
  smallInfeas_ = rhs.smallInfeas_;
  reasonableInfeas_ = rhs.reasonableInfeas_;
  exitDrop_ = rhs.exitDrop_;
  muAtExit_ = rhs.muAtExit_;
  exitFeasibility_ = rhs.exitFeasibility_;
  dropEnoughFeasibility_ = rhs.dropEnoughFeasibility_;
  dropEnoughWeighted_ = rhs.dropEnoughWeighted_;
  maxBigIts_ = rhs.maxBigIts_;
  maxIts_ = rhs.maxIts_;
  majorIterations_ = rhs.majorIterations_;
  logLevel_ = rhs.logLevel_;
  logFreq_ = rhs.logFreq_;
  checkFrequency_ = rhs.checkFrequency_;
  lambdaIterations_ = rhs.lambdaIterations_;
  maxIts2_ = rhs.maxIts2_;
  strategy_ = rhs.strategy_;
  lightWeight_=rhs.lightWeight_;
}
// Assignment operator. This copies the data
Idiot & 
Idiot::operator=(const Idiot & rhs)
{
  if (this != &rhs) {
    delete [] whenUsed_;
    model_ = rhs.model_;
    if (model_&&rhs.whenUsed_) {
      int numberColumns = model_->getNumCols();
      whenUsed_ = new int [numberColumns];
      memcpy(whenUsed_,rhs.whenUsed_,numberColumns*sizeof(int));
    } else {
      whenUsed_=NULL;
    }
    djTolerance_ = rhs.djTolerance_;
    mu_ = rhs.mu_;
    drop_ = rhs.drop_;
    muFactor_ = rhs.muFactor_;
    stopMu_ = rhs.stopMu_;
    smallInfeas_ = rhs.smallInfeas_;
    reasonableInfeas_ = rhs.reasonableInfeas_;
    exitDrop_ = rhs.exitDrop_;
    muAtExit_ = rhs.muAtExit_;
    exitFeasibility_ = rhs.exitFeasibility_;
    dropEnoughFeasibility_ = rhs.dropEnoughFeasibility_;
    dropEnoughWeighted_ = rhs.dropEnoughWeighted_;
    maxBigIts_ = rhs.maxBigIts_;
    maxIts_ = rhs.maxIts_;
    majorIterations_ = rhs.majorIterations_;
    logLevel_ = rhs.logLevel_;
    logFreq_ = rhs.logFreq_;
    checkFrequency_ = rhs.checkFrequency_;
    lambdaIterations_ = rhs.lambdaIterations_;
    maxIts2_ = rhs.maxIts2_;
    strategy_ = rhs.strategy_;
    lightWeight_=rhs.lightWeight_;
  }
  return *this;
}
Idiot::~Idiot()
{
  delete [] whenUsed_;
}