source: branches/devel/Cbc/src/Cbc_ampl.cpp @ 529

Last change on this file since 529 was 529, checked in by forrest, 13 years ago

for nonlinear

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 34.2 KB
Line 
1/****************************************************************
2Copyright (C) 1997-2000 Lucent Technologies
3Modifications for Coin -  Copyright (C) 2006, International Business Machines Corporation and others.
4All Rights Reserved
5
6Permission to use, copy, modify, and distribute this software and
7its documentation for any purpose and without fee is hereby
8granted, provided that the above copyright notice appear in all
9copies and that both that the copyright notice and this
10permission notice and warranty disclaimer appear in supporting
11documentation, and that the name of Lucent or any of its entities
12not be used in advertising or publicity pertaining to
13distribution of the software without specific, written prior
14permission.
15
16LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
17INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
18IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
19SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
20WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
21IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
22ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
23THIS SOFTWARE.
24****************************************************************/
25
26#include "CbcConfig.h"
27#ifdef HAVE_UNISTD_H
28# include "unistd.h"
29#endif
30#include "CoinUtilsConfig.h"
31#include "CoinHelperFunctions.hpp"
32#include "CoinModel.hpp"
33#include "CoinSort.hpp"
34#include "CoinPackedMatrix.hpp"
35#include "CoinMpsIO.hpp"
36#include "CoinFloatEqual.hpp"
37
38extern "C" {
39# include "getstub.h"
40}
41
42#include "Cbc_ampl.h"
43#include <string>
44#include <cassert>
45/* so decodePhrase and clpCheck can access */
46static ampl_info * saveInfo=NULL;
47// Set to 1 if algorithm found
48static char algFound[20]="";
49static char*
50checkPhrase(Option_Info *oi, keyword *kw, char *v)
51{
52  if (strlen(v))
53    printf("string %s\n",v);
54  // Say algorithm found
55  strcpy(algFound,kw->desc);;
56  return v;
57}
58static char*
59checkPhrase2(Option_Info *oi, keyword *kw, char *v)
60{
61  if (strlen(v))
62    printf("string %s\n",v);
63  // put out keyword
64  saveInfo->arguments=(char **) realloc(saveInfo->arguments,(saveInfo->numberArguments+1)*sizeof(char *));
65  saveInfo->arguments[saveInfo->numberArguments++]=strdup(kw->desc);
66  return v;
67}
68static fint
69decodePhrase(char * phrase,ftnlen length)
70{
71  char * blank = strchr(phrase,' ');
72  if (blank) {
73    /* split arguments */
74    *blank='\0';
75    saveInfo->arguments=(char **) realloc(saveInfo->arguments,(saveInfo->numberArguments+2)*sizeof(char *));
76    saveInfo->arguments[saveInfo->numberArguments++]=strdup(phrase);
77    *blank=' ';
78    phrase=blank+1; /* move on */
79    if (strlen(phrase))
80      saveInfo->arguments[saveInfo->numberArguments++]=strdup(phrase);
81  } else if (strlen(phrase)) {
82    saveInfo->arguments=(char **) realloc(saveInfo->arguments,(saveInfo->numberArguments+1)*sizeof(char *));
83    saveInfo->arguments[saveInfo->numberArguments++]=strdup(phrase);
84  }
85  return 0;
86}
87static void
88sos_kludge(int nsos, int *sosbeg, double *sosref,int * sosind)
89{
90  // Adjust sosref if necessary to make monotonic increasing
91  int i, j, k;
92  // first sort
93  for (i=0;i<nsos;i++) {
94    k = sosbeg[i];
95    int end=sosbeg[i+1];
96    CoinSort_2(sosref+k,sosref+end,sosind+k);
97  }
98  double t, t1;
99  for(i = j = 0; i++ < nsos; ) {
100    k = sosbeg[i];
101    t = sosref[j];
102    while(++j < k) {
103      t1 = sosref[j];
104      t += 1e-10;
105      if (t1 <= t)
106        sosref[j] = t1 = t + 1e-10;
107      t = t1;
108    }
109  }
110}
111static char xxxxxx[20];
112#define VP (char*)
113 static keyword keywds[] = { /* must be sorted */
114        { "barrier",    checkPhrase,            (char *) xxxxxx ,"-barrier" },
115        { "dual",       checkPhrase,            (char *) xxxxxx , "-dualsimplex"},
116        { "help",       checkPhrase2,           (char *) xxxxxx , "-?"},
117        { "initial",    checkPhrase,            (char *) xxxxxx , "-initialsolve"},
118        { "max",        checkPhrase2,           (char *) xxxxxx , "-maximize"},
119        { "maximize",   checkPhrase2,           (char *) xxxxxx , "-maximize"},
120        { "primal",     checkPhrase,            (char *) xxxxxx , "-primalsimplex"},
121        { "quit",       checkPhrase2,           (char *) xxxxxx , "-quit"},
122        { "wantsol",    WS_val,         NULL, "write .sol file (without -AMPL)" }
123        };
124static Option_Info Oinfo = {"cbc", "Cbc 1.01", "cbc_options", keywds, nkeywds, 0, "",
125                                0,decodePhrase,0,0,0, 20060130 };
126// strdup used to avoid g++ compiler warning
127 static SufDecl
128suftab[] = {
129#if 0
130        { "current", 0, ASL_Sufkind_con | ASL_Sufkind_outonly },
131        { "current", 0, ASL_Sufkind_var | ASL_Sufkind_outonly },
132        { "direction", 0, ASL_Sufkind_var },
133        { "down", 0, ASL_Sufkind_con | ASL_Sufkind_outonly },
134        { "down", 0, ASL_Sufkind_var | ASL_Sufkind_outonly },
135        { "priority", 0, ASL_Sufkind_var },
136#endif
137        { "cut", 0, ASL_Sufkind_con },
138        { "direction", 0, ASL_Sufkind_var },
139        { "downPseudocost", 0, ASL_Sufkind_var | ASL_Sufkind_real },
140        { "priority", 0, ASL_Sufkind_var },
141        { "ref", 0, ASL_Sufkind_var | ASL_Sufkind_real },
142        { "sos", 0, ASL_Sufkind_var },
143        { "sos", 0, ASL_Sufkind_con },
144        { "sosno", 0, ASL_Sufkind_var | ASL_Sufkind_real },
145        { "sosref", 0, ASL_Sufkind_var | ASL_Sufkind_real },
146        { strdup("sstatus"), 0, ASL_Sufkind_var, 0 },
147        { strdup("sstatus"), 0, ASL_Sufkind_con, 0 },
148        { "upPseudocost", 0, ASL_Sufkind_var | ASL_Sufkind_real }
149#if 0
150        { "unbdd", 0, ASL_Sufkind_var | ASL_Sufkind_outonly},
151        { "up", 0, ASL_Sufkind_con | ASL_Sufkind_outonly },
152        { "up", 0, ASL_Sufkind_var | ASL_Sufkind_outonly }
153#endif
154        };
155#include "float.h"
156#include "limits.h"
157static ASL *asl=NULL;
158static FILE *nl=NULL;
159
160static void
161mip_stuff(void)
162{
163  int i;
164  double *pseudoUp, *pseudoDown;
165  int *priority, *direction;
166  // To label cuts
167  int *cut;
168  SufDesc *dpup, *dpdown, *dpri, *ddir, *dcut;
169 
170  ddir = suf_get("direction", ASL_Sufkind_var);
171  direction = ddir->u.i;
172  dpri = suf_get("priority", ASL_Sufkind_var);
173  priority = dpri->u.i;
174  dcut = suf_get("cut", ASL_Sufkind_con);
175  cut = dcut->u.i;
176  dpdown = suf_get("downPseudocost", ASL_Sufkind_var);
177  pseudoDown = dpdown->u.r;
178  dpup = suf_get("upPseudocost", ASL_Sufkind_var);
179  pseudoUp = dpup->u.r;
180  assert(saveInfo);
181  int numberColumns = saveInfo->numberColumns;
182  if (direction) {
183    int baddir=0;
184    saveInfo->branchDirection = (int *) malloc(numberColumns*sizeof(int));
185    for (i=0;i<numberColumns;i++) {
186      int value = direction[i];
187      if (value<-1||value>1) {
188        baddir++;
189        value=0;
190      }
191      saveInfo->branchDirection[i]=value;
192    }
193    if (baddir)
194      fprintf(Stderr,
195              "Treating %d .direction values outside [-1, 1] as 0.\n",
196              baddir);
197  }
198  if (priority) {
199    int badpri=0;
200    saveInfo->priorities= (int *) malloc(numberColumns*sizeof(int));
201    for (i=0;i<numberColumns;i++) {
202      int value = priority[i];
203      if (value<0) {
204        badpri++;
205        value=0;
206      }
207      saveInfo->priorities[i]=value;
208    }
209    if (badpri)
210      fprintf(Stderr,
211              "Treating %d negative .priority values as 0\n",
212              badpri);
213  }
214  int numberRows = saveInfo->numberRows;
215  if (cut) {
216    int badcut=0;
217    saveInfo->cut= (int *) malloc(numberRows*sizeof(int));
218    for (i=0;i<numberRows;i++) {
219      int value = cut[i];
220      if (value<0) {
221        badcut++;
222        value=0;
223      }
224      saveInfo->cut[i]=value;
225    }
226    if (badcut)
227      fprintf(Stderr,
228              "Treating %d negative cut values as 0\n",
229              badcut);
230  }
231  if (pseudoDown||pseudoUp) {
232    int badpseudo=0;
233    if (!pseudoDown||!pseudoUp)
234      fprintf(Stderr,
235              "Only one set of pseudocosts - assumed same\n");
236    saveInfo->pseudoDown= (double *) malloc(numberColumns*sizeof(double));
237    saveInfo->pseudoUp = (double *) malloc(numberColumns*sizeof(double));
238    for (i=0;i<numberColumns;i++) {
239      double valueD=0.0, valueU=0.0;
240      if (pseudoDown) {
241        valueD = pseudoDown[i];
242        if (valueD<0) {
243          badpseudo++;
244          valueD=0.0;
245        }
246      }
247      if (pseudoUp) {
248        valueU = pseudoUp[i];
249        if (valueU<0) {
250          badpseudo++;
251          valueU=0.0;
252        }
253      }
254      if (!valueD)
255        valueD=valueU;
256      if (!valueU)
257        valueU=valueD;
258      saveInfo->pseudoDown[i]=valueD;
259      saveInfo->pseudoUp[i]=valueU;
260    }
261    if (badpseudo)
262      fprintf(Stderr,
263              "Treating %d negative pseudoCosts as 0.0\n",badpseudo);
264  }
265}
266static void
267stat_map(int *stat, int n, int *map, int mx, const char *what)
268{
269  int bad, i, i1=0, j, j1=0;
270  static char badfmt[] = "Coin driver: %s[%d] = %d\n";
271 
272  for(i = bad = 0; i < n; i++) {
273    if ((j = stat[i]) >= 0 && j <= mx)
274      stat[i] = map[j];
275    else {
276      stat[i] = 0;
277      i1 = i;
278      j1 = j;
279      if (!bad++)
280        fprintf(Stderr, badfmt, what, i, j);
281    }
282  }
283  if (bad > 1) {
284    if (bad == 2)
285      fprintf(Stderr, badfmt, what, i1, j1);
286    else
287      fprintf(Stderr,
288              "Coin driver: %d messages about bad %s values suppressed.\n",
289              bad-1, what);
290  }
291}
292int
293readAmpl(ampl_info * info, int argc, char **argv, void ** coinModel)
294{
295  char *stub;
296  ograd *og;
297  int i;
298  SufDesc *csd;
299  SufDesc *rsd;
300  /*bool *basis, *lower;*/
301  /*double *LU, *c, lb, objadj, *rshift, *shift, t, ub, *x, *x0, *x1;*/
302  char * environment = getenv("cbc_options");
303  char tempBuffer[20];
304  double * obj;
305  double * columnLower;
306  double * columnUpper;
307  double * rowLower;
308  double * rowUpper;
309  char ** saveArgv=argv;
310  char fileName[1000];
311  if (argc>1)
312    strcpy(fileName,argv[1]);
313  else
314    fileName[0]='\0';
315  int saveArgc = argc;
316  if (info->numberRows!=-1234567)
317    memset(info,0,sizeof(ampl_info)); // overwrite unless magic number set
318  /* save so can be accessed by decodePhrase */
319  saveInfo = info;
320  info->numberArguments=0;
321  info->arguments=(char **) malloc(2*sizeof(char *));
322  info->arguments[info->numberArguments++]=strdup("ampl");
323  info->arguments[info->numberArguments++]=strdup("cbc");
324  asl = ASL_alloc(ASL_read_f);
325  stub = getstub(&argv, &Oinfo);
326  if (!stub)
327    usage_ASL(&Oinfo, 1);
328  nl = jac0dim(stub, 0);
329  suf_declare(suftab, sizeof(suftab)/sizeof(SufDecl));
330 
331  /* set A_vals to get the constraints column-wise (malloc so can be freed) */
332  A_vals = (double *) malloc(nzc*sizeof(double));
333  if (!A_vals) {
334    printf("no memory\n");
335    return 1;
336  }
337  /* say we want primal solution */
338  want_xpi0=1;
339  /* for basis info */
340  info->columnStatus = (int *) malloc(n_var*sizeof(int));
341  info->rowStatus = (int *) malloc(n_con*sizeof(int));
342  csd = suf_iput("sstatus", ASL_Sufkind_var, info->columnStatus);
343  rsd = suf_iput("sstatus", ASL_Sufkind_con, info->rowStatus);
344  // testosi parameter - if >= 10 then go in through coinModel
345  int testOsi=-1;
346  for (i=0;i<saveArgc;i++) {
347    if (!strncmp(saveArgv[i],"testosi",7)) {
348      testOsi = atoi(saveArgv[i+1]);
349      break;
350    }
351  }
352  if (!(nlvc+nlvo)||testOsi>=10) {
353    /* read linear model*/
354    f_read(nl,0);
355    // see if any sos
356    if (true) {
357      char *sostype;
358      int nsosnz, *sosbeg, *sosind, * sospri;
359      double *sosref;
360      int nsos;
361      int i = ASL_suf_sos_explict_free;
362      int copri[2], **p_sospri;
363      copri[0] = 0;
364      copri[1] = 0;
365      p_sospri = &sospri;
366      nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
367                     &sosbeg, &sosind, &sosref);
368      if (nsos) {
369        info->numberSos=nsos;
370        info->sosType = (char *) malloc(nsos);
371        info->sosPriority = (int *) malloc(nsos*sizeof(int));
372        info->sosStart = (int *) malloc((nsos+1)*sizeof(int));
373        info->sosIndices = (int *) malloc(nsosnz*sizeof(int));
374        info->sosReference = (double *) malloc(nsosnz*sizeof(double));
375        sos_kludge(nsos, sosbeg, sosref,sosind);
376        for (int i=0;i<nsos;i++) {
377          int ichar = sostype[i];
378          assert (ichar=='1'||ichar=='2');
379          info->sosType[i]=ichar-'0';
380        }       
381        memcpy(info->sosPriority,sospri,nsos*sizeof(int));
382        memcpy(info->sosStart,sosbeg,(nsos+1)*sizeof(int));
383        memcpy(info->sosIndices,sosind,nsosnz*sizeof(int));
384        memcpy(info->sosReference,sosref,nsosnz*sizeof(double));
385      }
386    }
387   
388    /*sos_finish(&specialOrderedInfo, 0, &j, 0, 0, 0, 0, 0);*/
389    Oinfo.uinfo = tempBuffer;
390    if (getopts(argv, &Oinfo))
391      return 1;
392    /* objective*/
393    obj = (double *) malloc(n_var*sizeof(double));
394    for (i=0;i<n_var;i++)
395      obj[i]=0.0;;
396    if (n_obj) {
397      for (og = Ograd[0];og;og = og->next)
398        obj[og->varno] = og->coef;
399    }
400    if (objtype[0])
401      info->direction=-1.0;
402    else
403      info->direction=1.0;
404    info->offset=objconst(0);
405    /* Column bounds*/
406    columnLower = (double *) malloc(n_var*sizeof(double));
407    columnUpper = (double *) malloc(n_var*sizeof(double));
408#define COIN_DBL_MAX DBL_MAX
409    for (i=0;i<n_var;i++) {
410      columnLower[i]=LUv[2*i];
411      if (columnLower[i]<= negInfinity)
412        columnLower[i]=-COIN_DBL_MAX;
413      columnUpper[i]=LUv[2*i+1];
414      if (columnUpper[i]>= Infinity)
415        columnUpper[i]=COIN_DBL_MAX;
416    }
417    /* Row bounds*/
418    rowLower = (double *) malloc(n_con*sizeof(double));
419    rowUpper = (double *) malloc(n_con*sizeof(double));
420    for (i=0;i<n_con;i++) {
421      rowLower[i]=LUrhs[2*i];
422      if (rowLower[i]<= negInfinity)
423        rowLower[i]=-COIN_DBL_MAX;
424      rowUpper[i]=LUrhs[2*i+1];
425      if (rowUpper[i]>= Infinity)
426        rowUpper[i]=COIN_DBL_MAX;
427    }
428    info->numberRows=n_con;
429    info->numberColumns=n_var;
430    info->numberElements=nzc;;
431    info->numberBinary=nbv;
432    info->numberIntegers=niv;
433    info->objective=obj;
434    info->rowLower=rowLower;
435    info->rowUpper=rowUpper;
436    info->columnLower=columnLower;
437    info->columnUpper=columnUpper;
438    info->starts=A_colstarts;
439    /*A_colstarts=NULL;*/
440    info->rows=A_rownos;
441    /*A_rownos=NULL;*/
442    info->elements=A_vals;
443    /*A_vals=NULL;*/
444    info->primalSolution=NULL;
445    /* put in primalSolution if exists */
446    if (X0) {
447      info->primalSolution=(double *) malloc(n_var*sizeof(double));
448      memcpy(info->primalSolution,X0,n_var*sizeof(double));
449    }
450    info->dualSolution=NULL;
451    if (niv+nbv>0)
452      mip_stuff(); // get any extra info
453    if ((!(niv+nbv)&&(csd->kind & ASL_Sufkind_input))
454        ||(rsd->kind & ASL_Sufkind_input)) {
455      /* convert status - need info on map */
456      static int map[] = {1, 3, 1, 1, 2, 1, 1};
457      stat_map(info->columnStatus, n_var, map, 6, "incoming columnStatus");
458      stat_map(info->rowStatus, n_con, map, 6, "incoming rowStatus");
459    } else {
460      /* all slack basis */
461      // leave status for output */
462#if 0
463      free(info->rowStatus);
464      info->rowStatus=NULL;
465      free(info->columnStatus);
466      info->columnStatus=NULL;
467#endif
468    }
469  } else {
470    // QP
471    // Add .nl if not there
472    if (!strstr(fileName,".nl"))
473      strcat(fileName,".nl");
474    CoinModel * model = new CoinModel(1,fileName,info);
475    if (model->numberRows()>0)
476      *coinModel=(void *) model;
477    Oinfo.uinfo = tempBuffer;
478    if (getopts(argv, &Oinfo))
479      return 1;
480    Oinfo.wantsol=1;
481    if (objtype[0])
482      info->direction=-1.0;
483    else
484      info->direction=1.0;
485    info->offset=objconst(0);
486    info->numberRows=n_con;
487    info->numberColumns=n_var;
488    info->numberElements=nzc;;
489    info->numberBinary=nbv;
490    info->numberIntegers=niv;
491    if (nbv+niv+nlvbi+nlvci+nlvoi>0) {
492      mip_stuff(); // get any extra info
493      if (info->cut) 
494        model->setCutMarker(info->numberRows,info->cut);
495      if (info->priorities)
496        model->setPriorities(info->numberColumns,info->priorities);
497    }
498  }
499  /* add -solve - unless something there already
500   - also check for sleep=yes */
501  {
502    int found=0;
503    int foundLog=0;
504    int foundSleep=0;
505    const char * something[]={"solve","branch","duals","primals","user"};
506    for (i=0;i<info->numberArguments;i++) {
507      unsigned int j;
508      const char * argument = info->arguments[i];
509      for (j=0;j<sizeof(something)/sizeof(char *);j++) {
510        const char * check = something[j];
511        if (!strncmp(argument,check,sizeof(check))) {
512          found=(int)(j+1);
513        } else if (!strncmp(argument,"log",3)) {
514          foundLog=1;
515        } else if (!strncmp(argument,"sleep",5)) {
516          foundSleep=1;
517        }
518      }
519    }
520    if (foundLog) {
521      /* print options etc */
522      for (i=0;i<saveArgc;i++)
523        printf("%s ",saveArgv[i]);
524      printf("\n");
525      if (environment)
526        printf("env %s\n",environment);
527      /*printf("%d rows %d columns %d elements\n",n_con,n_var,nzc);*/
528    }
529    if (!found) {
530      if (!strlen(algFound)) {
531        info->arguments=(char **) realloc(info->arguments,(info->numberArguments+1)*sizeof(char *));
532        info->arguments[info->numberArguments++]=strdup("-solve");
533      } else {
534        // use algorithm from keyword
535        info->arguments=(char **) realloc(info->arguments,(info->numberArguments+1)*sizeof(char *));
536        info->arguments[info->numberArguments++]=strdup(algFound);
537      }
538    }
539    if (foundSleep) {
540      /* let user copy .nl file */
541      fprintf(stderr,"You can copy .nl file %s for debug purposes or attach debugger\n",saveArgv[1]);
542      fprintf(stderr,"Type q to quit, anything else to continue\n");
543      char getChar = getc(stdin);
544      if (getChar=='q'||getChar=='Q')
545        exit(1);
546    }
547  }
548  /* add -quit */
549  info->arguments=(char **) realloc(info->arguments,(info->numberArguments+1)*sizeof(char *));
550  info->arguments[info->numberArguments++]=strdup("-quit");
551  return 0;
552}
553void freeArrays1(ampl_info * info)
554{
555  free(info->objective);
556  info->objective=NULL;
557  free(info->rowLower);
558  info->rowLower=NULL;
559  free(info->rowUpper);
560  info->rowUpper=NULL;
561  free(info->columnLower);
562  info->columnLower=NULL;
563  free(info->columnUpper);
564  info->columnUpper=NULL;
565  /* this one not freed by ASL_free */
566  free(info->elements);
567  info->elements=NULL;
568  free(info->primalSolution);
569  info->primalSolution=NULL;
570  free(info->dualSolution);
571  info->dualSolution=NULL;
572  /*free(info->rowStatus);
573  info->rowStatus=NULL;
574  free(info->columnStatus);
575  info->columnStatus=NULL;*/
576}
577void freeArrays2(ampl_info * info)
578{
579  free(info->primalSolution);
580  info->primalSolution=NULL;
581  free(info->dualSolution);
582  info->dualSolution=NULL;
583  free(info->rowStatus);
584  info->rowStatus=NULL;
585  free(info->columnStatus);
586  info->columnStatus=NULL;
587  free(info->priorities);
588  info->priorities=NULL;
589  free(info->branchDirection);
590  info->branchDirection=NULL;
591  free(info->pseudoDown);
592  info->pseudoDown=NULL;
593  free(info->pseudoUp);
594  info->pseudoUp=NULL;
595  free(info->sosType);
596  info->sosType=NULL;
597  free(info->sosPriority);
598  info->sosPriority=NULL;
599  free(info->sosStart);
600  info->sosStart=NULL;
601  free(info->sosIndices);
602  info->sosIndices=NULL;
603  free(info->sosReference);
604  info->sosReference=NULL;
605  free(info->cut);
606  info->cut=NULL;
607  ASL_free(&asl);
608}
609void freeArgs(ampl_info * info)
610{
611  int i;
612  for ( i=0; i<info->numberArguments;i++)
613    free(info->arguments[i]);
614  free(info->arguments);
615}
616int ampl_obj_prec()
617{
618  return obj_prec();
619}
620void writeAmpl(ampl_info * info)
621{
622  char buf[1000];
623  typedef struct { const char *msg; int code; int wantObj; } Sol_info;
624  static Sol_info solinfo[] = {
625    { "optimal solution",                       000, 1 },
626    { "infeasible",                             200, 1 },
627    { "unbounded",                              300, 0 },
628    { "iteration limit etc",                    400, 1 },
629    { "solution limit",                         401, 1 },
630    { "ran out of space",                       500, 0 },
631    { "status unknown",                         501, 1 },
632    { "bug!",                                   502, 0 },
633    { "best MIP solution so far restored",      101, 1 },
634    { "failed to restore best MIP solution",    503, 1 },
635    { "optimal (?) solution",                   100, 1 }
636  };
637  /* convert status - need info on map */
638  static int map[] = {0, 3, 4, 1};
639  sprintf(buf,"%s %s",Oinfo.bsname,info->buffer);
640  solve_result_num = solinfo[info->problemStatus].code;
641  if (info->columnStatus) {
642    stat_map(info->columnStatus, n_var, map, 4, "outgoing columnStatus");
643    stat_map(info->rowStatus, n_con, map, 4, "outgoing rowStatus");
644    suf_iput("sstatus", ASL_Sufkind_var, info->columnStatus);
645    suf_iput("sstatus", ASL_Sufkind_con, info->rowStatus);
646  }
647  write_sol(buf,info->primalSolution,info->dualSolution,&Oinfo);
648}
649/* Read a problem from AMPL nl file
650 */
651CoinModel::CoinModel( int nonLinear, const char * fileName,const void * info)
652 :  numberRows_(0),
653    maximumRows_(0),
654    numberColumns_(0),
655    maximumColumns_(0),
656    numberElements_(0),
657    maximumElements_(0),
658    numberQuadraticElements_(0),
659    maximumQuadraticElements_(0),
660    optimizationDirection_(1.0),
661    objectiveOffset_(0.0),
662    rowLower_(NULL),
663    rowUpper_(NULL),
664    rowType_(NULL),
665    objective_(NULL),
666    columnLower_(NULL),
667    columnUpper_(NULL),
668    integerType_(NULL),
669    columnType_(NULL),
670    start_(NULL),
671    elements_(NULL),
672    quadraticElements_(NULL),
673    sortIndices_(NULL),
674    sortElements_(NULL),
675    sortSize_(0),
676    sizeAssociated_(0),
677    associated_(NULL),
678    numberSOS_(0),
679    startSOS_(NULL),
680    memberSOS_(NULL),
681    typeSOS_(NULL),
682    prioritySOS_(NULL),
683    referenceSOS_(NULL),
684    priority_(NULL),
685    cut_(NULL),
686    logLevel_(0),
687    type_(-1),
688    links_(0)
689{
690  problemName_ = strdup("");
691  int status=0;
692  if (!strcmp(fileName,"-")||!strcmp(fileName,"stdin")) {
693    // stdin
694  } else {
695    std::string name=fileName;
696    bool readable = fileCoinReadable(name);
697    if (!readable) {
698      std::cerr<<"Unable to open file "
699               <<fileName<<std::endl;
700      status = -1;
701    }
702  }
703  if (!status) {
704    gdb(nonLinear,fileName,info);
705  }
706}
707#if 0
708static real
709qterm(ASL *asl, fint *colq, fint *rowq, real *delsq)
710{
711  double t, t1, *x, *x0, *xe;
712  fint *rq0, *rqe;
713 
714  t = 0.;
715  x0 = x = X0;
716  xe = x + n_var;
717  rq0 = rowq;
718  while(x < xe) {
719    t1 = *x++;
720    rqe = rq0 + *++colq;
721    while(rowq < rqe)
722      t += t1*x0[*rowq++]**delsq++;
723  }
724  return 0.5 * t;
725}
726#endif
727void
728CoinModel::gdb( int nonLinear, const char * fileName,const void * info)
729{
730  const ampl_info * amplInfo = (const ampl_info *) info;
731  ograd *og=NULL;
732  int i;
733  SufDesc *csd=NULL;
734  SufDesc *rsd=NULL;
735  /*bool *basis, *lower;*/
736  /*double *LU, *c, lb, objadj, *rshift, *shift, t, ub, *x, *x0, *x1;*/
737  //char tempBuffer[20];
738  double * objective=NULL;
739  double * columnLower=NULL;
740  double * columnUpper=NULL;
741  double * rowLower=NULL;
742  double * rowUpper=NULL;
743  int * columnStatus=NULL;
744  int * rowStatus=NULL;
745  int numberRows=-1;
746  int numberColumns=-1;
747  int numberElements=-1;
748  int numberBinary=-1;
749  int numberIntegers=-1;
750  int numberAllNonLinearBoth=0;
751  int numberIntegerNonLinearBoth=0;
752  int numberAllNonLinearConstraints=0;
753  int numberIntegerNonLinearConstraints=0;
754  int numberAllNonLinearObjective=0;
755  int numberIntegerNonLinearObjective=0;
756  double * primalSolution=NULL;
757  double direction=1.0;
758  char * stub = strdup(fileName);
759  CoinPackedMatrix matrixByRow;
760  fint ** colqp=NULL;
761  int *z = NULL;
762  if (nonLinear==0) {
763    // linear
764    asl = ASL_alloc(ASL_read_f);
765    nl = jac0dim(stub, 0);
766    free(stub);
767    suf_declare(suftab, sizeof(suftab)/sizeof(SufDecl));
768   
769    /* set A_vals to get the constraints column-wise (malloc so can be freed) */
770    A_vals = (double *) malloc(nzc*sizeof(double));
771    if (!A_vals) {
772      printf("no memory\n");
773      return ;
774    }
775    /* say we want primal solution */
776    want_xpi0=1;
777    /* for basis info */
778    columnStatus = (int *) malloc(n_var*sizeof(int));
779    rowStatus = (int *) malloc(n_con*sizeof(int));
780    csd = suf_iput("sstatus", ASL_Sufkind_var, columnStatus);
781    rsd = suf_iput("sstatus", ASL_Sufkind_con, rowStatus);
782    /* read linear model*/
783    f_read(nl,0);
784    // see if any sos
785    if (true) {
786      char *sostype;
787      int nsosnz, *sosbeg, *sosind, * sospri;
788      double *sosref;
789      int nsos;
790      int i = ASL_suf_sos_explict_free;
791      int copri[2], **p_sospri;
792      copri[0] = 0;
793      copri[1] = 0;
794      p_sospri = &sospri;
795      nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
796                     &sosbeg, &sosind, &sosref);
797      if (nsos) {
798        abort();
799#if 0
800        info->numberSos=nsos;
801        info->sosType = (char *) malloc(nsos);
802        info->sosPriority = (int *) malloc(nsos*sizeof(int));
803        info->sosStart = (int *) malloc((nsos+1)*sizeof(int));
804        info->sosIndices = (int *) malloc(nsosnz*sizeof(int));
805        info->sosReference = (double *) malloc(nsosnz*sizeof(double));
806        sos_kludge(nsos, sosbeg, sosref,sosind);
807        for (int i=0;i<nsos;i++) {
808          int ichar = sostype[i];
809          assert (ichar=='1'||ichar=='2');
810          info->sosType[i]=ichar-'0';
811        }       
812        memcpy(info->sosPriority,sospri,nsos*sizeof(int));
813        memcpy(info->sosStart,sosbeg,(nsos+1)*sizeof(int));
814        memcpy(info->sosIndices,sosind,nsosnz*sizeof(int));
815        memcpy(info->sosReference,sosref,nsosnz*sizeof(double));
816#endif
817      }
818    }
819   
820    /*sos_finish(&specialOrderedInfo, 0, &j, 0, 0, 0, 0, 0);*/
821    //Oinfo.uinfo = tempBuffer;
822    //if (getopts(argv, &Oinfo))
823    //return 1;
824    /* objective*/
825    objective = (double *) malloc(n_var*sizeof(double));
826    for (i=0;i<n_var;i++)
827      objective[i]=0.0;;
828    if (n_obj) {
829      for (og = Ograd[0];og;og = og->next)
830        objective[og->varno] = og->coef;
831    }
832    if (objtype[0])
833      direction=-1.0;
834    else
835      direction=1.0;
836    objectiveOffset_=objconst(0);
837    /* Column bounds*/
838    columnLower = (double *) malloc(n_var*sizeof(double));
839    columnUpper = (double *) malloc(n_var*sizeof(double));
840#define COIN_DBL_MAX DBL_MAX
841    for (i=0;i<n_var;i++) {
842      columnLower[i]=LUv[2*i];
843      if (columnLower[i]<= negInfinity)
844        columnLower[i]=-COIN_DBL_MAX;
845      columnUpper[i]=LUv[2*i+1];
846      if (columnUpper[i]>= Infinity)
847        columnUpper[i]=COIN_DBL_MAX;
848    }
849    /* Row bounds*/
850    rowLower = (double *) malloc(n_con*sizeof(double));
851    rowUpper = (double *) malloc(n_con*sizeof(double));
852    for (i=0;i<n_con;i++) {
853      rowLower[i]=LUrhs[2*i];
854      if (rowLower[i]<= negInfinity)
855        rowLower[i]=-COIN_DBL_MAX;
856      rowUpper[i]=LUrhs[2*i+1];
857      if (rowUpper[i]>= Infinity)
858        rowUpper[i]=COIN_DBL_MAX;
859    }
860    numberRows=n_con;
861    numberColumns=n_var;
862    numberElements=nzc;;
863    numberBinary=nbv;
864    numberIntegers=niv;
865    /* put in primalSolution if exists */
866    if (X0) {
867      primalSolution=(double *) malloc(n_var*sizeof(double));
868      memcpy( primalSolution,X0,n_var*sizeof(double));
869    }
870    //double * dualSolution=NULL;
871    if (niv+nbv>0)
872      mip_stuff(); // get any extra info
873    if ((!(niv+nbv)&&(csd->kind & ASL_Sufkind_input))
874        ||(rsd->kind & ASL_Sufkind_input)) {
875      /* convert status - need info on map */
876      static int map[] = {1, 3, 1, 1, 2, 1, 1};
877      stat_map(columnStatus, n_var, map, 6, "incoming columnStatus");
878      stat_map(rowStatus, n_con, map, 6, "incoming rowStatus");
879    } else {
880      /* all slack basis */
881      // leave status for output */
882#if 0
883      free(rowStatus);
884      rowStatus=NULL;
885      free(columnStatus);
886      columnStatus=NULL;
887#endif
888    }
889    CoinPackedMatrix columnCopy(true,numberRows,numberColumns,numberElements,
890                                A_vals,A_rownos,A_colstarts,NULL);
891    matrixByRow.reverseOrderedCopyOf(columnCopy);
892  } else if (nonLinear==1) {
893    // quadratic
894    asl = ASL_alloc(ASL_read_fg);
895    nl = jac0dim(stub, (long) strlen(stub));
896    free(stub);
897    suf_declare(suftab, sizeof(suftab)/sizeof(SufDecl));
898    /* read  model*/
899    X0 = (double*) malloc(n_var*sizeof(double));
900    CoinZeroN(X0,n_var);
901    qp_read(nl,0);
902    assert (n_obj==1);
903    int nz = 1 + n_con;
904    colqp = (fint**) malloc(nz*(2*sizeof(int*)
905                                      + sizeof(double*)));
906    fint ** rowqp = colqp + nz;
907    double ** delsqp = (double **)(rowqp + nz);
908    z = (int*) malloc(nz*sizeof(int));
909    for (i=0;i<=n_con;i++) {
910      z[i] = nqpcheck(-i, rowqp+i, colqp+i, delsqp+i);
911    }
912    qp_opify();
913    /* objective*/
914    objective = (double *) malloc(n_var*sizeof(double));
915    for (i=0;i<n_var;i++)
916      objective[i]=0.0;;
917    if (n_obj) {
918      for (og = Ograd[0];og;og = og->next)
919        objective[og->varno] = og->coef;
920    }
921    if (objtype[0])
922      direction=-1.0;
923    else
924      direction=1.0;
925    objectiveOffset_=objconst(0);
926    /* Column bounds*/
927    columnLower = (double *) malloc(n_var*sizeof(double));
928    columnUpper = (double *) malloc(n_var*sizeof(double));
929    for (i=0;i<n_var;i++) {
930      columnLower[i]=LUv[2*i];
931      if (columnLower[i]<= negInfinity)
932        columnLower[i]=-COIN_DBL_MAX;
933      columnUpper[i]=LUv[2*i+1];
934      if (columnUpper[i]>= Infinity)
935        columnUpper[i]=COIN_DBL_MAX;
936    }
937    // Build by row from scratch
938    //matrixByRow.reserve(n_var,nzc,true);
939    // say row orderded
940    matrixByRow.transpose();
941    /* Row bounds*/
942    rowLower = (double *) malloc(n_con*sizeof(double));
943    rowUpper = (double *) malloc(n_con*sizeof(double));
944    int * column = new int [n_var];
945    double * element = new double [n_var];
946    for (i=0;i<n_con;i++) {
947      rowLower[i]=LUrhs[2*i];
948      if (rowLower[i]<= negInfinity)
949        rowLower[i]=-COIN_DBL_MAX;
950      rowUpper[i]=LUrhs[2*i+1];
951      if (rowUpper[i]>= Infinity)
952        rowUpper[i]=COIN_DBL_MAX;
953      int k=0;
954      for(cgrad * cg = Cgrad[i]; cg; cg = cg->next) {
955        column[k]=cg->varno;
956        element[k++]=cg->coef;
957      }
958      matrixByRow.appendRow(k,column,element);
959    }
960    delete [] column;
961    delete [] element;
962    numberRows=n_con;
963    numberColumns=n_var;
964    numberElements=nzc;;
965    numberBinary=nbv;
966    numberIntegers=niv;
967    numberAllNonLinearBoth=nlvb;
968    numberIntegerNonLinearBoth=nlvbi;
969    numberAllNonLinearConstraints=nlvc;
970    numberIntegerNonLinearConstraints=nlvci;
971    numberAllNonLinearObjective=nlvo;
972    numberIntegerNonLinearObjective=nlvoi;
973    /* say we want primal solution */
974    want_xpi0=1;
975    //double * dualSolution=NULL;
976  } else {
977    abort();
978  }
979  // set problem name
980  free(problemName_);
981  problemName_=strdup("???");
982
983  // Build by row from scratch
984  const double * element = matrixByRow.getElements();
985  const int * column = matrixByRow.getIndices();
986  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
987  const int * rowLength = matrixByRow.getVectorLengths();
988  for (i=0;i<numberRows;i++) {
989    addRow(rowLength[i],column+rowStart[i],
990           element+rowStart[i],rowLower[i],rowUpper[i]);
991  }
992  // Now do column part
993  for (i=0;i<numberColumns;i++) {
994    setColumnBounds(i,columnLower[i],columnUpper[i]);
995    setColumnObjective(i,objective[i]);
996  }
997  for ( i=numberColumns-numberBinary-numberIntegers;
998        i<numberColumns;i++) {
999      setColumnIsInteger(i,true);;
1000  }
1001  // and non linear
1002  for (i=numberAllNonLinearBoth-numberIntegerNonLinearBoth;
1003       i<numberAllNonLinearBoth;i++) {
1004    setColumnIsInteger(i,true);;
1005  }
1006  for (i=numberAllNonLinearConstraints-numberIntegerNonLinearConstraints;
1007       i<numberAllNonLinearConstraints;i++) {
1008    setColumnIsInteger(i,true);;
1009  }
1010  for (i=numberAllNonLinearObjective-numberIntegerNonLinearObjective;
1011       i<numberAllNonLinearObjective;i++) {
1012    setColumnIsInteger(i,true);;
1013  }
1014  free(columnLower);
1015  free(columnUpper);
1016  free(rowLower);
1017  free(rowUpper);
1018  free(objective);
1019  // do names
1020  int iRow;
1021  for (iRow=0;iRow<numberRows_;iRow++) {
1022    char name[9];
1023    sprintf(name,"r%7.7d",iRow);
1024    setRowName(iRow,name);
1025  }
1026  int iColumn;
1027  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1028    char name[9];
1029    sprintf(name,"c%7.7d",iColumn);
1030    setColumnName(iColumn,name);
1031  }
1032  if (colqp) {
1033    // add in quadratic
1034    int nz = 1 + n_con;
1035    int nOdd=0;
1036    fint ** rowqp = colqp + nz;
1037    double ** delsqp = (double **)(rowqp + nz);
1038    for (i=0;i<=n_con;i++) {
1039      int nels = z[i];
1040      if (nels) {
1041        double * element = delsqp[i];
1042        int * start = (int *) colqp[i];
1043        int * row = (int *) rowqp[i];
1044        if (!element) {
1045          // odd row - probably not quadratic
1046          nOdd++;
1047          continue;
1048        }
1049#if 0
1050        printf("%d quadratic els\n",nels);
1051        for (int j=0;j<n_var;j++) {
1052          for (int k=start[j];k<start[j+1];k++)
1053            printf("%d %d %g\n",j,row[k],element[k]);
1054        }
1055#endif
1056        if (i) {
1057          int iRow = i-1;
1058          for (int j=0;j<n_var;j++) {
1059            for (int k=start[j];k<start[j+1];k++) {
1060              int kColumn = row[k];
1061              double value = element[k];
1062              // ampl gives twice with assumed 0.5
1063              if (kColumn<j)
1064                continue;
1065              else if (kColumn==j)
1066                value *= 0.5;
1067              const char * expr = getElementAsString(iRow,j);
1068              double constant=0.0;
1069              bool linear;
1070              if (expr&&strcmp(expr,"Numeric")) {
1071                linear=false;
1072              } else {
1073                constant = getElement(iRow,j);
1074                linear=true;
1075              }
1076              char temp[1000];
1077              char temp2[30];
1078              if (value==1.0) 
1079                sprintf(temp2,"c%7.7d",kColumn);
1080              else
1081                sprintf(temp2,"%g*c%7.7d",value,kColumn);
1082              if (linear) {
1083                if (!constant)
1084                  strcpy(temp,temp2);
1085                else if (value>0.0) 
1086                  sprintf(temp,"%g+%s",constant,temp2);
1087                else
1088                  sprintf(temp,"%g%s",constant,temp2);
1089              } else {
1090                if (value>0.0) 
1091                  sprintf(temp,"%s+%s",expr,temp2);
1092                else
1093                  sprintf(temp,"%s%s",expr,temp2);
1094              }
1095              assert (strlen(temp)<1000);
1096              setElement(iRow,j,temp);
1097              if (amplInfo->logLevel>1)
1098                printf("el for row %d column c%7.7d is %s\n",iRow,j,temp);
1099            }
1100          }
1101        } else {
1102          // objective
1103          for (int j=0;j<n_var;j++) {
1104            for (int k=start[j];k<start[j+1];k++) {
1105              int kColumn = row[k];
1106              double value = element[k];
1107              // ampl gives twice with assumed 0.5
1108              if (kColumn<j)
1109                continue;
1110              else if (kColumn==j)
1111                value *= 0.5;
1112              const char * expr = getColumnObjectiveAsString(j);
1113              double constant=0.0;
1114              bool linear;
1115              if (expr&&strcmp(expr,"Numeric")) {
1116                linear=false;
1117              } else {
1118                constant = getColumnObjective(j);
1119                linear=true;
1120              }
1121              char temp[1000];
1122              char temp2[30];
1123              if (value==1.0) 
1124                sprintf(temp2,"c%7.7d",kColumn);
1125              else
1126                sprintf(temp2,"%g*c%7.7d",value,kColumn);
1127              if (linear) {
1128                if (!constant)
1129                  strcpy(temp,temp2);
1130                else if (value>0.0) 
1131                  sprintf(temp,"%g+%s",constant,temp2);
1132                else
1133                  sprintf(temp,"%g%s",constant,temp2);
1134              } else {
1135                if (value>0.0) 
1136                  sprintf(temp,"%s+%s",expr,temp2);
1137                else
1138                  sprintf(temp,"%s%s",expr,temp2);
1139              }
1140              assert (strlen(temp)<1000);
1141              setObjective(j,temp);
1142              if (amplInfo->logLevel>1)
1143                printf("el for objective column c%7.7d is %s\n",j,temp);
1144            }
1145          }
1146        }
1147      }
1148    }
1149    if (nOdd) {
1150      printf("%d non-linear constraints could not be converted to quadratic\n",nOdd);
1151      exit(77);
1152    }
1153  }
1154  free(colqp);
1155  free(z);
1156  // see if any sos
1157  {
1158    char *sostype;
1159    int nsosnz, *sosbeg, *sosind, * sospri;
1160    double *sosref;
1161    int nsos;
1162    int i = ASL_suf_sos_explict_free;
1163    int copri[2], **p_sospri;
1164    copri[0] = 0;
1165    copri[1] = 0;
1166    p_sospri = &sospri;
1167    nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
1168                   &sosbeg, &sosind, &sosref);
1169    if (nsos) {
1170      numberSOS_=nsos;
1171      typeSOS_ = new int [numberSOS_];
1172      prioritySOS_ = new int [numberSOS_];
1173      startSOS_ = new int [numberSOS_+1];
1174      memberSOS_ = new int[nsosnz];
1175      referenceSOS_ = new double [nsosnz];
1176      sos_kludge(nsos, sosbeg, sosref,sosind);
1177      for (int i=0;i<nsos;i++) {
1178        int ichar = sostype[i];
1179        assert (ichar=='1'||ichar=='2');
1180        typeSOS_[i]=ichar-'0';
1181      } 
1182      memcpy(prioritySOS_,sospri,nsos*sizeof(int));
1183      memcpy(startSOS_,sosbeg,(nsos+1)*sizeof(int));
1184      memcpy(memberSOS_,sosind,nsosnz*sizeof(int));
1185      memcpy(referenceSOS_,sosref,nsosnz*sizeof(double));
1186    }
1187  }
1188}
Note: See TracBrowser for help on using the repository browser.