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

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

update branches/devel for threads

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 41.9 KB
Line 
1/****************************************************************
2Copyright (C) 1997-2000 Lucent Technologies
3Modifications for Coin -  Copyright (C) 2006, International Business Machines Corporation and others.
4All Rights Reserved
5
6Permission to use, copy, modify, and distribute this software and
7its documentation for any purpose and without fee is hereby
8granted, provided that the above copyright notice appear in all
9copies and that both that the copyright notice and this
10permission notice and warranty disclaimer appear in supporting
11documentation, and that the name of Lucent or any of its entities
12not be used in advertising or publicity pertaining to
13distribution of the software without specific, written prior
14permission.
15
16LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
17INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
18IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
19SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
20WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
21IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
22ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
23THIS SOFTWARE.
24****************************************************************/
25#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    if (numberIntegers>0) {
535      mip_stuff(); // get any extra info
536      if (info->cut) 
537        model->setCutMarker(info->numberRows,info->cut);
538      if (info->priorities)
539        model->setPriorities(info->numberColumns,info->priorities);
540    }
541  }
542  /* add -solve - unless something there already
543   - also check for sleep=yes */
544  {
545    int found=0;
546    int foundLog=0;
547    int foundSleep=0;
548    const char * something[]={"solve","branch","duals","primals","user"};
549    for (i=0;i<info->numberArguments;i++) {
550      unsigned int j;
551      const char * argument = info->arguments[i];
552      for (j=0;j<sizeof(something)/sizeof(char *);j++) {
553        const char * check = something[j];
554        if (!strncmp(argument,check,sizeof(check))) {
555          found=(int)(j+1);
556        } else if (!strncmp(argument,"log",3)) {
557          foundLog=1;
558        } else if (!strncmp(argument,"sleep",5)) {
559          foundSleep=1;
560        }
561      }
562    }
563    if (foundLog) {
564      /* print options etc */
565      for (i=0;i<saveArgc;i++)
566        printf("%s ",saveArgv[i]);
567      printf("\n");
568      if (environment)
569        printf("env %s\n",environment);
570      /*printf("%d rows %d columns %d elements\n",n_con,n_var,nzc);*/
571    }
572    if (!found) {
573      if (!strlen(algFound)) {
574        info->arguments=(char **) realloc(info->arguments,(info->numberArguments+1)*sizeof(char *));
575        info->arguments[info->numberArguments++]=strdup("-solve");
576      } else {
577        // use algorithm from keyword
578        info->arguments=(char **) realloc(info->arguments,(info->numberArguments+1)*sizeof(char *));
579        info->arguments[info->numberArguments++]=strdup(algFound);
580      }
581    }
582    if (foundSleep) {
583      /* let user copy .nl file */
584      fprintf(stderr,"You can copy .nl file %s for debug purposes or attach debugger\n",saveArgv[1]);
585      fprintf(stderr,"Type q to quit, anything else to continue\n");
586      char getChar = getc(stdin);
587      if (getChar=='q'||getChar=='Q')
588        exit(1);
589    }
590  }
591  /* add -quit */
592  info->arguments=(char **) realloc(info->arguments,(info->numberArguments+1)*sizeof(char *));
593  info->arguments[info->numberArguments++]=strdup("-quit");
594  return 0;
595}
596void freeArrays1(ampl_info * info)
597{
598  free(info->objective);
599  info->objective=NULL;
600  free(info->rowLower);
601  info->rowLower=NULL;
602  free(info->rowUpper);
603  info->rowUpper=NULL;
604  free(info->columnLower);
605  info->columnLower=NULL;
606  free(info->columnUpper);
607  info->columnUpper=NULL;
608  /* this one not freed by ASL_free */
609  free(info->elements);
610  info->elements=NULL;
611  free(info->primalSolution);
612  info->primalSolution=NULL;
613  free(info->dualSolution);
614  info->dualSolution=NULL;
615  /*free(info->rowStatus);
616  info->rowStatus=NULL;
617  free(info->columnStatus);
618  info->columnStatus=NULL;*/
619}
620void freeArrays2(ampl_info * info)
621{
622  free(info->primalSolution);
623  info->primalSolution=NULL;
624  free(info->dualSolution);
625  info->dualSolution=NULL;
626  free(info->rowStatus);
627  info->rowStatus=NULL;
628  free(info->columnStatus);
629  info->columnStatus=NULL;
630  free(info->priorities);
631  info->priorities=NULL;
632  free(info->branchDirection);
633  info->branchDirection=NULL;
634  free(info->pseudoDown);
635  info->pseudoDown=NULL;
636  free(info->pseudoUp);
637  info->pseudoUp=NULL;
638  free(info->sosType);
639  info->sosType=NULL;
640  free(info->sosPriority);
641  info->sosPriority=NULL;
642  free(info->sosStart);
643  info->sosStart=NULL;
644  free(info->sosIndices);
645  info->sosIndices=NULL;
646  free(info->sosReference);
647  info->sosReference=NULL;
648  free(info->cut);
649  info->cut=NULL;
650  ASL_free(&asl);
651}
652void freeArgs(ampl_info * info)
653{
654  int i;
655  for ( i=0; i<info->numberArguments;i++)
656    free(info->arguments[i]);
657  free(info->arguments);
658}
659int ampl_obj_prec()
660{
661  return obj_prec();
662}
663void writeAmpl(ampl_info * info)
664{
665  char buf[1000];
666  typedef struct { const char *msg; int code; int wantObj; } Sol_info;
667  static Sol_info solinfo[] = {
668    { "optimal solution",                       000, 1 },
669    { "infeasible",                             200, 1 },
670    { "unbounded",                              300, 0 },
671    { "iteration limit etc",                    400, 1 },
672    { "solution limit",                         401, 1 },
673    { "ran out of space",                       500, 0 },
674    { "status unknown",                         501, 1 },
675    { "bug!",                                   502, 0 },
676    { "best MIP solution so far restored",      101, 1 },
677    { "failed to restore best MIP solution",    503, 1 },
678    { "optimal (?) solution",                   100, 1 }
679  };
680  /* convert status - need info on map */
681  static int map[] = {0, 3, 4, 1};
682  sprintf(buf,"%s %s",Oinfo.bsname,info->buffer);
683  solve_result_num = solinfo[info->problemStatus].code;
684  if (info->columnStatus) {
685    stat_map(info->columnStatus, n_var, map, 4, "outgoing columnStatus");
686    stat_map(info->rowStatus, n_con, map, 4, "outgoing rowStatus");
687    suf_iput("sstatus", ASL_Sufkind_var, info->columnStatus);
688    suf_iput("sstatus", ASL_Sufkind_con, info->rowStatus);
689  }
690  write_sol(buf,info->primalSolution,info->dualSolution,&Oinfo);
691}
692/* Read a problem from AMPL nl file
693 */
694CoinModel::CoinModel( int nonLinear, const char * fileName,const void * info)
695 :  numberRows_(0),
696    maximumRows_(0),
697    numberColumns_(0),
698    maximumColumns_(0),
699    numberElements_(0),
700    maximumElements_(0),
701    numberQuadraticElements_(0),
702    maximumQuadraticElements_(0),
703    optimizationDirection_(1.0),
704    objectiveOffset_(0.0),
705    rowLower_(NULL),
706    rowUpper_(NULL),
707    rowType_(NULL),
708    objective_(NULL),
709    columnLower_(NULL),
710    columnUpper_(NULL),
711    integerType_(NULL),
712    columnType_(NULL),
713    start_(NULL),
714    elements_(NULL),
715    quadraticElements_(NULL),
716    sortIndices_(NULL),
717    sortElements_(NULL),
718    sortSize_(0),
719    sizeAssociated_(0),
720    associated_(NULL),
721    numberSOS_(0),
722    startSOS_(NULL),
723    memberSOS_(NULL),
724    typeSOS_(NULL),
725    prioritySOS_(NULL),
726    referenceSOS_(NULL),
727    priority_(NULL),
728    cut_(NULL),
729    logLevel_(0),
730    type_(-1),
731    links_(0)
732{
733  problemName_ = strdup("");
734  int status=0;
735  if (!strcmp(fileName,"-")||!strcmp(fileName,"stdin")) {
736    // stdin
737  } else {
738    std::string name=fileName;
739    bool readable = fileCoinReadable(name);
740    if (!readable) {
741      std::cerr<<"Unable to open file "
742               <<fileName<<std::endl;
743      status = -1;
744    }
745  }
746  if (!status) {
747    gdb(nonLinear,fileName,info);
748  }
749}
750#if 0
751static real
752qterm(ASL *asl, fint *colq, fint *rowq, real *delsq)
753{
754  double t, t1, *x, *x0, *xe;
755  fint *rq0, *rqe;
756 
757  t = 0.;
758  x0 = x = X0;
759  xe = x + n_var;
760  rq0 = rowq;
761  while(x < xe) {
762    t1 = *x++;
763    rqe = rq0 + *++colq;
764    while(rowq < rqe)
765      t += t1*x0[*rowq++]**delsq++;
766  }
767  return 0.5 * t;
768}
769#endif
770// stolen from IPopt with changes
771typedef struct {
772  double obj_sign_;
773  ASL_pfgh * asl_;
774  double * non_const_x_;
775  int * column_; // for jacobian
776  int * rowStart_;
777  double * gradient_;
778  double * constraintValues_;
779  int nz_h_full_; // number of nonzeros in hessian
780  int nerror_;
781  bool objval_called_with_current_x_;
782  bool conval_called_with_current_x_;
783  bool jacval_called_with_current_x_;
784} CbcAmplInfo;
785
786void
787CoinModel::gdb( int nonLinear, const char * fileName,const void * info)
788{
789  const ampl_info * amplInfo = (const ampl_info *) info;
790  ograd *og=NULL;
791  int i;
792  SufDesc *csd=NULL;
793  SufDesc *rsd=NULL;
794  /*bool *basis, *lower;*/
795  /*double *LU, *c, lb, objadj, *rshift, *shift, t, ub, *x, *x0, *x1;*/
796  //char tempBuffer[20];
797  double * objective=NULL;
798  double * columnLower=NULL;
799  double * columnUpper=NULL;
800  double * rowLower=NULL;
801  double * rowUpper=NULL;
802  int * columnStatus=NULL;
803  int * rowStatus=NULL;
804  int numberRows=-1;
805  int numberColumns=-1;
806  int numberElements=-1;
807  int numberBinary=-1;
808  int numberIntegers=-1;
809  int numberAllNonLinearBoth=0;
810  int numberIntegerNonLinearBoth=0;
811  int numberAllNonLinearConstraints=0;
812  int numberIntegerNonLinearConstraints=0;
813  int numberAllNonLinearObjective=0;
814  int numberIntegerNonLinearObjective=0;
815  double * primalSolution=NULL;
816  double direction=1.0;
817  char * stub = strdup(fileName);
818  CoinPackedMatrix matrixByRow;
819  fint ** colqp=NULL;
820  int *z = NULL;
821  if (nonLinear==0) {
822    // linear
823    asl = ASL_alloc(ASL_read_f);
824    nl = jac0dim(stub, 0);
825    free(stub);
826    suf_declare(suftab, sizeof(suftab)/sizeof(SufDecl));
827   
828    /* set A_vals to get the constraints column-wise (malloc so can be freed) */
829    A_vals = (double *) malloc(nzc*sizeof(double));
830    if (!A_vals) {
831      printf("no memory\n");
832      return ;
833    }
834    /* say we want primal solution */
835    want_xpi0=1;
836    /* for basis info */
837    columnStatus = (int *) malloc(n_var*sizeof(int));
838    rowStatus = (int *) malloc(n_con*sizeof(int));
839    csd = suf_iput("sstatus", ASL_Sufkind_var, columnStatus);
840    rsd = suf_iput("sstatus", ASL_Sufkind_con, rowStatus);
841    /* read linear model*/
842    f_read(nl,0);
843    // see if any sos
844    if (true) {
845      char *sostype;
846      int nsosnz, *sosbeg, *sosind, * sospri;
847      double *sosref;
848      int nsos;
849      int i = ASL_suf_sos_explict_free;
850      int copri[2], **p_sospri;
851      copri[0] = 0;
852      copri[1] = 0;
853      p_sospri = &sospri;
854      nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
855                     &sosbeg, &sosind, &sosref);
856      if (nsos) {
857        abort();
858#if 0
859        info->numberSos=nsos;
860        info->sosType = (char *) malloc(nsos);
861        info->sosPriority = (int *) malloc(nsos*sizeof(int));
862        info->sosStart = (int *) malloc((nsos+1)*sizeof(int));
863        info->sosIndices = (int *) malloc(nsosnz*sizeof(int));
864        info->sosReference = (double *) malloc(nsosnz*sizeof(double));
865        sos_kludge(nsos, sosbeg, sosref,sosind);
866        for (int i=0;i<nsos;i++) {
867          int ichar = sostype[i];
868          assert (ichar=='1'||ichar=='2');
869          info->sosType[i]=ichar-'0';
870        }       
871        memcpy(info->sosPriority,sospri,nsos*sizeof(int));
872        memcpy(info->sosStart,sosbeg,(nsos+1)*sizeof(int));
873        memcpy(info->sosIndices,sosind,nsosnz*sizeof(int));
874        memcpy(info->sosReference,sosref,nsosnz*sizeof(double));
875#endif
876      }
877    }
878   
879    /*sos_finish(&specialOrderedInfo, 0, &j, 0, 0, 0, 0, 0);*/
880    //Oinfo.uinfo = tempBuffer;
881    //if (getopts(argv, &Oinfo))
882    //return 1;
883    /* objective*/
884    objective = (double *) malloc(n_var*sizeof(double));
885    for (i=0;i<n_var;i++)
886      objective[i]=0.0;
887    if (n_obj) {
888      for (og = Ograd[0];og;og = og->next)
889        objective[og->varno] = og->coef;
890    }
891    if (objtype[0])
892      direction=-1.0;
893    else
894      direction=1.0;
895    objectiveOffset_=objconst(0);
896    /* Column bounds*/
897    columnLower = (double *) malloc(n_var*sizeof(double));
898    columnUpper = (double *) malloc(n_var*sizeof(double));
899#define COIN_DBL_MAX DBL_MAX
900    for (i=0;i<n_var;i++) {
901      columnLower[i]=LUv[2*i];
902      if (columnLower[i]<= negInfinity)
903        columnLower[i]=-COIN_DBL_MAX;
904      columnUpper[i]=LUv[2*i+1];
905      if (columnUpper[i]>= Infinity)
906        columnUpper[i]=COIN_DBL_MAX;
907    }
908    /* Row bounds*/
909    rowLower = (double *) malloc(n_con*sizeof(double));
910    rowUpper = (double *) malloc(n_con*sizeof(double));
911    for (i=0;i<n_con;i++) {
912      rowLower[i]=LUrhs[2*i];
913      if (rowLower[i]<= negInfinity)
914        rowLower[i]=-COIN_DBL_MAX;
915      rowUpper[i]=LUrhs[2*i+1];
916      if (rowUpper[i]>= Infinity)
917        rowUpper[i]=COIN_DBL_MAX;
918    }
919    numberRows=n_con;
920    numberColumns=n_var;
921    numberElements=nzc;
922    numberBinary=nbv;
923    numberIntegers=niv;
924    /* put in primalSolution if exists */
925    if (X0) {
926      primalSolution=(double *) malloc(n_var*sizeof(double));
927      memcpy( primalSolution,X0,n_var*sizeof(double));
928    }
929    //double * dualSolution=NULL;
930    if (niv+nbv>0)
931      mip_stuff(); // get any extra info
932    if ((!(niv+nbv)&&(csd->kind & ASL_Sufkind_input))
933        ||(rsd->kind & ASL_Sufkind_input)) {
934      /* convert status - need info on map */
935      static int map[] = {1, 3, 1, 1, 2, 1, 1};
936      stat_map(columnStatus, n_var, map, 6, "incoming columnStatus");
937      stat_map(rowStatus, n_con, map, 6, "incoming rowStatus");
938    } else {
939      /* all slack basis */
940      // leave status for output */
941#if 0
942      free(rowStatus);
943      rowStatus=NULL;
944      free(columnStatus);
945      columnStatus=NULL;
946#endif
947    }
948    CoinPackedMatrix columnCopy(true,numberRows,numberColumns,numberElements,
949                                A_vals,A_rownos,A_colstarts,NULL);
950    matrixByRow.reverseOrderedCopyOf(columnCopy);
951  } else if (nonLinear==1) {
952    // quadratic
953    asl = ASL_alloc(ASL_read_fg);
954    nl = jac0dim(stub, (long) strlen(stub));
955    free(stub);
956    suf_declare(suftab, sizeof(suftab)/sizeof(SufDecl));
957    /* read  model*/
958    X0 = (double*) malloc(n_var*sizeof(double));
959    CoinZeroN(X0,n_var);
960    qp_read(nl,0);
961    assert (n_obj==1);
962    int nz = 1 + n_con;
963    colqp = (fint**) malloc(nz*(2*sizeof(int*)
964                                      + sizeof(double*)));
965    fint ** rowqp = colqp + nz;
966    double ** delsqp = (double **)(rowqp + nz);
967    z = (int*) malloc(nz*sizeof(int));
968    for (i=0;i<=n_con;i++) {
969      z[i] = nqpcheck(-i, rowqp+i, colqp+i, delsqp+i);
970    }
971    qp_opify();
972    /* objective*/
973    objective = (double *) malloc(n_var*sizeof(double));
974    for (i=0;i<n_var;i++)
975      objective[i]=0.0;
976    if (n_obj) {
977      for (og = Ograd[0];og;og = og->next)
978        objective[og->varno] = og->coef;
979    }
980    if (objtype[0])
981      direction=-1.0;
982    else
983      direction=1.0;
984    objectiveOffset_=objconst(0);
985    /* Column bounds*/
986    columnLower = (double *) malloc(n_var*sizeof(double));
987    columnUpper = (double *) malloc(n_var*sizeof(double));
988    for (i=0;i<n_var;i++) {
989      columnLower[i]=LUv[2*i];
990      if (columnLower[i]<= negInfinity)
991        columnLower[i]=-COIN_DBL_MAX;
992      columnUpper[i]=LUv[2*i+1];
993      if (columnUpper[i]>= Infinity)
994        columnUpper[i]=COIN_DBL_MAX;
995    }
996    // Build by row from scratch
997    //matrixByRow.reserve(n_var,nzc,true);
998    // say row orderded
999    matrixByRow.transpose();
1000    /* Row bounds*/
1001    rowLower = (double *) malloc(n_con*sizeof(double));
1002    rowUpper = (double *) malloc(n_con*sizeof(double));
1003    int * column = new int [n_var];
1004    double * element = new double [n_var];
1005    for (i=0;i<n_con;i++) {
1006      rowLower[i]=LUrhs[2*i];
1007      if (rowLower[i]<= negInfinity)
1008        rowLower[i]=-COIN_DBL_MAX;
1009      rowUpper[i]=LUrhs[2*i+1];
1010      if (rowUpper[i]>= Infinity)
1011        rowUpper[i]=COIN_DBL_MAX;
1012      int k=0;
1013      for(cgrad * cg = Cgrad[i]; cg; cg = cg->next) {
1014        column[k]=cg->varno;
1015        element[k++]=cg->coef;
1016      }
1017      matrixByRow.appendRow(k,column,element);
1018    }
1019    delete [] column;
1020    delete [] element;
1021    numberRows=n_con;
1022    numberColumns=n_var;
1023    numberElements=nzc;
1024    numberBinary=nbv;
1025    numberIntegers=niv;
1026    numberAllNonLinearBoth=nlvb;
1027    numberIntegerNonLinearBoth=nlvbi;
1028    numberAllNonLinearConstraints=nlvc;
1029    numberIntegerNonLinearConstraints=nlvci;
1030    numberAllNonLinearObjective=nlvo;
1031    numberIntegerNonLinearObjective=nlvoi;
1032    /* say we want primal solution */
1033    want_xpi0=1;
1034    //double * dualSolution=NULL;
1035    // save asl
1036    // Fix memory leak one day
1037    CbcAmplInfo * info = new CbcAmplInfo;
1038    //amplGamsData_ = info;
1039    info->asl_=NULL; // as wrong form asl;
1040    info->nz_h_full_=-1; // number of nonzeros in hessian
1041    info->objval_called_with_current_x_=false;
1042    info->nerror_=0;
1043    info->obj_sign_=direction;
1044    info->conval_called_with_current_x_=false;
1045    info->non_const_x_=NULL;
1046    info->jacval_called_with_current_x_=false;
1047    info->rowStart_ = NULL;
1048    info->column_ = NULL;
1049    info->gradient_ = NULL;
1050    info->constraintValues_ = NULL;
1051  } else if (nonLinear==2) {
1052    // General nonlinear!
1053    //ASL_pfgh* asl = (ASL_pfgh*)ASL_alloc(ASL_read_pfgh);
1054    asl = ASL_alloc(ASL_read_pfgh);
1055    nl = jac0dim(stub, (long) strlen(stub));
1056    free(stub);
1057    suf_declare(suftab, sizeof(suftab)/sizeof(SufDecl));
1058    /* read  model*/
1059    X0 = (double*) malloc(n_var*sizeof(double));
1060    CoinZeroN(X0,n_var);
1061    // code stolen from Ipopt
1062    int retcode = pfgh_read(nl, ASL_return_read_err | ASL_findgroups);
1063
1064    switch (retcode) {
1065      case ASL_readerr_none : {}
1066      break;
1067      case ASL_readerr_nofile : {
1068        printf( "Cannot open .nl file\n");
1069        exit(-1);
1070      }
1071      break;
1072      case ASL_readerr_nonlin : {
1073        assert(false); // this better not be an error!
1074        printf( "model involves nonlinearities (ed0read)\n");
1075        exit(-1);
1076      }
1077      break;
1078      case  ASL_readerr_argerr : {
1079        printf( "user-defined function with bad args\n");
1080        exit(-1);
1081      }
1082      break;
1083      case ASL_readerr_unavail : {
1084        printf( "user-defined function not available\n");
1085        exit(-1);
1086      }
1087      break;
1088      case ASL_readerr_corrupt : {
1089        printf( "corrupt .nl file\n");
1090        exit(-1);
1091      }
1092      break;
1093      case ASL_readerr_bug : {
1094        printf( "bug in .nl reader\n");
1095        exit(-1);
1096      }
1097      break;
1098      case ASL_readerr_CLP : {
1099        printf( "ASL error message: \"solver cannot handle CLP extensions\"\n");
1100        exit(-1);
1101      }
1102      break;
1103      default: {
1104        printf( "Unknown error in stub file read. retcode = %d\n", retcode);
1105        exit(-1);
1106      }
1107      break;
1108    }
1109
1110    // see "changes" in solvers directory of ampl code...
1111    hesset(1,0,1,0,nlc);
1112
1113    assert (n_obj==1);
1114    // find the nonzero structure for the hessian
1115    // parameters to sphsetup:
1116    int coeff_obj = 1; // coefficient of the objective fn ???
1117    int mult_supplied = 1; // multipliers will be supplied
1118    int uptri = 1; // only need the upper triangular part
1119    // save asl
1120    // Fix memory leak one day
1121    CbcAmplInfo * info = new CbcAmplInfo;
1122    moreInfo_ = (void *) info;
1123    //amplGamsData_ = info;
1124    info->asl_=(ASL_pfgh *) asl;
1125    // This is not easy to get from ampl so save
1126    info->nz_h_full_= sphsetup(-1, coeff_obj, mult_supplied, uptri);
1127    info->objval_called_with_current_x_=false;
1128    info->nerror_=0;
1129    info->obj_sign_=direction;
1130    info->conval_called_with_current_x_=false;
1131    info->non_const_x_=NULL;
1132    info->jacval_called_with_current_x_=false;
1133    // Look at nonlinear
1134    if (nzc) {
1135      n_conjac[1]=nlc; // just nonlinear
1136      int * rowStart = new int [nlc+1];
1137      info->rowStart_ = rowStart;
1138      // See how many
1139      int  current_nz = 0;
1140      for (int i=0; i<nlc; i++) {
1141        for (cgrad* cg=Cgrad[i]; cg; cg = cg->next) {
1142          current_nz++;
1143        }
1144      }
1145      // setup the structure
1146      int * column = new int [current_nz];
1147      info->column_ = column;
1148      current_nz = 0;
1149      rowStart[0]=0;
1150      for (int i=0; i<nlc; i++) {
1151        for (cgrad* cg=Cgrad[i]; cg; cg = cg->next) {
1152          cg->goff=current_nz;
1153          //iRow[cg->goff] = i ;
1154          //jCol[cg->goff] = cg->varno + 1;
1155          column[cg->goff] = cg->varno ;
1156          current_nz++;
1157        }
1158        rowStart[i+1]=current_nz;
1159      }
1160      info->gradient_ = new double [nzc];
1161      info->constraintValues_ = new double [nlc];
1162    }
1163    /* objective*/
1164    objective = (double *) malloc(n_var*sizeof(double));
1165    for (i=0;i<n_var;i++)
1166      objective[i]=0.0;
1167    if (n_obj) {
1168      for (og = Ograd[0];og;og = og->next)
1169        objective[og->varno] = og->coef;
1170    }
1171    if (objtype[0])
1172      direction=-1.0;
1173    else
1174      direction=1.0;
1175    objectiveOffset_=objconst(0);
1176    /* Column bounds*/
1177    columnLower = (double *) malloc(n_var*sizeof(double));
1178    columnUpper = (double *) malloc(n_var*sizeof(double));
1179    for (i=0;i<n_var;i++) {
1180      columnLower[i]=LUv[2*i];
1181      if (columnLower[i]<= negInfinity)
1182        columnLower[i]=-COIN_DBL_MAX;
1183      columnUpper[i]=LUv[2*i+1];
1184      if (columnUpper[i]>= Infinity)
1185        columnUpper[i]=COIN_DBL_MAX;
1186    }
1187    // Build by row from scratch
1188    //matrixByRow.reserve(n_var,nzc,true);
1189    // say row orderded
1190    matrixByRow.transpose();
1191    /* Row bounds*/
1192    rowLower = (double *) malloc(n_con*sizeof(double));
1193    rowUpper = (double *) malloc(n_con*sizeof(double));
1194    int * column = new int [n_var];
1195    double * element = new double [n_var];
1196    for (i=0;i<n_con;i++) {
1197      rowLower[i]=LUrhs[2*i];
1198      if (rowLower[i]<= negInfinity)
1199        rowLower[i]=-COIN_DBL_MAX;
1200      rowUpper[i]=LUrhs[2*i+1];
1201      if (rowUpper[i]>= Infinity)
1202        rowUpper[i]=COIN_DBL_MAX;
1203      int k=0;
1204      for(cgrad * cg = Cgrad[i]; cg; cg = cg->next) {
1205        column[k]=cg->varno;
1206        double value = cg->coef;
1207        if (!value)
1208          value = -1.2345e-29;
1209        element[k++]=value;
1210      }
1211      matrixByRow.appendRow(k,column,element);
1212    }
1213    delete [] column;
1214    delete [] element;
1215    numberRows=n_con;
1216    numberColumns=n_var;
1217    numberElements=nzc;
1218    numberBinary=nbv;
1219    numberIntegers=niv;
1220    numberAllNonLinearBoth=nlvb;
1221    numberIntegerNonLinearBoth=nlvbi;
1222    numberAllNonLinearConstraints=nlvc;
1223    numberIntegerNonLinearConstraints=nlvci;
1224    numberAllNonLinearObjective=nlvo;
1225    numberIntegerNonLinearObjective=nlvoi;
1226    /* say we want primal solution */
1227    want_xpi0=1;
1228    //double * dualSolution=NULL;
1229  } else {
1230    abort();
1231  }
1232  // set problem name
1233  free(problemName_);
1234  problemName_=strdup("???");
1235
1236  // Build by row from scratch
1237  const double * element = matrixByRow.getElements();
1238  const int * column = matrixByRow.getIndices();
1239  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
1240  const int * rowLength = matrixByRow.getVectorLengths();
1241  for (i=0;i<numberRows;i++) {
1242    addRow(rowLength[i],column+rowStart[i],
1243           element+rowStart[i],rowLower[i],rowUpper[i]);
1244  }
1245  // Now do column part
1246  for (i=0;i<numberColumns;i++) {
1247    setColumnBounds(i,columnLower[i],columnUpper[i]);
1248    setColumnObjective(i,objective[i]);
1249  }
1250  for ( i=numberColumns-numberBinary-numberIntegers;
1251        i<numberColumns;i++) {
1252    setColumnIsInteger(i,true);
1253  }
1254  // and non linear
1255  for (i=numberAllNonLinearBoth-numberIntegerNonLinearBoth;
1256       i<numberAllNonLinearBoth;i++) {
1257    setColumnIsInteger(i,true);
1258  }
1259  for (i=numberAllNonLinearConstraints-numberIntegerNonLinearConstraints;
1260       i<numberAllNonLinearConstraints;i++) {
1261    setColumnIsInteger(i,true);
1262  }
1263  for (i=numberAllNonLinearObjective-numberIntegerNonLinearObjective;
1264       i<numberAllNonLinearObjective;i++) {
1265    setColumnIsInteger(i,true);
1266  }
1267  free(columnLower);
1268  free(columnUpper);
1269  free(rowLower);
1270  free(rowUpper);
1271  free(objective);
1272  // do names
1273  int iRow;
1274  for (iRow=0;iRow<numberRows_;iRow++) {
1275    char name[9];
1276    sprintf(name,"r%7.7d",iRow);
1277    setRowName(iRow,name);
1278  }
1279  int iColumn;
1280  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1281    char name[9];
1282    sprintf(name,"c%7.7d",iColumn);
1283    setColumnName(iColumn,name);
1284  }
1285  if (colqp) {
1286    // add in quadratic
1287    int nz = 1 + n_con;
1288    int nOdd=0;
1289    fint ** rowqp = colqp + nz;
1290    double ** delsqp = (double **)(rowqp + nz);
1291    for (i=0;i<=n_con;i++) {
1292      int nels = z[i];
1293      if (nels) {
1294        double * element = delsqp[i];
1295        int * start = (int *) colqp[i];
1296        int * row = (int *) rowqp[i];
1297        if (!element) {
1298          // odd row - probably not quadratic
1299          nOdd++;
1300          continue;
1301        }
1302#if 0
1303        printf("%d quadratic els\n",nels);
1304        for (int j=0;j<n_var;j++) {
1305          for (int k=start[j];k<start[j+1];k++)
1306            printf("%d %d %g\n",j,row[k],element[k]);
1307        }
1308#endif
1309        if (i) {
1310          int iRow = i-1;
1311          for (int j=0;j<n_var;j++) {
1312            for (int k=start[j];k<start[j+1];k++) {
1313              int kColumn = row[k];
1314              double value = element[k];
1315              // ampl gives twice with assumed 0.5
1316              if (kColumn<j)
1317                continue;
1318              else if (kColumn==j)
1319                value *= 0.5;
1320              const char * expr = getElementAsString(iRow,j);
1321              double constant=0.0;
1322              bool linear;
1323              if (expr&&strcmp(expr,"Numeric")) {
1324                linear=false;
1325              } else {
1326                constant = getElement(iRow,j);
1327                linear=true;
1328              }
1329              char temp[1000];
1330              char temp2[30];
1331              if (value==1.0) 
1332                sprintf(temp2,"c%7.7d",kColumn);
1333              else
1334                sprintf(temp2,"%g*c%7.7d",value,kColumn);
1335              if (linear) {
1336                if (!constant)
1337                  strcpy(temp,temp2);
1338                else if (value>0.0) 
1339                  sprintf(temp,"%g+%s",constant,temp2);
1340                else
1341                  sprintf(temp,"%g%s",constant,temp2);
1342              } else {
1343                if (value>0.0) 
1344                  sprintf(temp,"%s+%s",expr,temp2);
1345                else
1346                  sprintf(temp,"%s%s",expr,temp2);
1347              }
1348              assert (strlen(temp)<1000);
1349              setElement(iRow,j,temp);
1350              if (amplInfo->logLevel>1)
1351                printf("el for row %d column c%7.7d is %s\n",iRow,j,temp);
1352            }
1353          }
1354        } else {
1355          // objective
1356          for (int j=0;j<n_var;j++) {
1357            for (int k=start[j];k<start[j+1];k++) {
1358              int kColumn = row[k];
1359              double value = element[k];
1360              // ampl gives twice with assumed 0.5
1361              if (kColumn<j)
1362                continue;
1363              else if (kColumn==j)
1364                value *= 0.5;
1365              const char * expr = getColumnObjectiveAsString(j);
1366              double constant=0.0;
1367              bool linear;
1368              if (expr&&strcmp(expr,"Numeric")) {
1369                linear=false;
1370              } else {
1371                constant = getColumnObjective(j);
1372                linear=true;
1373              }
1374              char temp[1000];
1375              char temp2[30];
1376              if (value==1.0) 
1377                sprintf(temp2,"c%7.7d",kColumn);
1378              else
1379                sprintf(temp2,"%g*c%7.7d",value,kColumn);
1380              if (linear) {
1381                if (!constant)
1382                  strcpy(temp,temp2);
1383                else if (value>0.0) 
1384                  sprintf(temp,"%g+%s",constant,temp2);
1385                else
1386                  sprintf(temp,"%g%s",constant,temp2);
1387              } else {
1388                if (value>0.0) 
1389                  sprintf(temp,"%s+%s",expr,temp2);
1390                else
1391                  sprintf(temp,"%s%s",expr,temp2);
1392              }
1393              assert (strlen(temp)<1000);
1394              setObjective(j,temp);
1395              if (amplInfo->logLevel>1)
1396                printf("el for objective column c%7.7d is %s\n",j,temp);
1397            }
1398          }
1399        }
1400      }
1401    }
1402    if (nOdd) {
1403      printf("%d non-linear constraints could not be converted to quadratic\n",nOdd);
1404      exit(77);
1405    }
1406  }
1407  free(colqp);
1408  free(z);
1409  // see if any sos
1410  {
1411    char *sostype;
1412    int nsosnz, *sosbeg, *sosind, * sospri;
1413    double *sosref;
1414    int nsos;
1415    int i = ASL_suf_sos_explict_free;
1416    int copri[2], **p_sospri;
1417    copri[0] = 0;
1418    copri[1] = 0;
1419    p_sospri = &sospri;
1420    nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
1421                   &sosbeg, &sosind, &sosref);
1422    if (nsos) {
1423      numberSOS_=nsos;
1424      typeSOS_ = new int [numberSOS_];
1425      prioritySOS_ = new int [numberSOS_];
1426      startSOS_ = new int [numberSOS_+1];
1427      memberSOS_ = new int[nsosnz];
1428      referenceSOS_ = new double [nsosnz];
1429      sos_kludge(nsos, sosbeg, sosref,sosind);
1430      for (int i=0;i<nsos;i++) {
1431        int ichar = sostype[i];
1432        assert (ichar=='1'||ichar=='2');
1433        typeSOS_[i]=ichar-'0';
1434      } 
1435      memcpy(prioritySOS_,sospri,nsos*sizeof(int));
1436      memcpy(startSOS_,sosbeg,(nsos+1)*sizeof(int));
1437      memcpy(memberSOS_,sosind,nsosnz*sizeof(int));
1438      memcpy(referenceSOS_,sosref,nsosnz*sizeof(double));
1439    }
1440  }
1441}
1442#else
1443#include "Cbc_ampl.h"
1444int
1445readAmpl(ampl_info * info, int argc, char **argv, void ** coinModel)
1446{
1447  return 0;
1448}
1449void freeArrays1(ampl_info * info)
1450{
1451}
1452void freeArrays2(ampl_info * info)
1453{
1454}
1455void freeArgs(ampl_info * info)
1456{
1457}
1458int ampl_obj_prec()
1459{
1460  return 0;
1461}
1462void writeAmpl(ampl_info * info)
1463{
1464}
1465#endif
Note: See TracBrowser for help on using the repository browser.