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

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

cuts or special

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 47.7 KB
Line 
1/****************************************************************
2Copyright (C) 1997-2000 Lucent Technologies
3Modifications for Coin -  Copyright (C) 2006, International Business Machines Corporation and others.
4All Rights Reserved
5
6Permission to use, copy, modify, and distribute this software and
7its documentation for any purpose and without fee is hereby
8granted, provided that the above copyright notice appear in all
9copies and that both that the copyright notice and this
10permission notice and warranty disclaimer appear in supporting
11documentation, and that the name of Lucent or any of its entities
12not be used in advertising or publicity pertaining to
13distribution of the software without specific, written prior
14permission.
15
16LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
17INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
18IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
19SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
20WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
21IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
22ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
23THIS SOFTWARE.
24****************************************************************/
25#ifdef COIN_HAS_ASL
26
27#include "CbcConfig.h"
28#ifdef HAVE_UNISTD_H
29# include "unistd.h"
30#endif
31#include "CoinUtilsConfig.h"
32#include "CoinHelperFunctions.hpp"
33#include "CoinModel.hpp"
34#include "CoinSort.hpp"
35#include "CoinPackedMatrix.hpp"
36#include "CoinMpsIO.hpp"
37#include "CoinFloatEqual.hpp"
38
39#include "Cbc_ampl.h"
40extern "C" {
41# include "getstub.h"
42# include "asl_pfgh.h"
43}
44
45#include <string>
46#include <cassert>
47/* so decodePhrase and clpCheck can access */
48static ampl_info * saveInfo=NULL;
49// Set to 1 if algorithm found
50static char algFound[20]="";
51static char*
52checkPhrase(Option_Info *oi, keyword *kw, char *v)
53{
54  if (strlen(v))
55    printf("string %s\n",v);
56  // Say algorithm found
57  strcpy(algFound,kw->desc);;
58  return v;
59}
60static char*
61checkPhrase2(Option_Info *oi, keyword *kw, char *v)
62{
63  if (strlen(v))
64    printf("string %s\n",v);
65  // put out keyword
66  saveInfo->arguments=(char **) realloc(saveInfo->arguments,(saveInfo->numberArguments+1)*sizeof(char *));
67  saveInfo->arguments[saveInfo->numberArguments++]=strdup(kw->desc);
68  return v;
69}
70static fint
71decodePhrase(char * phrase,ftnlen length)
72{
73  char * blank = strchr(phrase,' ');
74  if (blank) {
75    /* split arguments */
76    *blank='\0';
77    saveInfo->arguments=(char **) realloc(saveInfo->arguments,(saveInfo->numberArguments+2)*sizeof(char *));
78    saveInfo->arguments[saveInfo->numberArguments++]=strdup(phrase);
79    *blank=' ';
80    phrase=blank+1; /* move on */
81    if (strlen(phrase))
82      saveInfo->arguments[saveInfo->numberArguments++]=strdup(phrase);
83  } else if (strlen(phrase)) {
84    saveInfo->arguments=(char **) realloc(saveInfo->arguments,(saveInfo->numberArguments+1)*sizeof(char *));
85    saveInfo->arguments[saveInfo->numberArguments++]=strdup(phrase);
86  }
87  return 0;
88}
89static void
90sos_kludge(int nsos, int *sosbeg, double *sosref,int * sosind)
91{
92  // Adjust sosref if necessary to make monotonic increasing
93  int i, j, k;
94  // first sort
95  for (i=0;i<nsos;i++) {
96    k = sosbeg[i];
97    int end=sosbeg[i+1];
98    CoinSort_2(sosref+k,sosref+end,sosind+k);
99  }
100  double t, t1;
101  for(i = j = 0; i++ < nsos; ) {
102    k = sosbeg[i];
103    t = sosref[j];
104    while(++j < k) {
105      t1 = sosref[j];
106      t += 1e-10;
107      if (t1 <= t)
108        sosref[j] = t1 = t + 1e-10;
109      t = t1;
110    }
111  }
112}
113static char xxxxxx[20];
114#define VP (char*)
115 static keyword keywds[] = { /* must be sorted */
116        { "barrier",    checkPhrase,            (char *) xxxxxx ,"-barrier" },
117        { "dual",       checkPhrase,            (char *) xxxxxx , "-dualsimplex"},
118        { "help",       checkPhrase2,           (char *) xxxxxx , "-?"},
119        { "initial",    checkPhrase,            (char *) xxxxxx , "-initialsolve"},
120        { "max",        checkPhrase2,           (char *) xxxxxx , "-maximize"},
121        { "maximize",   checkPhrase2,           (char *) xxxxxx , "-maximize"},
122        { "primal",     checkPhrase,            (char *) xxxxxx , "-primalsimplex"},
123        { "quit",       checkPhrase2,           (char *) xxxxxx , "-quit"},
124        { "wantsol",    WS_val,         NULL, "write .sol file (without -AMPL)" }
125        };
126static Option_Info Oinfo = {"cbc", "Cbc 1.01", "cbc_options", keywds, nkeywds, 0, "",
127                                0,decodePhrase,0,0,0, 20060130 };
128// strdup used to avoid g++ compiler warning
129 static SufDecl
130suftab[] = {
131#if 0
132        { "current", 0, ASL_Sufkind_con | ASL_Sufkind_outonly },
133        { "current", 0, ASL_Sufkind_var | ASL_Sufkind_outonly },
134        { "direction", 0, ASL_Sufkind_var },
135        { "down", 0, ASL_Sufkind_con | ASL_Sufkind_outonly },
136        { "down", 0, ASL_Sufkind_var | ASL_Sufkind_outonly },
137        { "priority", 0, ASL_Sufkind_var },
138#endif
139        { "cut", 0, ASL_Sufkind_con },
140        { "direction", 0, ASL_Sufkind_var },
141        { "downPseudocost", 0, ASL_Sufkind_var | ASL_Sufkind_real },
142        { "priority", 0, ASL_Sufkind_var },
143        { "ref", 0, ASL_Sufkind_var | ASL_Sufkind_real },
144        { "sos", 0, ASL_Sufkind_var },
145        { "sos", 0, ASL_Sufkind_con },
146        { "sosno", 0, ASL_Sufkind_var | ASL_Sufkind_real },
147        { "sosref", 0, ASL_Sufkind_var | ASL_Sufkind_real },
148        { "special", 0, ASL_Sufkind_var },
149        { "special", 0, ASL_Sufkind_con },
150        /*{ "special", 0, ASL_Sufkind_con },*/
151        { strdup("sstatus"), 0, ASL_Sufkind_var, 0 },
152        { strdup("sstatus"), 0, ASL_Sufkind_con, 0 },
153        { "upPseudocost", 0, ASL_Sufkind_var | ASL_Sufkind_real }
154#if 0
155        { "unbdd", 0, ASL_Sufkind_var | ASL_Sufkind_outonly},
156        { "up", 0, ASL_Sufkind_con | ASL_Sufkind_outonly },
157        { "up", 0, ASL_Sufkind_var | ASL_Sufkind_outonly }
158#endif
159        };
160#include "float.h"
161#include "limits.h"
162static ASL *asl=NULL;
163static FILE *nl=NULL;
164
165static void
166mip_stuff(void)
167{
168  int i;
169  double *pseudoUp, *pseudoDown;
170  int *priority, *direction;
171  // To label cuts (there will be other uses for special)
172  int *cut;
173  // To label special variables - at present 1= must be >= 1 or <= -1
174  int * special;
175  SufDesc *dpup, *dpdown, *dpri, *ddir, *dcut, *dspecial;
176 
177  ddir = suf_get("direction", ASL_Sufkind_var);
178  direction = ddir->u.i;
179  dpri = suf_get("priority", ASL_Sufkind_var);
180  priority = dpri->u.i;
181  dspecial = suf_get("special", ASL_Sufkind_con);
182  dcut = suf_get("cut", ASL_Sufkind_con);
183  cut = dcut->u.i;
184  if (!cut) {
185    // try special
186    dcut = suf_get("special", ASL_Sufkind_con);
187    cut = dcut->u.i;
188  }
189  dspecial = suf_get("special", ASL_Sufkind_var);
190  special = dspecial->u.i;
191  dpdown = suf_get("downPseudocost", ASL_Sufkind_var);
192  pseudoDown = dpdown->u.r;
193  dpup = suf_get("upPseudocost", ASL_Sufkind_var);
194  pseudoUp = dpup->u.r;
195  assert(saveInfo);
196  int numberColumns = saveInfo->numberColumns;
197  if (direction) {
198    int baddir=0;
199    saveInfo->branchDirection = (int *) malloc(numberColumns*sizeof(int));
200    for (i=0;i<numberColumns;i++) {
201      int value = direction[i];
202      if (value<-1||value>1) {
203        baddir++;
204        value=0;
205      }
206      saveInfo->branchDirection[i]=value;
207    }
208    if (baddir)
209      fprintf(Stderr,
210              "Treating %d .direction values outside [-1, 1] as 0.\n",
211              baddir);
212  }
213  if (priority) {
214    int badpri=0;
215    saveInfo->priorities= (int *) malloc(numberColumns*sizeof(int));
216    for (i=0;i<numberColumns;i++) {
217      int value = priority[i];
218      if (value<0) {
219        badpri++;
220        value=0;
221      }
222      saveInfo->priorities[i]=value;
223    }
224    if (badpri)
225      fprintf(Stderr,
226              "Treating %d negative .priority values as 0\n",
227              badpri);
228  }
229  if (special) {
230    int badspecial=0;
231    saveInfo->special= (int *) malloc(numberColumns*sizeof(int));
232    for (i=0;i<numberColumns;i++) {
233      int value = special[i];
234      if (value<0) {
235        badspecial++;
236        value=0;
237      }
238      saveInfo->special[i]=value;
239    }
240    if (badspecial)
241      fprintf(Stderr,
242              "Treating %d negative special values as 0\n",
243              badspecial);
244  }
245  int numberRows = saveInfo->numberRows;
246  if (cut) {
247    int badcut=0;
248    saveInfo->cut= (int *) malloc(numberRows*sizeof(int));
249    for (i=0;i<numberRows;i++) {
250      int value = cut[i];
251      if (value<0) {
252        badcut++;
253        value=0;
254      }
255      saveInfo->cut[i]=value;
256    }
257    if (badcut)
258      fprintf(Stderr,
259              "Treating %d negative cut values as 0\n",
260              badcut);
261  }
262  if (pseudoDown||pseudoUp) {
263    int badpseudo=0;
264    if (!pseudoDown||!pseudoUp)
265      fprintf(Stderr,
266              "Only one set of pseudocosts - assumed same\n");
267    saveInfo->pseudoDown= (double *) malloc(numberColumns*sizeof(double));
268    saveInfo->pseudoUp = (double *) malloc(numberColumns*sizeof(double));
269    for (i=0;i<numberColumns;i++) {
270      double valueD=0.0, valueU=0.0;
271      if (pseudoDown) {
272        valueD = pseudoDown[i];
273        if (valueD<0) {
274          badpseudo++;
275          valueD=0.0;
276        }
277      }
278      if (pseudoUp) {
279        valueU = pseudoUp[i];
280        if (valueU<0) {
281          badpseudo++;
282          valueU=0.0;
283        }
284      }
285      if (!valueD)
286        valueD=valueU;
287      if (!valueU)
288        valueU=valueD;
289      saveInfo->pseudoDown[i]=valueD;
290      saveInfo->pseudoUp[i]=valueU;
291    }
292    if (badpseudo)
293      fprintf(Stderr,
294              "Treating %d negative pseudoCosts as 0.0\n",badpseudo);
295  }
296}
297static void
298stat_map(int *stat, int n, int *map, int mx, const char *what)
299{
300  int bad, i, i1=0, j, j1=0;
301  static char badfmt[] = "Coin driver: %s[%d] = %d\n";
302 
303  for(i = bad = 0; i < n; i++) {
304    if ((j = stat[i]) >= 0 && j <= mx)
305      stat[i] = map[j];
306    else {
307      stat[i] = 0;
308      i1 = i;
309      j1 = j;
310      if (!bad++)
311        fprintf(Stderr, badfmt, what, i, j);
312    }
313  }
314  if (bad > 1) {
315    if (bad == 2)
316      fprintf(Stderr, badfmt, what, i1, j1);
317    else
318      fprintf(Stderr,
319              "Coin driver: %d messages about bad %s values suppressed.\n",
320              bad-1, what);
321  }
322}
323int
324readAmpl(ampl_info * info, int argc, char **argv, void ** coinModel)
325{
326  char *stub;
327  ograd *og;
328  int i;
329  SufDesc *csd;
330  SufDesc *rsd;
331  /*bool *basis, *lower;*/
332  /*double *LU, *c, lb, objadj, *rshift, *shift, t, ub, *x, *x0, *x1;*/
333  char * environment = getenv("cbc_options");
334  char tempBuffer[20];
335  double * obj;
336  double * columnLower;
337  double * columnUpper;
338  double * rowLower;
339  double * rowUpper;
340  char ** saveArgv=argv;
341  char fileName[1000];
342  if (argc>1)
343    strcpy(fileName,argv[1]);
344  else
345    fileName[0]='\0';
346  int saveArgc = argc;
347  if (info->numberRows!=-1234567)
348    memset(info,0,sizeof(ampl_info)); // overwrite unless magic number set
349  /* save so can be accessed by decodePhrase */
350  saveInfo = info;
351  info->numberArguments=0;
352  info->arguments=(char **) malloc(2*sizeof(char *));
353  info->arguments[info->numberArguments++]=strdup("ampl");
354  info->arguments[info->numberArguments++]=strdup("cbc");
355  asl = ASL_alloc(ASL_read_f);
356  stub = getstub(&argv, &Oinfo);
357  if (!stub)
358    usage_ASL(&Oinfo, 1);
359  nl = jac0dim(stub, 0);
360  suf_declare(suftab, sizeof(suftab)/sizeof(SufDecl));
361 
362  /* set A_vals to get the constraints column-wise (malloc so can be freed) */
363  A_vals = (double *) malloc(nzc*sizeof(double));
364  if (!A_vals) {
365    printf("no memory\n");
366    return 1;
367  }
368  /* say we want primal solution */
369  want_xpi0=1;
370  /* for basis info */
371  info->columnStatus = (int *) malloc(n_var*sizeof(int));
372  info->rowStatus = (int *) malloc(n_con*sizeof(int));
373  csd = suf_iput("sstatus", ASL_Sufkind_var, info->columnStatus);
374  rsd = suf_iput("sstatus", ASL_Sufkind_con, info->rowStatus);
375  // testosi parameter - if >= 10 then go in through coinModel
376  int testOsi=-1;
377  for (i=0;i<saveArgc;i++) {
378    if (!strncmp(saveArgv[i],"testosi",7)) {
379      testOsi = atoi(saveArgv[i+1]);
380      break;
381    }
382  }
383  if (!(nlvc+nlvo)||testOsi>=10) {
384    /* read linear model*/
385    f_read(nl,0);
386    // see if any sos
387    if (true) {
388      char *sostype;
389      int nsosnz, *sosbeg, *sosind, * sospri;
390      double *sosref;
391      int nsos;
392      int i = ASL_suf_sos_explict_free;
393      int copri[2], **p_sospri;
394      copri[0] = 0;
395      copri[1] = 0;
396      p_sospri = &sospri;
397      nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
398                     &sosbeg, &sosind, &sosref);
399      if (nsos) {
400        info->numberSos=nsos;
401        info->sosType = (char *) malloc(nsos);
402        info->sosPriority = (int *) malloc(nsos*sizeof(int));
403        info->sosStart = (int *) malloc((nsos+1)*sizeof(int));
404        info->sosIndices = (int *) malloc(nsosnz*sizeof(int));
405        info->sosReference = (double *) malloc(nsosnz*sizeof(double));
406        sos_kludge(nsos, sosbeg, sosref,sosind);
407        for (int i=0;i<nsos;i++) {
408          int ichar = sostype[i];
409          assert (ichar=='1'||ichar=='2');
410          info->sosType[i]=ichar-'0';
411        }       
412        memcpy(info->sosPriority,sospri,nsos*sizeof(int));
413        memcpy(info->sosStart,sosbeg,(nsos+1)*sizeof(int));
414        memcpy(info->sosIndices,sosind,nsosnz*sizeof(int));
415        memcpy(info->sosReference,sosref,nsosnz*sizeof(double));
416      }
417    }
418   
419    /*sos_finish(&specialOrderedInfo, 0, &j, 0, 0, 0, 0, 0);*/
420    Oinfo.uinfo = tempBuffer;
421    if (getopts(argv, &Oinfo))
422      return 1;
423    /* objective*/
424    obj = (double *) malloc(n_var*sizeof(double));
425    for (i=0;i<n_var;i++)
426      obj[i]=0.0;;
427    if (n_obj) {
428      for (og = Ograd[0];og;og = og->next)
429        obj[og->varno] = og->coef;
430    }
431    if (objtype[0])
432      info->direction=-1.0;
433    else
434      info->direction=1.0;
435    info->offset=objconst(0);
436    /* Column bounds*/
437    columnLower = (double *) malloc(n_var*sizeof(double));
438    columnUpper = (double *) malloc(n_var*sizeof(double));
439#define COIN_DBL_MAX DBL_MAX
440    for (i=0;i<n_var;i++) {
441      columnLower[i]=LUv[2*i];
442      if (columnLower[i]<= negInfinity)
443        columnLower[i]=-COIN_DBL_MAX;
444      columnUpper[i]=LUv[2*i+1];
445      if (columnUpper[i]>= Infinity)
446        columnUpper[i]=COIN_DBL_MAX;
447    }
448    /* Row bounds*/
449    rowLower = (double *) malloc(n_con*sizeof(double));
450    rowUpper = (double *) malloc(n_con*sizeof(double));
451    for (i=0;i<n_con;i++) {
452      rowLower[i]=LUrhs[2*i];
453      if (rowLower[i]<= negInfinity)
454        rowLower[i]=-COIN_DBL_MAX;
455      rowUpper[i]=LUrhs[2*i+1];
456      if (rowUpper[i]>= Infinity)
457        rowUpper[i]=COIN_DBL_MAX;
458    }
459    info->numberRows=n_con;
460    info->numberColumns=n_var;
461    info->numberElements=nzc;;
462    info->numberBinary=nbv;
463    info->numberIntegers=niv;
464    info->objective=obj;
465    info->rowLower=rowLower;
466    info->rowUpper=rowUpper;
467    info->columnLower=columnLower;
468    info->columnUpper=columnUpper;
469    info->starts=A_colstarts;
470    /*A_colstarts=NULL;*/
471    info->rows=A_rownos;
472    /*A_rownos=NULL;*/
473    info->elements=A_vals;
474    /*A_vals=NULL;*/
475    info->primalSolution=NULL;
476    /* put in primalSolution if exists */
477    if (X0) {
478      info->primalSolution=(double *) malloc(n_var*sizeof(double));
479      memcpy(info->primalSolution,X0,n_var*sizeof(double));
480    }
481    info->dualSolution=NULL;
482    if (niv+nbv>0)
483      mip_stuff(); // get any extra info
484    if ((!(niv+nbv)&&(csd->kind & ASL_Sufkind_input))
485        ||(rsd->kind & ASL_Sufkind_input)) {
486      /* convert status - need info on map */
487      static int map[] = {1, 3, 1, 1, 2, 1, 1};
488      stat_map(info->columnStatus, n_var, map, 6, "incoming columnStatus");
489      stat_map(info->rowStatus, n_con, map, 6, "incoming rowStatus");
490    } else {
491      /* all slack basis */
492      // leave status for output */
493#if 0
494      free(info->rowStatus);
495      info->rowStatus=NULL;
496      free(info->columnStatus);
497      info->columnStatus=NULL;
498#endif
499    }
500  } else {
501    // QP
502    // Add .nl if not there
503    if (!strstr(fileName,".nl"))
504      strcat(fileName,".nl");
505    CoinModel * model = new CoinModel(1,fileName,info);
506    if (model->numberRows()>0)
507      *coinModel=(void *) model;
508    Oinfo.uinfo = tempBuffer;
509    if (getopts(argv, &Oinfo))
510      return 1;
511    Oinfo.wantsol=1;
512    if (objtype[0])
513      info->direction=-1.0;
514    else
515      info->direction=1.0;
516    info->offset=objconst(0);
517    info->numberRows=n_con;
518    info->numberColumns=n_var;
519    info->numberElements=nzc;;
520    info->numberBinary=nbv;
521    info->numberIntegers=niv;
522    if (nbv+niv+nlvbi+nlvci+nlvoi>0) {
523      mip_stuff(); // get any extra info
524      if (info->cut) 
525        model->setCutMarker(info->numberRows,info->cut);
526      if (info->priorities)
527        model->setPriorities(info->numberColumns,info->priorities);
528    }
529  }
530  /* add -solve - unless something there already
531   - also check for sleep=yes */
532  {
533    int found=0;
534    int foundLog=0;
535    int foundSleep=0;
536    const char * something[]={"solve","branch","duals","primals","user"};
537    for (i=0;i<info->numberArguments;i++) {
538      unsigned int j;
539      const char * argument = info->arguments[i];
540      for (j=0;j<sizeof(something)/sizeof(char *);j++) {
541        const char * check = something[j];
542        if (!strncmp(argument,check,sizeof(check))) {
543          found=(int)(j+1);
544        } else if (!strncmp(argument,"log",3)) {
545          foundLog=1;
546        } else if (!strncmp(argument,"sleep",5)) {
547          foundSleep=1;
548        }
549      }
550    }
551    if (foundLog) {
552      /* print options etc */
553      for (i=0;i<saveArgc;i++)
554        printf("%s ",saveArgv[i]);
555      printf("\n");
556      if (environment)
557        printf("env %s\n",environment);
558      /*printf("%d rows %d columns %d elements\n",n_con,n_var,nzc);*/
559    }
560    if (!found) {
561      if (!strlen(algFound)) {
562        info->arguments=(char **) realloc(info->arguments,(info->numberArguments+1)*sizeof(char *));
563        info->arguments[info->numberArguments++]=strdup("-solve");
564      } else {
565        // use algorithm from keyword
566        info->arguments=(char **) realloc(info->arguments,(info->numberArguments+1)*sizeof(char *));
567        info->arguments[info->numberArguments++]=strdup(algFound);
568      }
569    }
570    if (foundSleep) {
571      /* let user copy .nl file */
572      fprintf(stderr,"You can copy .nl file %s for debug purposes or attach debugger\n",saveArgv[1]);
573      fprintf(stderr,"Type q to quit, anything else to continue\n");
574      char getChar = getc(stdin);
575      if (getChar=='q'||getChar=='Q')
576        exit(1);
577    }
578  }
579  /* add -quit */
580  info->arguments=(char **) realloc(info->arguments,(info->numberArguments+1)*sizeof(char *));
581  info->arguments[info->numberArguments++]=strdup("-quit");
582  return 0;
583}
584void freeArrays1(ampl_info * info)
585{
586  free(info->objective);
587  info->objective=NULL;
588  free(info->rowLower);
589  info->rowLower=NULL;
590  free(info->rowUpper);
591  info->rowUpper=NULL;
592  free(info->columnLower);
593  info->columnLower=NULL;
594  free(info->columnUpper);
595  info->columnUpper=NULL;
596  /* this one not freed by ASL_free */
597  free(info->elements);
598  info->elements=NULL;
599  free(info->primalSolution);
600  info->primalSolution=NULL;
601  free(info->dualSolution);
602  info->dualSolution=NULL;
603  /*free(info->rowStatus);
604  info->rowStatus=NULL;
605  free(info->columnStatus);
606  info->columnStatus=NULL;*/
607}
608void freeArrays2(ampl_info * info)
609{
610  free(info->primalSolution);
611  info->primalSolution=NULL;
612  free(info->dualSolution);
613  info->dualSolution=NULL;
614  free(info->rowStatus);
615  info->rowStatus=NULL;
616  free(info->columnStatus);
617  info->columnStatus=NULL;
618  free(info->priorities);
619  info->priorities=NULL;
620  free(info->branchDirection);
621  info->branchDirection=NULL;
622  free(info->pseudoDown);
623  info->pseudoDown=NULL;
624  free(info->pseudoUp);
625  info->pseudoUp=NULL;
626  free(info->sosType);
627  info->sosType=NULL;
628  free(info->sosPriority);
629  info->sosPriority=NULL;
630  free(info->sosStart);
631  info->sosStart=NULL;
632  free(info->sosIndices);
633  info->sosIndices=NULL;
634  free(info->sosReference);
635  info->sosReference=NULL;
636  free(info->cut);
637  info->cut=NULL;
638  ASL_free(&asl);
639}
640void freeArgs(ampl_info * info)
641{
642  int i;
643  for ( i=0; i<info->numberArguments;i++)
644    free(info->arguments[i]);
645  free(info->arguments);
646}
647int ampl_obj_prec()
648{
649  return obj_prec();
650}
651void writeAmpl(ampl_info * info)
652{
653  char buf[1000];
654  typedef struct { const char *msg; int code; int wantObj; } Sol_info;
655  static Sol_info solinfo[] = {
656    { "optimal solution",                       000, 1 },
657    { "infeasible",                             200, 1 },
658    { "unbounded",                              300, 0 },
659    { "iteration limit etc",                    400, 1 },
660    { "solution limit",                         401, 1 },
661    { "ran out of space",                       500, 0 },
662    { "status unknown",                         501, 1 },
663    { "bug!",                                   502, 0 },
664    { "best MIP solution so far restored",      101, 1 },
665    { "failed to restore best MIP solution",    503, 1 },
666    { "optimal (?) solution",                   100, 1 }
667  };
668  /* convert status - need info on map */
669  static int map[] = {0, 3, 4, 1};
670  sprintf(buf,"%s %s",Oinfo.bsname,info->buffer);
671  solve_result_num = solinfo[info->problemStatus].code;
672  if (info->columnStatus) {
673    stat_map(info->columnStatus, n_var, map, 4, "outgoing columnStatus");
674    stat_map(info->rowStatus, n_con, map, 4, "outgoing rowStatus");
675    suf_iput("sstatus", ASL_Sufkind_var, info->columnStatus);
676    suf_iput("sstatus", ASL_Sufkind_con, info->rowStatus);
677  }
678  write_sol(buf,info->primalSolution,info->dualSolution,&Oinfo);
679}
680/* Read a problem from AMPL nl file
681 */
682CoinModel::CoinModel( int nonLinear, const char * fileName,const void * info)
683 :  numberRows_(0),
684    maximumRows_(0),
685    numberColumns_(0),
686    maximumColumns_(0),
687    numberElements_(0),
688    maximumElements_(0),
689    numberQuadraticElements_(0),
690    maximumQuadraticElements_(0),
691    optimizationDirection_(1.0),
692    objectiveOffset_(0.0),
693    rowLower_(NULL),
694    rowUpper_(NULL),
695    rowType_(NULL),
696    objective_(NULL),
697    columnLower_(NULL),
698    columnUpper_(NULL),
699    integerType_(NULL),
700    columnType_(NULL),
701    start_(NULL),
702    elements_(NULL),
703    quadraticElements_(NULL),
704    sortIndices_(NULL),
705    sortElements_(NULL),
706    sortSize_(0),
707    sizeAssociated_(0),
708    associated_(NULL),
709    numberSOS_(0),
710    startSOS_(NULL),
711    memberSOS_(NULL),
712    typeSOS_(NULL),
713    prioritySOS_(NULL),
714    referenceSOS_(NULL),
715    priority_(NULL),
716    cut_(NULL),
717    logLevel_(0),
718    type_(-1),
719    links_(0)
720{
721  problemName_ = strdup("");
722  int status=0;
723  if (!strcmp(fileName,"-")||!strcmp(fileName,"stdin")) {
724    // stdin
725  } else {
726    std::string name=fileName;
727    bool readable = fileCoinReadable(name);
728    if (!readable) {
729      std::cerr<<"Unable to open file "
730               <<fileName<<std::endl;
731      status = -1;
732    }
733  }
734  if (!status) {
735    gdb(nonLinear,fileName,info);
736  }
737}
738#if 0
739static real
740qterm(ASL *asl, fint *colq, fint *rowq, real *delsq)
741{
742  double t, t1, *x, *x0, *xe;
743  fint *rq0, *rqe;
744 
745  t = 0.;
746  x0 = x = X0;
747  xe = x + n_var;
748  rq0 = rowq;
749  while(x < xe) {
750    t1 = *x++;
751    rqe = rq0 + *++colq;
752    while(rowq < rqe)
753      t += t1*x0[*rowq++]**delsq++;
754  }
755  return 0.5 * t;
756}
757#endif
758// stolen from IPopt with changes
759typedef struct {
760  double obj_sign_;
761  ASL_pfgh * asl_;
762  double * non_const_x_;
763  int nz_h_full_; // number of nonzeros in hessian
764  int nerror_;
765  bool objval_called_with_current_x_;
766  bool conval_called_with_current_x_;
767} CbcAmplInfo;
768
769static bool get_nlp_info(void * amplInfo,int & n, int & m, int & nnz_jac_g,
770                              int & nnz_h_lag)
771{
772  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
773  ASL_pfgh* asl = info->asl_;
774 
775  n = n_var; // # of variables
776  m = n_con; // # of constraints
777  nnz_jac_g = nzc; // # of non-zeros in the jacobian
778  nnz_h_lag = info->nz_h_full_; // # of non-zeros in the hessian
779 
780  return true;
781}
782
783static bool get_bounds_info(void * amplInfo,int  n, double * x_l, 
784                            double * x_u, int  m, double * g_l, double * g_u)
785{
786  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
787  ASL_pfgh* asl = info->asl_;
788  assert(n == n_var);
789  assert(m == n_con);
790  int i;
791  for (i=0; i<n; i++) {
792    x_l[i] = LUv[2*i];
793    x_u[i] = LUv[2*i+1];
794  }
795 
796  for (i=0; i<m; i++) {
797    g_l[i] = LUrhs[2*i];
798    g_u[i] = LUrhs[2*i+1];
799  }
800  return true;
801}
802
803
804bool get_constraints_linearity(void * amplInfo,int  n,
805      int * const_types)
806{
807  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
808  ASL_pfgh* asl = info->asl_;
809  //check that n is good
810  assert(n == n_con);
811  // check that there are no network constraints
812  assert(nlnc == 0 && lnc == 0);
813  //the first nlc constraints are non linear the rest is linear
814  int i;
815  for (i=0; i<nlc; i++) {
816    const_types[i]=1;
817  }
818  // the rest is linear
819  for (i=nlc; i<n_con; i++)
820    const_types[i]=0;
821  return true;
822}
823#if 0
824bool get_starting_point(int  n, bool init_x, double * x, bool init_z,
825                        double * z_L, double * z_U, int  m, bool init_lambda, double * lambda)
826{
827  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
828  ASL_pfgh* asl = info->asl_;
829  assert(n == n_var);
830  assert(m == n_con);
831  int i;
832 
833  if (init_x) {
834    for (i=0; i<n; i++) {
835      if (havex0[i]) {
836        x[i] = X0[i];
837      }
838      else {
839        x[i] = 0.0;
840      }
841    }
842  }
843 
844  if (init_z) {
845    for (i=0; i<n; i++) {
846      z_L[i] = z_U[i] = 1.0;
847    }
848  }
849 
850  if (init_lambda) {
851    for (i=0; i<m; i++) {
852      if (havepi0[i]) {
853        lambda[i] = pi0[i];
854      }
855      else {
856        lambda[i] = 0.0;
857      }
858    }
859  }
860 
861  return true;
862}
863#endif
864static bool internal_objval(CbcAmplInfo * info ,double & obj_val)
865{
866  ASL_pfgh* asl = info->asl_;
867  info->objval_called_with_current_x_ = false; // in case the call below fails
868
869  if (n_obj==0) {
870    obj_val = 0;
871    info->objval_called_with_current_x_ = true;
872    return true;
873  }  else {
874    double  retval = objval(0, info->non_const_x_, (fint*)info->nerror_);
875    if (!info->nerror_) {
876      obj_val = info->obj_sign_*retval;
877      info->objval_called_with_current_x_ = true;
878      return true;
879    } else {
880      abort();
881    }
882  }
883 
884  return false;
885}
886
887static bool internal_conval(CbcAmplInfo * info ,int  m, double * g)
888{
889  ASL_pfgh* asl = info->asl_;
890  info->conval_called_with_current_x_ = false; // in case the call below fails
891 
892  bool allocated = false;
893  if (!g) {
894    g = new double[m];
895    allocated = true;
896  }
897 
898  conval(info->non_const_x_, g, (fint*)info->nerror_);
899
900  if (allocated) {
901    delete [] g;
902    g = NULL;
903  }
904 
905  if (!info->nerror_) {
906    info->conval_called_with_current_x_ = true;
907    return true;
908  } else {
909    abort();
910  }
911  return false;
912}
913
914
915static bool apply_new_x(CbcAmplInfo * info  ,bool new_x, int  n, const double * x)
916{
917  ASL_pfgh* asl = info->asl_;
918 
919  if (new_x) {
920    // update the flags so these methods are called
921    // before evaluating the hessian
922    info->conval_called_with_current_x_ = false;
923    info->objval_called_with_current_x_ = false;
924
925    //copy the data to the non_const_x_
926    if (!info->non_const_x_) {
927      info->non_const_x_ = new double [n];
928    }
929
930    for (int  i=0; i<n; i++) {
931      info->non_const_x_[i] = x[i];
932    }
933   
934    // tell ampl that we have a new x
935    xknowne(info->non_const_x_, (fint*)info->nerror_);
936    return info->nerror_ ? false : true;
937  }
938 
939  return true;
940}
941bool eval_f(void * amplInfo,int  n, const double * x, bool new_x, double & obj_value)
942{
943  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
944  if (!apply_new_x(info,new_x, n, x)) {
945    return false;
946  }
947 
948  return internal_objval(info,obj_value);
949}
950
951bool eval_grad_f(void * amplInfo,int  n, const double * x, bool new_x, double * grad_f)
952{
953  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
954  ASL_pfgh* asl = info->asl_;
955  if (!apply_new_x(info,new_x, n, x)) {
956    return false;
957  }
958  int i;
959 
960  if (n_obj==0) {
961    for (i=0; i<n; i++) {
962      grad_f[i] = 0.;
963    }
964  }
965  else {
966    objgrd(0, info->non_const_x_, grad_f, (fint*)info->nerror_);
967    if (info->nerror_) {
968      return false;
969    }
970   
971    if (info->obj_sign_==-1) {
972      for (i=0; i<n; i++) {
973        grad_f[i] = -grad_f[i];
974      }
975    }
976  }
977  return true;
978}
979
980bool eval_g(void * amplInfo,int  n, const double * x, bool new_x, int  m, double * g)
981{
982  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
983  ASL_pfgh* asl = info->asl_;
984  assert(n == n_var);
985  assert(m == n_con);
986 
987  if (!apply_new_x(info,new_x, n, x)) {
988    return false;
989  }
990 
991  return internal_conval(info,m, g);
992}
993
994bool eval_jac_g(void * amplInfo,int  n, const double * x, bool new_x,
995                int  m, int  nele_jac, int * iRow,
996                int  *jCol, double * values)
997{
998  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
999  ASL_pfgh* asl = info->asl_;
1000  assert(n == n_var);
1001  assert(m == n_con);
1002  int i;
1003 
1004  if (iRow && jCol && !values) {
1005    // setup the structure
1006    int  current_nz = 0;
1007    for (i=0; i<n_con; i++) {
1008      for (cgrad* cg=Cgrad[i]; cg; cg = cg->next) {
1009        //iRow[cg->goff] = i + 1;
1010        iRow[cg->goff] = i ;
1011        jCol[cg->goff] = cg->varno + 1;
1012        //                              iRow[current_nz] = i + 1;
1013        //                              jCol[current_nz] = cg->varno+1;
1014        current_nz++;
1015      }
1016    }
1017    assert(current_nz == nele_jac);
1018    return true;
1019  }
1020  else if (!iRow && !jCol && values) {
1021    if (!apply_new_x(info,new_x, n, x)) {
1022      return false;
1023    }
1024   
1025    jacval(info->non_const_x_, values, (fint*)info->nerror_);
1026    if (!info->nerror_) {
1027      return true;
1028    } else {
1029      abort();
1030    }
1031  }
1032  else {
1033    assert(false && "Invalid combination of iRow, jCol, and values pointers");
1034  }
1035 
1036  return false;
1037}
1038
1039bool eval_h(void * amplInfo,int  n, const double * x, bool new_x,
1040            double  obj_factor, int  m, const double * lambda,
1041            bool new_lambda, int  nele_hess, int * iRow,
1042            int * jCol, double * values)
1043{
1044  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
1045  ASL_pfgh* asl = info->asl_;
1046  assert(n == n_var);
1047  assert(m == n_con);
1048  int i;
1049 
1050  if (iRow && jCol && !values) {
1051    // setup the structure
1052    int k=0;
1053    for (int i=0; i<n; i++) {
1054      for (int j=sputinfo->hcolstarts[i]; j<sputinfo->hcolstarts[i+1]; j++) {
1055        //iRow[k] = i + 1;
1056        iRow[k] = i;
1057        jCol[k] = sputinfo->hrownos[j]+1;
1058        k++;
1059      }
1060    }
1061    assert(k==nele_hess);
1062    return true;
1063  }
1064  else if (!iRow & !jCol && values) {
1065    if (!apply_new_x(info,new_x, n, x)) {
1066      return false;
1067    }
1068    if (!info->objval_called_with_current_x_) {
1069      double  dummy;
1070      internal_objval(info,dummy);
1071      internal_conval(info,m,NULL);
1072    }
1073    if (!info->conval_called_with_current_x_) {
1074      internal_conval(info,m,NULL);
1075    }
1076    // copy lambda to non_const_lambda - note, we do not store a copy like
1077    // we do with x since lambda is only used here and not in other calls
1078    double * non_const_lambda = new double [m];
1079    for (i=0; i<m; i++) {
1080      non_const_lambda[i] = lambda[i];
1081    }
1082   
1083    real ow=info->obj_sign_*obj_factor;
1084    sphes(values, -1, &ow, non_const_lambda);
1085   
1086    delete [] non_const_lambda;
1087    return true;
1088  }
1089  else {
1090    assert(false && "Invalid combination of iRow, jCol, and values pointers");
1091  }
1092 
1093  return false;
1094}
1095void
1096CoinModel::gdb( int nonLinear, const char * fileName,const void * info)
1097{
1098  const ampl_info * amplInfo = (const ampl_info *) info;
1099  ograd *og=NULL;
1100  int i;
1101  SufDesc *csd=NULL;
1102  SufDesc *rsd=NULL;
1103  /*bool *basis, *lower;*/
1104  /*double *LU, *c, lb, objadj, *rshift, *shift, t, ub, *x, *x0, *x1;*/
1105  //char tempBuffer[20];
1106  double * objective=NULL;
1107  double * columnLower=NULL;
1108  double * columnUpper=NULL;
1109  double * rowLower=NULL;
1110  double * rowUpper=NULL;
1111  int * columnStatus=NULL;
1112  int * rowStatus=NULL;
1113  int numberRows=-1;
1114  int numberColumns=-1;
1115  int numberElements=-1;
1116  int numberBinary=-1;
1117  int numberIntegers=-1;
1118  int numberAllNonLinearBoth=0;
1119  int numberIntegerNonLinearBoth=0;
1120  int numberAllNonLinearConstraints=0;
1121  int numberIntegerNonLinearConstraints=0;
1122  int numberAllNonLinearObjective=0;
1123  int numberIntegerNonLinearObjective=0;
1124  double * primalSolution=NULL;
1125  double direction=1.0;
1126  char * stub = strdup(fileName);
1127  CoinPackedMatrix matrixByRow;
1128  fint ** colqp=NULL;
1129  int *z = NULL;
1130  if (nonLinear==0) {
1131    // linear
1132    asl = ASL_alloc(ASL_read_f);
1133    nl = jac0dim(stub, 0);
1134    free(stub);
1135    suf_declare(suftab, sizeof(suftab)/sizeof(SufDecl));
1136   
1137    /* set A_vals to get the constraints column-wise (malloc so can be freed) */
1138    A_vals = (double *) malloc(nzc*sizeof(double));
1139    if (!A_vals) {
1140      printf("no memory\n");
1141      return ;
1142    }
1143    /* say we want primal solution */
1144    want_xpi0=1;
1145    /* for basis info */
1146    columnStatus = (int *) malloc(n_var*sizeof(int));
1147    rowStatus = (int *) malloc(n_con*sizeof(int));
1148    csd = suf_iput("sstatus", ASL_Sufkind_var, columnStatus);
1149    rsd = suf_iput("sstatus", ASL_Sufkind_con, rowStatus);
1150    /* read linear model*/
1151    f_read(nl,0);
1152    // see if any sos
1153    if (true) {
1154      char *sostype;
1155      int nsosnz, *sosbeg, *sosind, * sospri;
1156      double *sosref;
1157      int nsos;
1158      int i = ASL_suf_sos_explict_free;
1159      int copri[2], **p_sospri;
1160      copri[0] = 0;
1161      copri[1] = 0;
1162      p_sospri = &sospri;
1163      nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
1164                     &sosbeg, &sosind, &sosref);
1165      if (nsos) {
1166        abort();
1167#if 0
1168        info->numberSos=nsos;
1169        info->sosType = (char *) malloc(nsos);
1170        info->sosPriority = (int *) malloc(nsos*sizeof(int));
1171        info->sosStart = (int *) malloc((nsos+1)*sizeof(int));
1172        info->sosIndices = (int *) malloc(nsosnz*sizeof(int));
1173        info->sosReference = (double *) malloc(nsosnz*sizeof(double));
1174        sos_kludge(nsos, sosbeg, sosref,sosind);
1175        for (int i=0;i<nsos;i++) {
1176          int ichar = sostype[i];
1177          assert (ichar=='1'||ichar=='2');
1178          info->sosType[i]=ichar-'0';
1179        }       
1180        memcpy(info->sosPriority,sospri,nsos*sizeof(int));
1181        memcpy(info->sosStart,sosbeg,(nsos+1)*sizeof(int));
1182        memcpy(info->sosIndices,sosind,nsosnz*sizeof(int));
1183        memcpy(info->sosReference,sosref,nsosnz*sizeof(double));
1184#endif
1185      }
1186    }
1187   
1188    /*sos_finish(&specialOrderedInfo, 0, &j, 0, 0, 0, 0, 0);*/
1189    //Oinfo.uinfo = tempBuffer;
1190    //if (getopts(argv, &Oinfo))
1191    //return 1;
1192    /* objective*/
1193    objective = (double *) malloc(n_var*sizeof(double));
1194    for (i=0;i<n_var;i++)
1195      objective[i]=0.0;;
1196    if (n_obj) {
1197      for (og = Ograd[0];og;og = og->next)
1198        objective[og->varno] = og->coef;
1199    }
1200    if (objtype[0])
1201      direction=-1.0;
1202    else
1203      direction=1.0;
1204    objectiveOffset_=objconst(0);
1205    /* Column bounds*/
1206    columnLower = (double *) malloc(n_var*sizeof(double));
1207    columnUpper = (double *) malloc(n_var*sizeof(double));
1208#define COIN_DBL_MAX DBL_MAX
1209    for (i=0;i<n_var;i++) {
1210      columnLower[i]=LUv[2*i];
1211      if (columnLower[i]<= negInfinity)
1212        columnLower[i]=-COIN_DBL_MAX;
1213      columnUpper[i]=LUv[2*i+1];
1214      if (columnUpper[i]>= Infinity)
1215        columnUpper[i]=COIN_DBL_MAX;
1216    }
1217    /* Row bounds*/
1218    rowLower = (double *) malloc(n_con*sizeof(double));
1219    rowUpper = (double *) malloc(n_con*sizeof(double));
1220    for (i=0;i<n_con;i++) {
1221      rowLower[i]=LUrhs[2*i];
1222      if (rowLower[i]<= negInfinity)
1223        rowLower[i]=-COIN_DBL_MAX;
1224      rowUpper[i]=LUrhs[2*i+1];
1225      if (rowUpper[i]>= Infinity)
1226        rowUpper[i]=COIN_DBL_MAX;
1227    }
1228    numberRows=n_con;
1229    numberColumns=n_var;
1230    numberElements=nzc;;
1231    numberBinary=nbv;
1232    numberIntegers=niv;
1233    /* put in primalSolution if exists */
1234    if (X0) {
1235      primalSolution=(double *) malloc(n_var*sizeof(double));
1236      memcpy( primalSolution,X0,n_var*sizeof(double));
1237    }
1238    //double * dualSolution=NULL;
1239    if (niv+nbv>0)
1240      mip_stuff(); // get any extra info
1241    if ((!(niv+nbv)&&(csd->kind & ASL_Sufkind_input))
1242        ||(rsd->kind & ASL_Sufkind_input)) {
1243      /* convert status - need info on map */
1244      static int map[] = {1, 3, 1, 1, 2, 1, 1};
1245      stat_map(columnStatus, n_var, map, 6, "incoming columnStatus");
1246      stat_map(rowStatus, n_con, map, 6, "incoming rowStatus");
1247    } else {
1248      /* all slack basis */
1249      // leave status for output */
1250#if 0
1251      free(rowStatus);
1252      rowStatus=NULL;
1253      free(columnStatus);
1254      columnStatus=NULL;
1255#endif
1256    }
1257    CoinPackedMatrix columnCopy(true,numberRows,numberColumns,numberElements,
1258                                A_vals,A_rownos,A_colstarts,NULL);
1259    matrixByRow.reverseOrderedCopyOf(columnCopy);
1260  } else if (nonLinear==1) {
1261    // quadratic
1262    asl = ASL_alloc(ASL_read_fg);
1263    nl = jac0dim(stub, (long) strlen(stub));
1264    free(stub);
1265    suf_declare(suftab, sizeof(suftab)/sizeof(SufDecl));
1266    /* read  model*/
1267    X0 = (double*) malloc(n_var*sizeof(double));
1268    CoinZeroN(X0,n_var);
1269    qp_read(nl,0);
1270    assert (n_obj==1);
1271    int nz = 1 + n_con;
1272    colqp = (fint**) malloc(nz*(2*sizeof(int*)
1273                                      + sizeof(double*)));
1274    fint ** rowqp = colqp + nz;
1275    double ** delsqp = (double **)(rowqp + nz);
1276    z = (int*) malloc(nz*sizeof(int));
1277    for (i=0;i<=n_con;i++) {
1278      z[i] = nqpcheck(-i, rowqp+i, colqp+i, delsqp+i);
1279    }
1280    qp_opify();
1281    /* objective*/
1282    objective = (double *) malloc(n_var*sizeof(double));
1283    for (i=0;i<n_var;i++)
1284      objective[i]=0.0;;
1285    if (n_obj) {
1286      for (og = Ograd[0];og;og = og->next)
1287        objective[og->varno] = og->coef;
1288    }
1289    if (objtype[0])
1290      direction=-1.0;
1291    else
1292      direction=1.0;
1293    objectiveOffset_=objconst(0);
1294    /* Column bounds*/
1295    columnLower = (double *) malloc(n_var*sizeof(double));
1296    columnUpper = (double *) malloc(n_var*sizeof(double));
1297    for (i=0;i<n_var;i++) {
1298      columnLower[i]=LUv[2*i];
1299      if (columnLower[i]<= negInfinity)
1300        columnLower[i]=-COIN_DBL_MAX;
1301      columnUpper[i]=LUv[2*i+1];
1302      if (columnUpper[i]>= Infinity)
1303        columnUpper[i]=COIN_DBL_MAX;
1304    }
1305    // Build by row from scratch
1306    //matrixByRow.reserve(n_var,nzc,true);
1307    // say row orderded
1308    matrixByRow.transpose();
1309    /* Row bounds*/
1310    rowLower = (double *) malloc(n_con*sizeof(double));
1311    rowUpper = (double *) malloc(n_con*sizeof(double));
1312    int * column = new int [n_var];
1313    double * element = new double [n_var];
1314    for (i=0;i<n_con;i++) {
1315      rowLower[i]=LUrhs[2*i];
1316      if (rowLower[i]<= negInfinity)
1317        rowLower[i]=-COIN_DBL_MAX;
1318      rowUpper[i]=LUrhs[2*i+1];
1319      if (rowUpper[i]>= Infinity)
1320        rowUpper[i]=COIN_DBL_MAX;
1321      int k=0;
1322      for(cgrad * cg = Cgrad[i]; cg; cg = cg->next) {
1323        column[k]=cg->varno;
1324        element[k++]=cg->coef;
1325      }
1326      matrixByRow.appendRow(k,column,element);
1327    }
1328    delete [] column;
1329    delete [] element;
1330    numberRows=n_con;
1331    numberColumns=n_var;
1332    numberElements=nzc;;
1333    numberBinary=nbv;
1334    numberIntegers=niv;
1335    numberAllNonLinearBoth=nlvb;
1336    numberIntegerNonLinearBoth=nlvbi;
1337    numberAllNonLinearConstraints=nlvc;
1338    numberIntegerNonLinearConstraints=nlvci;
1339    numberAllNonLinearObjective=nlvo;
1340    numberIntegerNonLinearObjective=nlvoi;
1341    /* say we want primal solution */
1342    want_xpi0=1;
1343    //double * dualSolution=NULL;
1344    // save asl
1345    // Fix memory leak one day
1346    CbcAmplInfo * info = new CbcAmplInfo;
1347    amplGamsData_ = info;
1348    info->asl_=NULL; // as wrong form asl;
1349    info->nz_h_full_=-1; // number of nonzeros in hessian
1350    info->objval_called_with_current_x_=false;
1351    info->nerror_=0;
1352    info->obj_sign_=direction;
1353    info->conval_called_with_current_x_=false;
1354    info->non_const_x_=NULL;
1355  } else if (nonLinear==2) {
1356    // General nonlinear!
1357    ASL_pfgh* asl = (ASL_pfgh*)ASL_alloc(ASL_read_pfgh);
1358    // was asl = ASL_alloc(ASL_read_fg);
1359    nl = jac0dim(stub, (long) strlen(stub));
1360    free(stub);
1361    suf_declare(suftab, sizeof(suftab)/sizeof(SufDecl));
1362    /* read  model*/
1363    X0 = (double*) malloc(n_var*sizeof(double));
1364    CoinZeroN(X0,n_var);
1365    // code stolen from Ipopt
1366    int retcode = pfgh_read(nl, ASL_return_read_err | ASL_findgroups);
1367
1368    switch (retcode) {
1369      case ASL_readerr_none : {}
1370      break;
1371      case ASL_readerr_nofile : {
1372        printf( "Cannot open .nl file\n");
1373        exit(-1);
1374      }
1375      break;
1376      case ASL_readerr_nonlin : {
1377        assert(false); // this better not be an error!
1378        printf( "model involves nonlinearities (ed0read)\n");
1379        exit(-1);
1380      }
1381      break;
1382      case  ASL_readerr_argerr : {
1383        printf( "user-defined function with bad args\n");
1384        exit(-1);
1385      }
1386      break;
1387      case ASL_readerr_unavail : {
1388        printf( "user-defined function not available\n");
1389        exit(-1);
1390      }
1391      break;
1392      case ASL_readerr_corrupt : {
1393        printf( "corrupt .nl file\n");
1394        exit(-1);
1395      }
1396      break;
1397      case ASL_readerr_bug : {
1398        printf( "bug in .nl reader\n");
1399        exit(-1);
1400      }
1401      break;
1402      case ASL_readerr_CLP : {
1403        printf( "ASL error message: \"solver cannot handle CLP extensions\"\n");
1404        exit(-1);
1405      }
1406      break;
1407      default: {
1408        printf( "Unknown error in stub file read. retcode = %d\n", retcode);
1409        exit(-1);
1410      }
1411      break;
1412    }
1413
1414    // see "changes" in solvers directory of ampl code...
1415    hesset(1,0,1,0,nlc);
1416
1417    assert (n_obj==1);
1418    // find the nonzero structure for the hessian
1419    // parameters to sphsetup:
1420    int coeff_obj = 1; // coefficient of the objective fn ???
1421    int mult_supplied = 1; // multipliers will be supplied
1422    int uptri = 1; // only need the upper triangular part
1423    // save asl
1424    // Fix memory leak one day
1425    CbcAmplInfo * info = new CbcAmplInfo;
1426    amplGamsData_ = info;
1427    info->asl_=asl;
1428    // This is not easy to get from ampl so save
1429    info->nz_h_full_= sphsetup(-1, coeff_obj, mult_supplied, uptri);
1430    info->objval_called_with_current_x_=false;
1431    info->nerror_=0;
1432    info->obj_sign_=direction;
1433    info->conval_called_with_current_x_=false;
1434    info->non_const_x_=NULL;
1435    /* objective*/
1436    objective = (double *) malloc(n_var*sizeof(double));
1437    for (i=0;i<n_var;i++)
1438      objective[i]=0.0;;
1439    if (n_obj) {
1440      for (og = Ograd[0];og;og = og->next)
1441        objective[og->varno] = og->coef;
1442    }
1443    if (objtype[0])
1444      direction=-1.0;
1445    else
1446      direction=1.0;
1447    objectiveOffset_=objconst(0);
1448    /* Column bounds*/
1449    columnLower = (double *) malloc(n_var*sizeof(double));
1450    columnUpper = (double *) malloc(n_var*sizeof(double));
1451    for (i=0;i<n_var;i++) {
1452      columnLower[i]=LUv[2*i];
1453      if (columnLower[i]<= negInfinity)
1454        columnLower[i]=-COIN_DBL_MAX;
1455      columnUpper[i]=LUv[2*i+1];
1456      if (columnUpper[i]>= Infinity)
1457        columnUpper[i]=COIN_DBL_MAX;
1458    }
1459    // Build by row from scratch
1460    //matrixByRow.reserve(n_var,nzc,true);
1461    // say row orderded
1462    matrixByRow.transpose();
1463    /* Row bounds*/
1464    rowLower = (double *) malloc(n_con*sizeof(double));
1465    rowUpper = (double *) malloc(n_con*sizeof(double));
1466    int * column = new int [n_var];
1467    double * element = new double [n_var];
1468    for (i=0;i<n_con;i++) {
1469      rowLower[i]=LUrhs[2*i];
1470      if (rowLower[i]<= negInfinity)
1471        rowLower[i]=-COIN_DBL_MAX;
1472      rowUpper[i]=LUrhs[2*i+1];
1473      if (rowUpper[i]>= Infinity)
1474        rowUpper[i]=COIN_DBL_MAX;
1475      int k=0;
1476      for(cgrad * cg = Cgrad[i]; cg; cg = cg->next) {
1477        column[k]=cg->varno;
1478        double value = cg->coef;
1479        if (!value)
1480          value = -1.2345e-29;
1481        element[k++]=value;
1482      }
1483      matrixByRow.appendRow(k,column,element);
1484    }
1485    delete [] column;
1486    delete [] element;
1487    numberRows=n_con;
1488    numberColumns=n_var;
1489    numberElements=nzc;;
1490    numberBinary=nbv;
1491    numberIntegers=niv;
1492    numberAllNonLinearBoth=nlvb;
1493    numberIntegerNonLinearBoth=nlvbi;
1494    numberAllNonLinearConstraints=nlvc;
1495    numberIntegerNonLinearConstraints=nlvci;
1496    numberAllNonLinearObjective=nlvo;
1497    numberIntegerNonLinearObjective=nlvoi;
1498    /* say we want primal solution */
1499    want_xpi0=1;
1500    //double * dualSolution=NULL;
1501  } else {
1502    abort();
1503  }
1504  // set problem name
1505  free(problemName_);
1506  problemName_=strdup("???");
1507
1508  // Build by row from scratch
1509  const double * element = matrixByRow.getElements();
1510  const int * column = matrixByRow.getIndices();
1511  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
1512  const int * rowLength = matrixByRow.getVectorLengths();
1513  for (i=0;i<numberRows;i++) {
1514    addRow(rowLength[i],column+rowStart[i],
1515           element+rowStart[i],rowLower[i],rowUpper[i]);
1516  }
1517  // Now do column part
1518  for (i=0;i<numberColumns;i++) {
1519    setColumnBounds(i,columnLower[i],columnUpper[i]);
1520    setColumnObjective(i,objective[i]);
1521  }
1522  for ( i=numberColumns-numberBinary-numberIntegers;
1523        i<numberColumns;i++) {
1524      setColumnIsInteger(i,true);;
1525  }
1526  // and non linear
1527  for (i=numberAllNonLinearBoth-numberIntegerNonLinearBoth;
1528       i<numberAllNonLinearBoth;i++) {
1529    setColumnIsInteger(i,true);;
1530  }
1531  for (i=numberAllNonLinearConstraints-numberIntegerNonLinearConstraints;
1532       i<numberAllNonLinearConstraints;i++) {
1533    setColumnIsInteger(i,true);;
1534  }
1535  for (i=numberAllNonLinearObjective-numberIntegerNonLinearObjective;
1536       i<numberAllNonLinearObjective;i++) {
1537    setColumnIsInteger(i,true);;
1538  }
1539  free(columnLower);
1540  free(columnUpper);
1541  free(rowLower);
1542  free(rowUpper);
1543  free(objective);
1544  // do names
1545  int iRow;
1546  for (iRow=0;iRow<numberRows_;iRow++) {
1547    char name[9];
1548    sprintf(name,"r%7.7d",iRow);
1549    setRowName(iRow,name);
1550  }
1551  int iColumn;
1552  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1553    char name[9];
1554    sprintf(name,"c%7.7d",iColumn);
1555    setColumnName(iColumn,name);
1556  }
1557  if (colqp) {
1558    // add in quadratic
1559    int nz = 1 + n_con;
1560    int nOdd=0;
1561    fint ** rowqp = colqp + nz;
1562    double ** delsqp = (double **)(rowqp + nz);
1563    for (i=0;i<=n_con;i++) {
1564      int nels = z[i];
1565      if (nels) {
1566        double * element = delsqp[i];
1567        int * start = (int *) colqp[i];
1568        int * row = (int *) rowqp[i];
1569        if (!element) {
1570          // odd row - probably not quadratic
1571          nOdd++;
1572          continue;
1573        }
1574#if 0
1575        printf("%d quadratic els\n",nels);
1576        for (int j=0;j<n_var;j++) {
1577          for (int k=start[j];k<start[j+1];k++)
1578            printf("%d %d %g\n",j,row[k],element[k]);
1579        }
1580#endif
1581        if (i) {
1582          int iRow = i-1;
1583          for (int j=0;j<n_var;j++) {
1584            for (int k=start[j];k<start[j+1];k++) {
1585              int kColumn = row[k];
1586              double value = element[k];
1587              // ampl gives twice with assumed 0.5
1588              if (kColumn<j)
1589                continue;
1590              else if (kColumn==j)
1591                value *= 0.5;
1592              const char * expr = getElementAsString(iRow,j);
1593              double constant=0.0;
1594              bool linear;
1595              if (expr&&strcmp(expr,"Numeric")) {
1596                linear=false;
1597              } else {
1598                constant = getElement(iRow,j);
1599                linear=true;
1600              }
1601              char temp[1000];
1602              char temp2[30];
1603              if (value==1.0) 
1604                sprintf(temp2,"c%7.7d",kColumn);
1605              else
1606                sprintf(temp2,"%g*c%7.7d",value,kColumn);
1607              if (linear) {
1608                if (!constant)
1609                  strcpy(temp,temp2);
1610                else if (value>0.0) 
1611                  sprintf(temp,"%g+%s",constant,temp2);
1612                else
1613                  sprintf(temp,"%g%s",constant,temp2);
1614              } else {
1615                if (value>0.0) 
1616                  sprintf(temp,"%s+%s",expr,temp2);
1617                else
1618                  sprintf(temp,"%s%s",expr,temp2);
1619              }
1620              assert (strlen(temp)<1000);
1621              setElement(iRow,j,temp);
1622              if (amplInfo->logLevel>1)
1623                printf("el for row %d column c%7.7d is %s\n",iRow,j,temp);
1624            }
1625          }
1626        } else {
1627          // objective
1628          for (int j=0;j<n_var;j++) {
1629            for (int k=start[j];k<start[j+1];k++) {
1630              int kColumn = row[k];
1631              double value = element[k];
1632              // ampl gives twice with assumed 0.5
1633              if (kColumn<j)
1634                continue;
1635              else if (kColumn==j)
1636                value *= 0.5;
1637              const char * expr = getColumnObjectiveAsString(j);
1638              double constant=0.0;
1639              bool linear;
1640              if (expr&&strcmp(expr,"Numeric")) {
1641                linear=false;
1642              } else {
1643                constant = getColumnObjective(j);
1644                linear=true;
1645              }
1646              char temp[1000];
1647              char temp2[30];
1648              if (value==1.0) 
1649                sprintf(temp2,"c%7.7d",kColumn);
1650              else
1651                sprintf(temp2,"%g*c%7.7d",value,kColumn);
1652              if (linear) {
1653                if (!constant)
1654                  strcpy(temp,temp2);
1655                else if (value>0.0) 
1656                  sprintf(temp,"%g+%s",constant,temp2);
1657                else
1658                  sprintf(temp,"%g%s",constant,temp2);
1659              } else {
1660                if (value>0.0) 
1661                  sprintf(temp,"%s+%s",expr,temp2);
1662                else
1663                  sprintf(temp,"%s%s",expr,temp2);
1664              }
1665              assert (strlen(temp)<1000);
1666              setObjective(j,temp);
1667              if (amplInfo->logLevel>1)
1668                printf("el for objective column c%7.7d is %s\n",j,temp);
1669            }
1670          }
1671        }
1672      }
1673    }
1674    if (nOdd) {
1675      printf("%d non-linear constraints could not be converted to quadratic\n",nOdd);
1676      exit(77);
1677    }
1678  }
1679  free(colqp);
1680  free(z);
1681  // see if any sos
1682  {
1683    char *sostype;
1684    int nsosnz, *sosbeg, *sosind, * sospri;
1685    double *sosref;
1686    int nsos;
1687    int i = ASL_suf_sos_explict_free;
1688    int copri[2], **p_sospri;
1689    copri[0] = 0;
1690    copri[1] = 0;
1691    p_sospri = &sospri;
1692    nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
1693                   &sosbeg, &sosind, &sosref);
1694    if (nsos) {
1695      numberSOS_=nsos;
1696      typeSOS_ = new int [numberSOS_];
1697      prioritySOS_ = new int [numberSOS_];
1698      startSOS_ = new int [numberSOS_+1];
1699      memberSOS_ = new int[nsosnz];
1700      referenceSOS_ = new double [nsosnz];
1701      sos_kludge(nsos, sosbeg, sosref,sosind);
1702      for (int i=0;i<nsos;i++) {
1703        int ichar = sostype[i];
1704        assert (ichar=='1'||ichar=='2');
1705        typeSOS_[i]=ichar-'0';
1706      } 
1707      memcpy(prioritySOS_,sospri,nsos*sizeof(int));
1708      memcpy(startSOS_,sosbeg,(nsos+1)*sizeof(int));
1709      memcpy(memberSOS_,sosind,nsosnz*sizeof(int));
1710      memcpy(referenceSOS_,sosref,nsosnz*sizeof(double));
1711    }
1712  }
1713}
1714#else
1715#include "Cbc_ampl.h"
1716int
1717readAmpl(ampl_info * info, int argc, char **argv, void ** coinModel)
1718{
1719  return 0;
1720}
1721void freeArrays1(ampl_info * info)
1722{
1723}
1724void freeArrays2(ampl_info * info)
1725{
1726}
1727void freeArgs(ampl_info * info)
1728{
1729}
1730int ampl_obj_prec()
1731{
1732  return 0;
1733}
1734void writeAmpl(ampl_info * info)
1735{
1736}
1737#endif
Note: See TracBrowser for help on using the repository browser.