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

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

sorting

  • 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    free(info->rowStatus);
412    info->rowStatus=NULL;
413    free(info->columnStatus);
414    info->columnStatus=NULL;
415  }
416  /* add -solve - unless something there already
417   - also check for sleep=yes */
418  {
419    int found=0;
420    int foundLog=0;
421    int foundSleep=0;
422    const char * something[]={"solve","branch","duals","primals","user"};
423    for (i=0;i<info->numberArguments;i++) {
424      unsigned int j;
425      const char * argument = info->arguments[i];
426      for (j=0;j<sizeof(something)/sizeof(char *);j++) {
427        const char * check = something[j];
428        if (!strncmp(argument,check,sizeof(check))) {
429          found=(int)(j+1);
430        } else if (!strncmp(argument,"log",3)) {
431          foundLog=1;
432        } else if (!strncmp(argument,"sleep",5)) {
433          foundSleep=1;
434        }
435      }
436    }
437    if (foundLog) {
438      /* print options etc */
439      for (i=0;i<saveArgc;i++)
440        printf("%s ",saveArgv[i]);
441      printf("\n");
442      if (environment)
443        printf("env %s\n",environment);
444      /*printf("%d rows %d columns %d elements\n",n_con,n_var,nzc);*/
445    }
446    if (!found) {
447      if (!strlen(algFound)) {
448        info->arguments=(char **) realloc(info->arguments,(info->numberArguments+1)*sizeof(char *));
449        info->arguments[info->numberArguments++]=strdup("-solve");
450      } else {
451        // use algorithm from keyword
452        info->arguments=(char **) realloc(info->arguments,(info->numberArguments+1)*sizeof(char *));
453        info->arguments[info->numberArguments++]=strdup(algFound);
454      }
455    }
456    if (foundSleep) {
457      /* let user copy .nl file */
458      fprintf(stderr,"You can copy .nl file %s for debug purposes or attach debugger\n",saveArgv[1]);
459      fprintf(stderr,"Type q to quit, anything else to continue\n");
460      char getChar = getc(stdin);
461      if (getChar=='q'||getChar=='Q')
462        exit(1);
463    }
464  }
465  /* add -quit */
466  info->arguments=(char **) realloc(info->arguments,(info->numberArguments+1)*sizeof(char *));
467  info->arguments[info->numberArguments++]=strdup("-quit");
468  return 0;
469}
470void freeArrays1(ampl_info * info)
471{
472  free(info->objective);
473  info->objective=NULL;
474  free(info->rowLower);
475  info->rowLower=NULL;
476  free(info->rowUpper);
477  info->rowUpper=NULL;
478  free(info->columnLower);
479  info->columnLower=NULL;
480  free(info->columnUpper);
481  info->columnUpper=NULL;
482  /* this one not freed by ASL_free */
483  free(info->elements);
484  info->elements=NULL;
485  free(info->primalSolution);
486  info->primalSolution=NULL;
487  free(info->dualSolution);
488  info->dualSolution=NULL;
489  /*free(info->rowStatus);
490  info->rowStatus=NULL;
491  free(info->columnStatus);
492  info->columnStatus=NULL;*/
493}
494void freeArrays2(ampl_info * info)
495{
496  free(info->primalSolution);
497  info->primalSolution=NULL;
498  free(info->dualSolution);
499  info->dualSolution=NULL;
500  free(info->rowStatus);
501  info->rowStatus=NULL;
502  free(info->columnStatus);
503  info->columnStatus=NULL;
504  free(info->priorities);
505  info->priorities=NULL;
506  free(info->branchDirection);
507  info->branchDirection=NULL;
508  free(info->pseudoDown);
509  info->pseudoDown=NULL;
510  free(info->pseudoUp);
511  info->pseudoUp=NULL;
512  free(info->sosType);
513  info->sosType=NULL;
514  free(info->sosPriority);
515  info->sosPriority=NULL;
516  free(info->sosStart);
517  info->sosStart=NULL;
518  free(info->sosIndices);
519  info->sosIndices=NULL;
520  free(info->sosReference);
521  info->sosReference=NULL;
522  ASL_free(&asl);
523}
524void freeArgs(ampl_info * info)
525{
526  int i;
527  for ( i=0; i<info->numberArguments;i++)
528    free(info->arguments[i]);
529  free(info->arguments);
530}
531int ampl_obj_prec()
532{
533  return obj_prec();
534}
535void writeAmpl(ampl_info * info)
536{
537  char buf[1000];
538  typedef struct { const char *msg; int code; int wantObj; } Sol_info;
539  static Sol_info solinfo[] = {
540    { "optimal solution",                       000, 1 },
541    { "infeasible",                             200, 1 },
542    { "unbounded",                              300, 0 },
543    { "iteration limit etc",                    400, 1 },
544    { "solution limit",                         401, 1 },
545    { "ran out of space",                       500, 0 },
546    { "status unknown",                         501, 1 },
547    { "bug!",                                   502, 0 },
548    { "best MIP solution so far restored",      101, 1 },
549    { "failed to restore best MIP solution",    503, 1 },
550    { "optimal (?) solution",                   100, 1 }
551  };
552  /* convert status - need info on map */
553  static int map[] = {0, 3, 4, 1};
554  sprintf(buf,"%s %s",Oinfo.bsname,info->buffer);
555  solve_result_num = solinfo[info->problemStatus].code;
556  if (info->columnStatus) {
557    stat_map(info->columnStatus, n_var, map, 4, "outgoing columnStatus");
558    stat_map(info->rowStatus, n_con, map, 4, "outgoing rowStatus");
559    suf_iput("sstatus", ASL_Sufkind_var, info->columnStatus);
560    suf_iput("sstatus", ASL_Sufkind_con, info->rowStatus);
561  }
562  write_sol(buf,info->primalSolution,info->dualSolution,&Oinfo);
563}
Note: See TracBrowser for help on using the repository browser.