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

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

priorities

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 32.1 KB
Line 
1/****************************************************************
2Copyright (C) 1997-2000 Lucent Technologies
3Modifications for Coin -  Copyright (C) 2006, International Business Machines Corporation and others.
4All Rights Reserved
5
6Permission to use, copy, modify, and distribute this software and
7its documentation for any purpose and without fee is hereby
8granted, provided that the above copyright notice appear in all
9copies and that both that the copyright notice and this
10permission notice and warranty disclaimer appear in supporting
11documentation, and that the name of Lucent or any of its entities
12not be used in advertising or publicity pertaining to
13distribution of the software without specific, written prior
14permission.
15
16LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
17INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
18IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
19SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
20WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
21IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
22ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
23THIS SOFTWARE.
24****************************************************************/
25
26#include "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  double * primalSolution=NULL;
718  double direction=1.0;
719  char * stub = strdup(fileName);
720  CoinPackedMatrix matrixByRow;
721  fint ** colqp=NULL;
722  int *z = NULL;
723  if (nonLinear==0) {
724    // linear
725    asl = ASL_alloc(ASL_read_f);
726    nl = jac0dim(stub, 0);
727    free(stub);
728    suf_declare(suftab, sizeof(suftab)/sizeof(SufDecl));
729   
730    /* set A_vals to get the constraints column-wise (malloc so can be freed) */
731    A_vals = (double *) malloc(nzc*sizeof(double));
732    if (!A_vals) {
733      printf("no memory\n");
734      return ;
735    }
736    /* say we want primal solution */
737    want_xpi0=1;
738    /* for basis info */
739    columnStatus = (int *) malloc(n_var*sizeof(int));
740    rowStatus = (int *) malloc(n_con*sizeof(int));
741    csd = suf_iput("sstatus", ASL_Sufkind_var, columnStatus);
742    rsd = suf_iput("sstatus", ASL_Sufkind_con, rowStatus);
743    /* read linear model*/
744    f_read(nl,0);
745    // see if any sos
746    if (true) {
747      char *sostype;
748      int nsosnz, *sosbeg, *sosind, * sospri;
749      double *sosref;
750      int nsos;
751      int i = ASL_suf_sos_explict_free;
752      int copri[2], **p_sospri;
753      copri[0] = 0;
754      copri[1] = 0;
755      p_sospri = &sospri;
756      nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
757                     &sosbeg, &sosind, &sosref);
758      if (nsos) {
759        abort();
760#if 0
761        info->numberSos=nsos;
762        info->sosType = (char *) malloc(nsos);
763        info->sosPriority = (int *) malloc(nsos*sizeof(int));
764        info->sosStart = (int *) malloc((nsos+1)*sizeof(int));
765        info->sosIndices = (int *) malloc(nsosnz*sizeof(int));
766        info->sosReference = (double *) malloc(nsosnz*sizeof(double));
767        sos_kludge(nsos, sosbeg, sosref,sosind);
768        for (int i=0;i<nsos;i++) {
769          int ichar = sostype[i];
770          assert (ichar=='1'||ichar=='2');
771          info->sosType[i]=ichar-'0';
772        }       
773        memcpy(info->sosPriority,sospri,nsos*sizeof(int));
774        memcpy(info->sosStart,sosbeg,(nsos+1)*sizeof(int));
775        memcpy(info->sosIndices,sosind,nsosnz*sizeof(int));
776        memcpy(info->sosReference,sosref,nsosnz*sizeof(double));
777#endif
778      }
779    }
780   
781    /*sos_finish(&specialOrderedInfo, 0, &j, 0, 0, 0, 0, 0);*/
782    //Oinfo.uinfo = tempBuffer;
783    //if (getopts(argv, &Oinfo))
784    //return 1;
785    /* objective*/
786    objective = (double *) malloc(n_var*sizeof(double));
787    for (i=0;i<n_var;i++)
788      objective[i]=0.0;;
789    if (n_obj) {
790      for (og = Ograd[0];og;og = og->next)
791        objective[og->varno] = og->coef;
792    }
793    if (objtype[0])
794      direction=-1.0;
795    else
796      direction=1.0;
797    objectiveOffset_=objconst(0);
798    /* Column bounds*/
799    columnLower = (double *) malloc(n_var*sizeof(double));
800    columnUpper = (double *) malloc(n_var*sizeof(double));
801#define COIN_DBL_MAX DBL_MAX
802    for (i=0;i<n_var;i++) {
803      columnLower[i]=LUv[2*i];
804      if (columnLower[i]<= negInfinity)
805        columnLower[i]=-COIN_DBL_MAX;
806      columnUpper[i]=LUv[2*i+1];
807      if (columnUpper[i]>= Infinity)
808        columnUpper[i]=COIN_DBL_MAX;
809    }
810    /* Row bounds*/
811    rowLower = (double *) malloc(n_con*sizeof(double));
812    rowUpper = (double *) malloc(n_con*sizeof(double));
813    for (i=0;i<n_con;i++) {
814      rowLower[i]=LUrhs[2*i];
815      if (rowLower[i]<= negInfinity)
816        rowLower[i]=-COIN_DBL_MAX;
817      rowUpper[i]=LUrhs[2*i+1];
818      if (rowUpper[i]>= Infinity)
819        rowUpper[i]=COIN_DBL_MAX;
820    }
821    numberRows=n_con;
822    numberColumns=n_var;
823    numberElements=nzc;;
824    numberBinary=nbv;
825    numberIntegers=niv;
826    /* put in primalSolution if exists */
827    if (X0) {
828      primalSolution=(double *) malloc(n_var*sizeof(double));
829      memcpy( primalSolution,X0,n_var*sizeof(double));
830    }
831    //double * dualSolution=NULL;
832    if (niv+nbv>0)
833      mip_stuff(); // get any extra info
834    if ((!(niv+nbv)&&(csd->kind & ASL_Sufkind_input))
835        ||(rsd->kind & ASL_Sufkind_input)) {
836      /* convert status - need info on map */
837      static int map[] = {1, 3, 1, 1, 2, 1, 1};
838      stat_map(columnStatus, n_var, map, 6, "incoming columnStatus");
839      stat_map(rowStatus, n_con, map, 6, "incoming rowStatus");
840    } else {
841      /* all slack basis */
842      // leave status for output */
843#if 0
844      free(rowStatus);
845      rowStatus=NULL;
846      free(columnStatus);
847      columnStatus=NULL;
848#endif
849    }
850    CoinPackedMatrix columnCopy(true,numberRows,numberColumns,numberElements,
851                                A_vals,A_rownos,A_colstarts,NULL);
852    matrixByRow.reverseOrderedCopyOf(columnCopy);
853  } else if (nonLinear==1) {
854    // quadratic
855    asl = ASL_alloc(ASL_read_fg);
856    nl = jac0dim(stub, (long) strlen(stub));
857    free(stub);
858    suf_declare(suftab, sizeof(suftab)/sizeof(SufDecl));
859    /* read  model*/
860    X0 = (double*) malloc(n_var*sizeof(double));
861    CoinZeroN(X0,n_var);
862    qp_read(nl,0);
863    assert (n_obj==1);
864    int nz = 1 + n_con;
865    colqp = (fint**) malloc(nz*(2*sizeof(int*)
866                                      + sizeof(double*)));
867    fint ** rowqp = colqp + nz;
868    double ** delsqp = (double **)(rowqp + nz);
869    z = (int*) malloc(nz*sizeof(int));
870    for (i=0;i<=n_con;i++) {
871      z[i] = nqpcheck(-i, rowqp+i, colqp+i, delsqp+i);
872    }
873    qp_opify();
874    /* objective*/
875    objective = (double *) malloc(n_var*sizeof(double));
876    for (i=0;i<n_var;i++)
877      objective[i]=0.0;;
878    if (n_obj) {
879      for (og = Ograd[0];og;og = og->next)
880        objective[og->varno] = og->coef;
881    }
882    if (objtype[0])
883      direction=-1.0;
884    else
885      direction=1.0;
886    objectiveOffset_=objconst(0);
887    /* Column bounds*/
888    columnLower = (double *) malloc(n_var*sizeof(double));
889    columnUpper = (double *) malloc(n_var*sizeof(double));
890    for (i=0;i<n_var;i++) {
891      columnLower[i]=LUv[2*i];
892      if (columnLower[i]<= negInfinity)
893        columnLower[i]=-COIN_DBL_MAX;
894      columnUpper[i]=LUv[2*i+1];
895      if (columnUpper[i]>= Infinity)
896        columnUpper[i]=COIN_DBL_MAX;
897    }
898    // Build by row from scratch
899    //matrixByRow.reserve(n_var,nzc,true);
900    // say row orderded
901    matrixByRow.transpose();
902    /* Row bounds*/
903    rowLower = (double *) malloc(n_con*sizeof(double));
904    rowUpper = (double *) malloc(n_con*sizeof(double));
905    int * column = new int [n_var];
906    double * element = new double [n_var];
907    for (i=0;i<n_con;i++) {
908      rowLower[i]=LUrhs[2*i];
909      if (rowLower[i]<= negInfinity)
910        rowLower[i]=-COIN_DBL_MAX;
911      rowUpper[i]=LUrhs[2*i+1];
912      if (rowUpper[i]>= Infinity)
913        rowUpper[i]=COIN_DBL_MAX;
914      int k=0;
915      for(cgrad * cg = Cgrad[i]; cg; cg = cg->next) {
916        column[k]=cg->varno;
917        element[k++]=cg->coef;
918      }
919      matrixByRow.appendRow(k,column,element);
920    }
921    numberRows=n_con;
922    numberColumns=n_var;
923    numberElements=nzc;;
924    numberBinary=nbv;
925    numberIntegers=niv+nlvci;
926    /* say we want primal solution */
927    want_xpi0=1;
928    //double * dualSolution=NULL;
929  } else {
930    abort();
931  }
932  // set problem name
933  free(problemName_);
934  problemName_=strdup("???");
935
936  // Build by row from scratch
937  const double * element = matrixByRow.getElements();
938  const int * column = matrixByRow.getIndices();
939  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
940  const int * rowLength = matrixByRow.getVectorLengths();
941  for (i=0;i<numberRows;i++) {
942    addRow(rowLength[i],column+rowStart[i],
943           element+rowStart[i],rowLower[i],rowUpper[i]);
944  }
945  // Now do column part
946  for (i=0;i<numberColumns;i++) {
947    setColumnBounds(i,columnLower[i],columnUpper[i]);
948    setColumnObjective(i,objective[i]);
949  }
950  for ( i=numberColumns-numberBinary-numberIntegers;
951        i<numberColumns;i++) {
952      setColumnIsInteger(i,true);;
953  }
954  // do names
955  int iRow;
956  for (iRow=0;iRow<numberRows_;iRow++) {
957    char name[9];
958    sprintf(name,"r%7.7d",iRow);
959    setRowName(iRow,name);
960  }
961   
962  int iColumn;
963  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
964    char name[9];
965    sprintf(name,"c%7.7d",iColumn);
966    setColumnName(iColumn,name);
967  }
968  if (colqp) {
969    // add in quadratic
970    int nz = 1 + n_con;
971    int nOdd=0;
972    fint ** rowqp = colqp + nz;
973    double ** delsqp = (double **)(rowqp + nz);
974    for (i=0;i<=n_con;i++) {
975      int nels = z[i];
976      if (nels) {
977        double * element = delsqp[i];
978        int * start = (int *) colqp[i];
979        int * row = (int *) rowqp[i];
980        if (!element) {
981          // odd row - probably not quadratic
982          nOdd++;
983          continue;
984        }
985#if 0
986        printf("%d quadratic els\n",nels);
987        for (int j=0;j<n_var;j++) {
988          for (int k=start[j];k<start[j+1];k++)
989            printf("%d %d %g\n",j,row[k],element[k]);
990        }
991#endif
992        if (i) {
993          int iRow = i-1;
994          for (int j=0;j<n_var;j++) {
995            for (int k=start[j];k<start[j+1];k++) {
996              int kColumn = row[k];
997              double value = element[k];
998              // ampl gives twice with assumed 0.5
999              if (kColumn>j)
1000                continue;
1001              else if (kColumn==j)
1002                value *= 0.5;
1003              const char * expr = getElementAsString(iRow,j);
1004              double constant=0.0;
1005              bool linear;
1006              if (expr&&strcmp(expr,"Numeric")) {
1007                linear=false;
1008              } else {
1009                constant = getElement(iRow,j);
1010                linear=true;
1011              }
1012              char temp[1000];
1013              char temp2[30];
1014              if (value==1.0) 
1015                sprintf(temp2,"c%7.7d",kColumn);
1016              else
1017                sprintf(temp2,"%g*c%7.7d",value,kColumn);
1018              if (linear) {
1019                if (!constant)
1020                  strcpy(temp,temp2);
1021                else if (value>0.0) 
1022                  sprintf(temp,"%g+%s",constant,temp2);
1023                else
1024                  sprintf(temp,"%g%s",constant,temp2);
1025              } else {
1026                if (value>0.0) 
1027                  sprintf(temp,"%s+%s",expr,temp2);
1028                else
1029                  sprintf(temp,"%s%s",expr,temp2);
1030              }
1031              assert (strlen(temp)<1000);
1032              setElement(iRow,j,temp);
1033              printf("el for row %d column c%7.7d is %s\n",iRow,j,temp);
1034            }
1035          }
1036        } else {
1037          // objective
1038          for (int j=0;j<n_var;j++) {
1039            for (int k=start[j];k<start[j+1];k++) {
1040              int kColumn = row[k];
1041              double value = element[k];
1042              // ampl gives twice with assumed 0.5
1043              if (kColumn>j)
1044                continue;
1045              else if (kColumn==j)
1046                value *= 0.5;
1047              const char * expr = getColumnObjectiveAsString(j);
1048              double constant=0.0;
1049              bool linear;
1050              if (expr&&strcmp(expr,"Numeric")) {
1051                linear=false;
1052              } else {
1053                constant = getColumnObjective(j);
1054                linear=true;
1055              }
1056              char temp[1000];
1057              char temp2[30];
1058              if (value==1.0) 
1059                sprintf(temp2,"c%7.7d",kColumn);
1060              else
1061                sprintf(temp2,"%g*c%7.7d",value,kColumn);
1062              if (linear) {
1063                if (!constant)
1064                  strcpy(temp,temp2);
1065                else if (value>0.0) 
1066                  sprintf(temp,"%g+%s",constant,temp2);
1067                else
1068                  sprintf(temp,"%g%s",constant,temp2);
1069              } else {
1070                if (value>0.0) 
1071                  sprintf(temp,"%s+%s",expr,temp2);
1072                else
1073                  sprintf(temp,"%s%s",expr,temp2);
1074              }
1075              assert (strlen(temp)<1000);
1076              setObjective(j,temp);
1077              printf("el for objective column c%7.7d is %s\n",j,temp);
1078            }
1079          }
1080        }
1081      }
1082    }
1083    if (nOdd) {
1084      printf("%d non-linear constraints could not be converted to quadratic\n",nOdd);
1085      exit(77);
1086    }
1087  }
1088  // see if any sos
1089  {
1090    char *sostype;
1091    int nsosnz, *sosbeg, *sosind, * sospri;
1092    double *sosref;
1093    int nsos;
1094    int i = ASL_suf_sos_explict_free;
1095    int copri[2], **p_sospri;
1096    copri[0] = 0;
1097    copri[1] = 0;
1098    p_sospri = &sospri;
1099    nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
1100                   &sosbeg, &sosind, &sosref);
1101    if (nsos) {
1102      numberSOS_=nsos;
1103      typeSOS_ = new int [numberSOS_];
1104      prioritySOS_ = new int [numberSOS_];
1105      startSOS_ = new int [numberSOS_+1];
1106      memberSOS_ = new int[nsosnz];
1107      referenceSOS_ = new double [nsosnz];
1108      sos_kludge(nsos, sosbeg, sosref,sosind);
1109      for (int i=0;i<nsos;i++) {
1110        int ichar = sostype[i];
1111        assert (ichar=='1'||ichar=='2');
1112        typeSOS_[i]=ichar-'0';
1113      } 
1114      memcpy(prioritySOS_,sospri,nsos*sizeof(int));
1115      memcpy(startSOS_,sosbeg,(nsos+1)*sizeof(int));
1116      memcpy(memberSOS_,sosind,nsosnz*sizeof(int));
1117      memcpy(referenceSOS_,sosref,nsosnz*sizeof(double));
1118    }
1119  }
1120}
Note: See TracBrowser for help on using the repository browser.