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

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

avoid segmentation fault

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