Changeset 162


Ignore:
Timestamp:
Apr 15, 2003 11:57:38 AM (16 years ago)
Author:
forrest
Message:

Piecewise linear

Location:
trunk
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • trunk/ClpNonLinearCost.cpp

    r159 r162  
    3131  infeasible_(NULL),
    3232  numberInfeasibilities_(-1),
    33   convex_(true)
     33  convex_(true),
     34  bothWays_(false)
    3435{
    3536
     
    4647  int numberTotal = numberRows_+numberColumns_;
    4748  convex_ = true;
     49  bothWays_ = false;
    4850  start_ = new int [numberTotal+1];
    4951  whichRange_ = new int [numberTotal];
     
    111113  int numberTotal = numberRows_+numberColumns_;
    112114  convex_ = true;
     115  bothWays_ = true;
    113116  start_ = new int [numberTotal+1];
    114117  whichRange_ = new int [numberTotal];
     
    171174      upperValue = rowUpper[iSequence-numberColumns_];
    172175      if (lowerValue>-1.0e30) {
     176        setInfeasible(put,true);
    173177        cost_[put++] = -infeasibilityCost;
    174178        lower_[put] = lowerValue;
     
    181185      upperValue = columnUpper[iSequence];
    182186      if (lowerValue>-1.0e30) {
     187        setInfeasible(put,true);
    183188        cost_[put++] = whichWay*cost[iSequence]-infeasibilityCost;
    184189        lower_[put] = lowerValue;
     
    202207    }
    203208    lower_[put] = upperValue;
     209    setInfeasible(put,true);
    204210    cost_[put++] = thisCost+infeasibilityCost;
    205211    if (upperValue<1.0e20) {
    206212      lower_[put] = COIN_DBL_MAX;
    207       setInfeasible(put-1,true);
    208213      cost_[put++] = 1.0e50;
    209214    }
     
    238243  infeasible_(NULL),
    239244  numberInfeasibilities_(-1),
    240   convex_(true)
     245  convex_(true),
     246  bothWays_(rhs.bothWays_)
    241247
    242248  if (numberRows_) {
     
    321327    sumInfeasibilities_ = rhs.sumInfeasibilities_;
    322328    convex_ = rhs.convex_;
     329    bothWays_ = rhs.bothWays_;
    323330  }
    324331  return *this;
     
    362369      if (value<lower_[iRange+1]+primalTolerance) {
    363370        // put in better range if infeasible
    364         if (value>=lower_[iRange+1]-primalTolerance&&iRange==start)
     371        if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start)
    365372          iRange++;
    366373        whichRange_[iSequence]=iRange;
     
    409416        }
    410417      }
    411       lower[iSequence] = lower_[iRange];
    412       upper[iSequence] = lower_[iRange+1];
    413       cost[iSequence] = cost_[iRange];
     418      //lower[iSequence] = lower_[iRange];
     419      //upper[iSequence] = lower_[iRange+1];
     420      //cost[iSequence] = cost_[iRange];
    414421      break;
    415422    case ClpSimplex::isFree:
     
    452459      break;
    453460    }
     461    lower[iSequence] = lower_[iRange];
     462    upper[iSequence] = lower_[iRange+1];
     463    cost[iSequence] = cost_[iRange];
    454464  }
    455465}
     
    562572      if (value<lower_[iRange+1]+primalTolerance) {
    563573        // put in better range
    564         if (value>=lower_[iRange+1]-primalTolerance&&iRange==start)
     574        if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start)
    565575          iRange++;
    566576        break;
     
    605615      if (value<lower_[iRange+1]+primalTolerance) {
    606616        // put in better range
    607         if (value>=lower_[iRange+1]-primalTolerance&&iRange==start)
     617        if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start)
    608618          iRange++;
    609619        break;
     
    638648  int start = start_[iPivot];
    639649  int end = start_[iPivot+1]-1;
    640   for (iRange=start; iRange<end;iRange++) {
    641     if (value<lower_[iRange+1]+primalTolerance) {
    642       // put in better range
    643       if (value>=lower_[iRange+1]-primalTolerance&&iRange==start)
    644         iRange++;
    645       break;
     650  if (!bothWays_) {
     651    for (iRange=start; iRange<end;iRange++) {
     652      if (value<lower_[iRange+1]+primalTolerance) {
     653        // put in better range
     654        if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start)
     655          iRange++;
     656        break;
     657      }
     658    }
     659  } else {
     660    // leave in current if possible
     661    iRange = whichRange_[iPivot];
     662    if (value<lower_[iRange]-primalTolerance||value>lower_[iRange+1]+primalTolerance) {
     663      for (iRange=start; iRange<end;iRange++) {
     664        if (value<lower_[iRange+1]+primalTolerance) {
     665          // put in better range
     666          if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start)
     667            iRange++;
     668          break;
     669        }
     670      }
    646671    }
    647672  }
  • trunk/ClpPrimalColumnSteepest.cpp

    r137 r162  
    88#include "CoinIndexedVector.hpp"
    99#include "ClpFactorization.hpp"
     10#include "ClpNonLinearCost.hpp"
    1011#include "ClpMessage.hpp"
    1112#include "CoinHelperFunctions.hpp"
     
    208209    model_->clpMatrix()->transposeTimes(model_,-1.0,
    209210                                        updates,spareColumn2,spareColumn1);
     211    // normal
    210212    for (iSection=0;iSection<2;iSection++) {
    211213     
     
    224226        addSequence = 0;
    225227      }
    226 
    227       for (j=0;j<number;j++) {
    228         int iSequence = index[j];
    229         double value = reducedCost[iSequence];
    230         value -= updateBy[iSequence];
    231         reducedCost[iSequence] = value;
    232         ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
    233 
    234         switch(status) {
     228      if (!model_->nonLinearCost()->lookBothWays()) {
     229       
     230        for (j=0;j<number;j++) {
     231          int iSequence = index[j];
     232          double value = reducedCost[iSequence];
     233          value -= updateBy[iSequence];
     234          reducedCost[iSequence] = value;
     235          ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
    235236         
    236         case ClpSimplex::basic:
    237           infeasible_->zero(iSequence+addSequence);
    238         case ClpSimplex::isFixed:
    239           break;
    240         case ClpSimplex::isFree:
    241         case ClpSimplex::superBasic:
    242           if (fabs(value)>1.0e2*tolerance) {
    243             // we are going to bias towards free (but only if reasonable)
    244             value *= FREE_BIAS;
    245             // store square in list
    246             if (infeas[iSequence+addSequence])
    247               infeas[iSequence+addSequence] = value*value; // already there
    248             else
    249               infeasible_->quickAdd(iSequence+addSequence,value*value);
    250           } else {
     237          switch(status) {
     238           
     239          case ClpSimplex::basic:
    251240            infeasible_->zero(iSequence+addSequence);
    252           }
    253           break;
    254         case ClpSimplex::atUpperBound:
    255           if (value>tolerance) {
    256             // store square in list
    257             if (infeas[iSequence+addSequence])
    258               infeas[iSequence+addSequence] = value*value; // already there
    259             else
    260               infeasible_->quickAdd(iSequence+addSequence,value*value);
    261           } else {
     241          case ClpSimplex::isFixed:
     242            break;
     243          case ClpSimplex::isFree:
     244          case ClpSimplex::superBasic:
     245            if (fabs(value)>1.0e2*tolerance) {
     246              // we are going to bias towards free (but only if reasonable)
     247              value *= FREE_BIAS;
     248              // store square in list
     249              if (infeas[iSequence+addSequence])
     250                infeas[iSequence+addSequence] = value*value; // already there
     251              else
     252                infeasible_->quickAdd(iSequence+addSequence,value*value);
     253            } else {
     254              infeasible_->zero(iSequence+addSequence);
     255            }
     256            break;
     257          case ClpSimplex::atUpperBound:
     258            if (value>tolerance) {
     259              // store square in list
     260              if (infeas[iSequence+addSequence])
     261                infeas[iSequence+addSequence] = value*value; // already there
     262              else
     263                infeasible_->quickAdd(iSequence+addSequence,value*value);
     264            } else {
     265              infeasible_->zero(iSequence+addSequence);
     266            }
     267            break;
     268          case ClpSimplex::atLowerBound:
     269            if (value<-tolerance) {
     270              // store square in list
     271              if (infeas[iSequence+addSequence])
     272                infeas[iSequence+addSequence] = value*value; // already there
     273              else
     274                infeasible_->quickAdd(iSequence+addSequence,value*value);
     275            } else {
     276              infeasible_->zero(iSequence+addSequence);
     277            }
     278          }
     279        }
     280      } else {
     281        ClpNonLinearCost * nonLinear = model_->nonLinearCost();
     282        // We can go up OR down
     283        for (j=0;j<number;j++) {
     284          int iSequence = index[j];
     285          double value = reducedCost[iSequence];
     286          value -= updateBy[iSequence];
     287          reducedCost[iSequence] = value;
     288          ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
     289         
     290          switch(status) {
     291           
     292          case ClpSimplex::basic:
    262293            infeasible_->zero(iSequence+addSequence);
    263           }
    264           break;
    265         case ClpSimplex::atLowerBound:
    266           if (value<-tolerance) {
    267             // store square in list
    268             if (infeas[iSequence+addSequence])
    269               infeas[iSequence+addSequence] = value*value; // already there
    270             else
    271               infeasible_->quickAdd(iSequence+addSequence,value*value);
    272           } else {
    273             infeasible_->zero(iSequence+addSequence);
     294          case ClpSimplex::isFixed:
     295            break;
     296          case ClpSimplex::isFree:
     297          case ClpSimplex::superBasic:
     298            if (fabs(value)>1.0e2*tolerance) {
     299              // we are going to bias towards free (but only if reasonable)
     300              value *= FREE_BIAS;
     301              // store square in list
     302              if (infeas[iSequence+addSequence])
     303                infeas[iSequence+addSequence] = value*value; // already there
     304              else
     305                infeasible_->quickAdd(iSequence+addSequence,value*value);
     306            } else {
     307              infeasible_->zero(iSequence+addSequence);
     308            }
     309            break;
     310          case ClpSimplex::atUpperBound:
     311            if (value>tolerance) {
     312              // store square in list
     313              if (infeas[iSequence+addSequence])
     314                infeas[iSequence+addSequence] = value*value; // already there
     315              else
     316                infeasible_->quickAdd(iSequence+addSequence,value*value);
     317            } else {
     318              // look other way - change up should be negative
     319              value -= nonLinear->changeUpInCost(iSequence+addSequence);
     320              if (value>-tolerance) {
     321                infeasible_->zero(iSequence+addSequence);
     322              } else {
     323                // store square in list
     324                if (infeas[iSequence+addSequence])
     325                  infeas[iSequence+addSequence] = value*value; // already there
     326                else
     327                  infeasible_->quickAdd(iSequence+addSequence,value*value);
     328              }
     329            }
     330            break;
     331          case ClpSimplex::atLowerBound:
     332            if (value<-tolerance) {
     333              // store square in list
     334              if (infeas[iSequence+addSequence])
     335                infeas[iSequence+addSequence] = value*value; // already there
     336              else
     337                infeasible_->quickAdd(iSequence+addSequence,value*value);
     338            } else {
     339              // look other way - change down should be positive
     340              value -= nonLinear->changeDownInCost(iSequence+addSequence);
     341              if (value<tolerance) {
     342                infeasible_->zero(iSequence+addSequence);
     343              } else {
     344                // store square in list
     345                if (infeas[iSequence+addSequence])
     346                  infeas[iSequence+addSequence] = value*value; // already there
     347                else
     348                  infeasible_->quickAdd(iSequence+addSequence,value*value);
     349              }
     350            }
    274351          }
    275352        }
     
    351428    }
    352429
    353     for (iSequence=0;iSequence<number;iSequence++) {
    354       double value = reducedCost[iSequence];
    355       ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
    356      
    357       switch(status) {
    358 
    359       case ClpSimplex::basic:
    360         infeasible_->zero(iSequence+addSequence);
    361       case ClpSimplex::isFixed:
    362         break;
    363       case ClpSimplex::isFree:
    364         if (fabs(value)>tolerance) {
    365           // we are going to bias towards free (but only if reasonable)
    366           value *= FREE_BIAS;
    367           // store square in list
    368           if (infeas[iSequence+addSequence])
    369             infeas[iSequence+addSequence] = value*value; // already there
    370           else
    371             infeasible_->quickAdd(iSequence+addSequence,value*value);
     430    if (!model_->nonLinearCost()->lookBothWays()) {
     431      for (iSequence=0;iSequence<number;iSequence++) {
     432        double value = reducedCost[iSequence];
     433        ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
     434       
     435        switch(status) {
     436         
     437        case ClpSimplex::basic:
     438          infeasible_->zero(iSequence+addSequence);
     439        case ClpSimplex::isFixed:
     440          break;
     441        case ClpSimplex::isFree:
     442          if (fabs(value)>tolerance) {
     443            // we are going to bias towards free (but only if reasonable)
     444            value *= FREE_BIAS;
     445            // store square in list
     446            if (infeas[iSequence+addSequence])
     447              infeas[iSequence+addSequence] = value*value; // already there
     448            else
     449              infeasible_->quickAdd(iSequence+addSequence,value*value);
    372450          } else {
    373451            infeasible_->zero(iSequence+addSequence);
    374452          }
    375         break;
    376       case ClpSimplex::atUpperBound:
    377         if (value>tolerance) {
     453          break;
     454        case ClpSimplex::atUpperBound:
     455          if (value>tolerance) {
    378456            // store square in list
    379           if (infeas[iSequence+addSequence])
    380             infeas[iSequence+addSequence] = value*value; // already there
    381           else
    382             infeasible_->quickAdd(iSequence+addSequence,value*value);
    383         } else {
     457            if (infeas[iSequence+addSequence])
     458              infeas[iSequence+addSequence] = value*value; // already there
     459            else
     460              infeasible_->quickAdd(iSequence+addSequence,value*value);
     461          } else {
     462            infeasible_->zero(iSequence+addSequence);
     463          }
     464          break;
     465        case ClpSimplex::atLowerBound:
     466          if (value<-tolerance) {
     467            // store square in list
     468            if (infeas[iSequence+addSequence])
     469              infeas[iSequence+addSequence] = value*value; // already there
     470            else
     471              infeasible_->quickAdd(iSequence+addSequence,value*value);
     472          } else {
     473            infeasible_->zero(iSequence+addSequence);
     474          }
     475        }
     476      }
     477    } else {
     478      // we can go both ways
     479      ClpNonLinearCost * nonLinear = model_->nonLinearCost();
     480      for (iSequence=0;iSequence<number;iSequence++) {
     481        double value = reducedCost[iSequence];
     482        ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
     483       
     484        switch(status) {
     485         
     486        case ClpSimplex::basic:
    384487          infeasible_->zero(iSequence+addSequence);
    385         }
    386         break;
    387       case ClpSimplex::atLowerBound:
    388         if (value<-tolerance) {
    389           // store square in list
    390           if (infeas[iSequence+addSequence])
    391             infeas[iSequence+addSequence] = value*value; // already there
    392           else
    393             infeasible_->quickAdd(iSequence+addSequence,value*value);
    394         } else {
    395           infeasible_->zero(iSequence+addSequence);
     488        case ClpSimplex::isFixed:
     489          break;
     490        case ClpSimplex::isFree:
     491          if (fabs(value)>tolerance) {
     492            // we are going to bias towards free (but only if reasonable)
     493            value *= FREE_BIAS;
     494            // store square in list
     495            if (infeas[iSequence+addSequence])
     496              infeas[iSequence+addSequence] = value*value; // already there
     497            else
     498              infeasible_->quickAdd(iSequence+addSequence,value*value);
     499          } else {
     500            infeasible_->zero(iSequence+addSequence);
     501          }
     502          break;
     503        case ClpSimplex::atUpperBound:
     504          if (value>tolerance) {
     505            // store square in list
     506            if (infeas[iSequence+addSequence])
     507              infeas[iSequence+addSequence] = value*value; // already there
     508            else
     509              infeasible_->quickAdd(iSequence+addSequence,value*value);
     510          } else {
     511            // look other way - change up should be negative
     512            value -= nonLinear->changeUpInCost(iSequence+addSequence);
     513            if (value>-tolerance) {
     514              infeasible_->zero(iSequence+addSequence);
     515            } else {
     516              // store square in list
     517              if (infeas[iSequence+addSequence])
     518                infeas[iSequence+addSequence] = value*value; // already there
     519              else
     520                infeasible_->quickAdd(iSequence+addSequence,value*value);
     521            }
     522          }
     523          break;
     524        case ClpSimplex::atLowerBound:
     525          if (value<-tolerance) {
     526            // store square in list
     527            if (infeas[iSequence+addSequence])
     528              infeas[iSequence+addSequence] = value*value; // already there
     529            else
     530              infeasible_->quickAdd(iSequence+addSequence,value*value);
     531          } else {
     532            // look other way - change down should be positive
     533            value -= nonLinear->changeDownInCost(iSequence+addSequence);
     534            if (value<tolerance) {
     535              infeasible_->zero(iSequence+addSequence);
     536            } else {
     537              // store square in list
     538              if (infeas[iSequence+addSequence])
     539                infeas[iSequence+addSequence] = value*value; // already there
     540              else
     541                infeasible_->quickAdd(iSequence+addSequence,value*value);
     542            }
     543          }
    396544        }
    397545      }
     
    683831    double * reducedCost = model_->djRegion();
    684832     
    685     for (iSequence=0;iSequence<number;iSequence++) {
    686       double value = reducedCost[iSequence];
    687       ClpSimplex::Status status = model_->getStatus(iSequence);
    688      
    689       switch(status) {
    690 
    691       case ClpSimplex::basic:
    692       case ClpSimplex::isFixed:
    693         break;
    694       case ClpSimplex::isFree:
    695       case ClpSimplex::superBasic:
    696         if (fabs(value)>1.0e2*tolerance) {
    697           // we are going to bias towards free (but only if reasonable)
    698           value *= FREE_BIAS;
    699           // store square in list
    700           infeasible_->quickAdd(iSequence,value*value);
     833    if (!model_->nonLinearCost()->lookBothWays()) {
     834      for (iSequence=0;iSequence<number;iSequence++) {
     835        double value = reducedCost[iSequence];
     836        ClpSimplex::Status status = model_->getStatus(iSequence);
     837       
     838        switch(status) {
     839         
     840        case ClpSimplex::basic:
     841        case ClpSimplex::isFixed:
     842          break;
     843        case ClpSimplex::isFree:
     844        case ClpSimplex::superBasic:
     845          if (fabs(value)>1.0e2*tolerance) {
     846            // we are going to bias towards free (but only if reasonable)
     847            value *= FREE_BIAS;
     848            // store square in list
     849            infeasible_->quickAdd(iSequence,value*value);
     850          }
     851          break;
     852        case ClpSimplex::atUpperBound:
     853          if (value>tolerance) {
     854            infeasible_->quickAdd(iSequence,value*value);
     855          }
     856          break;
     857        case ClpSimplex::atLowerBound:
     858          if (value<-tolerance) {
     859            infeasible_->quickAdd(iSequence,value*value);
     860          }
    701861        }
    702         break;
    703       case ClpSimplex::atUpperBound:
    704         if (value>tolerance) {
    705           infeasible_->quickAdd(iSequence,value*value);
    706         }
    707         break;
    708       case ClpSimplex::atLowerBound:
    709         if (value<-tolerance) {
    710           infeasible_->quickAdd(iSequence,value*value);
     862      }
     863    } else {
     864      ClpNonLinearCost * nonLinear = model_->nonLinearCost();
     865      // can go both ways
     866      for (iSequence=0;iSequence<number;iSequence++) {
     867        double value = reducedCost[iSequence];
     868        ClpSimplex::Status status = model_->getStatus(iSequence);
     869       
     870        switch(status) {
     871         
     872        case ClpSimplex::basic:
     873        case ClpSimplex::isFixed:
     874          break;
     875        case ClpSimplex::isFree:
     876        case ClpSimplex::superBasic:
     877          if (fabs(value)>1.0e2*tolerance) {
     878            // we are going to bias towards free (but only if reasonable)
     879            value *= FREE_BIAS;
     880            // store square in list
     881            infeasible_->quickAdd(iSequence,value*value);
     882          }
     883          break;
     884        case ClpSimplex::atUpperBound:
     885          if (value>tolerance) {
     886            infeasible_->quickAdd(iSequence,value*value);
     887          } else {
     888            // look other way - change up should be negative
     889            value -= nonLinear->changeUpInCost(iSequence);
     890            if (value<-tolerance) {
     891              // store square in list
     892              infeasible_->quickAdd(iSequence,value*value);
     893            }
     894          }
     895          break;
     896        case ClpSimplex::atLowerBound:
     897          if (value<-tolerance) {
     898            infeasible_->quickAdd(iSequence,value*value);
     899          } else {
     900            // look other way - change down should be positive
     901            value -= nonLinear->changeDownInCost(iSequence);
     902            if (value>tolerance) {
     903              // store square in list
     904              infeasible_->quickAdd(iSequence,value*value);
     905            }
     906          }
    711907        }
    712908      }
     
    9291125        double value=ADD_ONE;
    9301126        model_->factorization()->updateColumn(temp,alternateWeights_);
    931         int number = alternateWeights_->getNumElements();
    932         int j;
     1127        int number = alternateWeights_->getNumElements();       int j;
    9331128        for (j=0;j<number;j++) {
    9341129          int iRow=which[j];
  • trunk/ClpSimplexPrimal.cpp

    r156 r162  
    496496    sumPrimalInfeasibilities_ = 0.0;
    497497  }
    498   if (dualFeasible()||problemStatus_==-4||(type==3&&problemStatus_!=-5)) {
     498  // had ||(type==3&&problemStatus_!=-5) -- ??? why ????
     499  if (dualFeasible()||problemStatus_==-4) {
    499500   
    500501    if (nonLinearCost_->numberInfeasibilities()&&!alwaysOptimal) {
     
    10981099  if (sequenceIn_>=0) {
    10991100    valueIn_=solution_[sequenceIn_];
     1101    dualIn_=dj_[sequenceIn_];
     1102    if (nonLinearCost_->lookBothWays()) {
     1103      // double check
     1104      ClpSimplex::Status status = getStatus(sequenceIn_);
     1105     
     1106      switch(status) {
     1107      case ClpSimplex::atUpperBound:
     1108        if (dualIn_<0.0) {
     1109          // move to other side
     1110          printf("For %d U (%g, %g, %g) dj changed from %g",
     1111                 sequenceIn_,lower_[sequenceIn_],solution_[sequenceIn_],
     1112                 upper_[sequenceIn_],dualIn_);
     1113          dualIn_ -= nonLinearCost_->changeUpInCost(sequenceIn_);
     1114          printf(" to %g\n",dualIn_);
     1115          nonLinearCost_->setOne(sequenceIn_,upper_[sequenceIn_]+2.0*currentPrimalTolerance());
     1116          setStatus(sequenceIn_,ClpSimplex::atLowerBound);
     1117        }
     1118        break;
     1119      case ClpSimplex::atLowerBound:
     1120        if (dualIn_>0.0) {
     1121          // move to other side
     1122          printf("For %d L (%g, %g, %g) dj changed from %g",
     1123                 sequenceIn_,lower_[sequenceIn_],solution_[sequenceIn_],
     1124                 upper_[sequenceIn_],dualIn_);
     1125          dualIn_ -= nonLinearCost_->changeDownInCost(sequenceIn_);
     1126          printf(" to %g\n",dualIn_);
     1127          nonLinearCost_->setOne(sequenceIn_,lower_[sequenceIn_]-2.0*currentPrimalTolerance());
     1128          setStatus(sequenceIn_,ClpSimplex::atUpperBound);
     1129        }
     1130        break;
     1131      default:
     1132        break;
     1133      }
     1134    }
    11001135    lowerIn_=lower_[sequenceIn_];
    11011136    upperIn_=upper_[sequenceIn_];
    1102     dualIn_=dj_[sequenceIn_];
    11031137    if (dualIn_>0.0)
    11041138      directionIn_ = -1;
  • trunk/ClpSimplexPrimalQuadratic.cpp

    r156 r162  
    3030
    3131  int numberColumns = this->numberColumns();
     32  int numberRows = this->numberRows();
    3233  double * columnLower = this->columnLower();
    3334  double * columnUpper = this->columnUpper();
     
    9697  double lastObjective=1.0e31;
    9798  double * saveSolution = new double [numberColumns];
     99  double * savePi = new double [numberRows];
     100  unsigned char * saveStatus = new unsigned char[numberRows+numberColumns];
    98101  double targetDrop=1.0e31;
    99102  double objectiveOffset;
     
    106109  }
    107110  double goodMove=false;
     111  char * statusCheck = new char[numberColumns];
    108112  for (iPass=0;iPass<numberPasses;iPass++) {
    109113    // redo objective
     
    115119    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
    116120      iColumn=listNonLinearColumn[jNon];
     121      if (getColumnStatus(iColumn)==basic) {
     122        if (solution[iColumn]<columnLower[iColumn]+1.0e-8)
     123          statusCheck[iColumn]='l';
     124        else if (solution[iColumn]>columnUpper[iColumn]-1.0e-8)
     125          statusCheck[iColumn]='u';
     126        else
     127          statusCheck[iColumn]='B';
     128      } else {
     129        if (solution[iColumn]<columnLower[iColumn]+1.0e-8)
     130          statusCheck[iColumn]='L';
     131        else
     132          statusCheck[iColumn]='U';
     133      }
    117134      double valueI = solution[iColumn];
    118135      int j;
     
    185202      std::cout<<"True drop was "<<drop<<std::endl;
    186203      std::cout<<"largest delta is "<<maxDelta<<std::endl;
    187       if (maxDelta<deltaTolerance&&drop<1.0e-4) {
     204      if (maxDelta<deltaTolerance&&drop<1.0e-4&&goodMove) {
    188205        std::cout<<"Exiting"<<std::endl;
    189206        break;
    190207      }
    191208    } else {
    192       lastObjective += 1.0e-3; // to stop exit
     209      //lastObjective += 1.0e-3; // to stop exit
    193210    }
    194211    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     
    203220    if (goodMove) {
    204221      memcpy(saveSolution,solution,numberColumns*sizeof(double));
     222      memcpy(savePi,this->dualRowSolution(),numberRows*sizeof(double));
     223      memcpy(saveStatus,status_,numberRows+numberColumns);
    205224     
    206225      targetDrop=0.0;
    207226      if (iPass) {
    208227        // get reduced costs
    209         this->matrix()->transposeTimes(this->dualRowSolution(),
    210                                        this->dualColumnSolution());
     228        this->matrix()->transposeTimes(savePi,
     229                                     this->dualColumnSolution());
    211230        double * r = this->dualColumnSolution();
    212231        for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
    213232          iColumn=listNonLinearColumn[jNon];
    214233          double dj = objective[iColumn]-r[iColumn];
     234          r[iColumn]=dj;
    215235          if (dj<0.0)
    216236            targetDrop -= dj*(columnUpper[iColumn]-solution[iColumn]);
     
    223243               <<std::endl;
    224244      lastObjective = objValue;
     245      if (targetDrop<1.0e-5&&goodMove&&iPass) {
     246        printf("Exiting on target drop %g\n",targetDrop);
     247        break;
     248      }
     249      {
     250        double * r = this->dualColumnSolution();
     251        for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     252          iColumn=listNonLinearColumn[jNon];
     253          printf("Trust %d %g - solution %d %g obj %g dj %g state %c\n",
     254                 jNon,trust[jNon],iColumn,solution[iColumn],objective[iColumn],
     255                 r[iColumn],statusCheck[iColumn]);
     256        }
     257      }
     258      setLogLevel(63);
    225259      this->primal(1);
    226260    } else {
    227261      // bad pass - restore solution
     262      printf("Backtracking\n");
    228263      memcpy(solution,saveSolution,numberColumns*sizeof(double));
     264      memcpy(this->dualRowSolution(),savePi,numberRows*sizeof(double));
     265      memcpy(status_,saveStatus,numberRows+numberColumns);
    229266      iPass--;
    230267    }
     
    239276                             trueUpper[jNon]);
    240277  }
     278  delete [] statusCheck;
     279  delete [] savePi;
     280  delete [] saveStatus;
    241281  // redo objective
    242282  double offset=0.0;
     
    333373  int newNumberColumns = 3*numberColumns+ numberRows;
    334374  int numberElements = 2*matrix->getNumElements()
    335     +quadratic->getNumElements()
     375    +2*quadratic->getNumElements()
    336376    + 2*numberColumns;
    337377  double * elements2 = new double[numberElements];
  • trunk/Samples/Makefile.piece

    r155 r162  
    5959        LDFLAGS += -lhistory -lreadline -ltermcap
    6060endif
    61 LDFLAGS += -lefence
     61#LDFLAGS += -lefence
    6262###############################################################################
    6363
  • trunk/Samples/piece.cpp

    r159 r162  
    55   deletes every second column and solves the resulting problem */
    66
    7 void sortSegments(int, int*, double*, double*);
    87#include "ClpSimplex.hpp"
    98#include "ClpNonLinearCost.hpp"
     
    7473
    7574  int iColumn;
     75  for (iColumn=0;iColumn<numberColumns;iColumn++)
     76    printf("%g ",oldSolution[iColumn]);
     77  printf("\n");
     78   
    7679  numberColumns2=0;
    7780  numberElements=0;
    7881  start2[0]=0;
    7982  int segptr=0;
    80   int numseg=0;
    81 
    82   // treat first column separately
    83   iColumn=0;
     83
    8484  segstart[0]=0;
    85   columnLower2[numberColumns2] = columnLower1[iColumn];
    86   columnUpper2[numberColumns2] = columnUpper1[iColumn];
    87   objective2[numberColumns2] = objective1[iColumn];
    88   breakpt[segptr] = columnLower1[iColumn];
    89   slope[segptr++] = objective1[iColumn];
    90   double maxUpper = columnUpper1[iColumn];
    91   numseg=1;
    92   for (int j=start1[iColumn];j<start1[iColumn]+length1[iColumn];j++) {
    93     row2[numberElements]=row1[j];
    94     element2[numberElements++] = element1[j];
    95   }
    96   newSolution[0]=oldSolution[0];
    97   start2[++numberColumns2]=numberElements;
    9885
    9986  // Now check for duplicates
    100   for (iColumn=1;iColumn<numberColumns;iColumn++) {
    101     maxUpper = max(maxUpper, columnUpper1[iColumn]);
    102     // test if column identical to previous column
     87  for (iColumn=0;iColumn<numberColumns;iColumn++) {
     88    // test if column identical to next column
    10389    bool ifcopy=1;
    104     int  joff = length1[iColumn-1];
    105     for (int j=start1[iColumn];j<start1[iColumn]+length1[iColumn];j++) {
    106       if (row1[j] != row1[j-joff]){
    107         ifcopy=0;
    108         break;
    109       }
    110       if (element1[j] != element1[j-joff]){
    111         ifcopy=0;
    112         break;
    113       }
     90    if (iColumn<numberColumns-1) {
     91      int  joff = length1[iColumn];
     92      for (int j=start1[iColumn];j<start1[iColumn]+length1[iColumn];j++) {
     93        if (row1[j] != row1[j+joff]){
     94          ifcopy=0;
     95          break;
     96        }
     97        if (element1[j] != element1[j+joff]){
     98          ifcopy=0;
     99          break;
     100        }
     101      }
     102    } else {
     103      ifcopy=0;
    114104    }
     105    //if (iColumn>47||iColumn<45)
     106    //ifcopy=0;
    115107    if (ifcopy) {
     108      double lo1 = columnLower1[iColumn];
     109      double up1 = columnUpper1[iColumn];
     110      double obj1 = objective1[iColumn];
     111      double sol1 = oldSolution[iColumn];
     112      double lo2 = columnLower1[iColumn+1];
     113      double up2 = columnUpper1[iColumn+1];
     114      double obj2 = objective1[iColumn+1];
     115      double sol2 = oldSolution[iColumn+1];
     116      if(fabs(up1-lo2)>1.0e-8) {
     117        // try other way
     118        double temp;
     119        temp = lo1;
     120        lo1=lo2;
     121        lo2=temp;
     122        temp = up1;
     123        up1=up2;
     124        up2=temp;
     125        temp = obj1;
     126        obj1=obj2;
     127        obj2=temp;
     128        temp = sol1;
     129        sol1=sol2;
     130        sol2=temp;
     131        assert(fabs(up1-lo2)<1.0e-8);
     132      }
    116133      // subtract out from rhs
    117       double fixed = columnLower1[iColumn];
    118       if(fabs(fixed-columnUpper1[iColumn-1])>1.0e-8) {
    119         // try other way
    120         fixed = columnUpper1[iColumn];
    121         assert(fabs(fixed-columnLower1[iColumn-1])<1.0e-8);
    122       }
     134      double fixed = up1;
    123135      // do offset
    124       objectiveOffset += fixed*objective1[iColumn];
     136      objectiveOffset += fixed*obj2;
    125137      for (int j=start1[iColumn];j<start1[iColumn]+length1[iColumn];j++) {
    126138        int iRow = row1[j];
     
    131143          rowUpper2[iRow] -= value*fixed;
    132144      }
    133       if (fabs(oldSolution[iColumn]-fixed)>1.0e-8)
    134         newSolution[numberColumns2-1]=oldSolution[iColumn];
     145      newSolution[numberColumns2]=fixed;
     146      if (fabs(sol1-fixed)>1.0e-8)
     147        newSolution[numberColumns2]=sol1;
     148      if (fabs(sol2-fixed)>1.0e-8)
     149        newSolution[numberColumns2]=sol2;
     150      columnLower2[numberColumns2] = lo1;
     151      columnUpper2[numberColumns2] = up2;
     152      objective2[numberColumns2] = 0.0;
     153      breakpt[segptr] = lo1;
     154      slope[segptr++] = obj1;
     155      breakpt[segptr] = lo2;
     156      slope[segptr++] = obj2;
     157      for (int j=start1[iColumn];j<start1[iColumn]+length1[iColumn];j++) {
     158        row2[numberElements]=row1[j];
     159        element2[numberElements++] = element1[j];
     160      }
     161      start2[++numberColumns2]=numberElements;
     162      breakpt[segptr] = up2;
     163      slope[segptr++] = COIN_DBL_MAX;
     164      segstart[numberColumns2] = segptr;
     165      iColumn++; // skip next column
     166    } else {
     167      // normal column
     168      columnLower2[numberColumns2] = columnLower1[iColumn];
     169      columnUpper2[numberColumns2] = columnUpper1[iColumn];
     170      objective2[numberColumns2] = objective1[iColumn];
    135171      breakpt[segptr] = columnLower1[iColumn];
    136172      slope[segptr++] = objective1[iColumn];
    137       numseg++;
    138       continue;
     173      for (int j=start1[iColumn];j<start1[iColumn]+length1[iColumn];j++) {
     174        row2[numberElements]=row1[j];
     175        element2[numberElements++] = element1[j];
     176      }
     177      newSolution[numberColumns2]=oldSolution[iColumn];
     178      start2[++numberColumns2]=numberElements;
     179      breakpt[segptr] = columnUpper1[iColumn];
     180      slope[segptr++] = COIN_DBL_MAX;
     181      segstart[numberColumns2] = segptr;
    139182    }
    140     breakpt[segptr] = maxUpper;
    141     slope[segptr++] = COIN_DBL_MAX;
    142     segstart[numberColumns2] = segptr;
    143     // new column found
    144     columnLower2[numberColumns2] = columnLower1[iColumn];
    145     columnUpper2[numberColumns2] = columnUpper1[iColumn];
    146     objective2[numberColumns2] = objective1[iColumn];
    147     breakpt[segptr] = columnLower1[iColumn];
    148     slope[segptr++] = objective1[iColumn];
    149     maxUpper = columnUpper1[iColumn];
    150     for (int j=start1[iColumn];j<start1[iColumn]+length1[iColumn];j++) {
    151       row2[numberElements]=row1[j];
    152       element2[numberElements++] = element1[j];
    153     }
    154     newSolution[numberColumns2]=oldSolution[iColumn];
    155     start2[++numberColumns2]=numberElements;
    156   }
    157   breakpt[segptr] = maxUpper;
    158   slope[segptr++] = COIN_DBL_MAX;
    159   segstart[numberColumns2] = segptr;
     183  }
    160184
    161185  // print new number of columns, elements
     
    164188  printf("Objective offset is %g\n",objectiveOffset);
    165189
    166   /*  for (int k=0; k<20; k++)
    167     printf("%d  %e %e\n",segstart[k],breakpt[k],slope[k]);
    168   */
    169   sortSegments(numberColumns2, segstart, breakpt, slope);
    170 
    171   printf("After sorting\n");
    172   /*  for (int k=0; k<20; k++)
    173     printf("%d  %e %e\n",segstart[k],breakpt[k],slope[k]);
    174   */
    175190
    176191  ClpSimplex  model;
     
    183198                    rowLower2,rowUpper2);
    184199  model.scaling(0);
     200  model.setDblParam(ClpObjOffset,-objectiveOffset);
    185201  // Create nonlinear objective
    186202  int returnCode= model.createPiecewiseLinearCosts(segstart, breakpt, slope);
     
    207223  //memcpy(model.columnUpper(),newSolution,numberColumns2*sizeof(double));
    208224  delete [] newSolution;
     225  //model.setLogLevel(63);
    209226
    210227  const double * solution = model.primalColumnSolution();
     228  double * saveSol = new double[numberColumns2];
     229  memcpy(saveSol,solution,numberColumns2*sizeof(double));
    211230  for (iColumn=0;iColumn<numberColumns2;iColumn++)
    212231    printf("%g ",solution[iColumn]);
     
    214233  // solve
    215234  model.primal(1);
    216   for (iColumn=0;iColumn<numberColumns2;iColumn++)
    217     printf("%g ",solution[iColumn]);
     235  for (iColumn=0;iColumn<numberColumns2;iColumn++) {
     236    if (fabs(solution[iColumn]-saveSol[iColumn])>1.0e-3)
     237      printf(" ** was %g ",saveSol[iColumn]);
     238    printf("%g ",solution[iColumn]);
     239  }
    218240  printf("\n");
    219241  model.primal(1);
    220   for (iColumn=0;iColumn<numberColumns2;iColumn++)
    221     printf("%g ",solution[iColumn]);
     242  for (iColumn=0;iColumn<numberColumns2;iColumn++) {
     243    if (fabs(solution[iColumn]-saveSol[iColumn])>1.0e-3)
     244      printf(" ** was %g ",saveSol[iColumn]);
     245    printf("%g ",solution[iColumn]);
     246  }
    222247  printf("\n");
    223248  model.primal();
    224   for (iColumn=0;iColumn<numberColumns2;iColumn++)
    225     printf("%g ",solution[iColumn]);
     249  for (iColumn=0;iColumn<numberColumns2;iColumn++) {
     250    if (fabs(solution[iColumn]-saveSol[iColumn])>1.0e-3)
     251      printf(" ** was %g ",saveSol[iColumn]);
     252    printf("%g ",solution[iColumn]);
     253  }
    226254  printf("\n");
    227255  model.allSlackBasis();
    228   for (iColumn=0;iColumn<numberColumns2;iColumn++)
    229     printf("%g ",solution[iColumn]);
    230   printf("\n");
     256  for (iColumn=0;iColumn<numberColumns2;iColumn++) {
     257    if (fabs(solution[iColumn]-saveSol[iColumn])>1.0e-3)
     258      printf(" ** was %g ",saveSol[iColumn]);
     259    printf("%g ",solution[iColumn]);
     260  }
     261  printf("\n");
     262  model.setLogLevel(63);
    231263  model.primal();
    232   for (iColumn=0;iColumn<numberColumns2;iColumn++)
    233     printf("%g ",solution[iColumn]);
     264  for (iColumn=0;iColumn<numberColumns2;iColumn++) {
     265    if (fabs(solution[iColumn]-saveSol[iColumn])>1.0e-3)
     266      printf(" ** was %g ",saveSol[iColumn]);
     267    printf("%g ",solution[iColumn]);
     268  }
    234269  printf("\n");
    235270  return 0;
    236271}   
    237 
    238 // Quick and dirty sort for 2 segments only
    239 void sortSegments(int n_, int* ptr_, double* breakpt_, double* slope_){
    240   for (int i=0; i<n_; i++){
    241     if ( (ptr_[i+1] - ptr_[i]) <= 2)
    242       continue;
    243     if ( (ptr_[i+1] - ptr_[i]) > 3){
    244       printf("More than 2 real segments not implemented %d\n", i);
    245       exit(2);
    246     }
    247     if (breakpt_[ptr_[i]] < breakpt_[ptr_[i]+1])
    248       continue;
    249     double bsave = breakpt_[ptr_[i]];
    250     double ssave = slope_[ptr_[i]];
    251     breakpt_[ptr_[i]] = breakpt_[ptr_[i]+1];
    252     slope_[ptr_[i]] = slope_[ptr_[i]+1];
    253     breakpt_[ptr_[i]+1] = bsave;
    254     slope_[ptr_[i]+1] = ssave;
    255   }
    256   return;
    257 }
  • trunk/include/ClpNonLinearCost.hpp

    r114 r162  
    112112      return cost_[iRange]-cost_[iRange+1];
    113113  }
     114  inline double changeUpInCost(int sequence) const
     115  {
     116    int iRange = whichRange_[sequence]+offset_[sequence];
     117    if (iRange+1!=start_[sequence+1]&&!infeasible(iRange+1))
     118      return cost_[iRange]-cost_[iRange+1];
     119    else
     120      return -1.0e100;
     121  }
     122  inline double changeDownInCost(int sequence) const
     123  {
     124    int iRange = whichRange_[sequence]+offset_[sequence];
     125    if (iRange!=start_[sequence]&&!infeasible(iRange-1))
     126      return cost_[iRange]-cost_[iRange-1];
     127    else
     128      return 1.0e100;
     129  }
    114130  /// Returns current lower bound
    115131  inline double lower(int sequence) const
     
    140156  inline void setChangeInCost(double value)
    141157  {changeCost_ = value;};
     158  /// See if may want to look both ways
     159  inline bool lookBothWays() const
     160  { return bothWays_;};
    142161  //@}
    143162  ///@name Private functions to deal with infeasible regions
     
    188207  /// If all non-linear costs convex
    189208  bool convex_;
     209  /// If we should look both ways for djs
     210  bool bothWays_;
    190211  //@}
    191212};
  • trunk/include/ClpSimplex.hpp

    r159 r162  
    598598  void setRowScale(double * scale) { rowScale_ = scale;};
    599599  void setColumnScale(double * scale) { columnScale_ = scale;};
     600  /// Return pointer to details of costs
     601  inline ClpNonLinearCost * nonLinearCost() const
     602  { return nonLinearCost_;};
    600603  //@}
    601604  /**@name status methods */
Note: See TracChangeset for help on using the changeset viewer.