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

Last change on this file since 520 was 520, checked in by forrest, 12 years ago

fix bug in ampl interface and try slp

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 33.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 "CbcConfig.h"
27#ifdef HAVE_UNISTD_H
28# include "unistd.h"
29#endif
30#include "CoinUtilsConfig.h"
31#include "CoinHelperFunctions.hpp"
32#include "CoinModel.hpp"
33#include "CoinSort.hpp"
34#include "CoinPackedMatrix.hpp"
35#include "CoinMpsIO.hpp"
36#include "CoinFloatEqual.hpp"
37
38extern "C" {
39# include "getstub.h"
40}
41
42#include "Cbc_ampl.h"
43#include <string>
44#include <cassert>
45/* so decodePhrase and clpCheck can access */
46static ampl_info * saveInfo=NULL;
47// Set to 1 if algorithm found
48static char algFound[20]="";
49static char*
50checkPhrase(Option_Info *oi, keyword *kw, char *v)
51{
52  if (strlen(v))
53    printf("string %s\n",v);
54  // Say algorithm found
55  strcpy(algFound,kw->desc);;
56  return v;
57}
58static char*
59checkPhrase2(Option_Info *oi, keyword *kw, char *v)
60{
61  if (strlen(v))
62    printf("string %s\n",v);
63  // put out keyword
64  saveInfo->arguments=(char **) realloc(saveInfo->arguments,(saveInfo->numberArguments+1)*sizeof(char *));
65  saveInfo->arguments[saveInfo->numberArguments++]=strdup(kw->desc);
66  return v;
67}
68static fint
69decodePhrase(char * phrase,ftnlen length)
70{
71  char * blank = strchr(phrase,' ');
72  if (blank) {
73    /* split arguments */
74    *blank='\0';
75    saveInfo->arguments=(char **) realloc(saveInfo->arguments,(saveInfo->numberArguments+2)*sizeof(char *));
76    saveInfo->arguments[saveInfo->numberArguments++]=strdup(phrase);
77    *blank=' ';
78    phrase=blank+1; /* move on */
79    if (strlen(phrase))
80      saveInfo->arguments[saveInfo->numberArguments++]=strdup(phrase);
81  } else if (strlen(phrase)) {
82    saveInfo->arguments=(char **) realloc(saveInfo->arguments,(saveInfo->numberArguments+1)*sizeof(char *));
83    saveInfo->arguments[saveInfo->numberArguments++]=strdup(phrase);
84  }
85  return 0;
86}
87static void
88sos_kludge(int nsos, int *sosbeg, double *sosref,int * sosind)
89{
90  // Adjust sosref if necessary to make monotonic increasing
91  int i, j, k;
92  // first sort
93  for (i=0;i<nsos;i++) {
94    k = sosbeg[i];
95    int end=sosbeg[i+1];
96    CoinSort_2(sosref+k,sosref+end,sosind+k);
97  }
98  double t, t1;
99  for(i = j = 0; i++ < nsos; ) {
100    k = sosbeg[i];
101    t = sosref[j];
102    while(++j < k) {
103      t1 = sosref[j];
104      t += 1e-10;
105      if (t1 <= t)
106        sosref[j] = t1 = t + 1e-10;
107      t = t1;
108    }
109  }
110}
111static char xxxxxx[20];
112#define VP (char*)
113 static keyword keywds[] = { /* must be sorted */
114        { "barrier",    checkPhrase,            (char *) xxxxxx ,"-barrier" },
115        { "dual",       checkPhrase,            (char *) xxxxxx , "-dualsimplex"},
116        { "help",       checkPhrase2,           (char *) xxxxxx , "-?"},
117        { "initial",    checkPhrase,            (char *) xxxxxx , "-initialsolve"},
118        { "max",        checkPhrase2,           (char *) xxxxxx , "-maximize"},
119        { "maximize",   checkPhrase2,           (char *) xxxxxx , "-maximize"},
120        { "primal",     checkPhrase,            (char *) xxxxxx , "-primalsimplex"},
121        { "quit",       checkPhrase2,           (char *) xxxxxx , "-quit"},
122        { "wantsol",    WS_val,         NULL, "write .sol file (without -AMPL)" }
123        };
124static Option_Info Oinfo = {"cbc", "Cbc 1.01", "cbc_options", keywds, nkeywds, 0, "",
125                                0,decodePhrase,0,0,0, 20060130 };
126// strdup used to avoid g++ compiler warning
127 static SufDecl
128suftab[] = {
129#if 0
130        { "current", 0, ASL_Sufkind_con | ASL_Sufkind_outonly },
131        { "current", 0, ASL_Sufkind_var | ASL_Sufkind_outonly },
132        { "direction", 0, ASL_Sufkind_var },
133        { "down", 0, ASL_Sufkind_con | ASL_Sufkind_outonly },
134        { "down", 0, ASL_Sufkind_var | ASL_Sufkind_outonly },
135        { "priority", 0, ASL_Sufkind_var },
136#endif
137        { "direction", 0, ASL_Sufkind_var },
138        { "downPseudocost", 0, ASL_Sufkind_var | ASL_Sufkind_real },
139        { "priority", 0, ASL_Sufkind_var },
140        { "ref", 0, ASL_Sufkind_var | ASL_Sufkind_real },
141        { "sos", 0, ASL_Sufkind_var },
142        { "sos", 0, ASL_Sufkind_con },
143        { "sosno", 0, ASL_Sufkind_var | ASL_Sufkind_real },
144        { "sosref", 0, ASL_Sufkind_var | ASL_Sufkind_real },
145        { strdup("sstatus"), 0, ASL_Sufkind_var, 0 },
146        { strdup("sstatus"), 0, ASL_Sufkind_con, 0 },
147        { "upPseudocost", 0, ASL_Sufkind_var | ASL_Sufkind_real }
148#if 0
149        { "unbdd", 0, ASL_Sufkind_var | ASL_Sufkind_outonly},
150        { "up", 0, ASL_Sufkind_con | ASL_Sufkind_outonly },
151        { "up", 0, ASL_Sufkind_var | ASL_Sufkind_outonly }
152#endif
153        };
154#include "float.h"
155#include "limits.h"
156static ASL *asl=NULL;
157static FILE *nl=NULL;
158
159static void
160mip_stuff(void)
161{
162  int i;
163  double *pseudoUp, *pseudoDown;
164  int *priority, *direction;
165  SufDesc *dpup, *dpdown, *dpri, *ddir;
166 
167  ddir = suf_get("direction", ASL_Sufkind_var);
168  direction = ddir->u.i;
169  dpri = suf_get("priority", ASL_Sufkind_var);
170  priority = dpri->u.i;
171  dpdown = suf_get("downPseudocost", ASL_Sufkind_var);
172  pseudoDown = dpdown->u.r;
173  dpup = suf_get("upPseudocost", ASL_Sufkind_var);
174  pseudoUp = dpup->u.r;
175  assert(saveInfo);
176  int numberColumns = saveInfo->numberColumns;
177  if (direction) {
178    int baddir=0;
179    saveInfo->branchDirection = (int *) malloc(numberColumns*sizeof(int));
180    for (i=0;i<numberColumns;i++) {
181      int value = direction[i];
182      if (value<-1||value>1) {
183        baddir++;
184        value=0;
185      }
186      saveInfo->branchDirection[i]=value;
187    }
188    if (baddir)
189      fprintf(Stderr,
190              "Treating %d .direction values outside [-1, 1] as 0.\n",
191              baddir);
192  }
193  if (priority) {
194    int badpri=0;
195    saveInfo->priorities= (int *) malloc(numberColumns*sizeof(int));
196    for (i=0;i<numberColumns;i++) {
197      int value = priority[i];
198      if (value<0) {
199        badpri++;
200        value=0;
201      }
202      saveInfo->priorities[i]=value;
203    }
204    if (badpri)
205      fprintf(Stderr,
206              "Treating %d negative .priority values as 0\n",
207              badpri);
208  }
209  if (pseudoDown||pseudoUp) {
210    int badpseudo=0;
211    if (!pseudoDown||!pseudoUp)
212      fprintf(Stderr,
213              "Only one set of pseudocosts - assumed same\n");
214    saveInfo->pseudoDown= (double *) malloc(numberColumns*sizeof(double));
215    saveInfo->pseudoUp = (double *) malloc(numberColumns*sizeof(double));
216    for (i=0;i<numberColumns;i++) {
217      double valueD=0.0, valueU=0.0;
218      if (pseudoDown) {
219        valueD = pseudoDown[i];
220        if (valueD<0) {
221          badpseudo++;
222          valueD=0.0;
223        }
224      }
225      if (pseudoUp) {
226        valueU = pseudoUp[i];
227        if (valueU<0) {
228          badpseudo++;
229          valueU=0.0;
230        }
231      }
232      if (!valueD)
233        valueD=valueU;
234      if (!valueU)
235        valueU=valueD;
236      saveInfo->pseudoDown[i]=valueD;
237      saveInfo->pseudoUp[i]=valueU;
238    }
239    if (badpseudo)
240      fprintf(Stderr,
241              "Treating %d negative pseudoCosts as 0.0\n",badpseudo);
242  }
243}
244static void
245stat_map(int *stat, int n, int *map, int mx, const char *what)
246{
247  int bad, i, i1=0, j, j1=0;
248  static char badfmt[] = "Coin driver: %s[%d] = %d\n";
249 
250  for(i = bad = 0; i < n; i++) {
251    if ((j = stat[i]) >= 0 && j <= mx)
252      stat[i] = map[j];
253    else {
254      stat[i] = 0;
255      i1 = i;
256      j1 = j;
257      if (!bad++)
258        fprintf(Stderr, badfmt, what, i, j);
259    }
260  }
261  if (bad > 1) {
262    if (bad == 2)
263      fprintf(Stderr, badfmt, what, i1, j1);
264    else
265      fprintf(Stderr,
266              "Coin driver: %d messages about bad %s values suppressed.\n",
267              bad-1, what);
268  }
269}
270int
271readAmpl(ampl_info * info, int argc, char **argv, void ** coinModel)
272{
273  char *stub;
274  ograd *og;
275  int i;
276  SufDesc *csd;
277  SufDesc *rsd;
278  /*bool *basis, *lower;*/
279  /*double *LU, *c, lb, objadj, *rshift, *shift, t, ub, *x, *x0, *x1;*/
280  char * environment = getenv("cbc_options");
281  char tempBuffer[20];
282  double * obj;
283  double * columnLower;
284  double * columnUpper;
285  double * rowLower;
286  double * rowUpper;
287  char ** saveArgv=argv;
288  char fileName[1000];
289  if (argc>1)
290    strcpy(fileName,argv[1]);
291  else
292    fileName[0]='\0';
293  int saveArgc = argc;
294  memset(info,0,sizeof(ampl_info));
295  /* save so can be accessed by decodePhrase */
296  saveInfo = info;
297  info->numberArguments=0;
298  info->arguments=(char **) malloc(2*sizeof(char *));
299  info->arguments[info->numberArguments++]=strdup("ampl");
300  info->arguments[info->numberArguments++]=strdup("cbc");
301  asl = ASL_alloc(ASL_read_f);
302  stub = getstub(&argv, &Oinfo);
303  if (!stub)
304    usage_ASL(&Oinfo, 1);
305  nl = jac0dim(stub, 0);
306  suf_declare(suftab, sizeof(suftab)/sizeof(SufDecl));
307 
308  /* set A_vals to get the constraints column-wise (malloc so can be freed) */
309  A_vals = (double *) malloc(nzc*sizeof(double));
310  if (!A_vals) {
311    printf("no memory\n");
312    return 1;
313  }
314  /* say we want primal solution */
315  want_xpi0=1;
316  /* for basis info */
317  info->columnStatus = (int *) malloc(n_var*sizeof(int));
318  info->rowStatus = (int *) malloc(n_con*sizeof(int));
319  csd = suf_iput("sstatus", ASL_Sufkind_var, info->columnStatus);
320  rsd = suf_iput("sstatus", ASL_Sufkind_con, info->rowStatus);
321  // testosi parameter - if >= 10 then go in through coinModel
322  int testOsi=-1;
323  for (i=0;i<saveArgc;i++) {
324    if (!strncmp(saveArgv[i],"testosi",7)) {
325      testOsi = atoi(saveArgv[i+1]);
326      break;
327    }
328  }
329  if (!(nlvc+nlvo)||testOsi>=10) {
330    /* read linear model*/
331    f_read(nl,0);
332    // see if any sos
333    if (true) {
334      char *sostype;
335      int nsosnz, *sosbeg, *sosind, * sospri;
336      double *sosref;
337      int nsos;
338      int i = ASL_suf_sos_explict_free;
339      int copri[2], **p_sospri;
340      copri[0] = 0;
341      copri[1] = 0;
342      p_sospri = &sospri;
343      nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
344                     &sosbeg, &sosind, &sosref);
345      if (nsos) {
346        info->numberSos=nsos;
347        info->sosType = (char *) malloc(nsos);
348        info->sosPriority = (int *) malloc(nsos*sizeof(int));
349        info->sosStart = (int *) malloc((nsos+1)*sizeof(int));
350        info->sosIndices = (int *) malloc(nsosnz*sizeof(int));
351        info->sosReference = (double *) malloc(nsosnz*sizeof(double));
352        sos_kludge(nsos, sosbeg, sosref,sosind);
353        for (int i=0;i<nsos;i++) {
354          int ichar = sostype[i];
355          assert (ichar=='1'||ichar=='2');
356          info->sosType[i]=ichar-'0';
357        }       
358        memcpy(info->sosPriority,sospri,nsos*sizeof(int));
359        memcpy(info->sosStart,sosbeg,(nsos+1)*sizeof(int));
360        memcpy(info->sosIndices,sosind,nsosnz*sizeof(int));
361        memcpy(info->sosReference,sosref,nsosnz*sizeof(double));
362      }
363    }
364   
365    /*sos_finish(&specialOrderedInfo, 0, &j, 0, 0, 0, 0, 0);*/
366    Oinfo.uinfo = tempBuffer;
367    if (getopts(argv, &Oinfo))
368      return 1;
369    /* objective*/
370    obj = (double *) malloc(n_var*sizeof(double));
371    for (i=0;i<n_var;i++)
372      obj[i]=0.0;;
373    if (n_obj) {
374      for (og = Ograd[0];og;og = og->next)
375        obj[og->varno] = og->coef;
376    }
377    if (objtype[0])
378      info->direction=-1.0;
379    else
380      info->direction=1.0;
381    info->offset=objconst(0);
382    /* Column bounds*/
383    columnLower = (double *) malloc(n_var*sizeof(double));
384    columnUpper = (double *) malloc(n_var*sizeof(double));
385#define COIN_DBL_MAX DBL_MAX
386    for (i=0;i<n_var;i++) {
387      columnLower[i]=LUv[2*i];
388      if (columnLower[i]<= negInfinity)
389        columnLower[i]=-COIN_DBL_MAX;
390      columnUpper[i]=LUv[2*i+1];
391      if (columnUpper[i]>= Infinity)
392        columnUpper[i]=COIN_DBL_MAX;
393    }
394    /* Row bounds*/
395    rowLower = (double *) malloc(n_con*sizeof(double));
396    rowUpper = (double *) malloc(n_con*sizeof(double));
397    for (i=0;i<n_con;i++) {
398      rowLower[i]=LUrhs[2*i];
399      if (rowLower[i]<= negInfinity)
400        rowLower[i]=-COIN_DBL_MAX;
401      rowUpper[i]=LUrhs[2*i+1];
402      if (rowUpper[i]>= Infinity)
403        rowUpper[i]=COIN_DBL_MAX;
404    }
405    info->numberRows=n_con;
406    info->numberColumns=n_var;
407    info->numberElements=nzc;;
408    info->numberBinary=nbv;
409    info->numberIntegers=niv;
410    info->objective=obj;
411    info->rowLower=rowLower;
412    info->rowUpper=rowUpper;
413    info->columnLower=columnLower;
414    info->columnUpper=columnUpper;
415    info->starts=A_colstarts;
416    /*A_colstarts=NULL;*/
417    info->rows=A_rownos;
418    /*A_rownos=NULL;*/
419    info->elements=A_vals;
420    /*A_vals=NULL;*/
421    info->primalSolution=NULL;
422    /* put in primalSolution if exists */
423    if (X0) {
424      info->primalSolution=(double *) malloc(n_var*sizeof(double));
425      memcpy(info->primalSolution,X0,n_var*sizeof(double));
426    }
427    info->dualSolution=NULL;
428    if (niv+nbv>0)
429      mip_stuff(); // get any extra info
430    if ((!(niv+nbv)&&(csd->kind & ASL_Sufkind_input))
431        ||(rsd->kind & ASL_Sufkind_input)) {
432      /* convert status - need info on map */
433      static int map[] = {1, 3, 1, 1, 2, 1, 1};
434      stat_map(info->columnStatus, n_var, map, 6, "incoming columnStatus");
435      stat_map(info->rowStatus, n_con, map, 6, "incoming rowStatus");
436    } else {
437      /* all slack basis */
438      // leave status for output */
439#if 0
440      free(info->rowStatus);
441      info->rowStatus=NULL;
442      free(info->columnStatus);
443      info->columnStatus=NULL;
444#endif
445    }
446  } else {
447    // QP
448    // Add .nl if not there
449    if (!strstr(fileName,".nl"))
450      strcat(fileName,".nl");
451    CoinModel * model = new CoinModel(1,fileName);
452    if (model->numberRows()>0)
453      *coinModel=(void *) model;
454    Oinfo.uinfo = tempBuffer;
455    if (getopts(argv, &Oinfo))
456      return 1;
457    Oinfo.wantsol=1;
458    if (objtype[0])
459      info->direction=-1.0;
460    else
461      info->direction=1.0;
462    info->offset=objconst(0);
463    info->numberRows=n_con;
464    info->numberColumns=n_var;
465    info->numberElements=nzc;;
466    info->numberBinary=nbv;
467    info->numberIntegers=niv;
468    if (niv+nbv>0)
469      mip_stuff(); // get any extra info
470  }
471  /* add -solve - unless something there already
472   - also check for sleep=yes */
473  {
474    int found=0;
475    int foundLog=0;
476    int foundSleep=0;
477    const char * something[]={"solve","branch","duals","primals","user"};
478    for (i=0;i<info->numberArguments;i++) {
479      unsigned int j;
480      const char * argument = info->arguments[i];
481      for (j=0;j<sizeof(something)/sizeof(char *);j++) {
482        const char * check = something[j];
483        if (!strncmp(argument,check,sizeof(check))) {
484          found=(int)(j+1);
485        } else if (!strncmp(argument,"log",3)) {
486          foundLog=1;
487        } else if (!strncmp(argument,"sleep",5)) {
488          foundSleep=1;
489        }
490      }
491    }
492    if (foundLog) {
493      /* print options etc */
494      for (i=0;i<saveArgc;i++)
495        printf("%s ",saveArgv[i]);
496      printf("\n");
497      if (environment)
498        printf("env %s\n",environment);
499      /*printf("%d rows %d columns %d elements\n",n_con,n_var,nzc);*/
500    }
501    if (!found) {
502      if (!strlen(algFound)) {
503        info->arguments=(char **) realloc(info->arguments,(info->numberArguments+1)*sizeof(char *));
504        info->arguments[info->numberArguments++]=strdup("-solve");
505      } else {
506        // use algorithm from keyword
507        info->arguments=(char **) realloc(info->arguments,(info->numberArguments+1)*sizeof(char *));
508        info->arguments[info->numberArguments++]=strdup(algFound);
509      }
510    }
511    if (foundSleep) {
512      /* let user copy .nl file */
513      fprintf(stderr,"You can copy .nl file %s for debug purposes or attach debugger\n",saveArgv[1]);
514      fprintf(stderr,"Type q to quit, anything else to continue\n");
515      char getChar = getc(stdin);
516      if (getChar=='q'||getChar=='Q')
517        exit(1);
518    }
519  }
520  /* add -quit */
521  info->arguments=(char **) realloc(info->arguments,(info->numberArguments+1)*sizeof(char *));
522  info->arguments[info->numberArguments++]=strdup("-quit");
523  return 0;
524}
525void freeArrays1(ampl_info * info)
526{
527  free(info->objective);
528  info->objective=NULL;
529  free(info->rowLower);
530  info->rowLower=NULL;
531  free(info->rowUpper);
532  info->rowUpper=NULL;
533  free(info->columnLower);
534  info->columnLower=NULL;
535  free(info->columnUpper);
536  info->columnUpper=NULL;
537  /* this one not freed by ASL_free */
538  free(info->elements);
539  info->elements=NULL;
540  free(info->primalSolution);
541  info->primalSolution=NULL;
542  free(info->dualSolution);
543  info->dualSolution=NULL;
544  /*free(info->rowStatus);
545  info->rowStatus=NULL;
546  free(info->columnStatus);
547  info->columnStatus=NULL;*/
548}
549void freeArrays2(ampl_info * info)
550{
551  free(info->primalSolution);
552  info->primalSolution=NULL;
553  free(info->dualSolution);
554  info->dualSolution=NULL;
555  free(info->rowStatus);
556  info->rowStatus=NULL;
557  free(info->columnStatus);
558  info->columnStatus=NULL;
559  free(info->priorities);
560  info->priorities=NULL;
561  free(info->branchDirection);
562  info->branchDirection=NULL;
563  free(info->pseudoDown);
564  info->pseudoDown=NULL;
565  free(info->pseudoUp);
566  info->pseudoUp=NULL;
567  free(info->sosType);
568  info->sosType=NULL;
569  free(info->sosPriority);
570  info->sosPriority=NULL;
571  free(info->sosStart);
572  info->sosStart=NULL;
573  free(info->sosIndices);
574  info->sosIndices=NULL;
575  free(info->sosReference);
576  info->sosReference=NULL;
577  ASL_free(&asl);
578}
579void freeArgs(ampl_info * info)
580{
581  int i;
582  for ( i=0; i<info->numberArguments;i++)
583    free(info->arguments[i]);
584  free(info->arguments);
585}
586int ampl_obj_prec()
587{
588  return obj_prec();
589}
590void writeAmpl(ampl_info * info)
591{
592  char buf[1000];
593  typedef struct { const char *msg; int code; int wantObj; } Sol_info;
594  static Sol_info solinfo[] = {
595    { "optimal solution",                       000, 1 },
596    { "infeasible",                             200, 1 },
597    { "unbounded",                              300, 0 },
598    { "iteration limit etc",                    400, 1 },
599    { "solution limit",                         401, 1 },
600    { "ran out of space",                       500, 0 },
601    { "status unknown",                         501, 1 },
602    { "bug!",                                   502, 0 },
603    { "best MIP solution so far restored",      101, 1 },
604    { "failed to restore best MIP solution",    503, 1 },
605    { "optimal (?) solution",                   100, 1 }
606  };
607  /* convert status - need info on map */
608  static int map[] = {0, 3, 4, 1};
609  sprintf(buf,"%s %s",Oinfo.bsname,info->buffer);
610  solve_result_num = solinfo[info->problemStatus].code;
611  if (info->columnStatus) {
612    stat_map(info->columnStatus, n_var, map, 4, "outgoing columnStatus");
613    stat_map(info->rowStatus, n_con, map, 4, "outgoing rowStatus");
614    suf_iput("sstatus", ASL_Sufkind_var, info->columnStatus);
615    suf_iput("sstatus", ASL_Sufkind_con, info->rowStatus);
616  }
617  write_sol(buf,info->primalSolution,info->dualSolution,&Oinfo);
618}
619/* Read a problem from AMPL nl file
620 */
621CoinModel::CoinModel( int nonLinear, const char * fileName)
622 :  numberRows_(0),
623    maximumRows_(0),
624    numberColumns_(0),
625    maximumColumns_(0),
626    numberElements_(0),
627    maximumElements_(0),
628    numberQuadraticElements_(0),
629    maximumQuadraticElements_(0),
630    optimizationDirection_(1.0),
631    objectiveOffset_(0.0),
632    rowLower_(NULL),
633    rowUpper_(NULL),
634    rowType_(NULL),
635    objective_(NULL),
636    columnLower_(NULL),
637    columnUpper_(NULL),
638    integerType_(NULL),
639    columnType_(NULL),
640    start_(NULL),
641    elements_(NULL),
642    quadraticElements_(NULL),
643    sortIndices_(NULL),
644    sortElements_(NULL),
645    sortSize_(0),
646    sizeAssociated_(0),
647    associated_(NULL),
648    numberSOS_(0),
649    startSOS_(NULL),
650    memberSOS_(NULL),
651    typeSOS_(NULL),
652    prioritySOS_(NULL),
653    referenceSOS_(NULL),
654    priority_(NULL),
655    logLevel_(0),
656    type_(-1),
657    links_(0)
658{
659  problemName_ = strdup("");
660  int status=0;
661  if (!strcmp(fileName,"-")||!strcmp(fileName,"stdin")) {
662    // stdin
663  } else {
664    std::string name=fileName;
665    bool readable = fileCoinReadable(name);
666    if (!readable) {
667      std::cerr<<"Unable to open file "
668               <<fileName<<std::endl;
669      status = -1;
670    }
671  }
672  if (!status) {
673    gdb(nonLinear,fileName);
674  }
675}
676 static real
677qterm(ASL *asl, fint *colq, fint *rowq, real *delsq)
678{
679  double t, t1, *x, *x0, *xe;
680  fint *rq0, *rqe;
681 
682  t = 0.;
683  x0 = x = X0;
684  xe = x + n_var;
685  rq0 = rowq;
686  while(x < xe) {
687    t1 = *x++;
688    rqe = rq0 + *++colq;
689    while(rowq < rqe)
690      t += t1*x0[*rowq++]**delsq++;
691  }
692  return 0.5 * t;
693}
694
695void
696CoinModel::gdb( int nonLinear, const char * fileName)
697{
698  ograd *og=NULL;
699  int i;
700  SufDesc *csd=NULL;
701  SufDesc *rsd=NULL;
702  /*bool *basis, *lower;*/
703  /*double *LU, *c, lb, objadj, *rshift, *shift, t, ub, *x, *x0, *x1;*/
704  //char tempBuffer[20];
705  double * objective=NULL;
706  double * columnLower=NULL;
707  double * columnUpper=NULL;
708  double * rowLower=NULL;
709  double * rowUpper=NULL;
710  int * columnStatus=NULL;
711  int * rowStatus=NULL;
712  int numberRows=-1;
713  int numberColumns=-1;
714  int numberElements=-1;
715  int numberBinary=-1;
716  int numberIntegers=-1;
717  int numberAllNonLinearBoth=0;
718  int numberIntegerNonLinearBoth=0;
719  int numberAllNonLinearConstraints=0;
720  int numberIntegerNonLinearConstraints=0;
721  int numberAllNonLinearObjective=0;
722  int numberIntegerNonLinearObjective=0;
723  double * primalSolution=NULL;
724  double direction=1.0;
725  char * stub = strdup(fileName);
726  CoinPackedMatrix matrixByRow;
727  fint ** colqp=NULL;
728  int *z = NULL;
729  if (nonLinear==0) {
730    // linear
731    asl = ASL_alloc(ASL_read_f);
732    nl = jac0dim(stub, 0);
733    free(stub);
734    suf_declare(suftab, sizeof(suftab)/sizeof(SufDecl));
735   
736    /* set A_vals to get the constraints column-wise (malloc so can be freed) */
737    A_vals = (double *) malloc(nzc*sizeof(double));
738    if (!A_vals) {
739      printf("no memory\n");
740      return ;
741    }
742    /* say we want primal solution */
743    want_xpi0=1;
744    /* for basis info */
745    columnStatus = (int *) malloc(n_var*sizeof(int));
746    rowStatus = (int *) malloc(n_con*sizeof(int));
747    csd = suf_iput("sstatus", ASL_Sufkind_var, columnStatus);
748    rsd = suf_iput("sstatus", ASL_Sufkind_con, rowStatus);
749    /* read linear model*/
750    f_read(nl,0);
751    // see if any sos
752    if (true) {
753      char *sostype;
754      int nsosnz, *sosbeg, *sosind, * sospri;
755      double *sosref;
756      int nsos;
757      int i = ASL_suf_sos_explict_free;
758      int copri[2], **p_sospri;
759      copri[0] = 0;
760      copri[1] = 0;
761      p_sospri = &sospri;
762      nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
763                     &sosbeg, &sosind, &sosref);
764      if (nsos) {
765        abort();
766#if 0
767        info->numberSos=nsos;
768        info->sosType = (char *) malloc(nsos);
769        info->sosPriority = (int *) malloc(nsos*sizeof(int));
770        info->sosStart = (int *) malloc((nsos+1)*sizeof(int));
771        info->sosIndices = (int *) malloc(nsosnz*sizeof(int));
772        info->sosReference = (double *) malloc(nsosnz*sizeof(double));
773        sos_kludge(nsos, sosbeg, sosref,sosind);
774        for (int i=0;i<nsos;i++) {
775          int ichar = sostype[i];
776          assert (ichar=='1'||ichar=='2');
777          info->sosType[i]=ichar-'0';
778        }       
779        memcpy(info->sosPriority,sospri,nsos*sizeof(int));
780        memcpy(info->sosStart,sosbeg,(nsos+1)*sizeof(int));
781        memcpy(info->sosIndices,sosind,nsosnz*sizeof(int));
782        memcpy(info->sosReference,sosref,nsosnz*sizeof(double));
783#endif
784      }
785    }
786   
787    /*sos_finish(&specialOrderedInfo, 0, &j, 0, 0, 0, 0, 0);*/
788    //Oinfo.uinfo = tempBuffer;
789    //if (getopts(argv, &Oinfo))
790    //return 1;
791    /* objective*/
792    objective = (double *) malloc(n_var*sizeof(double));
793    for (i=0;i<n_var;i++)
794      objective[i]=0.0;;
795    if (n_obj) {
796      for (og = Ograd[0];og;og = og->next)
797        objective[og->varno] = og->coef;
798    }
799    if (objtype[0])
800      direction=-1.0;
801    else
802      direction=1.0;
803    objectiveOffset_=objconst(0);
804    /* Column bounds*/
805    columnLower = (double *) malloc(n_var*sizeof(double));
806    columnUpper = (double *) malloc(n_var*sizeof(double));
807#define COIN_DBL_MAX DBL_MAX
808    for (i=0;i<n_var;i++) {
809      columnLower[i]=LUv[2*i];
810      if (columnLower[i]<= negInfinity)
811        columnLower[i]=-COIN_DBL_MAX;
812      columnUpper[i]=LUv[2*i+1];
813      if (columnUpper[i]>= Infinity)
814        columnUpper[i]=COIN_DBL_MAX;
815    }
816    /* Row bounds*/
817    rowLower = (double *) malloc(n_con*sizeof(double));
818    rowUpper = (double *) malloc(n_con*sizeof(double));
819    for (i=0;i<n_con;i++) {
820      rowLower[i]=LUrhs[2*i];
821      if (rowLower[i]<= negInfinity)
822        rowLower[i]=-COIN_DBL_MAX;
823      rowUpper[i]=LUrhs[2*i+1];
824      if (rowUpper[i]>= Infinity)
825        rowUpper[i]=COIN_DBL_MAX;
826    }
827    numberRows=n_con;
828    numberColumns=n_var;
829    numberElements=nzc;;
830    numberBinary=nbv;
831    numberIntegers=niv;
832    /* put in primalSolution if exists */
833    if (X0) {
834      primalSolution=(double *) malloc(n_var*sizeof(double));
835      memcpy( primalSolution,X0,n_var*sizeof(double));
836    }
837    //double * dualSolution=NULL;
838    if (niv+nbv>0)
839      mip_stuff(); // get any extra info
840    if ((!(niv+nbv)&&(csd->kind & ASL_Sufkind_input))
841        ||(rsd->kind & ASL_Sufkind_input)) {
842      /* convert status - need info on map */
843      static int map[] = {1, 3, 1, 1, 2, 1, 1};
844      stat_map(columnStatus, n_var, map, 6, "incoming columnStatus");
845      stat_map(rowStatus, n_con, map, 6, "incoming rowStatus");
846    } else {
847      /* all slack basis */
848      // leave status for output */
849#if 0
850      free(rowStatus);
851      rowStatus=NULL;
852      free(columnStatus);
853      columnStatus=NULL;
854#endif
855    }
856    CoinPackedMatrix columnCopy(true,numberRows,numberColumns,numberElements,
857                                A_vals,A_rownos,A_colstarts,NULL);
858    matrixByRow.reverseOrderedCopyOf(columnCopy);
859  } else if (nonLinear==1) {
860    // quadratic
861    asl = ASL_alloc(ASL_read_fg);
862    nl = jac0dim(stub, (long) strlen(stub));
863    free(stub);
864    suf_declare(suftab, sizeof(suftab)/sizeof(SufDecl));
865    /* read  model*/
866    X0 = (double*) malloc(n_var*sizeof(double));
867    CoinZeroN(X0,n_var);
868    qp_read(nl,0);
869    assert (n_obj==1);
870    int nz = 1 + n_con;
871    colqp = (fint**) malloc(nz*(2*sizeof(int*)
872                                      + sizeof(double*)));
873    fint ** rowqp = colqp + nz;
874    double ** delsqp = (double **)(rowqp + nz);
875    z = (int*) malloc(nz*sizeof(int));
876    for (i=0;i<=n_con;i++) {
877      z[i] = nqpcheck(-i, rowqp+i, colqp+i, delsqp+i);
878    }
879    qp_opify();
880    /* objective*/
881    objective = (double *) malloc(n_var*sizeof(double));
882    for (i=0;i<n_var;i++)
883      objective[i]=0.0;;
884    if (n_obj) {
885      for (og = Ograd[0];og;og = og->next)
886        objective[og->varno] = og->coef;
887    }
888    if (objtype[0])
889      direction=-1.0;
890    else
891      direction=1.0;
892    objectiveOffset_=objconst(0);
893    /* Column bounds*/
894    columnLower = (double *) malloc(n_var*sizeof(double));
895    columnUpper = (double *) malloc(n_var*sizeof(double));
896    for (i=0;i<n_var;i++) {
897      columnLower[i]=LUv[2*i];
898      if (columnLower[i]<= negInfinity)
899        columnLower[i]=-COIN_DBL_MAX;
900      columnUpper[i]=LUv[2*i+1];
901      if (columnUpper[i]>= Infinity)
902        columnUpper[i]=COIN_DBL_MAX;
903    }
904    // Build by row from scratch
905    //matrixByRow.reserve(n_var,nzc,true);
906    // say row orderded
907    matrixByRow.transpose();
908    /* Row bounds*/
909    rowLower = (double *) malloc(n_con*sizeof(double));
910    rowUpper = (double *) malloc(n_con*sizeof(double));
911    int * column = new int [n_var];
912    double * element = new double [n_var];
913    for (i=0;i<n_con;i++) {
914      rowLower[i]=LUrhs[2*i];
915      if (rowLower[i]<= negInfinity)
916        rowLower[i]=-COIN_DBL_MAX;
917      rowUpper[i]=LUrhs[2*i+1];
918      if (rowUpper[i]>= Infinity)
919        rowUpper[i]=COIN_DBL_MAX;
920      int k=0;
921      for(cgrad * cg = Cgrad[i]; cg; cg = cg->next) {
922        column[k]=cg->varno;
923        element[k++]=cg->coef;
924      }
925      matrixByRow.appendRow(k,column,element);
926    }
927    numberRows=n_con;
928    numberColumns=n_var;
929    numberElements=nzc;;
930    numberBinary=nbv;
931    numberIntegers=niv;
932    numberAllNonLinearBoth=nlvb;
933    numberIntegerNonLinearBoth=nlvbi;
934    numberAllNonLinearConstraints=nlvc;
935    numberIntegerNonLinearConstraints=nlvci;
936    numberAllNonLinearObjective=nlvo;
937    numberIntegerNonLinearObjective=nlvoi;
938    /* say we want primal solution */
939    want_xpi0=1;
940    //double * dualSolution=NULL;
941  } else {
942    abort();
943  }
944  // set problem name
945  free(problemName_);
946  problemName_=strdup("???");
947
948  // Build by row from scratch
949  const double * element = matrixByRow.getElements();
950  const int * column = matrixByRow.getIndices();
951  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
952  const int * rowLength = matrixByRow.getVectorLengths();
953  for (i=0;i<numberRows;i++) {
954    addRow(rowLength[i],column+rowStart[i],
955           element+rowStart[i],rowLower[i],rowUpper[i]);
956  }
957  // Now do column part
958  for (i=0;i<numberColumns;i++) {
959    setColumnBounds(i,columnLower[i],columnUpper[i]);
960    setColumnObjective(i,objective[i]);
961  }
962  for ( i=numberColumns-numberBinary-numberIntegers;
963        i<numberColumns;i++) {
964      setColumnIsInteger(i,true);;
965  }
966  // and non linear
967  for (i=numberAllNonLinearBoth-numberIntegerNonLinearBoth;
968       i<numberAllNonLinearBoth;i++) {
969    setColumnIsInteger(i,true);;
970  }
971  for (i=numberAllNonLinearConstraints-numberIntegerNonLinearConstraints;
972       i<numberAllNonLinearConstraints;i++) {
973    setColumnIsInteger(i,true);;
974  }
975  for (i=numberAllNonLinearObjective-numberIntegerNonLinearObjective;
976       i<numberAllNonLinearObjective;i++) {
977    setColumnIsInteger(i,true);;
978  }
979  // do names
980  int iRow;
981  for (iRow=0;iRow<numberRows_;iRow++) {
982    char name[9];
983    sprintf(name,"r%7.7d",iRow);
984    setRowName(iRow,name);
985  }
986   
987  int iColumn;
988  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
989    char name[9];
990    sprintf(name,"c%7.7d",iColumn);
991    setColumnName(iColumn,name);
992  }
993  if (colqp) {
994    // add in quadratic
995    int nz = 1 + n_con;
996    int nOdd=0;
997    fint ** rowqp = colqp + nz;
998    double ** delsqp = (double **)(rowqp + nz);
999    for (i=0;i<=n_con;i++) {
1000      int nels = z[i];
1001      if (nels) {
1002        double * element = delsqp[i];
1003        int * start = (int *) colqp[i];
1004        int * row = (int *) rowqp[i];
1005        if (!element) {
1006          // odd row - probably not quadratic
1007          nOdd++;
1008          continue;
1009        }
1010#if 0
1011        printf("%d quadratic els\n",nels);
1012        for (int j=0;j<n_var;j++) {
1013          for (int k=start[j];k<start[j+1];k++)
1014            printf("%d %d %g\n",j,row[k],element[k]);
1015        }
1016#endif
1017        if (i) {
1018          int iRow = i-1;
1019          for (int j=0;j<n_var;j++) {
1020            for (int k=start[j];k<start[j+1];k++) {
1021              int kColumn = row[k];
1022              double value = element[k];
1023              // ampl gives twice with assumed 0.5
1024              if (kColumn>j)
1025                continue;
1026              else if (kColumn==j)
1027                value *= 0.5;
1028              const char * expr = getElementAsString(iRow,j);
1029              double constant=0.0;
1030              bool linear;
1031              if (expr&&strcmp(expr,"Numeric")) {
1032                linear=false;
1033              } else {
1034                constant = getElement(iRow,j);
1035                linear=true;
1036              }
1037              char temp[1000];
1038              char temp2[30];
1039              if (value==1.0) 
1040                sprintf(temp2,"c%7.7d",kColumn);
1041              else
1042                sprintf(temp2,"%g*c%7.7d",value,kColumn);
1043              if (linear) {
1044                if (!constant)
1045                  strcpy(temp,temp2);
1046                else if (value>0.0) 
1047                  sprintf(temp,"%g+%s",constant,temp2);
1048                else
1049                  sprintf(temp,"%g%s",constant,temp2);
1050              } else {
1051                if (value>0.0) 
1052                  sprintf(temp,"%s+%s",expr,temp2);
1053                else
1054                  sprintf(temp,"%s%s",expr,temp2);
1055              }
1056              assert (strlen(temp)<1000);
1057              setElement(iRow,j,temp);
1058              printf("el for row %d column c%7.7d is %s\n",iRow,j,temp);
1059            }
1060          }
1061        } else {
1062          // objective
1063          for (int j=0;j<n_var;j++) {
1064            for (int k=start[j];k<start[j+1];k++) {
1065              int kColumn = row[k];
1066              double value = element[k];
1067              // ampl gives twice with assumed 0.5
1068              if (kColumn>j)
1069                continue;
1070              else if (kColumn==j)
1071                value *= 0.5;
1072              const char * expr = getColumnObjectiveAsString(j);
1073              double constant=0.0;
1074              bool linear;
1075              if (expr&&strcmp(expr,"Numeric")) {
1076                linear=false;
1077              } else {
1078                constant = getColumnObjective(j);
1079                linear=true;
1080              }
1081              char temp[1000];
1082              char temp2[30];
1083              if (value==1.0) 
1084                sprintf(temp2,"c%7.7d",kColumn);
1085              else
1086                sprintf(temp2,"%g*c%7.7d",value,kColumn);
1087              if (linear) {
1088                if (!constant)
1089                  strcpy(temp,temp2);
1090                else if (value>0.0) 
1091                  sprintf(temp,"%g+%s",constant,temp2);
1092                else
1093                  sprintf(temp,"%g%s",constant,temp2);
1094              } else {
1095                if (value>0.0) 
1096                  sprintf(temp,"%s+%s",expr,temp2);
1097                else
1098                  sprintf(temp,"%s%s",expr,temp2);
1099              }
1100              assert (strlen(temp)<1000);
1101              setObjective(j,temp);
1102              printf("el for objective column c%7.7d is %s\n",j,temp);
1103            }
1104          }
1105        }
1106      }
1107    }
1108    if (nOdd) {
1109      printf("%d non-linear constraints could not be converted to quadratic\n",nOdd);
1110      exit(77);
1111    }
1112  }
1113  // see if any sos
1114  {
1115    char *sostype;
1116    int nsosnz, *sosbeg, *sosind, * sospri;
1117    double *sosref;
1118    int nsos;
1119    int i = ASL_suf_sos_explict_free;
1120    int copri[2], **p_sospri;
1121    copri[0] = 0;
1122    copri[1] = 0;
1123    p_sospri = &sospri;
1124    nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
1125                   &sosbeg, &sosind, &sosref);
1126    if (nsos) {
1127      numberSOS_=nsos;
1128      typeSOS_ = new int [numberSOS_];
1129      prioritySOS_ = new int [numberSOS_];
1130      startSOS_ = new int [numberSOS_+1];
1131      memberSOS_ = new int[nsosnz];
1132      referenceSOS_ = new double [nsosnz];
1133      sos_kludge(nsos, sosbeg, sosref,sosind);
1134      for (int i=0;i<nsos;i++) {
1135        int ichar = sostype[i];
1136        assert (ichar=='1'||ichar=='2');
1137        typeSOS_[i]=ichar-'0';
1138      } 
1139      memcpy(prioritySOS_,sospri,nsos*sizeof(int));
1140      memcpy(startSOS_,sosbeg,(nsos+1)*sizeof(int));
1141      memcpy(memberSOS_,sosind,nsosnz*sizeof(int));
1142      memcpy(referenceSOS_,sosref,nsosnz*sizeof(double));
1143    }
1144  }
1145}
Note: See TracBrowser for help on using the repository browser.