Changeset 491 for branches/devel/Cbc/src/Cbc_ampl.cpp
 Timestamp:
 Nov 14, 2006 4:49:44 PM (14 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

branches/devel/Cbc/src/Cbc_ampl.cpp
r468 r491 28 28 # include "unistd.h" 29 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" 30 37 31 38 extern "C" { … … 36 43 #include <string> 37 44 #include <cassert> 38 #include "CoinSort.hpp"39 45 /* so decodePhrase and clpCheck can access */ 40 46 static ampl_info * saveInfo=NULL; … … 572 578 write_sol(buf,info>primalSolution,info>dualSolution,&Oinfo); 573 579 } 580 /* Read a problem from AMPL nl file 581 */ 582 CoinModel::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 631 qterm(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 649 void 650 CoinModel::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 columnwise (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=numberColumnsnumberBinarynumberIntegers; 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 = i1; 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 TracChangeset
for help on using the changeset viewer.