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

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

nonlinear stuff

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 33.3 KB
Line 
1/****************************************************************
2Copyright (C) 1997-2000 Lucent Technologies
3Modifications for Coin -  Copyright (C) 2006, International Business Machines Corporation and others.
4All Rights Reserved
5
6Permission to use, copy, modify, and distribute this software and
7its documentation for any purpose and without fee is hereby
8granted, provided that the above copyright notice appear in all
9copies and that both that the copyright notice and this
10permission notice and warranty disclaimer appear in supporting
11documentation, and that the name of Lucent or any of its entities
12not be used in advertising or publicity pertaining to
13distribution of the software without specific, written prior
14permission.
15
16LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
17INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
18IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
19SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
20WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
21IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
22ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
23THIS SOFTWARE.
24****************************************************************/
25
26#include "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  if (info->numberRows!=-1234567)
295    memset(info,0,sizeof(ampl_info)); // overwrite unless magic number set
296  /* save so can be accessed by decodePhrase */
297  saveInfo = info;
298  info->numberArguments=0;
299  info->arguments=(char **) malloc(2*sizeof(char *));
300  info->arguments[info->numberArguments++]=strdup("ampl");
301  info->arguments[info->numberArguments++]=strdup("cbc");
302  asl = ASL_alloc(ASL_read_f);
303  stub = getstub(&argv, &Oinfo);
304  if (!stub)
305    usage_ASL(&Oinfo, 1);
306  nl = jac0dim(stub, 0);
307  suf_declare(suftab, sizeof(suftab)/sizeof(SufDecl));
308 
309  /* set A_vals to get the constraints column-wise (malloc so can be freed) */
310  A_vals = (double *) malloc(nzc*sizeof(double));
311  if (!A_vals) {
312    printf("no memory\n");
313    return 1;
314  }
315  /* say we want primal solution */
316  want_xpi0=1;
317  /* for basis info */
318  info->columnStatus = (int *) malloc(n_var*sizeof(int));
319  info->rowStatus = (int *) malloc(n_con*sizeof(int));
320  csd = suf_iput("sstatus", ASL_Sufkind_var, info->columnStatus);
321  rsd = suf_iput("sstatus", ASL_Sufkind_con, info->rowStatus);
322  // testosi parameter - if >= 10 then go in through coinModel
323  int testOsi=-1;
324  for (i=0;i<saveArgc;i++) {
325    if (!strncmp(saveArgv[i],"testosi",7)) {
326      testOsi = atoi(saveArgv[i+1]);
327      break;
328    }
329  }
330  if (!(nlvc+nlvo)||testOsi>=10) {
331    /* read linear model*/
332    f_read(nl,0);
333    // see if any sos
334    if (true) {
335      char *sostype;
336      int nsosnz, *sosbeg, *sosind, * sospri;
337      double *sosref;
338      int nsos;
339      int i = ASL_suf_sos_explict_free;
340      int copri[2], **p_sospri;
341      copri[0] = 0;
342      copri[1] = 0;
343      p_sospri = &sospri;
344      nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
345                     &sosbeg, &sosind, &sosref);
346      if (nsos) {
347        info->numberSos=nsos;
348        info->sosType = (char *) malloc(nsos);
349        info->sosPriority = (int *) malloc(nsos*sizeof(int));
350        info->sosStart = (int *) malloc((nsos+1)*sizeof(int));
351        info->sosIndices = (int *) malloc(nsosnz*sizeof(int));
352        info->sosReference = (double *) malloc(nsosnz*sizeof(double));
353        sos_kludge(nsos, sosbeg, sosref,sosind);
354        for (int i=0;i<nsos;i++) {
355          int ichar = sostype[i];
356          assert (ichar=='1'||ichar=='2');
357          info->sosType[i]=ichar-'0';
358        }       
359        memcpy(info->sosPriority,sospri,nsos*sizeof(int));
360        memcpy(info->sosStart,sosbeg,(nsos+1)*sizeof(int));
361        memcpy(info->sosIndices,sosind,nsosnz*sizeof(int));
362        memcpy(info->sosReference,sosref,nsosnz*sizeof(double));
363      }
364    }
365   
366    /*sos_finish(&specialOrderedInfo, 0, &j, 0, 0, 0, 0, 0);*/
367    Oinfo.uinfo = tempBuffer;
368    if (getopts(argv, &Oinfo))
369      return 1;
370    /* objective*/
371    obj = (double *) malloc(n_var*sizeof(double));
372    for (i=0;i<n_var;i++)
373      obj[i]=0.0;;
374    if (n_obj) {
375      for (og = Ograd[0];og;og = og->next)
376        obj[og->varno] = og->coef;
377    }
378    if (objtype[0])
379      info->direction=-1.0;
380    else
381      info->direction=1.0;
382    info->offset=objconst(0);
383    /* Column bounds*/
384    columnLower = (double *) malloc(n_var*sizeof(double));
385    columnUpper = (double *) malloc(n_var*sizeof(double));
386#define COIN_DBL_MAX DBL_MAX
387    for (i=0;i<n_var;i++) {
388      columnLower[i]=LUv[2*i];
389      if (columnLower[i]<= negInfinity)
390        columnLower[i]=-COIN_DBL_MAX;
391      columnUpper[i]=LUv[2*i+1];
392      if (columnUpper[i]>= Infinity)
393        columnUpper[i]=COIN_DBL_MAX;
394    }
395    /* Row bounds*/
396    rowLower = (double *) malloc(n_con*sizeof(double));
397    rowUpper = (double *) malloc(n_con*sizeof(double));
398    for (i=0;i<n_con;i++) {
399      rowLower[i]=LUrhs[2*i];
400      if (rowLower[i]<= negInfinity)
401        rowLower[i]=-COIN_DBL_MAX;
402      rowUpper[i]=LUrhs[2*i+1];
403      if (rowUpper[i]>= Infinity)
404        rowUpper[i]=COIN_DBL_MAX;
405    }
406    info->numberRows=n_con;
407    info->numberColumns=n_var;
408    info->numberElements=nzc;;
409    info->numberBinary=nbv;
410    info->numberIntegers=niv;
411    info->objective=obj;
412    info->rowLower=rowLower;
413    info->rowUpper=rowUpper;
414    info->columnLower=columnLower;
415    info->columnUpper=columnUpper;
416    info->starts=A_colstarts;
417    /*A_colstarts=NULL;*/
418    info->rows=A_rownos;
419    /*A_rownos=NULL;*/
420    info->elements=A_vals;
421    /*A_vals=NULL;*/
422    info->primalSolution=NULL;
423    /* put in primalSolution if exists */
424    if (X0) {
425      info->primalSolution=(double *) malloc(n_var*sizeof(double));
426      memcpy(info->primalSolution,X0,n_var*sizeof(double));
427    }
428    info->dualSolution=NULL;
429    if (niv+nbv>0)
430      mip_stuff(); // get any extra info
431    if ((!(niv+nbv)&&(csd->kind & ASL_Sufkind_input))
432        ||(rsd->kind & ASL_Sufkind_input)) {
433      /* convert status - need info on map */
434      static int map[] = {1, 3, 1, 1, 2, 1, 1};
435      stat_map(info->columnStatus, n_var, map, 6, "incoming columnStatus");
436      stat_map(info->rowStatus, n_con, map, 6, "incoming rowStatus");
437    } else {
438      /* all slack basis */
439      // leave status for output */
440#if 0
441      free(info->rowStatus);
442      info->rowStatus=NULL;
443      free(info->columnStatus);
444      info->columnStatus=NULL;
445#endif
446    }
447  } else {
448    // QP
449    // Add .nl if not there
450    if (!strstr(fileName,".nl"))
451      strcat(fileName,".nl");
452    CoinModel * model = new CoinModel(1,fileName,info);
453    if (model->numberRows()>0)
454      *coinModel=(void *) model;
455    Oinfo.uinfo = tempBuffer;
456    if (getopts(argv, &Oinfo))
457      return 1;
458    Oinfo.wantsol=1;
459    if (objtype[0])
460      info->direction=-1.0;
461    else
462      info->direction=1.0;
463    info->offset=objconst(0);
464    info->numberRows=n_con;
465    info->numberColumns=n_var;
466    info->numberElements=nzc;;
467    info->numberBinary=nbv;
468    info->numberIntegers=niv;
469    if (nbv+niv+nlvbi+nlvci+nlvoi>0)
470      mip_stuff(); // get any extra info
471  }
472  /* add -solve - unless something there already
473   - also check for sleep=yes */
474  {
475    int found=0;
476    int foundLog=0;
477    int foundSleep=0;
478    const char * something[]={"solve","branch","duals","primals","user"};
479    for (i=0;i<info->numberArguments;i++) {
480      unsigned int j;
481      const char * argument = info->arguments[i];
482      for (j=0;j<sizeof(something)/sizeof(char *);j++) {
483        const char * check = something[j];
484        if (!strncmp(argument,check,sizeof(check))) {
485          found=(int)(j+1);
486        } else if (!strncmp(argument,"log",3)) {
487          foundLog=1;
488        } else if (!strncmp(argument,"sleep",5)) {
489          foundSleep=1;
490        }
491      }
492    }
493    if (foundLog) {
494      /* print options etc */
495      for (i=0;i<saveArgc;i++)
496        printf("%s ",saveArgv[i]);
497      printf("\n");
498      if (environment)
499        printf("env %s\n",environment);
500      /*printf("%d rows %d columns %d elements\n",n_con,n_var,nzc);*/
501    }
502    if (!found) {
503      if (!strlen(algFound)) {
504        info->arguments=(char **) realloc(info->arguments,(info->numberArguments+1)*sizeof(char *));
505        info->arguments[info->numberArguments++]=strdup("-solve");
506      } else {
507        // use algorithm from keyword
508        info->arguments=(char **) realloc(info->arguments,(info->numberArguments+1)*sizeof(char *));
509        info->arguments[info->numberArguments++]=strdup(algFound);
510      }
511    }
512    if (foundSleep) {
513      /* let user copy .nl file */
514      fprintf(stderr,"You can copy .nl file %s for debug purposes or attach debugger\n",saveArgv[1]);
515      fprintf(stderr,"Type q to quit, anything else to continue\n");
516      char getChar = getc(stdin);
517      if (getChar=='q'||getChar=='Q')
518        exit(1);
519    }
520  }
521  /* add -quit */
522  info->arguments=(char **) realloc(info->arguments,(info->numberArguments+1)*sizeof(char *));
523  info->arguments[info->numberArguments++]=strdup("-quit");
524  return 0;
525}
526void freeArrays1(ampl_info * info)
527{
528  free(info->objective);
529  info->objective=NULL;
530  free(info->rowLower);
531  info->rowLower=NULL;
532  free(info->rowUpper);
533  info->rowUpper=NULL;
534  free(info->columnLower);
535  info->columnLower=NULL;
536  free(info->columnUpper);
537  info->columnUpper=NULL;
538  /* this one not freed by ASL_free */
539  free(info->elements);
540  info->elements=NULL;
541  free(info->primalSolution);
542  info->primalSolution=NULL;
543  free(info->dualSolution);
544  info->dualSolution=NULL;
545  /*free(info->rowStatus);
546  info->rowStatus=NULL;
547  free(info->columnStatus);
548  info->columnStatus=NULL;*/
549}
550void freeArrays2(ampl_info * info)
551{
552  free(info->primalSolution);
553  info->primalSolution=NULL;
554  free(info->dualSolution);
555  info->dualSolution=NULL;
556  free(info->rowStatus);
557  info->rowStatus=NULL;
558  free(info->columnStatus);
559  info->columnStatus=NULL;
560  free(info->priorities);
561  info->priorities=NULL;
562  free(info->branchDirection);
563  info->branchDirection=NULL;
564  free(info->pseudoDown);
565  info->pseudoDown=NULL;
566  free(info->pseudoUp);
567  info->pseudoUp=NULL;
568  free(info->sosType);
569  info->sosType=NULL;
570  free(info->sosPriority);
571  info->sosPriority=NULL;
572  free(info->sosStart);
573  info->sosStart=NULL;
574  free(info->sosIndices);
575  info->sosIndices=NULL;
576  free(info->sosReference);
577  info->sosReference=NULL;
578  ASL_free(&asl);
579}
580void freeArgs(ampl_info * info)
581{
582  int i;
583  for ( i=0; i<info->numberArguments;i++)
584    free(info->arguments[i]);
585  free(info->arguments);
586}
587int ampl_obj_prec()
588{
589  return obj_prec();
590}
591void writeAmpl(ampl_info * info)
592{
593  char buf[1000];
594  typedef struct { const char *msg; int code; int wantObj; } Sol_info;
595  static Sol_info solinfo[] = {
596    { "optimal solution",                       000, 1 },
597    { "infeasible",                             200, 1 },
598    { "unbounded",                              300, 0 },
599    { "iteration limit etc",                    400, 1 },
600    { "solution limit",                         401, 1 },
601    { "ran out of space",                       500, 0 },
602    { "status unknown",                         501, 1 },
603    { "bug!",                                   502, 0 },
604    { "best MIP solution so far restored",      101, 1 },
605    { "failed to restore best MIP solution",    503, 1 },
606    { "optimal (?) solution",                   100, 1 }
607  };
608  /* convert status - need info on map */
609  static int map[] = {0, 3, 4, 1};
610  sprintf(buf,"%s %s",Oinfo.bsname,info->buffer);
611  solve_result_num = solinfo[info->problemStatus].code;
612  if (info->columnStatus) {
613    stat_map(info->columnStatus, n_var, map, 4, "outgoing columnStatus");
614    stat_map(info->rowStatus, n_con, map, 4, "outgoing rowStatus");
615    suf_iput("sstatus", ASL_Sufkind_var, info->columnStatus);
616    suf_iput("sstatus", ASL_Sufkind_con, info->rowStatus);
617  }
618  write_sol(buf,info->primalSolution,info->dualSolution,&Oinfo);
619}
620/* Read a problem from AMPL nl file
621 */
622CoinModel::CoinModel( int nonLinear, const char * fileName,const void * info)
623 :  numberRows_(0),
624    maximumRows_(0),
625    numberColumns_(0),
626    maximumColumns_(0),
627    numberElements_(0),
628    maximumElements_(0),
629    numberQuadraticElements_(0),
630    maximumQuadraticElements_(0),
631    optimizationDirection_(1.0),
632    objectiveOffset_(0.0),
633    rowLower_(NULL),
634    rowUpper_(NULL),
635    rowType_(NULL),
636    objective_(NULL),
637    columnLower_(NULL),
638    columnUpper_(NULL),
639    integerType_(NULL),
640    columnType_(NULL),
641    start_(NULL),
642    elements_(NULL),
643    quadraticElements_(NULL),
644    sortIndices_(NULL),
645    sortElements_(NULL),
646    sortSize_(0),
647    sizeAssociated_(0),
648    associated_(NULL),
649    numberSOS_(0),
650    startSOS_(NULL),
651    memberSOS_(NULL),
652    typeSOS_(NULL),
653    prioritySOS_(NULL),
654    referenceSOS_(NULL),
655    priority_(NULL),
656    logLevel_(0),
657    type_(-1),
658    links_(0)
659{
660  problemName_ = strdup("");
661  int status=0;
662  if (!strcmp(fileName,"-")||!strcmp(fileName,"stdin")) {
663    // stdin
664  } else {
665    std::string name=fileName;
666    bool readable = fileCoinReadable(name);
667    if (!readable) {
668      std::cerr<<"Unable to open file "
669               <<fileName<<std::endl;
670      status = -1;
671    }
672  }
673  if (!status) {
674    gdb(nonLinear,fileName,info);
675  }
676}
677#if 0
678static real
679qterm(ASL *asl, fint *colq, fint *rowq, real *delsq)
680{
681  double t, t1, *x, *x0, *xe;
682  fint *rq0, *rqe;
683 
684  t = 0.;
685  x0 = x = X0;
686  xe = x + n_var;
687  rq0 = rowq;
688  while(x < xe) {
689    t1 = *x++;
690    rqe = rq0 + *++colq;
691    while(rowq < rqe)
692      t += t1*x0[*rowq++]**delsq++;
693  }
694  return 0.5 * t;
695}
696#endif
697void
698CoinModel::gdb( int nonLinear, const char * fileName,const void * info)
699{
700  const ampl_info * amplInfo = (const ampl_info *) info;
701  ograd *og=NULL;
702  int i;
703  SufDesc *csd=NULL;
704  SufDesc *rsd=NULL;
705  /*bool *basis, *lower;*/
706  /*double *LU, *c, lb, objadj, *rshift, *shift, t, ub, *x, *x0, *x1;*/
707  //char tempBuffer[20];
708  double * objective=NULL;
709  double * columnLower=NULL;
710  double * columnUpper=NULL;
711  double * rowLower=NULL;
712  double * rowUpper=NULL;
713  int * columnStatus=NULL;
714  int * rowStatus=NULL;
715  int numberRows=-1;
716  int numberColumns=-1;
717  int numberElements=-1;
718  int numberBinary=-1;
719  int numberIntegers=-1;
720  int numberAllNonLinearBoth=0;
721  int numberIntegerNonLinearBoth=0;
722  int numberAllNonLinearConstraints=0;
723  int numberIntegerNonLinearConstraints=0;
724  int numberAllNonLinearObjective=0;
725  int numberIntegerNonLinearObjective=0;
726  double * primalSolution=NULL;
727  double direction=1.0;
728  char * stub = strdup(fileName);
729  CoinPackedMatrix matrixByRow;
730  fint ** colqp=NULL;
731  int *z = NULL;
732  if (nonLinear==0) {
733    // linear
734    asl = ASL_alloc(ASL_read_f);
735    nl = jac0dim(stub, 0);
736    free(stub);
737    suf_declare(suftab, sizeof(suftab)/sizeof(SufDecl));
738   
739    /* set A_vals to get the constraints column-wise (malloc so can be freed) */
740    A_vals = (double *) malloc(nzc*sizeof(double));
741    if (!A_vals) {
742      printf("no memory\n");
743      return ;
744    }
745    /* say we want primal solution */
746    want_xpi0=1;
747    /* for basis info */
748    columnStatus = (int *) malloc(n_var*sizeof(int));
749    rowStatus = (int *) malloc(n_con*sizeof(int));
750    csd = suf_iput("sstatus", ASL_Sufkind_var, columnStatus);
751    rsd = suf_iput("sstatus", ASL_Sufkind_con, rowStatus);
752    /* read linear model*/
753    f_read(nl,0);
754    // see if any sos
755    if (true) {
756      char *sostype;
757      int nsosnz, *sosbeg, *sosind, * sospri;
758      double *sosref;
759      int nsos;
760      int i = ASL_suf_sos_explict_free;
761      int copri[2], **p_sospri;
762      copri[0] = 0;
763      copri[1] = 0;
764      p_sospri = &sospri;
765      nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
766                     &sosbeg, &sosind, &sosref);
767      if (nsos) {
768        abort();
769#if 0
770        info->numberSos=nsos;
771        info->sosType = (char *) malloc(nsos);
772        info->sosPriority = (int *) malloc(nsos*sizeof(int));
773        info->sosStart = (int *) malloc((nsos+1)*sizeof(int));
774        info->sosIndices = (int *) malloc(nsosnz*sizeof(int));
775        info->sosReference = (double *) malloc(nsosnz*sizeof(double));
776        sos_kludge(nsos, sosbeg, sosref,sosind);
777        for (int i=0;i<nsos;i++) {
778          int ichar = sostype[i];
779          assert (ichar=='1'||ichar=='2');
780          info->sosType[i]=ichar-'0';
781        }       
782        memcpy(info->sosPriority,sospri,nsos*sizeof(int));
783        memcpy(info->sosStart,sosbeg,(nsos+1)*sizeof(int));
784        memcpy(info->sosIndices,sosind,nsosnz*sizeof(int));
785        memcpy(info->sosReference,sosref,nsosnz*sizeof(double));
786#endif
787      }
788    }
789   
790    /*sos_finish(&specialOrderedInfo, 0, &j, 0, 0, 0, 0, 0);*/
791    //Oinfo.uinfo = tempBuffer;
792    //if (getopts(argv, &Oinfo))
793    //return 1;
794    /* objective*/
795    objective = (double *) malloc(n_var*sizeof(double));
796    for (i=0;i<n_var;i++)
797      objective[i]=0.0;;
798    if (n_obj) {
799      for (og = Ograd[0];og;og = og->next)
800        objective[og->varno] = og->coef;
801    }
802    if (objtype[0])
803      direction=-1.0;
804    else
805      direction=1.0;
806    objectiveOffset_=objconst(0);
807    /* Column bounds*/
808    columnLower = (double *) malloc(n_var*sizeof(double));
809    columnUpper = (double *) malloc(n_var*sizeof(double));
810#define COIN_DBL_MAX DBL_MAX
811    for (i=0;i<n_var;i++) {
812      columnLower[i]=LUv[2*i];
813      if (columnLower[i]<= negInfinity)
814        columnLower[i]=-COIN_DBL_MAX;
815      columnUpper[i]=LUv[2*i+1];
816      if (columnUpper[i]>= Infinity)
817        columnUpper[i]=COIN_DBL_MAX;
818    }
819    /* Row bounds*/
820    rowLower = (double *) malloc(n_con*sizeof(double));
821    rowUpper = (double *) malloc(n_con*sizeof(double));
822    for (i=0;i<n_con;i++) {
823      rowLower[i]=LUrhs[2*i];
824      if (rowLower[i]<= negInfinity)
825        rowLower[i]=-COIN_DBL_MAX;
826      rowUpper[i]=LUrhs[2*i+1];
827      if (rowUpper[i]>= Infinity)
828        rowUpper[i]=COIN_DBL_MAX;
829    }
830    numberRows=n_con;
831    numberColumns=n_var;
832    numberElements=nzc;;
833    numberBinary=nbv;
834    numberIntegers=niv;
835    /* put in primalSolution if exists */
836    if (X0) {
837      primalSolution=(double *) malloc(n_var*sizeof(double));
838      memcpy( primalSolution,X0,n_var*sizeof(double));
839    }
840    //double * dualSolution=NULL;
841    if (niv+nbv>0)
842      mip_stuff(); // get any extra info
843    if ((!(niv+nbv)&&(csd->kind & ASL_Sufkind_input))
844        ||(rsd->kind & ASL_Sufkind_input)) {
845      /* convert status - need info on map */
846      static int map[] = {1, 3, 1, 1, 2, 1, 1};
847      stat_map(columnStatus, n_var, map, 6, "incoming columnStatus");
848      stat_map(rowStatus, n_con, map, 6, "incoming rowStatus");
849    } else {
850      /* all slack basis */
851      // leave status for output */
852#if 0
853      free(rowStatus);
854      rowStatus=NULL;
855      free(columnStatus);
856      columnStatus=NULL;
857#endif
858    }
859    CoinPackedMatrix columnCopy(true,numberRows,numberColumns,numberElements,
860                                A_vals,A_rownos,A_colstarts,NULL);
861    matrixByRow.reverseOrderedCopyOf(columnCopy);
862  } else if (nonLinear==1) {
863    // quadratic
864    asl = ASL_alloc(ASL_read_fg);
865    nl = jac0dim(stub, (long) strlen(stub));
866    free(stub);
867    suf_declare(suftab, sizeof(suftab)/sizeof(SufDecl));
868    /* read  model*/
869    X0 = (double*) malloc(n_var*sizeof(double));
870    CoinZeroN(X0,n_var);
871    qp_read(nl,0);
872    assert (n_obj==1);
873    int nz = 1 + n_con;
874    colqp = (fint**) malloc(nz*(2*sizeof(int*)
875                                      + sizeof(double*)));
876    fint ** rowqp = colqp + nz;
877    double ** delsqp = (double **)(rowqp + nz);
878    z = (int*) malloc(nz*sizeof(int));
879    for (i=0;i<=n_con;i++) {
880      z[i] = nqpcheck(-i, rowqp+i, colqp+i, delsqp+i);
881    }
882    qp_opify();
883    /* objective*/
884    objective = (double *) malloc(n_var*sizeof(double));
885    for (i=0;i<n_var;i++)
886      objective[i]=0.0;;
887    if (n_obj) {
888      for (og = Ograd[0];og;og = og->next)
889        objective[og->varno] = og->coef;
890    }
891    if (objtype[0])
892      direction=-1.0;
893    else
894      direction=1.0;
895    objectiveOffset_=objconst(0);
896    /* Column bounds*/
897    columnLower = (double *) malloc(n_var*sizeof(double));
898    columnUpper = (double *) malloc(n_var*sizeof(double));
899    for (i=0;i<n_var;i++) {
900      columnLower[i]=LUv[2*i];
901      if (columnLower[i]<= negInfinity)
902        columnLower[i]=-COIN_DBL_MAX;
903      columnUpper[i]=LUv[2*i+1];
904      if (columnUpper[i]>= Infinity)
905        columnUpper[i]=COIN_DBL_MAX;
906    }
907    // Build by row from scratch
908    //matrixByRow.reserve(n_var,nzc,true);
909    // say row orderded
910    matrixByRow.transpose();
911    /* Row bounds*/
912    rowLower = (double *) malloc(n_con*sizeof(double));
913    rowUpper = (double *) malloc(n_con*sizeof(double));
914    int * column = new int [n_var];
915    double * element = new double [n_var];
916    for (i=0;i<n_con;i++) {
917      rowLower[i]=LUrhs[2*i];
918      if (rowLower[i]<= negInfinity)
919        rowLower[i]=-COIN_DBL_MAX;
920      rowUpper[i]=LUrhs[2*i+1];
921      if (rowUpper[i]>= Infinity)
922        rowUpper[i]=COIN_DBL_MAX;
923      int k=0;
924      for(cgrad * cg = Cgrad[i]; cg; cg = cg->next) {
925        column[k]=cg->varno;
926        element[k++]=cg->coef;
927      }
928      matrixByRow.appendRow(k,column,element);
929    }
930    numberRows=n_con;
931    numberColumns=n_var;
932    numberElements=nzc;;
933    numberBinary=nbv;
934    numberIntegers=niv;
935    numberAllNonLinearBoth=nlvb;
936    numberIntegerNonLinearBoth=nlvbi;
937    numberAllNonLinearConstraints=nlvc;
938    numberIntegerNonLinearConstraints=nlvci;
939    numberAllNonLinearObjective=nlvo;
940    numberIntegerNonLinearObjective=nlvoi;
941    /* say we want primal solution */
942    want_xpi0=1;
943    //double * dualSolution=NULL;
944  } else {
945    abort();
946  }
947  // set problem name
948  free(problemName_);
949  problemName_=strdup("???");
950
951  // Build by row from scratch
952  const double * element = matrixByRow.getElements();
953  const int * column = matrixByRow.getIndices();
954  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
955  const int * rowLength = matrixByRow.getVectorLengths();
956  for (i=0;i<numberRows;i++) {
957    addRow(rowLength[i],column+rowStart[i],
958           element+rowStart[i],rowLower[i],rowUpper[i]);
959  }
960  // Now do column part
961  for (i=0;i<numberColumns;i++) {
962    setColumnBounds(i,columnLower[i],columnUpper[i]);
963    setColumnObjective(i,objective[i]);
964  }
965  for ( i=numberColumns-numberBinary-numberIntegers;
966        i<numberColumns;i++) {
967      setColumnIsInteger(i,true);;
968  }
969  // and non linear
970  for (i=numberAllNonLinearBoth-numberIntegerNonLinearBoth;
971       i<numberAllNonLinearBoth;i++) {
972    setColumnIsInteger(i,true);;
973  }
974  for (i=numberAllNonLinearConstraints-numberIntegerNonLinearConstraints;
975       i<numberAllNonLinearConstraints;i++) {
976    setColumnIsInteger(i,true);;
977  }
978  for (i=numberAllNonLinearObjective-numberIntegerNonLinearObjective;
979       i<numberAllNonLinearObjective;i++) {
980    setColumnIsInteger(i,true);;
981  }
982  // do names
983  int iRow;
984  for (iRow=0;iRow<numberRows_;iRow++) {
985    char name[9];
986    sprintf(name,"r%7.7d",iRow);
987    setRowName(iRow,name);
988  }
989   
990  int iColumn;
991  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
992    char name[9];
993    sprintf(name,"c%7.7d",iColumn);
994    setColumnName(iColumn,name);
995  }
996  if (colqp) {
997    // add in quadratic
998    int nz = 1 + n_con;
999    int nOdd=0;
1000    fint ** rowqp = colqp + nz;
1001    double ** delsqp = (double **)(rowqp + nz);
1002    for (i=0;i<=n_con;i++) {
1003      int nels = z[i];
1004      if (nels) {
1005        double * element = delsqp[i];
1006        int * start = (int *) colqp[i];
1007        int * row = (int *) rowqp[i];
1008        if (!element) {
1009          // odd row - probably not quadratic
1010          nOdd++;
1011          continue;
1012        }
1013#if 0
1014        printf("%d quadratic els\n",nels);
1015        for (int j=0;j<n_var;j++) {
1016          for (int k=start[j];k<start[j+1];k++)
1017            printf("%d %d %g\n",j,row[k],element[k]);
1018        }
1019#endif
1020        if (i) {
1021          int iRow = i-1;
1022          for (int j=0;j<n_var;j++) {
1023            for (int k=start[j];k<start[j+1];k++) {
1024              int kColumn = row[k];
1025              double value = element[k];
1026              // ampl gives twice with assumed 0.5
1027              if (kColumn<j)
1028                continue;
1029              else if (kColumn==j)
1030                value *= 0.5;
1031              const char * expr = getElementAsString(iRow,j);
1032              double constant=0.0;
1033              bool linear;
1034              if (expr&&strcmp(expr,"Numeric")) {
1035                linear=false;
1036              } else {
1037                constant = getElement(iRow,j);
1038                linear=true;
1039              }
1040              char temp[1000];
1041              char temp2[30];
1042              if (value==1.0) 
1043                sprintf(temp2,"c%7.7d",kColumn);
1044              else
1045                sprintf(temp2,"%g*c%7.7d",value,kColumn);
1046              if (linear) {
1047                if (!constant)
1048                  strcpy(temp,temp2);
1049                else if (value>0.0) 
1050                  sprintf(temp,"%g+%s",constant,temp2);
1051                else
1052                  sprintf(temp,"%g%s",constant,temp2);
1053              } else {
1054                if (value>0.0) 
1055                  sprintf(temp,"%s+%s",expr,temp2);
1056                else
1057                  sprintf(temp,"%s%s",expr,temp2);
1058              }
1059              assert (strlen(temp)<1000);
1060              setElement(iRow,j,temp);
1061              if (amplInfo->logLevel>1)
1062                printf("el for row %d column c%7.7d is %s\n",iRow,j,temp);
1063            }
1064          }
1065        } else {
1066          // objective
1067          for (int j=0;j<n_var;j++) {
1068            for (int k=start[j];k<start[j+1];k++) {
1069              int kColumn = row[k];
1070              double value = element[k];
1071              // ampl gives twice with assumed 0.5
1072              if (kColumn<j)
1073                continue;
1074              else if (kColumn==j)
1075                value *= 0.5;
1076              const char * expr = getColumnObjectiveAsString(j);
1077              double constant=0.0;
1078              bool linear;
1079              if (expr&&strcmp(expr,"Numeric")) {
1080                linear=false;
1081              } else {
1082                constant = getColumnObjective(j);
1083                linear=true;
1084              }
1085              char temp[1000];
1086              char temp2[30];
1087              if (value==1.0) 
1088                sprintf(temp2,"c%7.7d",kColumn);
1089              else
1090                sprintf(temp2,"%g*c%7.7d",value,kColumn);
1091              if (linear) {
1092                if (!constant)
1093                  strcpy(temp,temp2);
1094                else if (value>0.0) 
1095                  sprintf(temp,"%g+%s",constant,temp2);
1096                else
1097                  sprintf(temp,"%g%s",constant,temp2);
1098              } else {
1099                if (value>0.0) 
1100                  sprintf(temp,"%s+%s",expr,temp2);
1101                else
1102                  sprintf(temp,"%s%s",expr,temp2);
1103              }
1104              assert (strlen(temp)<1000);
1105              setObjective(j,temp);
1106              if (amplInfo->logLevel>1)
1107                printf("el for objective column c%7.7d is %s\n",j,temp);
1108            }
1109          }
1110        }
1111      }
1112    }
1113    if (nOdd) {
1114      printf("%d non-linear constraints could not be converted to quadratic\n",nOdd);
1115      exit(77);
1116    }
1117  }
1118  // see if any sos
1119  {
1120    char *sostype;
1121    int nsosnz, *sosbeg, *sosind, * sospri;
1122    double *sosref;
1123    int nsos;
1124    int i = ASL_suf_sos_explict_free;
1125    int copri[2], **p_sospri;
1126    copri[0] = 0;
1127    copri[1] = 0;
1128    p_sospri = &sospri;
1129    nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
1130                   &sosbeg, &sosind, &sosref);
1131    if (nsos) {
1132      numberSOS_=nsos;
1133      typeSOS_ = new int [numberSOS_];
1134      prioritySOS_ = new int [numberSOS_];
1135      startSOS_ = new int [numberSOS_+1];
1136      memberSOS_ = new int[nsosnz];
1137      referenceSOS_ = new double [nsosnz];
1138      sos_kludge(nsos, sosbeg, sosref,sosind);
1139      for (int i=0;i<nsos;i++) {
1140        int ichar = sostype[i];
1141        assert (ichar=='1'||ichar=='2');
1142        typeSOS_[i]=ichar-'0';
1143      } 
1144      memcpy(prioritySOS_,sospri,nsos*sizeof(int));
1145      memcpy(startSOS_,sosbeg,(nsos+1)*sizeof(int));
1146      memcpy(memberSOS_,sosind,nsosnz*sizeof(int));
1147      memcpy(referenceSOS_,sosref,nsosnz*sizeof(double));
1148    }
1149  }
1150}
Note: See TracBrowser for help on using the repository browser.