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

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

for owning solution in OsiBranchingObject?

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 30.5 KB
Line 
1/****************************************************************
2Copyright (C) 1997-2000 Lucent Technologies
3Modifications for Coin -  Copyright (C) 2006, International Business Machines Corporation and others.
4All Rights Reserved
5
6Permission to use, copy, modify, and distribute this software and
7its documentation for any purpose and without fee is hereby
8granted, provided that the above copyright notice appear in all
9copies and that both that the copyright notice and this
10permission notice and warranty disclaimer appear in supporting
11documentation, and that the name of Lucent or any of its entities
12not be used in advertising or publicity pertaining to
13distribution of the software without specific, written prior
14permission.
15
16LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
17INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
18IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
19SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
20WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
21IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
22ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
23THIS SOFTWARE.
24****************************************************************/
25
26#include "CbcConfig.h"
27#ifdef HAVE_UNISTD_H
28# include "unistd.h"
29#endif
30#include "CoinUtilsConfig.h"
31#include "CoinHelperFunctions.hpp"
32#include "CoinModel.hpp"
33#include "CoinSort.hpp"
34#include "CoinPackedMatrix.hpp"
35#include "CoinMpsIO.hpp"
36#include "CoinFloatEqual.hpp"
37
38extern "C" {
39# include "getstub.h"
40}
41
42#include "Cbc_ampl.h"
43#include <string>
44#include <cassert>
45/* so decodePhrase and clpCheck can access */
46static ampl_info * saveInfo=NULL;
47// Set to 1 if algorithm found
48static char algFound[20]="";
49static char*
50checkPhrase(Option_Info *oi, keyword *kw, char *v)
51{
52  if (strlen(v))
53    printf("string %s\n",v);
54  // Say algorithm found
55  strcpy(algFound,kw->desc);;
56  return v;
57}
58static char*
59checkPhrase2(Option_Info *oi, keyword *kw, char *v)
60{
61  if (strlen(v))
62    printf("string %s\n",v);
63  // put out keyword
64  saveInfo->arguments=(char **) realloc(saveInfo->arguments,(saveInfo->numberArguments+1)*sizeof(char *));
65  saveInfo->arguments[saveInfo->numberArguments++]=strdup(kw->desc);
66  return v;
67}
68static fint
69decodePhrase(char * phrase,ftnlen length)
70{
71  char * blank = strchr(phrase,' ');
72  if (blank) {
73    /* split arguments */
74    *blank='\0';
75    saveInfo->arguments=(char **) realloc(saveInfo->arguments,(saveInfo->numberArguments+2)*sizeof(char *));
76    saveInfo->arguments[saveInfo->numberArguments++]=strdup(phrase);
77    *blank=' ';
78    phrase=blank+1; /* move on */
79    if (strlen(phrase))
80      saveInfo->arguments[saveInfo->numberArguments++]=strdup(phrase);
81  } else if (strlen(phrase)) {
82    saveInfo->arguments=(char **) realloc(saveInfo->arguments,(saveInfo->numberArguments+1)*sizeof(char *));
83    saveInfo->arguments[saveInfo->numberArguments++]=strdup(phrase);
84  }
85  return 0;
86}
87static void
88sos_kludge(int nsos, int *sosbeg, double *sosref,int * sosind)
89{
90  // Adjust sosref if necessary to make monotonic increasing
91  int i, j, k;
92  // first sort
93  for (i=0;i<nsos;i++) {
94    k = sosbeg[i];
95    int end=sosbeg[i+1];
96    CoinSort_2(sosref+k,sosref+end,sosind+k);
97  }
98  double t, t1;
99  for(i = j = 0; i++ < nsos; ) {
100    k = sosbeg[i];
101    t = sosref[j];
102    while(++j < k) {
103      t1 = sosref[j];
104      t += 1e-10;
105      if (t1 <= t)
106        sosref[j] = t1 = t + 1e-10;
107      t = t1;
108    }
109  }
110}
111static char xxxxxx[20];
112#define VP (char*)
113 static keyword keywds[] = { /* must be sorted */
114        { "barrier",    checkPhrase,            (char *) xxxxxx ,"-barrier" },
115        { "dual",       checkPhrase,            (char *) xxxxxx , "-dualsimplex"},
116        { "help",       checkPhrase2,           (char *) xxxxxx , "-?"},
117        { "initial",    checkPhrase,            (char *) xxxxxx , "-initialsolve"},
118        { "max",        checkPhrase2,           (char *) xxxxxx , "-maximize"},
119        { "maximize",   checkPhrase2,           (char *) xxxxxx , "-maximize"},
120        { "primal",     checkPhrase,            (char *) xxxxxx , "-primalsimplex"},
121        { "quit",       checkPhrase2,           (char *) xxxxxx , "-quit"},
122        { "wantsol",    WS_val,         NULL, "write .sol file (without -AMPL)" }
123        };
124static Option_Info Oinfo = {"cbc", "Cbc 1.01", "cbc_options", keywds, nkeywds, 0, "",
125                                0,decodePhrase,0,0,0, 20060130 };
126// strdup used to avoid g++ compiler warning
127 static SufDecl
128suftab[] = {
129#if 0
130        { "current", 0, ASL_Sufkind_con | ASL_Sufkind_outonly },
131        { "current", 0, ASL_Sufkind_var | ASL_Sufkind_outonly },
132        { "direction", 0, ASL_Sufkind_var },
133        { "down", 0, ASL_Sufkind_con | ASL_Sufkind_outonly },
134        { "down", 0, ASL_Sufkind_var | ASL_Sufkind_outonly },
135        { "priority", 0, ASL_Sufkind_var },
136#endif
137        { "direction", 0, ASL_Sufkind_var },
138        { "downPseudocost", 0, ASL_Sufkind_var | ASL_Sufkind_real },
139        { "priority", 0, ASL_Sufkind_var },
140        { "ref", 0, ASL_Sufkind_var | ASL_Sufkind_real },
141        { "sos", 0, ASL_Sufkind_var },
142        { "sos", 0, ASL_Sufkind_con },
143        { "sosno", 0, ASL_Sufkind_var | ASL_Sufkind_real },
144        { "sosref", 0, ASL_Sufkind_var | ASL_Sufkind_real },
145        { strdup("sstatus"), 0, ASL_Sufkind_var, 0 },
146        { strdup("sstatus"), 0, ASL_Sufkind_con, 0 },
147        { "upPseudocost", 0, ASL_Sufkind_var | ASL_Sufkind_real }
148#if 0
149        { "unbdd", 0, ASL_Sufkind_var | ASL_Sufkind_outonly},
150        { "up", 0, ASL_Sufkind_con | ASL_Sufkind_outonly },
151        { "up", 0, ASL_Sufkind_var | ASL_Sufkind_outonly }
152#endif
153        };
154#include "float.h"
155#include "limits.h"
156static ASL *asl=NULL;
157static FILE *nl=NULL;
158
159static void
160mip_stuff(void)
161{
162  int i;
163  double *pseudoUp, *pseudoDown;
164  int *priority, *direction;
165  SufDesc *dpup, *dpdown, *dpri, *ddir;
166 
167  ddir = suf_get("direction", ASL_Sufkind_var);
168  direction = ddir->u.i;
169  dpri = suf_get("priority", ASL_Sufkind_var);
170  priority = dpri->u.i;
171  dpdown = suf_get("downPseudocost", ASL_Sufkind_var);
172  pseudoDown = dpdown->u.r;
173  dpup = suf_get("upPseudocost", ASL_Sufkind_var);
174  pseudoUp = dpup->u.r;
175  assert(saveInfo);
176  int numberColumns = saveInfo->numberColumns;
177  if (direction) {
178    int baddir=0;
179    saveInfo->branchDirection = (int *) malloc(numberColumns*sizeof(int));
180    for (i=0;i<numberColumns;i++) {
181      int value = direction[i];
182      if (value<-1||value>1) {
183        baddir++;
184        value=0;
185      }
186      saveInfo->branchDirection[i]=value;
187    }
188    if (baddir)
189      fprintf(Stderr,
190              "Treating %d .direction values outside [-1, 1] as 0.\n",
191              baddir);
192  }
193  if (priority) {
194    int badpri=0;
195    saveInfo->priorities= (int *) malloc(numberColumns*sizeof(int));
196    for (i=0;i<numberColumns;i++) {
197      int value = priority[i];
198      if (value<0) {
199        badpri++;
200        value=0;
201      }
202      saveInfo->priorities[i]=value;
203    }
204    if (badpri)
205      fprintf(Stderr,
206              "Treating %d negative .priority values as 0\n",
207              badpri);
208  }
209  if (pseudoDown||pseudoUp) {
210    int badpseudo=0;
211    if (!pseudoDown||!pseudoUp)
212      fprintf(Stderr,
213              "Only one set of pseudocosts - assumed same\n");
214    saveInfo->pseudoDown= (double *) malloc(numberColumns*sizeof(double));
215    saveInfo->pseudoUp = (double *) malloc(numberColumns*sizeof(double));
216    for (i=0;i<numberColumns;i++) {
217      double valueD=0.0, valueU=0.0;
218      if (pseudoDown) {
219        valueD = pseudoDown[i];
220        if (valueD<0) {
221          badpseudo++;
222          valueD=0.0;
223        }
224      }
225      if (pseudoUp) {
226        valueU = pseudoUp[i];
227        if (valueU<0) {
228          badpseudo++;
229          valueU=0.0;
230        }
231      }
232      if (!valueD)
233        valueD=valueU;
234      if (!valueU)
235        valueU=valueD;
236      saveInfo->pseudoDown[i]=valueD;
237      saveInfo->pseudoUp[i]=valueU;
238    }
239    if (badpseudo)
240      fprintf(Stderr,
241              "Treating %d negative pseudoCosts as 0.0\n",badpseudo);
242  }
243}
244static void
245stat_map(int *stat, int n, int *map, int mx, const char *what)
246{
247  int bad, i, i1=0, j, j1=0;
248  static char badfmt[] = "Coin driver: %s[%d] = %d\n";
249 
250  for(i = bad = 0; i < n; i++) {
251    if ((j = stat[i]) >= 0 && j <= mx)
252      stat[i] = map[j];
253    else {
254      stat[i] = 0;
255      i1 = i;
256      j1 = j;
257      if (!bad++)
258        fprintf(Stderr, badfmt, what, i, j);
259    }
260  }
261  if (bad > 1) {
262    if (bad == 2)
263      fprintf(Stderr, badfmt, what, i1, j1);
264    else
265      fprintf(Stderr,
266              "Coin driver: %d messages about bad %s values suppressed.\n",
267              bad-1, what);
268  }
269}
270int
271readAmpl(ampl_info * info, int argc, char **argv)
272{
273  char *stub;
274  ograd *og;
275  int i;
276  SufDesc *csd;
277  SufDesc *rsd;
278  /*bool *basis, *lower;*/
279  /*double *LU, *c, lb, objadj, *rshift, *shift, t, ub, *x, *x0, *x1;*/
280  char * environment = getenv("cbc_options");
281  char tempBuffer[20];
282  double * obj;
283  double * columnLower;
284  double * columnUpper;
285  double * rowLower;
286  double * rowUpper;
287  char ** saveArgv=argv;
288  int saveArgc = argc;
289  memset(info,0,sizeof(ampl_info));
290  /* save so can be accessed by decodePhrase */
291  saveInfo = info;
292  info->numberArguments=0;
293  info->arguments=(char **) malloc(2*sizeof(char *));
294  info->arguments[info->numberArguments++]=strdup("ampl");
295  info->arguments[info->numberArguments++]=strdup("cbc");
296  asl = ASL_alloc(ASL_read_f);
297  stub = getstub(&argv, &Oinfo);
298  if (!stub)
299    usage_ASL(&Oinfo, 1);
300  nl = jac0dim(stub, 0);
301  suf_declare(suftab, sizeof(suftab)/sizeof(SufDecl));
302 
303  /* set A_vals to get the constraints column-wise (malloc so can be freed) */
304  A_vals = (double *) malloc(nzc*sizeof(double));
305  if (!A_vals) {
306    printf("no memory\n");
307    return 1;
308  }
309  /* say we want primal solution */
310  want_xpi0=1;
311  /* for basis info */
312  info->columnStatus = (int *) malloc(n_var*sizeof(int));
313  info->rowStatus = (int *) malloc(n_con*sizeof(int));
314  csd = suf_iput("sstatus", ASL_Sufkind_var, info->columnStatus);
315  rsd = suf_iput("sstatus", ASL_Sufkind_con, info->rowStatus);
316  /* read linear model*/
317  f_read(nl,0);
318  // see if any sos
319  if (true) {
320    char *sostype;
321    int nsosnz, *sosbeg, *sosind, * sospri;
322    double *sosref;
323    int nsos;
324    int i = ASL_suf_sos_explict_free;
325    int copri[2], **p_sospri;
326    copri[0] = 0;
327    copri[1] = 0;
328    p_sospri = &sospri;
329    nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
330                                &sosbeg, &sosind, &sosref);
331    if (nsos) {
332      info->numberSos=nsos;
333      info->sosType = (char *) malloc(nsos);
334      info->sosPriority = (int *) malloc(nsos*sizeof(int));
335      info->sosStart = (int *) malloc((nsos+1)*sizeof(int));
336      info->sosIndices = (int *) malloc(nsosnz*sizeof(int));
337      info->sosReference = (double *) malloc(nsosnz*sizeof(double));
338      sos_kludge(nsos, sosbeg, sosref,sosind);
339      for (int i=0;i<nsos;i++) {
340        int ichar = sostype[i];
341        assert (ichar=='1'||ichar=='2');
342        info->sosType[i]=ichar-'0';
343      } 
344      memcpy(info->sosPriority,sospri,nsos*sizeof(int));
345      memcpy(info->sosStart,sosbeg,(nsos+1)*sizeof(int));
346      memcpy(info->sosIndices,sosind,nsosnz*sizeof(int));
347      memcpy(info->sosReference,sosref,nsosnz*sizeof(double));
348    }
349  }
350
351  /*sos_finish(&specialOrderedInfo, 0, &j, 0, 0, 0, 0, 0);*/
352  Oinfo.uinfo = tempBuffer;
353  if (getopts(argv, &Oinfo))
354    return 1;
355  /* objective*/
356  obj = (double *) malloc(n_var*sizeof(double));
357  for (i=0;i<n_var;i++)
358    obj[i]=0.0;;
359  if (n_obj) {
360    for (og = Ograd[0];og;og = og->next)
361      obj[og->varno] = og->coef;
362  }
363  if (objtype[0])
364    info->direction=-1.0;
365  else
366    info->direction=1.0;
367  info->offset=objconst(0);
368  /* Column bounds*/
369  columnLower = (double *) malloc(n_var*sizeof(double));
370  columnUpper = (double *) malloc(n_var*sizeof(double));
371#define COIN_DBL_MAX DBL_MAX
372  for (i=0;i<n_var;i++) {
373    columnLower[i]=LUv[2*i];
374    if (columnLower[i]<= negInfinity)
375      columnLower[i]=-COIN_DBL_MAX;
376    columnUpper[i]=LUv[2*i+1];
377    if (columnUpper[i]>= Infinity)
378      columnUpper[i]=COIN_DBL_MAX;
379  }
380  /* Row bounds*/
381  rowLower = (double *) malloc(n_con*sizeof(double));
382  rowUpper = (double *) malloc(n_con*sizeof(double));
383  for (i=0;i<n_con;i++) {
384    rowLower[i]=LUrhs[2*i];
385    if (rowLower[i]<= negInfinity)
386      rowLower[i]=-COIN_DBL_MAX;
387    rowUpper[i]=LUrhs[2*i+1];
388    if (rowUpper[i]>= Infinity)
389      rowUpper[i]=COIN_DBL_MAX;
390  }
391  info->numberRows=n_con;
392  info->numberColumns=n_var;
393  info->numberElements=nzc;;
394  info->numberBinary=nbv;
395  info->numberIntegers=niv;
396  info->objective=obj;
397  info->rowLower=rowLower;
398  info->rowUpper=rowUpper;
399  info->columnLower=columnLower;
400  info->columnUpper=columnUpper;
401  info->starts=A_colstarts;
402  /*A_colstarts=NULL;*/
403  info->rows=A_rownos;
404  /*A_rownos=NULL;*/
405  info->elements=A_vals;
406  /*A_vals=NULL;*/
407  info->primalSolution=NULL;
408  /* put in primalSolution if exists */
409  if (X0) {
410    info->primalSolution=(double *) malloc(n_var*sizeof(double));
411    memcpy(info->primalSolution,X0,n_var*sizeof(double));
412  }
413  info->dualSolution=NULL;
414  if (niv+nbv>0)
415    mip_stuff(); // get any extra info
416  if ((!(niv+nbv)&&(csd->kind & ASL_Sufkind_input))
417      ||(rsd->kind & ASL_Sufkind_input)) {
418    /* convert status - need info on map */
419    static int map[] = {1, 3, 1, 1, 2, 1, 1};
420    stat_map(info->columnStatus, n_var, map, 6, "incoming columnStatus");
421    stat_map(info->rowStatus, n_con, map, 6, "incoming rowStatus");
422  } else {
423    /* all slack basis */
424    // leave status for output */
425#if 0
426    free(info->rowStatus);
427    info->rowStatus=NULL;
428    free(info->columnStatus);
429    info->columnStatus=NULL;
430#endif
431  }
432  /* add -solve - unless something there already
433   - also check for sleep=yes */
434  {
435    int found=0;
436    int foundLog=0;
437    int foundSleep=0;
438    const char * something[]={"solve","branch","duals","primals","user"};
439    for (i=0;i<info->numberArguments;i++) {
440      unsigned int j;
441      const char * argument = info->arguments[i];
442      for (j=0;j<sizeof(something)/sizeof(char *);j++) {
443        const char * check = something[j];
444        if (!strncmp(argument,check,sizeof(check))) {
445          found=(int)(j+1);
446        } else if (!strncmp(argument,"log",3)) {
447          foundLog=1;
448        } else if (!strncmp(argument,"sleep",5)) {
449          foundSleep=1;
450        }
451      }
452    }
453    if (foundLog) {
454      /* print options etc */
455      for (i=0;i<saveArgc;i++)
456        printf("%s ",saveArgv[i]);
457      printf("\n");
458      if (environment)
459        printf("env %s\n",environment);
460      /*printf("%d rows %d columns %d elements\n",n_con,n_var,nzc);*/
461    }
462    if (!found) {
463      if (!strlen(algFound)) {
464        info->arguments=(char **) realloc(info->arguments,(info->numberArguments+1)*sizeof(char *));
465        info->arguments[info->numberArguments++]=strdup("-solve");
466      } else {
467        // use algorithm from keyword
468        info->arguments=(char **) realloc(info->arguments,(info->numberArguments+1)*sizeof(char *));
469        info->arguments[info->numberArguments++]=strdup(algFound);
470      }
471    }
472    if (foundSleep) {
473      /* let user copy .nl file */
474      fprintf(stderr,"You can copy .nl file %s for debug purposes or attach debugger\n",saveArgv[1]);
475      fprintf(stderr,"Type q to quit, anything else to continue\n");
476      char getChar = getc(stdin);
477      if (getChar=='q'||getChar=='Q')
478        exit(1);
479    }
480  }
481  /* add -quit */
482  info->arguments=(char **) realloc(info->arguments,(info->numberArguments+1)*sizeof(char *));
483  info->arguments[info->numberArguments++]=strdup("-quit");
484  return 0;
485}
486void freeArrays1(ampl_info * info)
487{
488  free(info->objective);
489  info->objective=NULL;
490  free(info->rowLower);
491  info->rowLower=NULL;
492  free(info->rowUpper);
493  info->rowUpper=NULL;
494  free(info->columnLower);
495  info->columnLower=NULL;
496  free(info->columnUpper);
497  info->columnUpper=NULL;
498  /* this one not freed by ASL_free */
499  free(info->elements);
500  info->elements=NULL;
501  free(info->primalSolution);
502  info->primalSolution=NULL;
503  free(info->dualSolution);
504  info->dualSolution=NULL;
505  /*free(info->rowStatus);
506  info->rowStatus=NULL;
507  free(info->columnStatus);
508  info->columnStatus=NULL;*/
509}
510void freeArrays2(ampl_info * info)
511{
512  free(info->primalSolution);
513  info->primalSolution=NULL;
514  free(info->dualSolution);
515  info->dualSolution=NULL;
516  free(info->rowStatus);
517  info->rowStatus=NULL;
518  free(info->columnStatus);
519  info->columnStatus=NULL;
520  free(info->priorities);
521  info->priorities=NULL;
522  free(info->branchDirection);
523  info->branchDirection=NULL;
524  free(info->pseudoDown);
525  info->pseudoDown=NULL;
526  free(info->pseudoUp);
527  info->pseudoUp=NULL;
528  free(info->sosType);
529  info->sosType=NULL;
530  free(info->sosPriority);
531  info->sosPriority=NULL;
532  free(info->sosStart);
533  info->sosStart=NULL;
534  free(info->sosIndices);
535  info->sosIndices=NULL;
536  free(info->sosReference);
537  info->sosReference=NULL;
538  ASL_free(&asl);
539}
540void freeArgs(ampl_info * info)
541{
542  int i;
543  for ( i=0; i<info->numberArguments;i++)
544    free(info->arguments[i]);
545  free(info->arguments);
546}
547int ampl_obj_prec()
548{
549  return obj_prec();
550}
551void writeAmpl(ampl_info * info)
552{
553  char buf[1000];
554  typedef struct { const char *msg; int code; int wantObj; } Sol_info;
555  static Sol_info solinfo[] = {
556    { "optimal solution",                       000, 1 },
557    { "infeasible",                             200, 1 },
558    { "unbounded",                              300, 0 },
559    { "iteration limit etc",                    400, 1 },
560    { "solution limit",                         401, 1 },
561    { "ran out of space",                       500, 0 },
562    { "status unknown",                         501, 1 },
563    { "bug!",                                   502, 0 },
564    { "best MIP solution so far restored",      101, 1 },
565    { "failed to restore best MIP solution",    503, 1 },
566    { "optimal (?) solution",                   100, 1 }
567  };
568  /* convert status - need info on map */
569  static int map[] = {0, 3, 4, 1};
570  sprintf(buf,"%s %s",Oinfo.bsname,info->buffer);
571  solve_result_num = solinfo[info->problemStatus].code;
572  if (info->columnStatus) {
573    stat_map(info->columnStatus, n_var, map, 4, "outgoing columnStatus");
574    stat_map(info->rowStatus, n_con, map, 4, "outgoing rowStatus");
575    suf_iput("sstatus", ASL_Sufkind_var, info->columnStatus);
576    suf_iput("sstatus", ASL_Sufkind_con, info->rowStatus);
577  }
578  write_sol(buf,info->primalSolution,info->dualSolution,&Oinfo);
579}
580/* Read a problem from AMPL nl file
581 */
582CoinModel::CoinModel( int nonLinear, const char * fileName)
583 :  numberRows_(0),
584    maximumRows_(0),
585    numberColumns_(0),
586    maximumColumns_(0),
587    numberElements_(0),
588    maximumElements_(0),
589    numberQuadraticElements_(0),
590    maximumQuadraticElements_(0),
591    optimizationDirection_(1.0),
592    objectiveOffset_(0.0),
593    rowLower_(NULL),
594    rowUpper_(NULL),
595    rowType_(NULL),
596    objective_(NULL),
597    columnLower_(NULL),
598    columnUpper_(NULL),
599    integerType_(NULL),
600    columnType_(NULL),
601    start_(NULL),
602    elements_(NULL),
603    quadraticElements_(NULL),
604    sortIndices_(NULL),
605    sortElements_(NULL),
606    sortSize_(0),
607    sizeAssociated_(0),
608    associated_(NULL),
609    logLevel_(0),
610    type_(-1),
611    links_(0)
612{
613  problemName_ = strdup("");
614  int status=0;
615  if (!strcmp(fileName,"-")||!strcmp(fileName,"stdin")) {
616    // stdin
617  } else {
618    std::string name=fileName;
619    bool readable = fileCoinReadable(name);
620    if (!readable) {
621      std::cerr<<"Unable to open file "
622               <<fileName<<std::endl;
623      status = -1;
624    }
625  }
626  if (!status) {
627    gdb(nonLinear,fileName);
628  }
629}
630 static real
631qterm(ASL *asl, fint *colq, fint *rowq, real *delsq)
632{
633  double t, t1, *x, *x0, *xe;
634  fint *rq0, *rqe;
635 
636  t = 0.;
637  x0 = x = X0;
638  xe = x + n_var;
639  rq0 = rowq;
640  while(x < xe) {
641    t1 = *x++;
642    rqe = rq0 + *++colq;
643    while(rowq < rqe)
644      t += t1*x0[*rowq++]**delsq++;
645  }
646  return 0.5 * t;
647}
648
649void
650CoinModel::gdb( int nonLinear, const char * fileName)
651{
652  ograd *og=NULL;
653  int i;
654  SufDesc *csd=NULL;
655  SufDesc *rsd=NULL;
656  /*bool *basis, *lower;*/
657  /*double *LU, *c, lb, objadj, *rshift, *shift, t, ub, *x, *x0, *x1;*/
658  //char tempBuffer[20];
659  double * objective=NULL;
660  double * columnLower=NULL;
661  double * columnUpper=NULL;
662  double * rowLower=NULL;
663  double * rowUpper=NULL;
664  int * columnStatus=NULL;
665  int * rowStatus=NULL;
666  int numberRows=-1;
667  int numberColumns=-1;
668  int numberElements=-1;
669  int numberBinary=-1;
670  int numberIntegers=-1;
671  double * primalSolution=NULL;
672  double direction=1.0;
673  char * stub = strdup(fileName);
674  CoinPackedMatrix matrixByRow;
675  fint ** colqp=NULL;
676  int *z = NULL;
677  if (nonLinear==0) {
678    // linear
679    asl = ASL_alloc(ASL_read_f);
680    nl = jac0dim(stub, 0);
681    free(stub);
682    suf_declare(suftab, sizeof(suftab)/sizeof(SufDecl));
683   
684    /* set A_vals to get the constraints column-wise (malloc so can be freed) */
685    A_vals = (double *) malloc(nzc*sizeof(double));
686    if (!A_vals) {
687      printf("no memory\n");
688      return ;
689    }
690    /* say we want primal solution */
691    want_xpi0=1;
692    /* for basis info */
693    columnStatus = (int *) malloc(n_var*sizeof(int));
694    rowStatus = (int *) malloc(n_con*sizeof(int));
695    csd = suf_iput("sstatus", ASL_Sufkind_var, columnStatus);
696    rsd = suf_iput("sstatus", ASL_Sufkind_con, rowStatus);
697    /* read linear model*/
698    f_read(nl,0);
699    // see if any sos
700    if (true) {
701      char *sostype;
702      int nsosnz, *sosbeg, *sosind, * sospri;
703      double *sosref;
704      int nsos;
705      int i = ASL_suf_sos_explict_free;
706      int copri[2], **p_sospri;
707      copri[0] = 0;
708      copri[1] = 0;
709      p_sospri = &sospri;
710      nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
711                     &sosbeg, &sosind, &sosref);
712      if (nsos) {
713        abort();
714#if 0
715        info->numberSos=nsos;
716        info->sosType = (char *) malloc(nsos);
717        info->sosPriority = (int *) malloc(nsos*sizeof(int));
718        info->sosStart = (int *) malloc((nsos+1)*sizeof(int));
719        info->sosIndices = (int *) malloc(nsosnz*sizeof(int));
720        info->sosReference = (double *) malloc(nsosnz*sizeof(double));
721        sos_kludge(nsos, sosbeg, sosref,sosind);
722        for (int i=0;i<nsos;i++) {
723          int ichar = sostype[i];
724          assert (ichar=='1'||ichar=='2');
725          info->sosType[i]=ichar-'0';
726        }       
727        memcpy(info->sosPriority,sospri,nsos*sizeof(int));
728        memcpy(info->sosStart,sosbeg,(nsos+1)*sizeof(int));
729        memcpy(info->sosIndices,sosind,nsosnz*sizeof(int));
730        memcpy(info->sosReference,sosref,nsosnz*sizeof(double));
731#endif
732      }
733    }
734   
735    /*sos_finish(&specialOrderedInfo, 0, &j, 0, 0, 0, 0, 0);*/
736    //Oinfo.uinfo = tempBuffer;
737    //if (getopts(argv, &Oinfo))
738    //return 1;
739    /* objective*/
740    objective = (double *) malloc(n_var*sizeof(double));
741    for (i=0;i<n_var;i++)
742      objective[i]=0.0;;
743    if (n_obj) {
744      for (og = Ograd[0];og;og = og->next)
745        objective[og->varno] = og->coef;
746    }
747    if (objtype[0])
748      direction=-1.0;
749    else
750      direction=1.0;
751    objectiveOffset_=objconst(0);
752    /* Column bounds*/
753    columnLower = (double *) malloc(n_var*sizeof(double));
754    columnUpper = (double *) malloc(n_var*sizeof(double));
755#define COIN_DBL_MAX DBL_MAX
756    for (i=0;i<n_var;i++) {
757      columnLower[i]=LUv[2*i];
758      if (columnLower[i]<= negInfinity)
759        columnLower[i]=-COIN_DBL_MAX;
760      columnUpper[i]=LUv[2*i+1];
761      if (columnUpper[i]>= Infinity)
762        columnUpper[i]=COIN_DBL_MAX;
763    }
764    /* Row bounds*/
765    rowLower = (double *) malloc(n_con*sizeof(double));
766    rowUpper = (double *) malloc(n_con*sizeof(double));
767    for (i=0;i<n_con;i++) {
768      rowLower[i]=LUrhs[2*i];
769      if (rowLower[i]<= negInfinity)
770        rowLower[i]=-COIN_DBL_MAX;
771      rowUpper[i]=LUrhs[2*i+1];
772      if (rowUpper[i]>= Infinity)
773        rowUpper[i]=COIN_DBL_MAX;
774    }
775    numberRows=n_con;
776    numberColumns=n_var;
777    numberElements=nzc;;
778    numberBinary=nbv;
779    numberIntegers=niv;
780    /* put in primalSolution if exists */
781    if (X0) {
782      primalSolution=(double *) malloc(n_var*sizeof(double));
783      memcpy( primalSolution,X0,n_var*sizeof(double));
784    }
785    //double * dualSolution=NULL;
786    if (niv+nbv>0)
787      mip_stuff(); // get any extra info
788    if ((!(niv+nbv)&&(csd->kind & ASL_Sufkind_input))
789        ||(rsd->kind & ASL_Sufkind_input)) {
790      /* convert status - need info on map */
791      static int map[] = {1, 3, 1, 1, 2, 1, 1};
792      stat_map(columnStatus, n_var, map, 6, "incoming columnStatus");
793      stat_map(rowStatus, n_con, map, 6, "incoming rowStatus");
794    } else {
795      /* all slack basis */
796      // leave status for output */
797#if 0
798      free(rowStatus);
799      rowStatus=NULL;
800      free(columnStatus);
801      columnStatus=NULL;
802#endif
803    }
804    CoinPackedMatrix columnCopy(true,numberRows,numberColumns,numberElements,
805                                A_vals,A_rownos,A_colstarts,NULL);
806    matrixByRow.reverseOrderedCopyOf(columnCopy);
807  } else if (nonLinear==1) {
808    // quadratic
809    asl = ASL_alloc(ASL_read_fg);
810    nl = jac0dim(stub, (long) strlen(stub));
811    free(stub);
812    /* read  model*/
813    X0 = (double*) malloc(n_var*sizeof(double));
814    CoinZeroN(X0,n_var);
815    qp_read(nl,0);
816    assert (n_obj==1);
817    int nz = 1 + n_con;
818    colqp = (fint**) malloc(nz*(2*sizeof(int*)
819                                      + sizeof(double*)));
820    fint ** rowqp = colqp + nz;
821    double ** delsqp = (double **)(rowqp + nz);
822    z = (int*) malloc(nz*sizeof(int));
823    for (i=0;i<=n_con;i++) {
824      z[i] = nqpcheck(-i, rowqp+i, colqp+i, delsqp+i);
825    }
826    qp_opify();
827    /* objective*/
828    objective = (double *) malloc(n_var*sizeof(double));
829    for (i=0;i<n_var;i++)
830      objective[i]=0.0;;
831    if (n_obj) {
832      for (og = Ograd[0];og;og = og->next)
833        objective[og->varno] = og->coef;
834    }
835    if (objtype[0])
836      direction=-1.0;
837    else
838      direction=1.0;
839    objectiveOffset_=objconst(0);
840    /* Column bounds*/
841    columnLower = (double *) malloc(n_var*sizeof(double));
842    columnUpper = (double *) malloc(n_var*sizeof(double));
843    for (i=0;i<n_var;i++) {
844      columnLower[i]=LUv[2*i];
845      if (columnLower[i]<= negInfinity)
846        columnLower[i]=-COIN_DBL_MAX;
847      columnUpper[i]=LUv[2*i+1];
848      if (columnUpper[i]>= Infinity)
849        columnUpper[i]=COIN_DBL_MAX;
850    }
851    // Build by row from scratch
852    //matrixByRow.reserve(n_var,nzc,true);
853    // say row orderded
854    matrixByRow.transpose();
855    /* Row bounds*/
856    rowLower = (double *) malloc(n_con*sizeof(double));
857    rowUpper = (double *) malloc(n_con*sizeof(double));
858    int * column = new int [n_var];
859    double * element = new double [n_var];
860    for (i=0;i<n_con;i++) {
861      rowLower[i]=LUrhs[2*i];
862      if (rowLower[i]<= negInfinity)
863        rowLower[i]=-COIN_DBL_MAX;
864      rowUpper[i]=LUrhs[2*i+1];
865      if (rowUpper[i]>= Infinity)
866        rowUpper[i]=COIN_DBL_MAX;
867      int k=0;
868      for(cgrad * cg = Cgrad[i]; cg; cg = cg->next) {
869        column[k]=cg->varno;
870        element[k++]=cg->coef;
871      }
872      matrixByRow.appendRow(k,column,element);
873    }
874    numberRows=n_con;
875    numberColumns=n_var;
876    numberElements=nzc;;
877    numberBinary=nbv;
878    numberIntegers=niv+nlvci;
879    //double * dualSolution=NULL;
880  } else {
881    abort();
882  }
883  // set problem name
884  free(problemName_);
885  problemName_=strdup("???");
886
887  // Build by row from scratch
888  const double * element = matrixByRow.getElements();
889  const int * column = matrixByRow.getIndices();
890  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
891  const int * rowLength = matrixByRow.getVectorLengths();
892  for (i=0;i<numberRows;i++) {
893    addRow(rowLength[i],column+rowStart[i],
894           element+rowStart[i],rowLower[i],rowUpper[i]);
895  }
896  // Now do column part
897  for (i=0;i<numberColumns;i++) {
898    setColumnBounds(i,columnLower[i],columnUpper[i]);
899    setColumnObjective(i,objective[i]);
900  }
901  for ( i=numberColumns-numberBinary-numberIntegers;
902        i<numberColumns;i++) {
903      setColumnIsInteger(i,true);;
904  }
905  // do names
906  int iRow;
907  for (iRow=0;iRow<numberRows_;iRow++) {
908    char name[9];
909    sprintf(name,"r%7.7d",iRow);
910    setRowName(iRow,name);
911  }
912   
913  int iColumn;
914  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
915    char name[9];
916    sprintf(name,"c%7.7d",iColumn);
917    setColumnName(iColumn,name);
918  }
919  if (colqp) {
920    // add in quadratic
921    int nz = 1 + n_con;
922    fint ** rowqp = colqp + nz;
923    double ** delsqp = (double **)(rowqp + nz);
924    for (i=0;i<=n_con;i++) {
925      int nels = z[i];
926      if (nels) {
927        double * element = delsqp[i];
928        int * start = (int *) colqp[i];
929        int * row = (int *) rowqp[i];
930#if 0
931        printf("%d quadratic els\n",nels);
932        for (int j=0;j<n_var;j++) {
933          for (int k=start[j];k<start[j+1];k++)
934            printf("%d %d %g\n",j,row[k],element[k]);
935        }
936#endif
937        if (i) {
938          int iRow = i-1;
939          for (int j=0;j<n_var;j++) {
940            for (int k=start[j];k<start[j+1];k++) {
941              int kColumn = row[k];
942              double value = element[k];
943              // ampl gives twice with assumed 0.5
944              if (kColumn>j)
945                continue;
946              else if (kColumn==j)
947                value *= 0.5;
948              const char * expr = getElementAsString(iRow,j);
949              double constant=0.0;
950              bool linear;
951              if (expr&&strcmp(expr,"Numeric")) {
952                linear=false;
953              } else {
954                constant = getElement(iRow,j);
955                linear=true;
956              }
957              char temp[100];
958              char temp2[30];
959              if (value==1.0) 
960                sprintf(temp2,"c%7.7d",kColumn);
961              else
962                sprintf(temp2,"%g*c%7.7d",value,kColumn);
963              if (linear) {
964                if (!constant)
965                  strcpy(temp,temp2);
966                else if (value>0.0) 
967                  sprintf(temp,"%g+%s",constant,temp2);
968                else
969                  sprintf(temp,"%g%s",constant,temp2);
970              } else {
971                if (value>0.0) 
972                  sprintf(temp,"%s+%s",expr,temp2);
973                else
974                  sprintf(temp,"%s%s",expr,temp2);
975              }
976              setElement(iRow,j,temp);
977              printf("el for row %d column c%7.7d is %s\n",iRow,j,temp);
978            }
979          }
980        } else {
981          // objective
982          for (int j=0;j<n_var;j++) {
983            for (int k=start[j];k<start[j+1];k++) {
984              int kColumn = row[k];
985              double value = element[k];
986              // ampl gives twice with assumed 0.5
987              if (kColumn>j)
988                continue;
989              else if (kColumn==j)
990                value *= 0.5;
991              const char * expr = getColumnObjectiveAsString(j);
992              double constant=0.0;
993              bool linear;
994              if (expr&&strcmp(expr,"Numeric")) {
995                linear=false;
996              } else {
997                constant = getColumnObjective(j);
998                linear=true;
999              }
1000              char temp[100];
1001              char temp2[30];
1002              if (value==1.0) 
1003                sprintf(temp2,"c%7.7d",kColumn);
1004              else
1005                sprintf(temp2,"%g*c%7.7d",value,kColumn);
1006              if (linear) {
1007                if (!constant)
1008                  strcpy(temp,temp2);
1009                else if (value>0.0) 
1010                  sprintf(temp,"%g+%s",constant,temp2);
1011                else
1012                  sprintf(temp,"%g%s",constant,temp2);
1013              } else {
1014                if (value>0.0) 
1015                  sprintf(temp,"%s+%s",expr,temp2);
1016                else
1017                  sprintf(temp,"%s%s",expr,temp2);
1018              }
1019              setObjective(j,temp);
1020              printf("el for objective column c%7.7d is %s\n",j,temp);
1021            }
1022          }
1023        }
1024      }
1025    }
1026  }
1027  // see if any sos
1028  {
1029    char *sostype;
1030    int nsosnz, *sosbeg, *sosind, * sospri;
1031    double *sosref;
1032    int nsos;
1033    int i = ASL_suf_sos_explict_free;
1034    int copri[2], **p_sospri;
1035    copri[0] = 0;
1036    copri[1] = 0;
1037    p_sospri = &sospri;
1038    nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
1039                   &sosbeg, &sosind, &sosref);
1040    if (nsos) {
1041      numberSOS_=nsos;
1042      typeSOS_ = new int [numberSOS_];
1043      prioritySOS_ = new int [numberSOS_];
1044      startSOS_ = new int [numberSOS_+1];
1045      memberSOS_ = new int[nsosnz];
1046      referenceSOS_ = new double [nsosnz];
1047      sos_kludge(nsos, sosbeg, sosref,sosind);
1048      for (int i=0;i<nsos;i++) {
1049        int ichar = sostype[i];
1050        assert (ichar=='1'||ichar=='2');
1051        typeSOS_[i]=ichar-'0';
1052      } 
1053      memcpy(prioritySOS_,sospri,nsos*sizeof(int));
1054      memcpy(startSOS_,sosbeg,(nsos+1)*sizeof(int));
1055      memcpy(memberSOS_,sosind,nsosnz*sizeof(int));
1056      memcpy(referenceSOS_,sosref,nsosnz*sizeof(double));
1057    }
1058  }
1059}
Note: See TracBrowser for help on using the repository browser.