Changeset 1643 for trunk/Clp


Ignore:
Timestamp:
Dec 2, 2010 7:14:40 AM (9 years ago)
Author:
forrest
Message:

stab at parametrics and sparser matrix multiply

Location:
trunk/Clp/src
Files:
9 edited

Legend:

Unmodified
Added
Removed
  • trunk/Clp/src/CbcOrClpParam.cpp

    r1605 r1643  
    24702470values, 2 saves values, 3 with greater accuracy and 4 in IEEE."
    24712471     );
     2472#ifdef COIN_HAS_CLP
     2473     parameters[numberParameters++] =
     2474          CbcOrClpParam("para!metrics", "Import data from file and do parametrics",
     2475                        CLP_PARAM_ACTION_PARAMETRICS, 3);
     2476     parameters[numberParameters-1].setLonghelp
     2477     (
     2478          "This will read a file with parametric data from the given file name \
     2479and then do parametrics.  It will use the default\
     2480 directory given by 'directory'.  A name of '$' will use the previous value for the name.  This\
     2481 is initialized to '', i.e. it must be set.  This can not read from compressed files. \
     2482File is in modified csv format - a line ROWS will be followed by rows data \
     2483while a line COLUMNS will be followed by column data.  The last line \
     2484should be ENDATA. The ROWS line must exist and is in the format \
     2485ROWS, inital theta, final theta, interval theta, n where n is 0 to get \
     2486CLPI0062 message at interval or at each change of theta \
     2487and 1 to get CLPI0063 message at each iteration.  If interval theta is 0.0 \
     2488or >= final theta then no interval reporting.  n may be missed out when it is \
     2489taken as 0.  If there is Row data then \
     2490there is a headings line with allowed headings - name, number, \
     2491lower(rhs change), upper(rhs change), rhs(change).  Either the lower and upper \
     2492fields should be given or the rhs field. \
     2493The optional COLUMNS line is followed by a headings line with allowed \
     2494headings - name, number, objective(change), lower(change), upper(change). \
     2495 Exactly one of name and number must be given for either section and \
     2496missing ones have value 0.0."
     2497     );
     2498#endif
    24722499#ifdef COIN_HAS_CBC
    24732500     parameters[numberParameters++] =
  • trunk/Clp/src/CbcOrClpParam.hpp

    r1525 r1643  
    246246     CLP_PARAM_ACTION_STOREDFILE,
    247247     CLP_PARAM_ACTION_ENVIRONMENT,
     248     CLP_PARAM_ACTION_PARAMETRICS,
    248249
    249250     CBC_PARAM_ACTION_BAB = 351,
  • trunk/Clp/src/ClpMain.cpp

    r1642 r1643  
    12931293                                        time1 = time2;
    12941294                                   }
     1295                              } else {
     1296                                   std::cout << "** Current model not valid" << std::endl;
     1297                              }
     1298                              break;
     1299                         case CLP_PARAM_ACTION_PARAMETRICS:
     1300                              if (goodModels[iModel]) {
     1301                                   // get next field
     1302                                   field = CoinReadGetString(argc, argv);
     1303                                   if (field == "$") {
     1304                                        field = parameters[iParam].stringValue();
     1305                                   } else if (field == "EOL") {
     1306                                        parameters[iParam].printString();
     1307                                        break;
     1308                                   } else {
     1309                                        parameters[iParam].setStringValue(field);
     1310                                   }
     1311                                   std::string fileName;
     1312                                   //bool canOpen = false;
     1313                                   if (field[0] == '/' || field[0] == '\\') {
     1314                                        fileName = field;
     1315                                   } else if (field[0] == '~') {
     1316                                        char * environVar = getenv("HOME");
     1317                                        if (environVar) {
     1318                                             std::string home(environVar);
     1319                                             field = field.erase(0, 1);
     1320                                             fileName = home + field;
     1321                                        } else {
     1322                                             fileName = field;
     1323                                        }
     1324                                   } else {
     1325                                        fileName = directory + field;
     1326                                   }
     1327                                   ClpSimplex * model2 = models + iModel;
     1328                                   static_cast<ClpSimplexOther *> (model2)->parametrics(fileName.c_str());
     1329                                   time2 = CoinCpuTime();
     1330                                   totalTime += time2 - time1;
     1331                                   time1 = time2;
    12951332                              } else {
    12961333                                   std::cout << "** Current model not valid" << std::endl;
  • trunk/Clp/src/ClpMessage.cpp

    r1525 r1643  
    100100     {CLP_BAD_STRING_VALUES, 3005, 1, "%d string elements had no values associated with them"},
    101101     {CLP_CRUNCH_STATS, 61, 2, "Crunch %d (%d) rows, %d (%d) columns and %d (%d) elements"},
     102     {CLP_PARAMETRICS_STATS, 62, 1, "Theta %g - objective %g"},
     103     {CLP_PARAMETRICS_STATS2, 63, 2, "Theta %g - objective %g, %s in, %s out"},
    102104     {CLP_GENERAL, 1000, 1, "%s"},
    103105     {CLP_DUMMY_END, 999999, 0, ""}
  • trunk/Clp/src/ClpMessage.hpp

    r1525 r1643  
    100100     CLP_BAD_STRING_VALUES,
    101101     CLP_CRUNCH_STATS,
     102     CLP_PARAMETRICS_STATS,
     103     CLP_PARAMETRICS_STATS2,
    102104     CLP_GENERAL,
    103105     CLP_DUMMY_END
  • trunk/Clp/src/ClpPackedMatrix.cpp

    r1525 r1643  
    11491149               int * index = columnArray->getIndices();
    11501150               double * array = columnArray->denseVector();
     1151#ifdef COIN_SPARSE_MATRIX
    11511152               assert (!y->getNumElements());
    11521153               // and set up mark as char array
    1153                //char * marked = (char *) (index+columnArray->capacity());
    1154                //int * lookup = y->getIndices();
    1155                int numberColumns = matrix_->getNumCols();
     1154               char * marked = reinterpret_cast<char *> (index+columnArray->capacity());
     1155               int * lookup = y->getIndices();
    11561156               //int numberRows = matrix_->getNumRows();
    11571157#ifndef NDEBUG
     
    11591159               //assert(!marked[i]);
    11601160#endif
    1161                //if (numberInRowArray>0)
     1161               numberNonZero=gutsOfTransposeTimesByRowGE3(rowArray,index,array,
     1162                                           lookup,marked,zeroTolerance,scalar);
     1163#else
     1164               int numberColumns = matrix_->getNumCols();
    11621165               numberNonZero = gutsOfTransposeTimesByRowGEK(rowArray, index, array,
    11631166                               numberColumns, zeroTolerance, scalar);
    1164                //else
    1165                //numberNonZero=gutsOfTransposeTimesByRowGE3(rowArray,index,array,
    1166                //                                  lookup,marked,zeroTolerance,scalar);
     1167#endif
    11671168               columnArray->setNumElements(numberNonZero);
    11681169          } else {
  • trunk/Clp/src/ClpPresolve.cpp

    r1641 r1643  
    9898{
    9999     // Check matrix
     100     int checkType = ((si.specialOptions() & 128) != 0) ? 14 : 15;
    100101     if (!si.clpMatrix()->allElementsInRange(&si, si.getSmallElementValue(),
    101                                              1.0e20))
     102                                             1.0e20,checkType))
    102103          return NULL;
    103104     else
  • trunk/Clp/src/ClpSimplexOther.cpp

    r1525 r1643  
    2424#include <string>
    2525#include <stdio.h>
    26 #include <iostream>
     26#include <iostream> 
    2727/* Dual ranging.
    2828   This computes increase/decrease in cost for each given variable and corresponding
     
    19541954                    returnCode = -1;
    19551955               }
    1956                endingTheta = CoinMin(endingTheta, maxTheta);
     1956               if (maxTheta < endingTheta) {
     1957                 char line[100];
     1958                 sprintf(line,"Crossover considerations reduce ending  theta from %g to %g\n",
     1959                         endingTheta,maxTheta);
     1960                 handler_->message(CLP_GENERAL,messages_)
     1961                   << line << CoinMessageEol;
     1962                 endingTheta = maxTheta;
     1963               }
    19571964               if (endingTheta < startingTheta) {
    19581965                    // bad initial
     
    19721979               assert (!problemStatus_);
    19731980               // Now do parametrics
    1974                printf("at starting theta of %g objective value is %g\n", startingTheta,
    1975                       objectiveValue());
     1981               handler_->message(CLP_PARAMETRICS_STATS, messages_)
     1982                 << startingTheta << objectiveValue() << CoinMessageEol;
    19761983               while (!returnCode) {
    1977                     assert (reportIncrement);
     1984                    //assert (reportIncrement);
    19781985                    returnCode = parametricsLoop(startingTheta, endingTheta, reportIncrement,
    19791986                                                 chgLower, chgUpper, chgObjective, data,
     
    19881995                         //cost_[i] += change*chgObjective[i];
    19891996                         //}
    1990                          printf("at theta of %g objective value is %g\n", startingTheta,
    1991                                 objectiveValue());
     1997                         handler_->message(CLP_PARAMETRICS_STATS, messages_)
     1998                           << startingTheta << objectiveValue() << CoinMessageEol;
    19921999                         if (startingTheta >= endingTheta)
    19932000                              break;
     
    19952002                         // trouble - do external solve
    19962003                         needToDoSomething = true;
     2004                    } else if (problemStatus_==1) {
     2005                      // can't move any further
     2006                      if (!canTryQuick) {
     2007                        handler_->message(CLP_PARAMETRICS_STATS, messages_)
     2008                          << endingTheta << objectiveValue() << CoinMessageEol;
     2009                        problemStatus_=0;
     2010                      }
    19972011                    } else {
    19982012                         abort();
     
    20432057                    copyModel.dual();
    20442058                    if (copyModel.problemStatus()) {
    2045                          printf("Can not get to theta of %g\n", startingTheta);
     2059                      char line[100];
     2060                      sprintf(line,"Can not get to theta of %g\n", startingTheta);
     2061                      handler_->message(CLP_GENERAL,messages_)
     2062                        << line << CoinMessageEol;
    20462063                         canTryQuick = false; // do slowly to get exact amount
    20472064                         // back to last known good
     
    20652082     }
    20662083     perturbation_ = savePerturbation;
     2084     char line[100];
     2085     sprintf(line,"Ending theta %g\n", endingTheta);
     2086     handler_->message(CLP_GENERAL,messages_)
     2087       << line << CoinMessageEol;
    20672088     return problemStatus_;
     2089}
     2090/* Version of parametrics which reads from file
     2091   See CbcClpParam.cpp for details of format
     2092   Returns -2 if unable to open file */
     2093int
     2094ClpSimplexOther::parametrics(const char * dataFile)
     2095{
     2096  int returnCode=-2;
     2097  FILE *fp = fopen(dataFile, "r");
     2098  char line[200];
     2099  if (!fp) {
     2100    handler_->message(CLP_UNABLE_OPEN, messages_)
     2101      << dataFile << CoinMessageEol;
     2102    return -2;
     2103  }
     2104
     2105  if (!fgets(line, 200, fp)) {
     2106    sprintf(line,"Empty parametrics file %s?",dataFile);
     2107    handler_->message(CLP_GENERAL,messages_)
     2108      << line << CoinMessageEol;
     2109    fclose(fp);
     2110    return -2;
     2111  }
     2112  char * pos = line;
     2113  char * put = line;
     2114  while (*pos >= ' ' && *pos != '\n') {
     2115    if (*pos != ' ' && *pos != '\t') {
     2116      *put = static_cast<char>(tolower(*pos));
     2117      put++;
     2118    }
     2119    pos++;
     2120  }
     2121  *put = '\0';
     2122  pos = line;
     2123  double startTheta=0.0;
     2124  double endTheta=0.0;
     2125  double intervalTheta=COIN_DBL_MAX;
     2126  int detail=0;
     2127  bool good = true;
     2128  while (good) {
     2129    good=false;
     2130    // check ROWS
     2131    char * comma = strchr(pos, ',');
     2132    if (!comma)
     2133      break;
     2134    *comma = '\0';
     2135    if (strcmp(pos,"rows"))
     2136      break;
     2137    *comma = ',';
     2138    pos = comma+1;
     2139    // check lower theta
     2140    comma = strchr(pos, ',');
     2141    if (!comma)
     2142      break;
     2143    *comma = '\0';
     2144    startTheta = atof(pos);
     2145    *comma = ',';
     2146    pos = comma+1;
     2147    // check upper theta
     2148    comma = strchr(pos, ',');
     2149    good=true;
     2150    if (comma)
     2151      *comma = '\0';
     2152    endTheta = atof(pos);
     2153    if (comma) {
     2154      *comma = ',';
     2155      pos = comma+1;
     2156      comma = strchr(pos, ',');
     2157      if (comma)
     2158        *comma = '\0';
     2159      intervalTheta = atof(pos);
     2160      if (comma) {
     2161        *comma = ',';
     2162        pos = comma+1;
     2163        comma = strchr(pos, ',');
     2164        if (comma)
     2165          *comma = '\0';
     2166        detail = atoi(pos);
     2167        if (comma)
     2168        *comma = ',';
     2169      }
     2170    }
     2171    break;
     2172  }
     2173  if (good) {
     2174    if (startTheta<0.0||
     2175        startTheta>endTheta||
     2176        intervalTheta<0.0)
     2177      good=false;
     2178    if (detail<0||detail>1)
     2179      good=false;
     2180  }
     2181  if (intervalTheta>=endTheta)
     2182    intervalTheta=0.0;
     2183  if (!good) {
     2184    sprintf(line,"Odd first line %s on file %s?",line,dataFile);
     2185    handler_->message(CLP_GENERAL,messages_)
     2186      << line << CoinMessageEol;
     2187    fclose(fp);
     2188    return -2;
     2189  }
     2190  if (!fgets(line, 200, fp)) {
     2191    sprintf(line,"Not enough records on parametrics file %s?",dataFile);
     2192    handler_->message(CLP_GENERAL,messages_)
     2193      << line << CoinMessageEol;
     2194    fclose(fp);
     2195    return -2;
     2196  }
     2197  double * lowerRowMove = NULL;
     2198  double * upperRowMove = NULL;
     2199  double * lowerColumnMove = NULL;
     2200  double * upperColumnMove = NULL;
     2201  double * objectiveMove = NULL;
     2202  char saveLine[200];
     2203  saveLine[0]='\0';
     2204  std::string headingsRow[] = {"name", "number", "lower", "upper", "rhs"};
     2205  int gotRow[] = { -1, -1, -1, -1, -1};
     2206  int orderRow[5];
     2207  assert(sizeof(gotRow) == sizeof(orderRow));
     2208  int nAcross = 0;
     2209  pos = line;
     2210  put = line;
     2211  while (*pos >= ' ' && *pos != '\n') {
     2212    if (*pos != ' ' && *pos != '\t') {
     2213      *put = static_cast<char>(tolower(*pos));
     2214      put++;
     2215    }
     2216    pos++;
     2217  }
     2218  *put = '\0';
     2219  pos = line;
     2220  int i;
     2221  good = true;
     2222  if (strncmp(line,"column",6)) {
     2223    while (pos) {
     2224      char * comma = strchr(pos, ',');
     2225      if (comma)
     2226        *comma = '\0';
     2227      for (i = 0; i < static_cast<int> (sizeof(gotRow) / sizeof(int)); i++) {
     2228        if (headingsRow[i] == pos) {
     2229          if (gotRow[i] < 0) {
     2230            orderRow[nAcross] = i;
     2231            gotRow[i] = nAcross++;
     2232          } else {
     2233            // duplicate
     2234            good = false;
     2235          }
     2236          break;
     2237        }
     2238      }
     2239      if (i == static_cast<int> (sizeof(gotRow) / sizeof(int)))
     2240        good = false;
     2241      if (comma) {
     2242        *comma = ',';
     2243        pos = comma + 1;
     2244      } else {
     2245        break;
     2246      }
     2247    }
     2248    if (gotRow[0] < 0 && gotRow[1] < 0)
     2249      good = false;
     2250    if (gotRow[0] >= 0 && gotRow[1] >= 0)
     2251      good = false;
     2252    if (gotRow[0] >= 0 && !lengthNames())
     2253      good = false;
     2254    if (gotRow[4]<0) {
     2255      if (gotRow[2] < 0 && gotRow[3] >= 0)
     2256        good = false;
     2257      else if (gotRow[3] < 0 && gotRow[2] >= 0)
     2258        good = false;
     2259    } else if (gotRow[2]>=0||gotRow[3]>=0) {
     2260      good = false;
     2261    }
     2262    if (good) {
     2263      char ** rowNames = new char * [numberRows_];
     2264      int iRow;
     2265      for (iRow = 0; iRow < numberRows_; iRow++) {
     2266        rowNames[iRow] =
     2267          CoinStrdup(rowName(iRow).c_str());
     2268      }
     2269      lowerRowMove = new double [numberRows_];
     2270      memset(lowerRowMove,0,numberRows_*sizeof(double));
     2271      upperRowMove = new double [numberRows_];
     2272      memset(upperRowMove,0,numberRows_*sizeof(double));
     2273      int nLine = 0;
     2274      int nBadLine = 0;
     2275      int nBadName = 0;
     2276      bool goodLine=false;
     2277      while (fgets(line, 200, fp)) {
     2278        goodLine=true;
     2279        if (!strncmp(line, "ENDATA", 6)||
     2280            !strncmp(line, "COLUMN",6))
     2281          break;
     2282        goodLine=false;
     2283        nLine++;
     2284        iRow = -1;
     2285        double upper = 0.0;
     2286        double lower = 0.0;
     2287        char * pos = line;
     2288        char * put = line;
     2289        while (*pos >= ' ' && *pos != '\n') {
     2290          if (*pos != ' ' && *pos != '\t') {
     2291            *put = *pos;
     2292            put++;
     2293          }
     2294          pos++;
     2295        }
     2296        *put = '\0';
     2297        pos = line;
     2298        for (int i = 0; i < nAcross; i++) {
     2299          char * comma = strchr(pos, ',');
     2300          if (comma) {
     2301            *comma = '\0';
     2302          } else if (i < nAcross - 1) {
     2303            nBadLine++;
     2304            break;
     2305          }
     2306          switch (orderRow[i]) {
     2307            // name
     2308          case 0:
     2309            // For large problems this could be slow
     2310            for (iRow = 0; iRow < numberRows_; iRow++) {
     2311              if (!strcmp(rowNames[iRow], pos))
     2312                break;
     2313            }
     2314            if (iRow == numberRows_)
     2315              iRow = -1;
     2316            break;
     2317            // number
     2318          case 1:
     2319            iRow = atoi(pos);
     2320            if (iRow < 0 || iRow >= numberRows_)
     2321              iRow = -1;
     2322            break;
     2323            // lower
     2324          case 2:
     2325            upper = atof(pos);
     2326            break;
     2327            // upper
     2328          case 3:
     2329            lower = atof(pos);
     2330            break;
     2331            // rhs
     2332          case 4:
     2333            lower = atof(pos);
     2334            upper = lower;
     2335            break;
     2336          }
     2337          if (comma) {
     2338            *comma = ',';
     2339            pos = comma + 1;
     2340          }
     2341        }
     2342        if (iRow >= 0) {
     2343          if (rowLower_[iRow]>-1.0e20)
     2344            lowerRowMove[iRow] = lower;
     2345          else
     2346            lowerRowMove[iRow]=0.0;
     2347          if (rowUpper_[iRow]<1.0e20)
     2348            upperRowMove[iRow] = upper;
     2349          else
     2350            upperRowMove[iRow] = lower;
     2351        } else {
     2352          nBadName++;
     2353          if(saveLine[0]=='\0')
     2354            strcpy(saveLine,line);
     2355        }
     2356      }
     2357      sprintf(line,"%d Row fields and %d records", nAcross, nLine);
     2358      handler_->message(CLP_GENERAL,messages_)
     2359        << line << CoinMessageEol;
     2360      if (nBadName) {
     2361        sprintf(line," ** %d records did not match on name/sequence, first bad %s", nBadName,saveLine);
     2362        handler_->message(CLP_GENERAL,messages_)
     2363          << line << CoinMessageEol;
     2364        returnCode=-1;
     2365        good=false;
     2366      }
     2367      for (iRow = 0; iRow < numberRows_; iRow++) {
     2368        free(rowNames[iRow]);
     2369      }
     2370      delete [] rowNames;
     2371    } else {
     2372      sprintf(line,"Duplicate or unknown keyword - or name/number fields wrong");
     2373      handler_->message(CLP_GENERAL,messages_)
     2374        << line << CoinMessageEol;
     2375      returnCode=-1;
     2376      good=false;
     2377    }
     2378  }
     2379  if (good&&(!strncmp(line, "COLUMN",6)||!strncmp(line, "column",6))) {
     2380    if (!fgets(line, 200, fp)) {
     2381      sprintf(line,"Not enough records on parametrics file %s after COLUMNS?",dataFile);
     2382      handler_->message(CLP_GENERAL,messages_)
     2383        << line << CoinMessageEol;
     2384      fclose(fp);
     2385      return -2;
     2386    }
     2387    std::string headingsColumn[] = {"name", "number", "lower", "upper", "objective"};
     2388    saveLine[0]='\0';
     2389    int gotColumn[] = { -1, -1, -1, -1, -1};
     2390    int orderColumn[5];
     2391    assert(sizeof(gotColumn) == sizeof(orderColumn));
     2392    nAcross = 0;
     2393    pos = line;
     2394    put = line;
     2395    while (*pos >= ' ' && *pos != '\n') {
     2396      if (*pos != ' ' && *pos != '\t') {
     2397        *put = static_cast<char>(tolower(*pos));
     2398        put++;
     2399      }
     2400      pos++;
     2401    }
     2402    *put = '\0';
     2403    pos = line;
     2404    int i;
     2405    if (strncmp(line,"endata",6)&&good) {
     2406      while (pos) {
     2407        char * comma = strchr(pos, ',');
     2408        if (comma)
     2409          *comma = '\0';
     2410        for (i = 0; i < static_cast<int> (sizeof(gotColumn) / sizeof(int)); i++) {
     2411          if (headingsColumn[i] == pos) {
     2412            if (gotColumn[i] < 0) {
     2413              orderColumn[nAcross] = i;
     2414              gotColumn[i] = nAcross++;
     2415            } else {
     2416              // duplicate
     2417              good = false;
     2418            }
     2419            break;
     2420          }
     2421        }
     2422        if (i == static_cast<int> (sizeof(gotColumn) / sizeof(int)))
     2423          good = false;
     2424        if (comma) {
     2425          *comma = ',';
     2426          pos = comma + 1;
     2427        } else {
     2428          break;
     2429        }
     2430      }
     2431      if (gotColumn[0] < 0 && gotColumn[1] < 0)
     2432        good = false;
     2433      if (gotColumn[0] >= 0 && gotColumn[1] >= 0)
     2434        good = false;
     2435      if (gotColumn[0] >= 0 && !lengthNames())
     2436        good = false;
     2437      if (good) {
     2438        char ** columnNames = new char * [numberColumns_];
     2439        int iColumn;
     2440        for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     2441          columnNames[iColumn] =
     2442            CoinStrdup(columnName(iColumn).c_str());
     2443        }
     2444        lowerColumnMove = reinterpret_cast<double *> (malloc(numberColumns_ * sizeof(double)));
     2445        memset(lowerColumnMove,0,numberColumns_*sizeof(double));
     2446        upperColumnMove = reinterpret_cast<double *> (malloc(numberColumns_ * sizeof(double)));
     2447        memset(upperColumnMove,0,numberColumns_*sizeof(double));
     2448        objectiveMove = reinterpret_cast<double *> (malloc(numberColumns_ * sizeof(double)));
     2449        memset(objectiveMove,0,numberColumns_*sizeof(double));
     2450        int nLine = 0;
     2451        int nBadLine = 0;
     2452        int nBadName = 0;
     2453        bool goodLine=false;
     2454        while (fgets(line, 200, fp)) {
     2455          goodLine=true;
     2456          if (!strncmp(line, "ENDATA", 6))
     2457            break;
     2458          goodLine=false;
     2459          nLine++;
     2460          iColumn = -1;
     2461          double upper = 0.0;
     2462          double lower = 0.0;
     2463          double obj =0.0;
     2464          char * pos = line;
     2465          char * put = line;
     2466          while (*pos >= ' ' && *pos != '\n') {
     2467            if (*pos != ' ' && *pos != '\t') {
     2468              *put = *pos;
     2469              put++;
     2470            }
     2471            pos++;
     2472          }
     2473          *put = '\0';
     2474          pos = line;
     2475          for (int i = 0; i < nAcross; i++) {
     2476            char * comma = strchr(pos, ',');
     2477            if (comma) {
     2478              *comma = '\0';
     2479            } else if (i < nAcross - 1) {
     2480              nBadLine++;
     2481              break;
     2482            }
     2483            switch (orderColumn[i]) {
     2484              // name
     2485            case 0:
     2486              // For large problems this could be slow
     2487              for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     2488                if (!strcmp(columnNames[iColumn], pos))
     2489                  break;
     2490              }
     2491              if (iColumn == numberColumns_)
     2492                iColumn = -1;
     2493              break;
     2494              // number
     2495            case 1:
     2496              iColumn = atoi(pos);
     2497              if (iColumn < 0 || iColumn >= numberColumns_)
     2498                iColumn = -1;
     2499              break;
     2500              // lower
     2501            case 2:
     2502              upper = atof(pos);
     2503              break;
     2504              // upper
     2505            case 3:
     2506              lower = atof(pos);
     2507              break;
     2508              // objective
     2509            case 4:
     2510              obj = atof(pos);
     2511              upper = lower;
     2512              break;
     2513            }
     2514            if (comma) {
     2515              *comma = ',';
     2516              pos = comma + 1;
     2517            }
     2518          }
     2519          if (iColumn >= 0) {
     2520            if (columnLower_[iColumn]>-1.0e20)
     2521              lowerColumnMove[iColumn] = lower;
     2522            else
     2523              lowerColumnMove[iColumn]=0.0;
     2524            if (columnUpper_[iColumn]<1.0e20)
     2525              upperColumnMove[iColumn] = upper;
     2526            else
     2527              upperColumnMove[iColumn] = lower;
     2528            objectiveMove[iColumn] = obj;
     2529          } else {
     2530            nBadName++;
     2531            if(saveLine[0]=='\0')
     2532              strcpy(saveLine,line);
     2533          }
     2534        }
     2535        sprintf(line,"%d Column fields and %d records", nAcross, nLine);
     2536        handler_->message(CLP_GENERAL,messages_)
     2537          << line << CoinMessageEol;
     2538        if (nBadName) {
     2539          sprintf(line," ** %d records did not match on name/sequence, first bad %s", nBadName,saveLine);
     2540          handler_->message(CLP_GENERAL,messages_)
     2541            << line << CoinMessageEol;
     2542          returnCode=-1;
     2543          good=false;
     2544        }
     2545        for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     2546          free(columnNames[iColumn]);
     2547        }
     2548        delete [] columnNames;
     2549      } else {
     2550        sprintf(line,"Duplicate or unknown keyword - or name/number fields wrong");
     2551        handler_->message(CLP_GENERAL,messages_)
     2552          << line << CoinMessageEol;
     2553        returnCode=-1;
     2554        good=false;
     2555      }
     2556    }
     2557  }
     2558  returnCode=-1;
     2559  if (good) {
     2560    // clean arrays
     2561    if (lowerRowMove) {
     2562      bool empty=true;
     2563      for (int i=0;i<numberRows_;i++) {
     2564        if (lowerRowMove[i]) {
     2565          empty=false;
     2566        break;
     2567        }
     2568      }
     2569      if (empty) {
     2570        delete [] lowerRowMove;
     2571        lowerRowMove=NULL;
     2572      }
     2573    }
     2574    if (upperRowMove) {
     2575      bool empty=true;
     2576      for (int i=0;i<numberRows_;i++) {
     2577        if (upperRowMove[i]) {
     2578          empty=false;
     2579        break;
     2580        }
     2581      }
     2582      if (empty) {
     2583        delete [] upperRowMove;
     2584        upperRowMove=NULL;
     2585      }
     2586    }
     2587    if (lowerColumnMove) {
     2588      bool empty=true;
     2589      for (int i=0;i<numberColumns_;i++) {
     2590        if (lowerColumnMove[i]) {
     2591          empty=false;
     2592        break;
     2593        }
     2594      }
     2595      if (empty) {
     2596        delete [] lowerColumnMove;
     2597        lowerColumnMove=NULL;
     2598      }
     2599    }
     2600    if (upperColumnMove) {
     2601      bool empty=true;
     2602      for (int i=0;i<numberColumns_;i++) {
     2603        if (upperColumnMove[i]) {
     2604          empty=false;
     2605        break;
     2606        }
     2607      }
     2608      if (empty) {
     2609        delete [] upperColumnMove;
     2610        upperColumnMove=NULL;
     2611      }
     2612    }
     2613    if (objectiveMove) {
     2614      bool empty=true;
     2615      for (int i=0;i<numberColumns_;i++) {
     2616        if (objectiveMove[i]) {
     2617          empty=false;
     2618        break;
     2619        }
     2620      }
     2621      if (empty) {
     2622        delete [] objectiveMove;
     2623        objectiveMove=NULL;
     2624      }
     2625    }
     2626    int saveScaling = scalingFlag_;
     2627    scalingFlag_ = 0;
     2628    int saveLogLevel = handler_->logLevel();
     2629    if (detail>0&&!intervalTheta)
     2630      handler_->setLogLevel(3);
     2631    else
     2632      handler_->setLogLevel(1);
     2633    returnCode = parametrics(startTheta,endTheta,intervalTheta,
     2634                             lowerColumnMove,upperColumnMove,
     2635                             lowerRowMove,upperRowMove,
     2636                             objectiveMove);
     2637    scalingFlag_ = saveScaling;
     2638    handler_->setLogLevel(saveLogLevel);
     2639  }
     2640  delete [] lowerRowMove;
     2641  delete [] upperRowMove;
     2642  delete [] lowerColumnMove;
     2643  delete [] upperColumnMove;
     2644  delete [] objectiveMove;
     2645  fclose(fp);
     2646  return returnCode;
    20682647}
    20692648int
     
    21452724
    21462725          // exit if victory declared
    2147           if (problemStatus_ >= 0)
     2726          if (problemStatus_ >= 0 &&
     2727              (canTryQuick || startingTheta>=endingTheta-1.0e-7) )
    21482728               break;
    21492729
     
    21632743          }
    21642744          // Do iterations
     2745          problemStatus_=-1;
    21652746          if (canTryQuick) {
    21662747               double * saveDuals = NULL;
     
    21702751                              changeLower, changeUpper,
    21712752                              changeObjective);
     2753               startingTheta = endingTheta;
    21722754          }
    21732755     }
     
    23032885   +1 looks infeasible
    23042886   +3 max iterations
     2887   +4 accuracy problems
    23052888*/
    23062889int
     
    23272910     int returnCode = -1;
    23282911     double saveSumDual = sumDualInfeasibilities_; // so we know to be careful
     2912     double lastTheta = startingTheta;
    23292913     double useTheta = startingTheta;
    2330      double * primalChange = new double[numberRows_];
    2331      double * dualChange = new double[numberColumns_];
    23322914     int numberTotal = numberColumns_ + numberRows_;
     2915     double * primalChange = new double[numberTotal];
     2916     double * dualChange = new double[numberTotal];
    23332917     int iSequence;
    23342918     // See if bounds
     
    23492933     assert (type);
    23502934     while (problemStatus_ == -1) {
    2351           double increaseTheta = CoinMin(endingTheta - useTheta, 1.0e50);
     2935          double increaseTheta = CoinMin(endingTheta - lastTheta, 1.0e50);
    23522936
    23532937          // Get theta for bounds - we know can't crossover
    23542938          int pivotType = nextTheta(type, increaseTheta, primalChange, dualChange,
    23552939                                    changeLower, changeUpper, changeObjective);
    2356           if (pivotType)
    2357                abort();
     2940          useTheta += theta_;
     2941          double change = useTheta - lastTheta;
     2942          for (int i = 0; i < numberTotal; i++) {
     2943            lower_[i] += change * changeLower[i];
     2944            upper_[i] += change * changeUpper[i];
     2945            switch(getStatus(i)) {
     2946             
     2947            case basic:
     2948            case isFree:
     2949            case superBasic:
     2950              break;
     2951            case isFixed:
     2952            case atUpperBound:
     2953              solution_[i] = upper_[i];
     2954              break;
     2955            case atLowerBound:
     2956              solution_[i] = lower_[i];
     2957              break;
     2958            }
     2959            cost_[i] += change * changeObjective[i];
     2960            assert (solution_[i]>lower_[i]-1.0e-5&&
     2961                    solution_[i]<upper_[i]+1.0e-5);
     2962          }
     2963          sequenceIn_=-1;
     2964          if (pivotType) {
     2965              problemStatus_ = -2;
     2966              endingTheta = useTheta;
     2967              return 4;
     2968          }
    23582969          // choose row to go out
    2359           // dualRow will go to virtual row pivot choice algorithm
    2360           reinterpret_cast<ClpSimplexDual *> ( this)->dualRow(-1);
     2970          //reinterpret_cast<ClpSimplexDual *> ( this)->dualRow(-1);
    23612971          if (pivotRow_ >= 0) {
    23622972               // we found a pivot row
     
    25783188                    //setStatus(sequenceIn_,basic);
    25793189                    dj_[sequenceIn_] = 0.0;
    2580                     double oldValue = valueIn_;
     3190                    //double oldValue = valueIn_;
    25813191                    if (directionIn_ == -1) {
    25823192                         // as if from upper bound
     
    25863196                         valueIn_ = lowerIn_ + dualOut_;
    25873197                    }
    2588                     objectiveChange += cost_[sequenceIn_] * (valueIn_ - oldValue);
     3198                    objectiveChange = 0.0;
     3199                    for (int i=0;i<numberTotal;i++)
     3200                      objectiveChange += solution_[i]*cost_[i];
     3201                    objectiveChange -= objectiveValue_;
    25893202                    // outgoing
    25903203                    // set dj to zero unless values pass
     
    25983211                    solution_[sequenceOut_] = valueOut_;
    25993212                    int whatNext = housekeeping(objectiveChange);
     3213                    {
     3214                      char in[200],out[200];
     3215                      int iSequence=sequenceIn_;
     3216                      if (iSequence<numberColumns_) {
     3217                        if (lengthNames_)
     3218                          strcpy(in,columnNames_[iSequence].c_str());
     3219                         else
     3220                          sprintf(in,"C%7.7d",iSequence);
     3221                      } else {
     3222                        iSequence -= numberColumns_;
     3223                        if (lengthNames_)
     3224                          strcpy(in,rowNames_[iSequence].c_str());
     3225                         else
     3226                          sprintf(in,"R%7.7d",iSequence);
     3227                      }
     3228                      iSequence=sequenceOut_;
     3229                      if (iSequence<numberColumns_) {
     3230                        if (lengthNames_)
     3231                          strcpy(out,columnNames_[iSequence].c_str());
     3232                         else
     3233                          sprintf(out,"C%7.7d",iSequence);
     3234                      } else {
     3235                        iSequence -= numberColumns_;
     3236                        if (lengthNames_)
     3237                          strcpy(out,rowNames_[iSequence].c_str());
     3238                         else
     3239                          sprintf(out,"R%7.7d",iSequence);
     3240                      }
     3241                      handler_->message(CLP_PARAMETRICS_STATS2, messages_)
     3242                        << useTheta << objectiveValue()
     3243                        << in << out << CoinMessageEol;
     3244                    }
     3245                    if (useTheta>lastTheta+1.0e-9) {
     3246                      handler_->message(CLP_PARAMETRICS_STATS, messages_)
     3247                        << useTheta << objectiveValue() << CoinMessageEol;
     3248                      lastTheta = useTheta;
     3249                    }
    26003250                    // and set bounds correctly
    26013251                    reinterpret_cast<ClpSimplexDual *> ( this)->originalBound(sequenceIn_);
     
    27983448     delete [] primalChange;
    27993449     delete [] dualChange;
     3450     endingTheta = lastTheta;
    28003451     return returnCode;
    28013452}
     
    28323483          // use array
    28333484          double * array = rowArray_[1]->denseVector();
     3485          // put slacks in
     3486          for (int i=0;i<numberRows_;i++)
     3487            array[i] = - primalChange[i+numberColumns_];
    28343488          times(1.0, primalChange, array);
    28353489          int * index = rowArray_[1]->getIndices();
    28363490          int number = 0;
     3491          pivotRow_ = -1;
    28373492          for (iRow = 0; iRow < numberRows_; iRow++) {
    28383493               double value = array[iRow];
    28393494               if (value) {
    2840                     array[iRow] = value;
    28413495                    index[number++] = iRow;
    28423496               }
     
    28453499          rowArray_[1]->setNumElements(number);
    28463500          factorization_->updateColumn(rowArray_[0], rowArray_[1]);
    2847           number = rowArray_[1]->getNumElements();
    2848           pivotRow_ = -1;
    2849           for (iRow = 0; iRow < number; iRow++) {
    2850                int iPivot = index[iRow];
     3501          //number = rowArray_[1]->getNumElements();
     3502          for (int iPivot = 0; iPivot < numberRows_; iPivot++) {
     3503            //int iPivot = index[iRow];
    28513504               iSequence = pivotVariable_[iPivot];
    28523505               // solution value will be sol - theta*alpha
     
    28613514               double hitsLower = COIN_DBL_MAX;
    28623515               thetaCoefficient = changeLower[iSequence] + alpha;
    2863                if (fabs(thetaCoefficient) > 1.0e-8)
    2864                     hitsLower = (currentSolution - currentLower) / thetaCoefficient;
    2865                if (hitsLower < 0.0) {
     3516               if (thetaCoefficient > 1.0e-8)
     3517                hitsLower = (currentSolution - currentLower) / thetaCoefficient;
     3518               //if (hitsLower < 0.0) {
    28663519                    // does not hit - but should we check further
    2867                     hitsLower = COIN_DBL_MAX;
    2868                }
     3520               //   hitsLower = COIN_DBL_MAX;
     3521               //}
    28693522               double hitsUpper = COIN_DBL_MAX;
    28703523               thetaCoefficient = changeUpper[iSequence] + alpha;
    2871                if (fabs(thetaCoefficient) > 1.0e-8)
     3524               if (thetaCoefficient < -1.0e-8)
    28723525                    hitsUpper = (currentSolution - currentUpper) / thetaCoefficient;
    2873                if (hitsUpper < 0.0) {
     3526               //if (hitsUpper < 0.0) {
    28743527                    // does not hit - but should we check further
    2875                     hitsUpper = COIN_DBL_MAX;
    2876                }
     3528               //   hitsUpper = COIN_DBL_MAX;
     3529               //}
    28773530               if (CoinMin(hitsLower, hitsUpper) < theta_) {
    28783531                    theta_ = CoinMin(hitsLower, hitsUpper);
     
    28853538          abort();
    28863539     }
     3540     theta_ = CoinMax(theta_,0.0);
     3541     // update solution
     3542     double * array = rowArray_[1]->denseVector();
     3543     int * index = rowArray_[1]->getIndices();
     3544     int number = rowArray_[1]->getNumElements();
     3545     for (int iRow = 0; iRow < number; iRow++) {
     3546       int iPivot = index[iRow];
     3547       iSequence = pivotVariable_[iPivot];
     3548       // solution value will be sol - theta*alpha
     3549       double alpha = array[iPivot];
     3550       solution_[iSequence] -= theta_ * alpha;
     3551     }
    28873552     if (pivotRow_ >= 0) {
    28883553          sequenceOut_ = pivotVariable_[pivotRow_];
    28893554          valueOut_ = solution_[sequenceOut_];
    2890           lowerOut_ = lower_[sequenceOut_];
    2891           upperOut_ = upper_[sequenceOut_];
     3555          lowerOut_ = lower_[sequenceOut_]+theta_*changeLower[sequenceOut_];
     3556          upperOut_ = upper_[sequenceOut_]+theta_*changeUpper[sequenceOut_];
    28923557          if (!toLower) {
    28933558               directionOut_ = -1;
    28943559               dualOut_ = valueOut_ - upperOut_;
    2895           } else if (valueOut_ < lowerOut_) {
     3560          } else {
    28963561               directionOut_ = 1;
    28973562               dualOut_ = lowerOut_ - valueOut_;
  • trunk/Clp/src/ClpSimplexOther.hpp

    r1525 r1643  
    8585                     const double * changeLowerRhs, const double * changeUpperRhs,
    8686                     const double * changeObjective);
     87     /** Version of parametrics which reads from file
     88         See CbcClpParam.cpp for details of format
     89         Returns -2 if unable to open file */
     90     int parametrics(const char * dataFile);
     91
    8792private:
    8893     /** Parametrics - inner loop
Note: See TracChangeset for help on using the changeset viewer.