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

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

for ampl and version

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