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

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

take out misleading setting of integer variables because of ampl

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