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

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

for lou and callable cbcmain

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