source: trunk/Cbc/src/Cbc_ampl.cpp @ 369

Last change on this file since 369 was 369, checked in by andreasw, 14 years ago

Enabled Ampl executable for Cbc
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 17.0 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 "getstub.h"
27#include "Cbc_ampl.h"
28#include "unistd.h"
29#include <string>
30#include <cassert>
31/* so decodePhrase and clpCheck can access */
32static ampl_info * saveInfo=NULL;
33// Set to 1 if algorithm found
34static char algFound[20]="";
35static char*
36checkPhrase(Option_Info *oi, keyword *kw, char *v)
37{
38  if (strlen(v))
39    printf("string %s\n",v);
40  // Say algorithm found
41  strcpy(algFound,kw->desc);;
42  return v;
43}
44static char*
45checkPhrase2(Option_Info *oi, keyword *kw, char *v)
46{
47  if (strlen(v))
48    printf("string %s\n",v);
49  // put out keyword
50  saveInfo->arguments=(char **) realloc(saveInfo->arguments,(saveInfo->numberArguments+1)*sizeof(char *));
51  saveInfo->arguments[saveInfo->numberArguments++]=strdup(kw->desc);
52  return v;
53}
54static fint
55decodePhrase(char * phrase,ftnlen length)
56{
57  char * blank = strchr(phrase,' ');
58  if (blank) {
59    /* split arguments */
60    *blank='\0';
61    saveInfo->arguments=(char **) realloc(saveInfo->arguments,(saveInfo->numberArguments+2)*sizeof(char *));
62    saveInfo->arguments[saveInfo->numberArguments++]=strdup(phrase);
63    *blank=' ';
64    phrase=blank+1; /* move on */
65    saveInfo->arguments[saveInfo->numberArguments++]=strdup(phrase);
66  } else {
67    saveInfo->arguments=(char **) realloc(saveInfo->arguments,(saveInfo->numberArguments+1)*sizeof(char *));
68    saveInfo->arguments[saveInfo->numberArguments++]=strdup(phrase);
69  }
70  return 0;
71}
72static void
73sos_kludge(int nsos, int *sosbeg, double *sosref)
74{
75  // Adjust sosref if necessary to make monotonic increasing
76  int i, j, k;
77  double t, t1;
78  for(i = j = 0; i++ < nsos; ) {
79    k = sosbeg[i];
80    t = sosref[j];
81    while(++j < k) {
82      t1 = sosref[j];
83      t += 1e-10;
84      if (t1 <= t)
85        sosref[j] = t1 = t + 1e-10;
86      t = t1;
87    }
88  }
89}
90static char xxxxxx[20];
91#define VP (char*)
92 static keyword keywds[] = { /* must be sorted */
93        { "barrier",    checkPhrase,            (char *) xxxxxx ,"-barrier" },
94        { "dual",       checkPhrase,            (char *) xxxxxx , "-dualsimplex"},
95        { "help",       checkPhrase2,           (char *) xxxxxx , "-?"},
96        { "initial",    checkPhrase,            (char *) xxxxxx , "-initialsolve"},
97        { "max",        checkPhrase2,           (char *) xxxxxx , "-maximize"},
98        { "maximize",   checkPhrase2,           (char *) xxxxxx , "-maximize"},
99        { "primal",     checkPhrase,            (char *) xxxxxx , "-primalsimplex"},
100        { "quit",       checkPhrase2,           (char *) xxxxxx , "-quit"},
101        { "wantsol",    WS_val,         NULL, "write .sol file (without -AMPL)" }
102        };
103static Option_Info Oinfo = {"cbc", "Cbc 1.01", "cbc_options", keywds, nkeywds, 0, "",
104                                0,decodePhrase,0,0,0, 20060130 };
105// strdup used to avoid g++ compiler warning
106 static SufDecl
107suftab[] = {
108#if 0
109        { "current", 0, ASL_Sufkind_con | ASL_Sufkind_outonly },
110        { "current", 0, ASL_Sufkind_var | ASL_Sufkind_outonly },
111        { "direction", 0, ASL_Sufkind_var },
112        { "down", 0, ASL_Sufkind_con | ASL_Sufkind_outonly },
113        { "down", 0, ASL_Sufkind_var | ASL_Sufkind_outonly },
114        { "priority", 0, ASL_Sufkind_var },
115#endif
116        { "direction", 0, ASL_Sufkind_var },
117        { "downPseudocost", 0, ASL_Sufkind_var | ASL_Sufkind_real },
118        { "priority", 0, ASL_Sufkind_var },
119        { "ref", 0, ASL_Sufkind_var | ASL_Sufkind_real },
120        { "sos", 0, ASL_Sufkind_var },
121        { "sos", 0, ASL_Sufkind_con },
122        { "sosno", 0, ASL_Sufkind_var | ASL_Sufkind_real },
123        { "sosref", 0, ASL_Sufkind_var | ASL_Sufkind_real },
124        { strdup("sstatus"), 0, ASL_Sufkind_var, 0 },
125        { strdup("sstatus"), 0, ASL_Sufkind_con, 0 },
126        { "upPseudocost", 0, ASL_Sufkind_var | ASL_Sufkind_real }
127#if 0
128        { "unbdd", 0, ASL_Sufkind_var | ASL_Sufkind_outonly},
129        { "up", 0, ASL_Sufkind_con | ASL_Sufkind_outonly },
130        { "up", 0, ASL_Sufkind_var | ASL_Sufkind_outonly }
131#endif
132        };
133#include "float.h"
134#include "limits.h"
135static ASL *asl=NULL;
136static FILE *nl=NULL;
137
138static void
139mip_stuff(void)
140{
141  int i;
142  double *pseudoUp, *pseudoDown;
143  int *priority, *direction;
144  SufDesc *dpup, *dpdown, *dpri, *ddir;
145 
146  ddir = suf_get("direction", ASL_Sufkind_var);
147  direction = ddir->u.i;
148  dpri = suf_get("priority", ASL_Sufkind_var);
149  priority = dpri->u.i;
150  dpdown = suf_get("downPseudocost", ASL_Sufkind_var);
151  pseudoDown = dpdown->u.r;
152  dpup = suf_get("upPseudocost", ASL_Sufkind_var);
153  pseudoUp = dpup->u.r;
154  assert(saveInfo);
155  int numberColumns = saveInfo->numberColumns;
156  if (direction) {
157    int baddir=0;
158    saveInfo->branchDirection = (int *) malloc(numberColumns*sizeof(int));
159    for (i=0;i<numberColumns;i++) {
160      int value = direction[i];
161      if (value<-1||value>1) {
162        baddir++;
163        value=0;
164      }
165      saveInfo->branchDirection[i]=value;
166    }
167    if (baddir)
168      fprintf(Stderr,
169              "Treating %d .direction values outside [-1, 1] as 0.\n",
170              baddir);
171  }
172  if (priority) {
173    int badpri=0;
174    saveInfo->priorities= (int *) malloc(numberColumns*sizeof(int));
175    for (i=0;i<numberColumns;i++) {
176      int value = priority[i];
177      if (value<0) {
178        badpri++;
179        value=0;
180      }
181      saveInfo->priorities[i]=value;
182    }
183    if (badpri)
184      fprintf(Stderr,
185              "Treating %d negative .priority values as 0\n",
186              badpri);
187  }
188  if (pseudoDown||pseudoUp) {
189    int badpseudo=0;
190    if (!pseudoDown||!pseudoUp)
191      fprintf(Stderr,
192              "Only one set of pseudocosts - assumed same\n");
193    saveInfo->pseudoDown= (double *) malloc(numberColumns*sizeof(double));
194    saveInfo->pseudoUp = (double *) malloc(numberColumns*sizeof(double));
195    for (i=0;i<numberColumns;i++) {
196      double valueD=0.0, valueU=0.0;
197      if (pseudoDown) {
198        valueD = pseudoDown[i];
199        if (valueD<0) {
200          badpseudo++;
201          valueD=0.0;
202        }
203      }
204      if (pseudoUp) {
205        valueU = pseudoUp[i];
206        if (valueU<0) {
207          badpseudo++;
208          valueU=0.0;
209        }
210      }
211      if (!valueD)
212        valueD=valueU;
213      if (!valueU)
214        valueU=valueD;
215      saveInfo->pseudoDown[i]=valueD;
216      saveInfo->pseudoUp[i]=valueU;
217    }
218    if (badpseudo)
219      fprintf(Stderr,
220              "Treating %d negative pseudoCosts as 0.0\n",badpseudo);
221  }
222}
223static void
224stat_map(int *stat, int n, int *map, int mx, const char *what)
225{
226  int bad, i, i1=0, j, j1=0;
227  static char badfmt[] = "Coin driver: %s[%d] = %d\n";
228 
229  for(i = bad = 0; i < n; i++) {
230    if ((j = stat[i]) >= 0 && j <= mx)
231      stat[i] = map[j];
232    else {
233      stat[i] = 0;
234      i1 = i;
235      j1 = j;
236      if (!bad++)
237        fprintf(Stderr, badfmt, what, i, j);
238    }
239  }
240  if (bad > 1) {
241    if (bad == 2)
242      fprintf(Stderr, badfmt, what, i1, j1);
243    else
244      fprintf(Stderr,
245              "Coin driver: %d messages about bad %s values suppressed.\n",
246              bad-1, what);
247  }
248}
249int
250readAmpl(ampl_info * info, int argc, char **argv)
251{
252  char *stub;
253  ograd *og;
254  int i;
255  SufDesc *csd;
256  SufDesc *rsd;
257  /*bool *basis, *lower;*/
258  /*double *LU, *c, lb, objadj, *rshift, *shift, t, ub, *x, *x0, *x1;*/
259  char * environment = getenv("cbc_options");
260  char tempBuffer[20];
261  double * obj;
262  double * columnLower;
263  double * columnUpper;
264  double * rowLower;
265  double * rowUpper;
266  char ** saveArgv=argv;
267  int saveArgc = argc;
268  memset(info,0,sizeof(ampl_info));
269  /* save so can be accessed by decodePhrase */
270  saveInfo = info;
271  info->numberArguments=0;
272  info->arguments=(char **) malloc(2*sizeof(char *));
273  info->arguments[info->numberArguments++]=strdup("ampl");
274  info->arguments[info->numberArguments++]=strdup("cbc");
275  asl = ASL_alloc(ASL_read_f);
276  stub = getstub(&argv, &Oinfo);
277  if (!stub)
278    usage_ASL(&Oinfo, 1);
279  nl = jac0dim(stub, 0);
280  suf_declare(suftab, sizeof(suftab)/sizeof(SufDecl));
281 
282  /* set A_vals to get the constraints column-wise (malloc so can be freed) */
283  A_vals = (double *) malloc(nzc*sizeof(double));
284  if (!A_vals) {
285    printf("no memory\n");
286    return 1;
287  }
288  /* say we want primal solution */
289  want_xpi0=1;
290  /* for basis info */
291  info->columnStatus = (int *) malloc(n_var*sizeof(int));
292  info->rowStatus = (int *) malloc(n_con*sizeof(int));
293  csd = suf_iput("sstatus", ASL_Sufkind_var, info->columnStatus);
294  rsd = suf_iput("sstatus", ASL_Sufkind_con, info->rowStatus);
295  /* read linear model*/
296  f_read(nl,0);
297  // see if any sos
298  if (true) {
299    char *sostype;
300    int nsosnz, *sosbeg, *sosind, * sospri;
301    double *sosref;
302    int nsos;
303    int i = ASL_suf_sos_explict_free;
304    int copri[2], **p_sospri;
305    copri[0] = 0;
306    copri[1] = 0;
307    p_sospri = &sospri;
308    nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
309                                &sosbeg, &sosind, &sosref);
310    if (nsos) {
311      info->numberSos=nsos;
312      info->sosType = (char *) malloc(nsos);
313      info->sosPriority = (int *) malloc(nsos*sizeof(int));
314      info->sosStart = (int *) malloc((nsos+1)*sizeof(int));
315      info->sosIndices = (int *) malloc(nsosnz*sizeof(int));
316      info->sosReference = (double *) malloc(nsosnz*sizeof(double));
317      sos_kludge(nsos, sosbeg, sosref);
318      for (int i=0;i<nsos;i++) {
319        int ichar = sostype[i];
320        assert (ichar=='1'||ichar=='2');
321        info->sosType[i]=ichar-'0';
322      } 
323      memcpy(info->sosPriority,sospri,nsos*sizeof(int));
324      memcpy(info->sosStart,sosbeg,(nsos+1)*sizeof(int));
325      memcpy(info->sosIndices,sosind,nsosnz*sizeof(int));
326      memcpy(info->sosReference,sosref,nsosnz*sizeof(double));
327    }
328  }
329
330  /*sos_finish(&specialOrderedInfo, 0, &j, 0, 0, 0, 0, 0);*/
331  Oinfo.uinfo = tempBuffer;
332  if (getopts(argv, &Oinfo))
333    return 1;
334  /* objective*/
335  obj = (double *) malloc(n_var*sizeof(double));
336  for (i=0;i<n_var;i++)
337    obj[i]=0.0;;
338  if (n_obj) {
339    for (og = Ograd[0];og;og = og->next)
340      obj[og->varno] = og->coef;
341  }
342  if (objtype[0])
343    info->direction=-1.0;
344  else
345    info->direction=1.0;
346  info->offset=objconst(0);
347  /* Column bounds*/
348  columnLower = (double *) malloc(n_var*sizeof(double));
349  columnUpper = (double *) malloc(n_var*sizeof(double));
350#define COIN_DBL_MAX DBL_MAX
351  for (i=0;i<n_var;i++) {
352    columnLower[i]=LUv[2*i];
353    if (columnLower[i]<= negInfinity)
354      columnLower[i]=-COIN_DBL_MAX;
355    columnUpper[i]=LUv[2*i+1];
356    if (columnUpper[i]>= Infinity)
357      columnUpper[i]=COIN_DBL_MAX;
358  }
359  /* Row bounds*/
360  rowLower = (double *) malloc(n_con*sizeof(double));
361  rowUpper = (double *) malloc(n_con*sizeof(double));
362  for (i=0;i<n_con;i++) {
363    rowLower[i]=LUrhs[2*i];
364    if (rowLower[i]<= negInfinity)
365      rowLower[i]=-COIN_DBL_MAX;
366    rowUpper[i]=LUrhs[2*i+1];
367    if (rowUpper[i]>= Infinity)
368      rowUpper[i]=COIN_DBL_MAX;
369  }
370  info->numberRows=n_con;
371  info->numberColumns=n_var;
372  info->numberElements=nzc;;
373  info->numberBinary=nbv;
374  info->numberIntegers=niv;
375  info->objective=obj;
376  info->rowLower=rowLower;
377  info->rowUpper=rowUpper;
378  info->columnLower=columnLower;
379  info->columnUpper=columnUpper;
380  info->starts=A_colstarts;
381  /*A_colstarts=NULL;*/
382  info->rows=A_rownos;
383  /*A_rownos=NULL;*/
384  info->elements=A_vals;
385  /*A_vals=NULL;*/
386  info->primalSolution=NULL;
387  /* put in primalSolution if exists */
388  if (X0) {
389    info->primalSolution=(double *) malloc(n_var*sizeof(double));
390    memcpy(info->primalSolution,X0,n_var*sizeof(double));
391  }
392  info->dualSolution=NULL;
393  if (niv+nbv>0)
394    mip_stuff(); // get any extra info
395  if ((!(niv+nbv)&&(csd->kind & ASL_Sufkind_input))
396      ||(rsd->kind & ASL_Sufkind_input)) {
397    /* convert status - need info on map */
398    static int map[] = {1, 3, 1, 1, 2, 1, 1};
399    stat_map(info->columnStatus, n_var, map, 6, "incoming columnStatus");
400    stat_map(info->rowStatus, n_con, map, 6, "incoming rowStatus");
401  } else {
402    /* all slack basis */
403    free(info->rowStatus);
404    info->rowStatus=NULL;
405    free(info->columnStatus);
406    info->columnStatus=NULL;
407  }
408  /* add -solve - unless something there already
409   - also check for sleep=yes */
410  {
411    int found=0;
412    int foundLog=0;
413    int foundSleep=0;
414    const char * something[]={"solve","branch","duals","primals","user"};
415    for (i=0;i<info->numberArguments;i++) {
416      unsigned int j;
417      const char * argument = info->arguments[i];
418      for (j=0;j<sizeof(something)/sizeof(char *);j++) {
419        const char * check = something[j];
420        if (!strncmp(argument,check,sizeof(check))) {
421          found=(int)(j+1);
422        } else if (!strncmp(argument,"log",3)) {
423          foundLog=1;
424        } else if (!strncmp(argument,"sleep",5)) {
425          foundSleep=1;
426        }
427      }
428    }
429    if (foundLog) {
430      /* print options etc */
431      for (i=0;i<saveArgc;i++)
432        printf("%s ",saveArgv[i]);
433      printf("\n");
434      if (environment)
435        printf("env %s\n",environment);
436      /*printf("%d rows %d columns %d elements\n",n_con,n_var,nzc);*/
437    }
438    if (!found) {
439      if (!strlen(algFound)) {
440        info->arguments=(char **) realloc(info->arguments,(info->numberArguments+1)*sizeof(char *));
441        info->arguments[info->numberArguments++]=strdup("-solve");
442      } else {
443        // use algorithm from keyword
444        info->arguments=(char **) realloc(info->arguments,(info->numberArguments+1)*sizeof(char *));
445        info->arguments[info->numberArguments++]=strdup(algFound);
446      }
447    }
448    if (foundSleep) {
449      /* let user copy .nl file */
450      fprintf(stderr,"You can copy .nl file %s for debug purposes or attach debugger\n",saveArgv[1]);
451      fprintf(stderr,"Type q to quit, anything else to continue\n");
452      char getChar = getc(stdin);
453      if (getChar=='q'||getChar=='Q')
454        exit(1);
455    }
456  }
457  /* add -quit */
458  info->arguments=(char **) realloc(info->arguments,(info->numberArguments+1)*sizeof(char *));
459  info->arguments[info->numberArguments++]=strdup("-quit");
460  return 0;
461}
462void freeArrays1(ampl_info * info)
463{
464  free(info->objective);
465  info->objective=NULL;
466  free(info->rowLower);
467  info->rowLower=NULL;
468  free(info->rowUpper);
469  info->rowUpper=NULL;
470  free(info->columnLower);
471  info->columnLower=NULL;
472  free(info->columnUpper);
473  info->columnUpper=NULL;
474  /* this one not freed by ASL_free */
475  free(info->elements);
476  info->elements=NULL;
477  free(info->primalSolution);
478  info->primalSolution=NULL;
479  free(info->dualSolution);
480  info->dualSolution=NULL;
481  /*free(info->rowStatus);
482  info->rowStatus=NULL;
483  free(info->columnStatus);
484  info->columnStatus=NULL;*/
485}
486void freeArrays2(ampl_info * info)
487{
488  free(info->primalSolution);
489  info->primalSolution=NULL;
490  free(info->dualSolution);
491  info->dualSolution=NULL;
492  free(info->rowStatus);
493  info->rowStatus=NULL;
494  free(info->columnStatus);
495  info->columnStatus=NULL;
496  free(info->priorities);
497  info->priorities=NULL;
498  free(info->branchDirection);
499  info->branchDirection=NULL;
500  free(info->pseudoDown);
501  info->pseudoDown=NULL;
502  free(info->pseudoUp);
503  info->pseudoUp=NULL;
504  free(info->sosType);
505  info->sosType=NULL;
506  free(info->sosPriority);
507  info->sosPriority=NULL;
508  free(info->sosStart);
509  info->sosStart=NULL;
510  free(info->sosIndices);
511  info->sosIndices=NULL;
512  free(info->sosReference);
513  info->sosReference=NULL;
514  ASL_free(&asl);
515}
516void freeArgs(ampl_info * info)
517{
518  int i;
519  for ( i=0; i<info->numberArguments;i++)
520    free(info->arguments[i]);
521  free(info->arguments);
522}
523int ampl_obj_prec()
524{
525  return obj_prec();
526}
527void writeAmpl(ampl_info * info)
528{
529  char buf[1000];
530  typedef struct { const char *msg; int code; int wantObj; } Sol_info;
531  static Sol_info solinfo[] = {
532    { "optimal solution",                       000, 1 },
533    { "infeasible",                             200, 1 },
534    { "unbounded",                              300, 0 },
535    { "iteration limit etc",                    400, 1 },
536    { "solution limit",                         401, 1 },
537    { "ran out of space",                       500, 0 },
538    { "status unknown",                         501, 1 },
539    { "bug!",                                   502, 0 },
540    { "best MIP solution so far restored",      101, 1 },
541    { "failed to restore best MIP solution",    503, 1 },
542    { "optimal (?) solution",                   100, 1 }
543  };
544  /* convert status - need info on map */
545  static int map[] = {0, 3, 4, 1};
546  sprintf(buf,"%s %s",Oinfo.bsname,info->buffer);
547  solve_result_num = solinfo[info->problemStatus].code;
548  if (info->columnStatus) {
549    stat_map(info->columnStatus, n_var, map, 4, "outgoing columnStatus");
550    stat_map(info->rowStatus, n_con, map, 4, "outgoing rowStatus");
551    suf_iput("sstatus", ASL_Sufkind_var, info->columnStatus);
552    suf_iput("sstatus", ASL_Sufkind_con, info->rowStatus);
553  }
554  write_sol(buf,info->primalSolution,info->dualSolution,&Oinfo);
555}
Note: See TracBrowser for help on using the repository browser.