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

Last change on this file since 492 was 492, checked in by forrest, 14 years ago

for nonlinear

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 30.8 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    numberSOS_(0),
610    startSOS_(NULL),
611    memberSOS_(NULL),
612    typeSOS_(NULL),
613    prioritySOS_(NULL),
614    referenceSOS_(NULL),
615    priority_(NULL),
616    logLevel_(0),
617    type_(-1),
618    links_(0)
619{
620  problemName_ = strdup("");
621  int status=0;
622  if (!strcmp(fileName,"-")||!strcmp(fileName,"stdin")) {
623    // stdin
624  } else {
625    std::string name=fileName;
626    bool readable = fileCoinReadable(name);
627    if (!readable) {
628      std::cerr<<"Unable to open file "
629               <<fileName<<std::endl;
630      status = -1;
631    }
632  }
633  if (!status) {
634    gdb(nonLinear,fileName);
635  }
636}
637 static real
638qterm(ASL *asl, fint *colq, fint *rowq, real *delsq)
639{
640  double t, t1, *x, *x0, *xe;
641  fint *rq0, *rqe;
642 
643  t = 0.;
644  x0 = x = X0;
645  xe = x + n_var;
646  rq0 = rowq;
647  while(x < xe) {
648    t1 = *x++;
649    rqe = rq0 + *++colq;
650    while(rowq < rqe)
651      t += t1*x0[*rowq++]**delsq++;
652  }
653  return 0.5 * t;
654}
655
656void
657CoinModel::gdb( int nonLinear, const char * fileName)
658{
659  ograd *og=NULL;
660  int i;
661  SufDesc *csd=NULL;
662  SufDesc *rsd=NULL;
663  /*bool *basis, *lower;*/
664  /*double *LU, *c, lb, objadj, *rshift, *shift, t, ub, *x, *x0, *x1;*/
665  //char tempBuffer[20];
666  double * objective=NULL;
667  double * columnLower=NULL;
668  double * columnUpper=NULL;
669  double * rowLower=NULL;
670  double * rowUpper=NULL;
671  int * columnStatus=NULL;
672  int * rowStatus=NULL;
673  int numberRows=-1;
674  int numberColumns=-1;
675  int numberElements=-1;
676  int numberBinary=-1;
677  int numberIntegers=-1;
678  double * primalSolution=NULL;
679  double direction=1.0;
680  char * stub = strdup(fileName);
681  CoinPackedMatrix matrixByRow;
682  fint ** colqp=NULL;
683  int *z = NULL;
684  if (nonLinear==0) {
685    // linear
686    asl = ASL_alloc(ASL_read_f);
687    nl = jac0dim(stub, 0);
688    free(stub);
689    suf_declare(suftab, sizeof(suftab)/sizeof(SufDecl));
690   
691    /* set A_vals to get the constraints column-wise (malloc so can be freed) */
692    A_vals = (double *) malloc(nzc*sizeof(double));
693    if (!A_vals) {
694      printf("no memory\n");
695      return ;
696    }
697    /* say we want primal solution */
698    want_xpi0=1;
699    /* for basis info */
700    columnStatus = (int *) malloc(n_var*sizeof(int));
701    rowStatus = (int *) malloc(n_con*sizeof(int));
702    csd = suf_iput("sstatus", ASL_Sufkind_var, columnStatus);
703    rsd = suf_iput("sstatus", ASL_Sufkind_con, rowStatus);
704    /* read linear model*/
705    f_read(nl,0);
706    // see if any sos
707    if (true) {
708      char *sostype;
709      int nsosnz, *sosbeg, *sosind, * sospri;
710      double *sosref;
711      int nsos;
712      int i = ASL_suf_sos_explict_free;
713      int copri[2], **p_sospri;
714      copri[0] = 0;
715      copri[1] = 0;
716      p_sospri = &sospri;
717      nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
718                     &sosbeg, &sosind, &sosref);
719      if (nsos) {
720        abort();
721#if 0
722        info->numberSos=nsos;
723        info->sosType = (char *) malloc(nsos);
724        info->sosPriority = (int *) malloc(nsos*sizeof(int));
725        info->sosStart = (int *) malloc((nsos+1)*sizeof(int));
726        info->sosIndices = (int *) malloc(nsosnz*sizeof(int));
727        info->sosReference = (double *) malloc(nsosnz*sizeof(double));
728        sos_kludge(nsos, sosbeg, sosref,sosind);
729        for (int i=0;i<nsos;i++) {
730          int ichar = sostype[i];
731          assert (ichar=='1'||ichar=='2');
732          info->sosType[i]=ichar-'0';
733        }       
734        memcpy(info->sosPriority,sospri,nsos*sizeof(int));
735        memcpy(info->sosStart,sosbeg,(nsos+1)*sizeof(int));
736        memcpy(info->sosIndices,sosind,nsosnz*sizeof(int));
737        memcpy(info->sosReference,sosref,nsosnz*sizeof(double));
738#endif
739      }
740    }
741   
742    /*sos_finish(&specialOrderedInfo, 0, &j, 0, 0, 0, 0, 0);*/
743    //Oinfo.uinfo = tempBuffer;
744    //if (getopts(argv, &Oinfo))
745    //return 1;
746    /* objective*/
747    objective = (double *) malloc(n_var*sizeof(double));
748    for (i=0;i<n_var;i++)
749      objective[i]=0.0;;
750    if (n_obj) {
751      for (og = Ograd[0];og;og = og->next)
752        objective[og->varno] = og->coef;
753    }
754    if (objtype[0])
755      direction=-1.0;
756    else
757      direction=1.0;
758    objectiveOffset_=objconst(0);
759    /* Column bounds*/
760    columnLower = (double *) malloc(n_var*sizeof(double));
761    columnUpper = (double *) malloc(n_var*sizeof(double));
762#define COIN_DBL_MAX DBL_MAX
763    for (i=0;i<n_var;i++) {
764      columnLower[i]=LUv[2*i];
765      if (columnLower[i]<= negInfinity)
766        columnLower[i]=-COIN_DBL_MAX;
767      columnUpper[i]=LUv[2*i+1];
768      if (columnUpper[i]>= Infinity)
769        columnUpper[i]=COIN_DBL_MAX;
770    }
771    /* Row bounds*/
772    rowLower = (double *) malloc(n_con*sizeof(double));
773    rowUpper = (double *) malloc(n_con*sizeof(double));
774    for (i=0;i<n_con;i++) {
775      rowLower[i]=LUrhs[2*i];
776      if (rowLower[i]<= negInfinity)
777        rowLower[i]=-COIN_DBL_MAX;
778      rowUpper[i]=LUrhs[2*i+1];
779      if (rowUpper[i]>= Infinity)
780        rowUpper[i]=COIN_DBL_MAX;
781    }
782    numberRows=n_con;
783    numberColumns=n_var;
784    numberElements=nzc;;
785    numberBinary=nbv;
786    numberIntegers=niv;
787    /* put in primalSolution if exists */
788    if (X0) {
789      primalSolution=(double *) malloc(n_var*sizeof(double));
790      memcpy( primalSolution,X0,n_var*sizeof(double));
791    }
792    //double * dualSolution=NULL;
793    if (niv+nbv>0)
794      mip_stuff(); // get any extra info
795    if ((!(niv+nbv)&&(csd->kind & ASL_Sufkind_input))
796        ||(rsd->kind & ASL_Sufkind_input)) {
797      /* convert status - need info on map */
798      static int map[] = {1, 3, 1, 1, 2, 1, 1};
799      stat_map(columnStatus, n_var, map, 6, "incoming columnStatus");
800      stat_map(rowStatus, n_con, map, 6, "incoming rowStatus");
801    } else {
802      /* all slack basis */
803      // leave status for output */
804#if 0
805      free(rowStatus);
806      rowStatus=NULL;
807      free(columnStatus);
808      columnStatus=NULL;
809#endif
810    }
811    CoinPackedMatrix columnCopy(true,numberRows,numberColumns,numberElements,
812                                A_vals,A_rownos,A_colstarts,NULL);
813    matrixByRow.reverseOrderedCopyOf(columnCopy);
814  } else if (nonLinear==1) {
815    // quadratic
816    asl = ASL_alloc(ASL_read_fg);
817    nl = jac0dim(stub, (long) strlen(stub));
818    free(stub);
819    suf_declare(suftab, sizeof(suftab)/sizeof(SufDecl));
820    /* read  model*/
821    X0 = (double*) malloc(n_var*sizeof(double));
822    CoinZeroN(X0,n_var);
823    qp_read(nl,0);
824    assert (n_obj==1);
825    int nz = 1 + n_con;
826    colqp = (fint**) malloc(nz*(2*sizeof(int*)
827                                      + sizeof(double*)));
828    fint ** rowqp = colqp + nz;
829    double ** delsqp = (double **)(rowqp + nz);
830    z = (int*) malloc(nz*sizeof(int));
831    for (i=0;i<=n_con;i++) {
832      z[i] = nqpcheck(-i, rowqp+i, colqp+i, delsqp+i);
833    }
834    qp_opify();
835    /* objective*/
836    objective = (double *) malloc(n_var*sizeof(double));
837    for (i=0;i<n_var;i++)
838      objective[i]=0.0;;
839    if (n_obj) {
840      for (og = Ograd[0];og;og = og->next)
841        objective[og->varno] = og->coef;
842    }
843    if (objtype[0])
844      direction=-1.0;
845    else
846      direction=1.0;
847    objectiveOffset_=objconst(0);
848    /* Column bounds*/
849    columnLower = (double *) malloc(n_var*sizeof(double));
850    columnUpper = (double *) malloc(n_var*sizeof(double));
851    for (i=0;i<n_var;i++) {
852      columnLower[i]=LUv[2*i];
853      if (columnLower[i]<= negInfinity)
854        columnLower[i]=-COIN_DBL_MAX;
855      columnUpper[i]=LUv[2*i+1];
856      if (columnUpper[i]>= Infinity)
857        columnUpper[i]=COIN_DBL_MAX;
858    }
859    // Build by row from scratch
860    //matrixByRow.reserve(n_var,nzc,true);
861    // say row orderded
862    matrixByRow.transpose();
863    /* Row bounds*/
864    rowLower = (double *) malloc(n_con*sizeof(double));
865    rowUpper = (double *) malloc(n_con*sizeof(double));
866    int * column = new int [n_var];
867    double * element = new double [n_var];
868    for (i=0;i<n_con;i++) {
869      rowLower[i]=LUrhs[2*i];
870      if (rowLower[i]<= negInfinity)
871        rowLower[i]=-COIN_DBL_MAX;
872      rowUpper[i]=LUrhs[2*i+1];
873      if (rowUpper[i]>= Infinity)
874        rowUpper[i]=COIN_DBL_MAX;
875      int k=0;
876      for(cgrad * cg = Cgrad[i]; cg; cg = cg->next) {
877        column[k]=cg->varno;
878        element[k++]=cg->coef;
879      }
880      matrixByRow.appendRow(k,column,element);
881    }
882    numberRows=n_con;
883    numberColumns=n_var;
884    numberElements=nzc;;
885    numberBinary=nbv;
886    numberIntegers=niv+nlvci;
887    //double * dualSolution=NULL;
888  } else {
889    abort();
890  }
891  // set problem name
892  free(problemName_);
893  problemName_=strdup("???");
894
895  // Build by row from scratch
896  const double * element = matrixByRow.getElements();
897  const int * column = matrixByRow.getIndices();
898  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
899  const int * rowLength = matrixByRow.getVectorLengths();
900  for (i=0;i<numberRows;i++) {
901    addRow(rowLength[i],column+rowStart[i],
902           element+rowStart[i],rowLower[i],rowUpper[i]);
903  }
904  // Now do column part
905  for (i=0;i<numberColumns;i++) {
906    setColumnBounds(i,columnLower[i],columnUpper[i]);
907    setColumnObjective(i,objective[i]);
908  }
909  for ( i=numberColumns-numberBinary-numberIntegers;
910        i<numberColumns;i++) {
911      setColumnIsInteger(i,true);;
912  }
913  // do names
914  int iRow;
915  for (iRow=0;iRow<numberRows_;iRow++) {
916    char name[9];
917    sprintf(name,"r%7.7d",iRow);
918    setRowName(iRow,name);
919  }
920   
921  int iColumn;
922  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
923    char name[9];
924    sprintf(name,"c%7.7d",iColumn);
925    setColumnName(iColumn,name);
926  }
927  if (colqp) {
928    // add in quadratic
929    int nz = 1 + n_con;
930    fint ** rowqp = colqp + nz;
931    double ** delsqp = (double **)(rowqp + nz);
932    for (i=0;i<=n_con;i++) {
933      int nels = z[i];
934      if (nels) {
935        double * element = delsqp[i];
936        int * start = (int *) colqp[i];
937        int * row = (int *) rowqp[i];
938#if 0
939        printf("%d quadratic els\n",nels);
940        for (int j=0;j<n_var;j++) {
941          for (int k=start[j];k<start[j+1];k++)
942            printf("%d %d %g\n",j,row[k],element[k]);
943        }
944#endif
945        if (i) {
946          int iRow = i-1;
947          for (int j=0;j<n_var;j++) {
948            for (int k=start[j];k<start[j+1];k++) {
949              int kColumn = row[k];
950              double value = element[k];
951              // ampl gives twice with assumed 0.5
952              if (kColumn>j)
953                continue;
954              else if (kColumn==j)
955                value *= 0.5;
956              const char * expr = getElementAsString(iRow,j);
957              double constant=0.0;
958              bool linear;
959              if (expr&&strcmp(expr,"Numeric")) {
960                linear=false;
961              } else {
962                constant = getElement(iRow,j);
963                linear=true;
964              }
965              char temp[1000];
966              char temp2[30];
967              if (value==1.0) 
968                sprintf(temp2,"c%7.7d",kColumn);
969              else
970                sprintf(temp2,"%g*c%7.7d",value,kColumn);
971              if (linear) {
972                if (!constant)
973                  strcpy(temp,temp2);
974                else if (value>0.0) 
975                  sprintf(temp,"%g+%s",constant,temp2);
976                else
977                  sprintf(temp,"%g%s",constant,temp2);
978              } else {
979                if (value>0.0) 
980                  sprintf(temp,"%s+%s",expr,temp2);
981                else
982                  sprintf(temp,"%s%s",expr,temp2);
983              }
984              assert (strlen(temp)<1000);
985              setElement(iRow,j,temp);
986              printf("el for row %d column c%7.7d is %s\n",iRow,j,temp);
987            }
988          }
989        } else {
990          // objective
991          for (int j=0;j<n_var;j++) {
992            for (int k=start[j];k<start[j+1];k++) {
993              int kColumn = row[k];
994              double value = element[k];
995              // ampl gives twice with assumed 0.5
996              if (kColumn>j)
997                continue;
998              else if (kColumn==j)
999                value *= 0.5;
1000              const char * expr = getColumnObjectiveAsString(j);
1001              double constant=0.0;
1002              bool linear;
1003              if (expr&&strcmp(expr,"Numeric")) {
1004                linear=false;
1005              } else {
1006                constant = getColumnObjective(j);
1007                linear=true;
1008              }
1009              char temp[1000];
1010              char temp2[30];
1011              if (value==1.0) 
1012                sprintf(temp2,"c%7.7d",kColumn);
1013              else
1014                sprintf(temp2,"%g*c%7.7d",value,kColumn);
1015              if (linear) {
1016                if (!constant)
1017                  strcpy(temp,temp2);
1018                else if (value>0.0) 
1019                  sprintf(temp,"%g+%s",constant,temp2);
1020                else
1021                  sprintf(temp,"%g%s",constant,temp2);
1022              } else {
1023                if (value>0.0) 
1024                  sprintf(temp,"%s+%s",expr,temp2);
1025                else
1026                  sprintf(temp,"%s%s",expr,temp2);
1027              }
1028              assert (strlen(temp)<1000);
1029              setObjective(j,temp);
1030              printf("el for objective column c%7.7d is %s\n",j,temp);
1031            }
1032          }
1033        }
1034      }
1035    }
1036  }
1037  // see if any sos
1038  {
1039    char *sostype;
1040    int nsosnz, *sosbeg, *sosind, * sospri;
1041    double *sosref;
1042    int nsos;
1043    int i = ASL_suf_sos_explict_free;
1044    int copri[2], **p_sospri;
1045    copri[0] = 0;
1046    copri[1] = 0;
1047    p_sospri = &sospri;
1048    nsos = suf_sos(i, &nsosnz, &sostype, p_sospri, copri,
1049                   &sosbeg, &sosind, &sosref);
1050    if (nsos) {
1051      numberSOS_=nsos;
1052      typeSOS_ = new int [numberSOS_];
1053      prioritySOS_ = new int [numberSOS_];
1054      startSOS_ = new int [numberSOS_+1];
1055      memberSOS_ = new int[nsosnz];
1056      referenceSOS_ = new double [nsosnz];
1057      sos_kludge(nsos, sosbeg, sosref,sosind);
1058      for (int i=0;i<nsos;i++) {
1059        int ichar = sostype[i];
1060        assert (ichar=='1'||ichar=='2');
1061        typeSOS_[i]=ichar-'0';
1062      } 
1063      memcpy(prioritySOS_,sospri,nsos*sizeof(int));
1064      memcpy(startSOS_,sosbeg,(nsos+1)*sizeof(int));
1065      memcpy(memberSOS_,sosind,nsosnz*sizeof(int));
1066      memcpy(referenceSOS_,sosref,nsosnz*sizeof(double));
1067    }
1068  }
1069}
Note: See TracBrowser for help on using the repository browser.