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

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

for output

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