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

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

to compile on cygwin

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