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

Last change on this file since 496 was 496, checked in by forrest, 14 years ago

nonlinear

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 31.7 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    CoinModel * model = new CoinModel(1,fileName);
449    if (model->numberRows()>0)
450      *coinModel=(void *) model;
451    Oinfo.uinfo = tempBuffer;
452    if (getopts(argv, &Oinfo))
453      return 1;
454    if (objtype[0])
455      info->direction=-1.0;
456    else
457      info->direction=1.0;
458    info->offset=objconst(0);
459    info->numberRows=n_con;
460    info->numberColumns=n_var;
461    info->numberElements=nzc;;
462    info->numberBinary=nbv;
463    info->numberIntegers=niv;
464  }
465  /* add -solve - unless something there already
466   - also check for sleep=yes */
467  {
468    int found=0;
469    int foundLog=0;
470    int foundSleep=0;
471    const char * something[]={"solve","branch","duals","primals","user"};
472    for (i=0;i<info->numberArguments;i++) {
473      unsigned int j;
474      const char * argument = info->arguments[i];
475      for (j=0;j<sizeof(something)/sizeof(char *);j++) {
476        const char * check = something[j];
477        if (!strncmp(argument,check,sizeof(check))) {
478          found=(int)(j+1);
479        } else if (!strncmp(argument,"log",3)) {
480          foundLog=1;
481        } else if (!strncmp(argument,"sleep",5)) {
482          foundSleep=1;
483        }
484      }
485    }
486    if (foundLog) {
487      /* print options etc */
488      for (i=0;i<saveArgc;i++)
489        printf("%s ",saveArgv[i]);
490      printf("\n");
491      if (environment)
492        printf("env %s\n",environment);
493      /*printf("%d rows %d columns %d elements\n",n_con,n_var,nzc);*/
494    }
495    if (!found) {
496      if (!strlen(algFound)) {
497        info->arguments=(char **) realloc(info->arguments,(info->numberArguments+1)*sizeof(char *));
498        info->arguments[info->numberArguments++]=strdup("-solve");
499      } else {
500        // use algorithm from keyword
501        info->arguments=(char **) realloc(info->arguments,(info->numberArguments+1)*sizeof(char *));
502        info->arguments[info->numberArguments++]=strdup(algFound);
503      }
504    }
505    if (foundSleep) {
506      /* let user copy .nl file */
507      fprintf(stderr,"You can copy .nl file %s for debug purposes or attach debugger\n",saveArgv[1]);
508      fprintf(stderr,"Type q to quit, anything else to continue\n");
509      char getChar = getc(stdin);
510      if (getChar=='q'||getChar=='Q')
511        exit(1);
512    }
513  }
514  /* add -quit */
515  info->arguments=(char **) realloc(info->arguments,(info->numberArguments+1)*sizeof(char *));
516  info->arguments[info->numberArguments++]=strdup("-quit");
517  return 0;
518}
519void freeArrays1(ampl_info * info)
520{
521  free(info->objective);
522  info->objective=NULL;
523  free(info->rowLower);
524  info->rowLower=NULL;
525  free(info->rowUpper);
526  info->rowUpper=NULL;
527  free(info->columnLower);
528  info->columnLower=NULL;
529  free(info->columnUpper);
530  info->columnUpper=NULL;
531  /* this one not freed by ASL_free */
532  free(info->elements);
533  info->elements=NULL;
534  free(info->primalSolution);
535  info->primalSolution=NULL;
536  free(info->dualSolution);
537  info->dualSolution=NULL;
538  /*free(info->rowStatus);
539  info->rowStatus=NULL;
540  free(info->columnStatus);
541  info->columnStatus=NULL;*/
542}
543void freeArrays2(ampl_info * info)
544{
545  free(info->primalSolution);
546  info->primalSolution=NULL;
547  free(info->dualSolution);
548  info->dualSolution=NULL;
549  free(info->rowStatus);
550  info->rowStatus=NULL;
551  free(info->columnStatus);
552  info->columnStatus=NULL;
553  free(info->priorities);
554  info->priorities=NULL;
555  free(info->branchDirection);
556  info->branchDirection=NULL;
557  free(info->pseudoDown);
558  info->pseudoDown=NULL;
559  free(info->pseudoUp);
560  info->pseudoUp=NULL;
561  free(info->sosType);
562  info->sosType=NULL;
563  free(info->sosPriority);
564  info->sosPriority=NULL;
565  free(info->sosStart);
566  info->sosStart=NULL;
567  free(info->sosIndices);
568  info->sosIndices=NULL;
569  free(info->sosReference);
570  info->sosReference=NULL;
571  ASL_free(&asl);
572}
573void freeArgs(ampl_info * info)
574{
575  int i;
576  for ( i=0; i<info->numberArguments;i++)
577    free(info->arguments[i]);
578  free(info->arguments);
579}
580int ampl_obj_prec()
581{
582  return obj_prec();
583}
584void writeAmpl(ampl_info * info)
585{
586  char buf[1000];
587  typedef struct { const char *msg; int code; int wantObj; } Sol_info;
588  static Sol_info solinfo[] = {
589    { "optimal solution",                       000, 1 },
590    { "infeasible",                             200, 1 },
591    { "unbounded",                              300, 0 },
592    { "iteration limit etc",                    400, 1 },
593    { "solution limit",                         401, 1 },
594    { "ran out of space",                       500, 0 },
595    { "status unknown",                         501, 1 },
596    { "bug!",                                   502, 0 },
597    { "best MIP solution so far restored",      101, 1 },
598    { "failed to restore best MIP solution",    503, 1 },
599    { "optimal (?) solution",                   100, 1 }
600  };
601  /* convert status - need info on map */
602  static int map[] = {0, 3, 4, 1};
603  sprintf(buf,"%s %s",Oinfo.bsname,info->buffer);
604  solve_result_num = solinfo[info->problemStatus].code;
605  if (info->columnStatus) {
606    stat_map(info->columnStatus, n_var, map, 4, "outgoing columnStatus");
607    stat_map(info->rowStatus, n_con, map, 4, "outgoing rowStatus");
608    suf_iput("sstatus", ASL_Sufkind_var, info->columnStatus);
609    suf_iput("sstatus", ASL_Sufkind_con, info->rowStatus);
610  }
611  write_sol(buf,info->primalSolution,info->dualSolution,&Oinfo);
612}
613/* Read a problem from AMPL nl file
614 */
615CoinModel::CoinModel( int nonLinear, const char * fileName)
616 :  numberRows_(0),
617    maximumRows_(0),
618    numberColumns_(0),
619    maximumColumns_(0),
620    numberElements_(0),
621    maximumElements_(0),
622    numberQuadraticElements_(0),
623    maximumQuadraticElements_(0),
624    optimizationDirection_(1.0),
625    objectiveOffset_(0.0),
626    rowLower_(NULL),
627    rowUpper_(NULL),
628    rowType_(NULL),
629    objective_(NULL),
630    columnLower_(NULL),
631    columnUpper_(NULL),
632    integerType_(NULL),
633    columnType_(NULL),
634    start_(NULL),
635    elements_(NULL),
636    quadraticElements_(NULL),
637    sortIndices_(NULL),
638    sortElements_(NULL),
639    sortSize_(0),
640    sizeAssociated_(0),
641    associated_(NULL),
642    numberSOS_(0),
643    startSOS_(NULL),
644    memberSOS_(NULL),
645    typeSOS_(NULL),
646    prioritySOS_(NULL),
647    referenceSOS_(NULL),
648    priority_(NULL),
649    logLevel_(0),
650    type_(-1),
651    links_(0)
652{
653  problemName_ = strdup("");
654  int status=0;
655  if (!strcmp(fileName,"-")||!strcmp(fileName,"stdin")) {
656    // stdin
657  } else {
658    std::string name=fileName;
659    bool readable = fileCoinReadable(name);
660    if (!readable) {
661      std::cerr<<"Unable to open file "
662               <<fileName<<std::endl;
663      status = -1;
664    }
665  }
666  if (!status) {
667    gdb(nonLinear,fileName);
668  }
669}
670 static real
671qterm(ASL *asl, fint *colq, fint *rowq, real *delsq)
672{
673  double t, t1, *x, *x0, *xe;
674  fint *rq0, *rqe;
675 
676  t = 0.;
677  x0 = x = X0;
678  xe = x + n_var;
679  rq0 = rowq;
680  while(x < xe) {
681    t1 = *x++;
682    rqe = rq0 + *++colq;
683    while(rowq < rqe)
684      t += t1*x0[*rowq++]**delsq++;
685  }
686  return 0.5 * t;
687}
688
689void
690CoinModel::gdb( int nonLinear, const char * fileName)
691{
692  ograd *og=NULL;
693  int i;
694  SufDesc *csd=NULL;
695  SufDesc *rsd=NULL;
696  /*bool *basis, *lower;*/
697  /*double *LU, *c, lb, objadj, *rshift, *shift, t, ub, *x, *x0, *x1;*/
698  //char tempBuffer[20];
699  double * objective=NULL;
700  double * columnLower=NULL;
701  double * columnUpper=NULL;
702  double * rowLower=NULL;
703  double * rowUpper=NULL;
704  int * columnStatus=NULL;
705  int * rowStatus=NULL;
706  int numberRows=-1;
707  int numberColumns=-1;
708  int numberElements=-1;
709  int numberBinary=-1;
710  int numberIntegers=-1;
711  double * primalSolution=NULL;
712  double direction=1.0;
713  char * stub = strdup(fileName);
714  CoinPackedMatrix matrixByRow;
715  fint ** colqp=NULL;
716  int *z = NULL;
717  if (nonLinear==0) {
718    // linear
719    asl = ASL_alloc(ASL_read_f);
720    nl = jac0dim(stub, 0);
721    free(stub);
722    suf_declare(suftab, sizeof(suftab)/sizeof(SufDecl));
723   
724    /* set A_vals to get the constraints column-wise (malloc so can be freed) */
725    A_vals = (double *) malloc(nzc*sizeof(double));
726    if (!A_vals) {
727      printf("no memory\n");
728      return ;
729    }
730    /* say we want primal solution */
731    want_xpi0=1;
732    /* for basis info */
733    columnStatus = (int *) malloc(n_var*sizeof(int));
734    rowStatus = (int *) malloc(n_con*sizeof(int));
735    csd = suf_iput("sstatus", ASL_Sufkind_var, columnStatus);
736    rsd = suf_iput("sstatus", ASL_Sufkind_con, rowStatus);
737    /* read linear model*/
738    f_read(nl,0);
739    // see if any sos
740    if (true) {
741      char *sostype;
742      int nsosnz, *sosbeg, *sosind, * sospri;
743      double *sosref;
744      int nsos;
745      int i = ASL_suf_sos_explict_free;
746      int copri[2], **p_sospri;
747      copri[0] = 0;
748      copri[1] = 0;
749      p_sospri = &sospri;
750      nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
751                     &sosbeg, &sosind, &sosref);
752      if (nsos) {
753        abort();
754#if 0
755        info->numberSos=nsos;
756        info->sosType = (char *) malloc(nsos);
757        info->sosPriority = (int *) malloc(nsos*sizeof(int));
758        info->sosStart = (int *) malloc((nsos+1)*sizeof(int));
759        info->sosIndices = (int *) malloc(nsosnz*sizeof(int));
760        info->sosReference = (double *) malloc(nsosnz*sizeof(double));
761        sos_kludge(nsos, sosbeg, sosref,sosind);
762        for (int i=0;i<nsos;i++) {
763          int ichar = sostype[i];
764          assert (ichar=='1'||ichar=='2');
765          info->sosType[i]=ichar-'0';
766        }       
767        memcpy(info->sosPriority,sospri,nsos*sizeof(int));
768        memcpy(info->sosStart,sosbeg,(nsos+1)*sizeof(int));
769        memcpy(info->sosIndices,sosind,nsosnz*sizeof(int));
770        memcpy(info->sosReference,sosref,nsosnz*sizeof(double));
771#endif
772      }
773    }
774   
775    /*sos_finish(&specialOrderedInfo, 0, &j, 0, 0, 0, 0, 0);*/
776    //Oinfo.uinfo = tempBuffer;
777    //if (getopts(argv, &Oinfo))
778    //return 1;
779    /* objective*/
780    objective = (double *) malloc(n_var*sizeof(double));
781    for (i=0;i<n_var;i++)
782      objective[i]=0.0;;
783    if (n_obj) {
784      for (og = Ograd[0];og;og = og->next)
785        objective[og->varno] = og->coef;
786    }
787    if (objtype[0])
788      direction=-1.0;
789    else
790      direction=1.0;
791    objectiveOffset_=objconst(0);
792    /* Column bounds*/
793    columnLower = (double *) malloc(n_var*sizeof(double));
794    columnUpper = (double *) malloc(n_var*sizeof(double));
795#define COIN_DBL_MAX DBL_MAX
796    for (i=0;i<n_var;i++) {
797      columnLower[i]=LUv[2*i];
798      if (columnLower[i]<= negInfinity)
799        columnLower[i]=-COIN_DBL_MAX;
800      columnUpper[i]=LUv[2*i+1];
801      if (columnUpper[i]>= Infinity)
802        columnUpper[i]=COIN_DBL_MAX;
803    }
804    /* Row bounds*/
805    rowLower = (double *) malloc(n_con*sizeof(double));
806    rowUpper = (double *) malloc(n_con*sizeof(double));
807    for (i=0;i<n_con;i++) {
808      rowLower[i]=LUrhs[2*i];
809      if (rowLower[i]<= negInfinity)
810        rowLower[i]=-COIN_DBL_MAX;
811      rowUpper[i]=LUrhs[2*i+1];
812      if (rowUpper[i]>= Infinity)
813        rowUpper[i]=COIN_DBL_MAX;
814    }
815    numberRows=n_con;
816    numberColumns=n_var;
817    numberElements=nzc;;
818    numberBinary=nbv;
819    numberIntegers=niv;
820    /* put in primalSolution if exists */
821    if (X0) {
822      primalSolution=(double *) malloc(n_var*sizeof(double));
823      memcpy( primalSolution,X0,n_var*sizeof(double));
824    }
825    //double * dualSolution=NULL;
826    if (niv+nbv>0)
827      mip_stuff(); // get any extra info
828    if ((!(niv+nbv)&&(csd->kind & ASL_Sufkind_input))
829        ||(rsd->kind & ASL_Sufkind_input)) {
830      /* convert status - need info on map */
831      static int map[] = {1, 3, 1, 1, 2, 1, 1};
832      stat_map(columnStatus, n_var, map, 6, "incoming columnStatus");
833      stat_map(rowStatus, n_con, map, 6, "incoming rowStatus");
834    } else {
835      /* all slack basis */
836      // leave status for output */
837#if 0
838      free(rowStatus);
839      rowStatus=NULL;
840      free(columnStatus);
841      columnStatus=NULL;
842#endif
843    }
844    CoinPackedMatrix columnCopy(true,numberRows,numberColumns,numberElements,
845                                A_vals,A_rownos,A_colstarts,NULL);
846    matrixByRow.reverseOrderedCopyOf(columnCopy);
847  } else if (nonLinear==1) {
848    // quadratic
849    asl = ASL_alloc(ASL_read_fg);
850    nl = jac0dim(stub, (long) strlen(stub));
851    free(stub);
852    suf_declare(suftab, sizeof(suftab)/sizeof(SufDecl));
853    /* read  model*/
854    X0 = (double*) malloc(n_var*sizeof(double));
855    CoinZeroN(X0,n_var);
856    qp_read(nl,0);
857    assert (n_obj==1);
858    int nz = 1 + n_con;
859    colqp = (fint**) malloc(nz*(2*sizeof(int*)
860                                      + sizeof(double*)));
861    fint ** rowqp = colqp + nz;
862    double ** delsqp = (double **)(rowqp + nz);
863    z = (int*) malloc(nz*sizeof(int));
864    for (i=0;i<=n_con;i++) {
865      z[i] = nqpcheck(-i, rowqp+i, colqp+i, delsqp+i);
866    }
867    qp_opify();
868    /* objective*/
869    objective = (double *) malloc(n_var*sizeof(double));
870    for (i=0;i<n_var;i++)
871      objective[i]=0.0;;
872    if (n_obj) {
873      for (og = Ograd[0];og;og = og->next)
874        objective[og->varno] = og->coef;
875    }
876    if (objtype[0])
877      direction=-1.0;
878    else
879      direction=1.0;
880    objectiveOffset_=objconst(0);
881    /* Column bounds*/
882    columnLower = (double *) malloc(n_var*sizeof(double));
883    columnUpper = (double *) malloc(n_var*sizeof(double));
884    for (i=0;i<n_var;i++) {
885      columnLower[i]=LUv[2*i];
886      if (columnLower[i]<= negInfinity)
887        columnLower[i]=-COIN_DBL_MAX;
888      columnUpper[i]=LUv[2*i+1];
889      if (columnUpper[i]>= Infinity)
890        columnUpper[i]=COIN_DBL_MAX;
891    }
892    // Build by row from scratch
893    //matrixByRow.reserve(n_var,nzc,true);
894    // say row orderded
895    matrixByRow.transpose();
896    /* Row bounds*/
897    rowLower = (double *) malloc(n_con*sizeof(double));
898    rowUpper = (double *) malloc(n_con*sizeof(double));
899    int * column = new int [n_var];
900    double * element = new double [n_var];
901    for (i=0;i<n_con;i++) {
902      rowLower[i]=LUrhs[2*i];
903      if (rowLower[i]<= negInfinity)
904        rowLower[i]=-COIN_DBL_MAX;
905      rowUpper[i]=LUrhs[2*i+1];
906      if (rowUpper[i]>= Infinity)
907        rowUpper[i]=COIN_DBL_MAX;
908      int k=0;
909      for(cgrad * cg = Cgrad[i]; cg; cg = cg->next) {
910        column[k]=cg->varno;
911        element[k++]=cg->coef;
912      }
913      matrixByRow.appendRow(k,column,element);
914    }
915    numberRows=n_con;
916    numberColumns=n_var;
917    numberElements=nzc;;
918    numberBinary=nbv;
919    numberIntegers=niv+nlvci;
920    //double * dualSolution=NULL;
921  } else {
922    abort();
923  }
924  // set problem name
925  free(problemName_);
926  problemName_=strdup("???");
927
928  // Build by row from scratch
929  const double * element = matrixByRow.getElements();
930  const int * column = matrixByRow.getIndices();
931  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
932  const int * rowLength = matrixByRow.getVectorLengths();
933  for (i=0;i<numberRows;i++) {
934    addRow(rowLength[i],column+rowStart[i],
935           element+rowStart[i],rowLower[i],rowUpper[i]);
936  }
937  // Now do column part
938  for (i=0;i<numberColumns;i++) {
939    setColumnBounds(i,columnLower[i],columnUpper[i]);
940    setColumnObjective(i,objective[i]);
941  }
942  for ( i=numberColumns-numberBinary-numberIntegers;
943        i<numberColumns;i++) {
944      setColumnIsInteger(i,true);;
945  }
946  // do names
947  int iRow;
948  for (iRow=0;iRow<numberRows_;iRow++) {
949    char name[9];
950    sprintf(name,"r%7.7d",iRow);
951    setRowName(iRow,name);
952  }
953   
954  int iColumn;
955  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
956    char name[9];
957    sprintf(name,"c%7.7d",iColumn);
958    setColumnName(iColumn,name);
959  }
960  if (colqp) {
961    // add in quadratic
962    int nz = 1 + n_con;
963    fint ** rowqp = colqp + nz;
964    double ** delsqp = (double **)(rowqp + nz);
965    for (i=0;i<=n_con;i++) {
966      int nels = z[i];
967      if (nels) {
968        double * element = delsqp[i];
969        int * start = (int *) colqp[i];
970        int * row = (int *) rowqp[i];
971#if 0
972        printf("%d quadratic els\n",nels);
973        for (int j=0;j<n_var;j++) {
974          for (int k=start[j];k<start[j+1];k++)
975            printf("%d %d %g\n",j,row[k],element[k]);
976        }
977#endif
978        if (i) {
979          int iRow = i-1;
980          for (int j=0;j<n_var;j++) {
981            for (int k=start[j];k<start[j+1];k++) {
982              int kColumn = row[k];
983              double value = element[k];
984              // ampl gives twice with assumed 0.5
985              if (kColumn>j)
986                continue;
987              else if (kColumn==j)
988                value *= 0.5;
989              const char * expr = getElementAsString(iRow,j);
990              double constant=0.0;
991              bool linear;
992              if (expr&&strcmp(expr,"Numeric")) {
993                linear=false;
994              } else {
995                constant = getElement(iRow,j);
996                linear=true;
997              }
998              char temp[1000];
999              char temp2[30];
1000              if (value==1.0) 
1001                sprintf(temp2,"c%7.7d",kColumn);
1002              else
1003                sprintf(temp2,"%g*c%7.7d",value,kColumn);
1004              if (linear) {
1005                if (!constant)
1006                  strcpy(temp,temp2);
1007                else if (value>0.0) 
1008                  sprintf(temp,"%g+%s",constant,temp2);
1009                else
1010                  sprintf(temp,"%g%s",constant,temp2);
1011              } else {
1012                if (value>0.0) 
1013                  sprintf(temp,"%s+%s",expr,temp2);
1014                else
1015                  sprintf(temp,"%s%s",expr,temp2);
1016              }
1017              assert (strlen(temp)<1000);
1018              setElement(iRow,j,temp);
1019              printf("el for row %d column c%7.7d is %s\n",iRow,j,temp);
1020            }
1021          }
1022        } else {
1023          // objective
1024          for (int j=0;j<n_var;j++) {
1025            for (int k=start[j];k<start[j+1];k++) {
1026              int kColumn = row[k];
1027              double value = element[k];
1028              // ampl gives twice with assumed 0.5
1029              if (kColumn>j)
1030                continue;
1031              else if (kColumn==j)
1032                value *= 0.5;
1033              const char * expr = getColumnObjectiveAsString(j);
1034              double constant=0.0;
1035              bool linear;
1036              if (expr&&strcmp(expr,"Numeric")) {
1037                linear=false;
1038              } else {
1039                constant = getColumnObjective(j);
1040                linear=true;
1041              }
1042              char temp[1000];
1043              char temp2[30];
1044              if (value==1.0) 
1045                sprintf(temp2,"c%7.7d",kColumn);
1046              else
1047                sprintf(temp2,"%g*c%7.7d",value,kColumn);
1048              if (linear) {
1049                if (!constant)
1050                  strcpy(temp,temp2);
1051                else if (value>0.0) 
1052                  sprintf(temp,"%g+%s",constant,temp2);
1053                else
1054                  sprintf(temp,"%g%s",constant,temp2);
1055              } else {
1056                if (value>0.0) 
1057                  sprintf(temp,"%s+%s",expr,temp2);
1058                else
1059                  sprintf(temp,"%s%s",expr,temp2);
1060              }
1061              assert (strlen(temp)<1000);
1062              setObjective(j,temp);
1063              printf("el for objective column c%7.7d is %s\n",j,temp);
1064            }
1065          }
1066        }
1067      }
1068    }
1069  }
1070  // see if any sos
1071  {
1072    char *sostype;
1073    int nsosnz, *sosbeg, *sosind, * sospri;
1074    double *sosref;
1075    int nsos;
1076    int i = ASL_suf_sos_explict_free;
1077    int copri[2], **p_sospri;
1078    copri[0] = 0;
1079    copri[1] = 0;
1080    p_sospri = &sospri;
1081    nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
1082                   &sosbeg, &sosind, &sosref);
1083    if (nsos) {
1084      numberSOS_=nsos;
1085      typeSOS_ = new int [numberSOS_];
1086      prioritySOS_ = new int [numberSOS_];
1087      startSOS_ = new int [numberSOS_+1];
1088      memberSOS_ = new int[nsosnz];
1089      referenceSOS_ = new double [nsosnz];
1090      sos_kludge(nsos, sosbeg, sosref,sosind);
1091      for (int i=0;i<nsos;i++) {
1092        int ichar = sostype[i];
1093        assert (ichar=='1'||ichar=='2');
1094        typeSOS_[i]=ichar-'0';
1095      } 
1096      memcpy(prioritySOS_,sospri,nsos*sizeof(int));
1097      memcpy(startSOS_,sosbeg,(nsos+1)*sizeof(int));
1098      memcpy(memberSOS_,sosind,nsosnz*sizeof(int));
1099      memcpy(referenceSOS_,sosref,nsosnz*sizeof(double));
1100    }
1101  }
1102}
Note: See TracBrowser for help on using the repository browser.