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

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

quadratic again

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 31.9 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    int nOdd=0;
964    fint ** rowqp = colqp + nz;
965    double ** delsqp = (double **)(rowqp + nz);
966    for (i=0;i<=n_con;i++) {
967      int nels = z[i];
968      if (nels) {
969        double * element = delsqp[i];
970        int * start = (int *) colqp[i];
971        int * row = (int *) rowqp[i];
972        if (!element) {
973          // odd row - probably not quadratic
974          nOdd++;
975          continue;
976        }
977#if 0
978        printf("%d quadratic els\n",nels);
979        for (int j=0;j<n_var;j++) {
980          for (int k=start[j];k<start[j+1];k++)
981            printf("%d %d %g\n",j,row[k],element[k]);
982        }
983#endif
984        if (i) {
985          int iRow = i-1;
986          for (int j=0;j<n_var;j++) {
987            for (int k=start[j];k<start[j+1];k++) {
988              int kColumn = row[k];
989              double value = element[k];
990              // ampl gives twice with assumed 0.5
991              if (kColumn>j)
992                continue;
993              else if (kColumn==j)
994                value *= 0.5;
995              const char * expr = getElementAsString(iRow,j);
996              double constant=0.0;
997              bool linear;
998              if (expr&&strcmp(expr,"Numeric")) {
999                linear=false;
1000              } else {
1001                constant = getElement(iRow,j);
1002                linear=true;
1003              }
1004              char temp[1000];
1005              char temp2[30];
1006              if (value==1.0) 
1007                sprintf(temp2,"c%7.7d",kColumn);
1008              else
1009                sprintf(temp2,"%g*c%7.7d",value,kColumn);
1010              if (linear) {
1011                if (!constant)
1012                  strcpy(temp,temp2);
1013                else if (value>0.0) 
1014                  sprintf(temp,"%g+%s",constant,temp2);
1015                else
1016                  sprintf(temp,"%g%s",constant,temp2);
1017              } else {
1018                if (value>0.0) 
1019                  sprintf(temp,"%s+%s",expr,temp2);
1020                else
1021                  sprintf(temp,"%s%s",expr,temp2);
1022              }
1023              assert (strlen(temp)<1000);
1024              setElement(iRow,j,temp);
1025              printf("el for row %d column c%7.7d is %s\n",iRow,j,temp);
1026            }
1027          }
1028        } else {
1029          // objective
1030          for (int j=0;j<n_var;j++) {
1031            for (int k=start[j];k<start[j+1];k++) {
1032              int kColumn = row[k];
1033              double value = element[k];
1034              // ampl gives twice with assumed 0.5
1035              if (kColumn>j)
1036                continue;
1037              else if (kColumn==j)
1038                value *= 0.5;
1039              const char * expr = getColumnObjectiveAsString(j);
1040              double constant=0.0;
1041              bool linear;
1042              if (expr&&strcmp(expr,"Numeric")) {
1043                linear=false;
1044              } else {
1045                constant = getColumnObjective(j);
1046                linear=true;
1047              }
1048              char temp[1000];
1049              char temp2[30];
1050              if (value==1.0) 
1051                sprintf(temp2,"c%7.7d",kColumn);
1052              else
1053                sprintf(temp2,"%g*c%7.7d",value,kColumn);
1054              if (linear) {
1055                if (!constant)
1056                  strcpy(temp,temp2);
1057                else if (value>0.0) 
1058                  sprintf(temp,"%g+%s",constant,temp2);
1059                else
1060                  sprintf(temp,"%g%s",constant,temp2);
1061              } else {
1062                if (value>0.0) 
1063                  sprintf(temp,"%s+%s",expr,temp2);
1064                else
1065                  sprintf(temp,"%s%s",expr,temp2);
1066              }
1067              assert (strlen(temp)<1000);
1068              setObjective(j,temp);
1069              printf("el for objective column c%7.7d is %s\n",j,temp);
1070            }
1071          }
1072        }
1073      }
1074    }
1075    if (nOdd) {
1076      printf("%d non-linear constraints could not be converted to quadratic\n",nOdd);
1077      exit(77);
1078    }
1079  }
1080  // see if any sos
1081  {
1082    char *sostype;
1083    int nsosnz, *sosbeg, *sosind, * sospri;
1084    double *sosref;
1085    int nsos;
1086    int i = ASL_suf_sos_explict_free;
1087    int copri[2], **p_sospri;
1088    copri[0] = 0;
1089    copri[1] = 0;
1090    p_sospri = &sospri;
1091    nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
1092                   &sosbeg, &sosind, &sosref);
1093    if (nsos) {
1094      numberSOS_=nsos;
1095      typeSOS_ = new int [numberSOS_];
1096      prioritySOS_ = new int [numberSOS_];
1097      startSOS_ = new int [numberSOS_+1];
1098      memberSOS_ = new int[nsosnz];
1099      referenceSOS_ = new double [nsosnz];
1100      sos_kludge(nsos, sosbeg, sosref,sosind);
1101      for (int i=0;i<nsos;i++) {
1102        int ichar = sostype[i];
1103        assert (ichar=='1'||ichar=='2');
1104        typeSOS_[i]=ichar-'0';
1105      } 
1106      memcpy(prioritySOS_,sospri,nsos*sizeof(int));
1107      memcpy(startSOS_,sosbeg,(nsos+1)*sizeof(int));
1108      memcpy(memberSOS_,sosind,nsosnz*sizeof(int));
1109      memcpy(referenceSOS_,sosref,nsosnz*sizeof(double));
1110    }
1111  }
1112}
Note: See TracBrowser for help on using the repository browser.