#include /* Bring in the declarations for the string functions */ #include #include /* Include declaration for functions at end of program */ static int populatebyrow (CPXENVptr env, CPXLPptr lp, double **norm, double **tum, int nnorm, int ntum, int nge, int *active, int numactive); static void free_and_null (char **ptr), usage (char *progname); static void fillvec (int node, int *bitvector); static void compute_margin (int l, int r, int *bitvector, double *margin); #define NR 300 #define NC 15000 #define NT 5000000 #define EPS 0.000001 #define N 1000 #define N2 100000 #define LEN 600 #define MAX 1000 #define YSTEP 12 #define XSTEP 35 #define ENDX 500 int enc = 1; char s[N][256]; int NUMROWS; int NUMCOLS; int NUMNZ; int status = 0; double obj[NC]; double lb[NC]; double ub[NC]; char *colname[NC]; int rmatbeg[NR]; int rmatind[NT]; double rmatval[NT]; double rhs[NR]; char sense[NR]; char *rowname[NR]; int *leftchild, *rightchild; double *margin; double **pair; int nnorm; int ntum; int nge; int *ypos; void yplace(int node) { if(node < ntum) { ypos[node] = enc * YSTEP; enc++; } else { yplace(leftchild[node]); yplace(rightchild[node]); ypos[node] = (ypos[leftchild[node]] + ypos[rightchild[node]] )/2; } } void draw(int node, int depth) { printf("newpath\n%d %d moveto\n",depth*XSTEP,ypos[node]); printf("%d %d lineto\nstroke\n",depth*XSTEP,ypos[leftchild[node]]); printf("newpath\n%d %d moveto\n",depth*XSTEP,ypos[node]); printf("%d %d lineto\nstroke\n",depth*XSTEP,ypos[rightchild[node]]); printf("newpath\n%d %d moveto\n",depth*XSTEP+1,ypos[node]-2); printf("(%2d:%1.2f) show\n",node/1000-ntum,margin[node]); printf("newpath\n%d %d moveto\n",depth*XSTEP,ypos[leftchild[node]]); if(leftchild[node] < ntum) { printf("%d %d lineto\nstroke\n",ENDX,ypos[leftchild[node]]); printf("newpath\n%d %d moveto\n",ENDX,ypos[leftchild[node]]); printf("(%d:%s) show\n",leftchild[node],s[leftchild[node]]); } else { printf("%d %d lineto\nstroke\n",(depth+1)*XSTEP,ypos[leftchild[node]]); draw(leftchild[node],depth+1); } printf("newpath\n%d %d moveto\n",depth*XSTEP,ypos[rightchild[node]]); if(rightchild[node] < ntum) { printf("%d %d lineto\nstroke\n",ENDX,ypos[rightchild[node]]); printf("newpath\n%d %d moveto\n",ENDX,ypos[rightchild[node]]); printf("(%d:%s) show\n",rightchild[node],s[rightchild[node]]); } else { printf("%d %d lineto\nstroke\n",(depth+1)*XSTEP,ypos[rightchild[node]]); draw(rightchild[node],depth+1); } } int main (argc, argv) int argc; char **argv; { /* Declare and allocate space for the variables and arrays where we will store the optimization results including the status, objective value, variable values, dual values, row slacks and variable reduced costs. */ int solstat; double objval; double *x = NULL; double *pi = NULL; double *slack = NULL; double *dj = NULL; double **norm; double **tum; double tempfloat; int i,j; CPXENVptr env = NULL; CPXLPptr lp = NULL; int status = 0; int cur_numrows, cur_numcols; int *activetum; int numactive; int *parent; double **weightvector; FILE *fp; int root; int *bitvector, *maxbitvector; double m, max_margin; int left, right; char str; int k = 0; int pos = 0; fp = fopen("tumornames.txt","r"); while((str = fgetc(fp)) != EOF) { if(str == '\n') { s[k++][pos] = '\0'; pos = 0; } else s[k][pos++] = str; } fclose(fp); /* Check the command line arguments */ fprintf(stderr,"Input number of normals, number of tumors and number of genes\n"); scanf("%d%d%d",&nnorm,&ntum,&nge); norm = (double **)calloc(nnorm,sizeof(double *)); for(i = 0; i < nnorm; i++) norm[i] = (double *)calloc(nge,sizeof(double)); tum = (double **)calloc(ntum,sizeof(double *)); for( i = 0; i < ntum; i++) tum[i] = (double *)calloc(nge,sizeof(double)); activetum = (int *)calloc(ntum,sizeof(int)); leftchild = (int *)calloc(2*ntum-1, sizeof(int)); rightchild = (int *)calloc(2*ntum-1, sizeof(int)); margin = (double *) calloc(2*ntum-1, sizeof(double)); parent = (int *) calloc(2*ntum-1, sizeof(int)); ypos = (int *) calloc(2*ntum-1, sizeof(int)); for( i = 0; i < ntum; i++) { leftchild[i] = i; rightchild[i] = i; parent[i] = i; } pair = (double **) calloc(ntum, sizeof(double *)); for( i = 0; i < ntum; i++) pair[i] = (double *) calloc(ntum, sizeof(double)); weightvector = (double **)calloc(2*ntum-1, sizeof(double *)); for( i = 0; i < 2*ntum-1; i++) weightvector[i] = (double *) calloc(nge+2,sizeof(double)); fp = fopen("normal.txt","r"); for( j = 0; j < nge; j++) for( i = 0; i < nnorm; i++) { fscanf(fp,"%lf",&tempfloat); norm[i][j] = tempfloat; } fclose(fp); fp = fopen("tumor.txt","r"); for( j = 0; j < nge; j++) for( i = 0; i < ntum; i++) { fscanf(fp,"%lf",&tempfloat); tum[i][j] = tempfloat; } fclose(fp); fp = fopen("cluster.txt","w"); NUMCOLS = nge+2; for( j = 0; j < nge; j++) for( i = 1; i <= nnorm; i++) { rmatval[i*NUMCOLS+j] = norm[i-1][j]; rmatind[i*NUMCOLS+j] = j; } /* Initialize the CPLEX environment */ env = CPXopenCPLEX (&status); /* If an error occurs, the status value indicates the reason for failure. A call to CPXgeterrorstring will produce the text of the error message. Note that CPXopenCPLEX produces no output, so the only way to see the cause of the error is to use CPXgeterrorstring. For other CPLEX routines, the errors will be seen if the CPX_PARAM_SCRIND indicator is set to CPX_ON. */ if ( env == NULL ) { char errmsg[1024]; fprintf (stderr, "Could not open CPLEX environment.\n"); CPXgeterrorstring (env, status, errmsg); fprintf (stderr, "%s", errmsg); goto TERMINATE; } /* Turn on data checking */ status = CPXsetintparam (env, CPX_PARAM_DATACHECK, CPX_ON); if ( status ) { fprintf (stderr, "Failure to turn on data checking, error %d.\n", status); goto TERMINATE; } for( i = 0; i < ntum; i++) for( j = i; j < ntum; j++) { activetum[0] = i; activetum[1] = j; numactive = 2; /* Create the problem. */ lp = CPXcreateprob (env, &status, "lpex1"); /* A returned pointer of NULL may mean that not enough memory was available or there was some other problem. In the case of failure, an error message will have been written to the error channel from inside CPLEX. In this example, the setting of the parameter CPX_PARAM_SCRIND causes the error message to appear on stdout. */ if ( lp == NULL ) { fprintf (stderr, "Failed to create LP.\n"); goto TERMINATE; } /* Now populate the problem with the data. For building large problems, consider setting the row, column and nonzero growth parameters before performing this task. */ status = populatebyrow (env, lp, norm,tum, nnorm,ntum,nge, activetum,numactive); if ( status ) { fprintf (stderr, "Failed to populate problem.\n"); goto TERMINATE; } /* Optimize the problem and obtain solution. */ status = CPXlpopt (env, lp); if ( status ) { fprintf (stderr, "Failed to optimize LP.\n"); goto TERMINATE; } /* The size of the problem should be obtained by asking CPLEX what the actual size is, rather than using sizes from when the problem was built. cur_numrows and cur_numcols store the current number of rows and columns, respectively. */ cur_numrows = CPXgetnumrows (env, lp); cur_numcols = CPXgetnumcols (env, lp); x = (double *) malloc (cur_numcols * sizeof(double)); slack = (double *) malloc (cur_numrows * sizeof(double)); dj = (double *) malloc (cur_numcols * sizeof(double)); pi = (double *) malloc (cur_numrows * sizeof(double)); if ( x == NULL || slack == NULL || dj == NULL || pi == NULL ) { status = CPXERR_NO_MEMORY; fprintf (stderr, "Could not allocate memory for solution.\n"); goto TERMINATE; } status = CPXsolution (env, lp, &solstat, &objval, x, pi, slack, dj); if ( status ) { fprintf (stderr, "Failed to obtain solution.\n"); goto TERMINATE; } /* Write the output to the screen. */ /* fprintf (stderr,"%d %d %f\n\n", i,j,objval); */ pair[i][j] = pair[j][i] = objval; if(i == j) { margin[i] = objval; memcpy(weightvector[i], x, cur_numcols * sizeof(double)); } free_and_null ((char **) &x); free_and_null ((char **) &slack); free_and_null ((char **) &dj); free_and_null ((char **) &pi); if ( lp != NULL ) { status = CPXfreeprob (env, &lp); if ( status ) { fprintf (stderr, "CPXfreeprob failed, error code %d.\n", status); } } } root = ntum; bitvector = (int *) calloc(ntum, sizeof(int)); maxbitvector = (int *) calloc(ntum, sizeof(int)); while(root < 2 * ntum-1) { max_margin = -100; for( i = 0; i < root; i++) if(parent[i] == i) for( j = i+1; j < root; j++) if(parent[j] == j) { compute_margin(i,j,bitvector,&m); if(m > max_margin) { max_margin = m; left = i; right = j; memcpy(maxbitvector, bitvector, ntum * sizeof(int)); } } leftchild[root] = left; rightchild[root] = right; parent[root] = root; parent[left] = parent[right] = root; numactive = 0; for( i = 0; i < ntum; i++) if(maxbitvector[i] == 1) activetum[numactive++] = i; /* Create the problem. */ lp = CPXcreateprob (env, &status, "lpex1"); /* A returned pointer of NULL may mean that not enough memory was available or there was some other problem. In the case of failure, an error message will have been written to the error channel from inside CPLEX. In this example, the setting of the parameter CPX_PARAM_SCRIND causes the error message to appear on stdout. */ if ( lp == NULL ) { fprintf (stderr, "Failed to create LP.\n"); goto TERMINATE; } /* Now populate the problem with the data. For building large problems, consider setting the row, column and nonzero growth parameters before performing this task. */ status = populatebyrow (env, lp, norm,tum, nnorm,ntum,nge, activetum,numactive); if ( status ) { fprintf (stderr, "Failed to populate problem.\n"); goto TERMINATE; } /* Optimize the problem and obtain solution. */ status = CPXlpopt (env, lp); if ( status ) { fprintf (stderr, "Failed to optimize LP.\n"); goto TERMINATE; } /* The size of the problem should be obtained by asking CPLEX what the actual size is, rather than using sizes from when the problem was built. cur_numrows and cur_numcols store the current number of rows and columns, respectively. */ cur_numrows = CPXgetnumrows (env, lp); cur_numcols = CPXgetnumcols (env, lp); x = (double *) malloc (cur_numcols * sizeof(double)); slack = (double *) malloc (cur_numrows * sizeof(double)); dj = (double *) malloc (cur_numcols * sizeof(double)); pi = (double *) malloc (cur_numrows * sizeof(double)); if ( x == NULL || slack == NULL || dj == NULL || pi == NULL ) { status = CPXERR_NO_MEMORY; fprintf (stderr, "Could not allocate memory for solution.\n"); goto TERMINATE; } status = CPXsolution (env, lp, &solstat, &objval, x, pi, slack, dj); if ( status ) { fprintf (stderr, "Failed to obtain solution.\n"); goto TERMINATE; } /* Write the output to the screen. */ margin[root] = objval; memcpy(weightvector[root], x, cur_numcols * sizeof(double)); fprintf(fp,"====================================\n"); for(i = 0; i < ntum; i++) if(maxbitvector[i] == 1) fprintf(fp,"%d ",i); fprintf(fp,"\n====================================\n"); fprintf(fp,"%2.3f\n=====================================\n",margin[root]); for(j = 0; j < nge; j++) if(x[j] > EPS) fprintf(fp,"%d %f\n",j,x[j]); root++; free_and_null ((char **) &x); free_and_null ((char **) &slack); free_and_null ((char **) &dj); free_and_null ((char **) &pi); if ( lp != NULL ) { status = CPXfreeprob (env, &lp); if ( status ) { fprintf (stderr, "CPXfreeprob failed, error code %d.\n", status); } } } root--; fclose(fp); yplace(root); printf("/Times-Roman findfont\n 9 scalefont\n setfont\n"); draw(root,1); TERMINATE: free_and_null ((char **) &x); free_and_null ((char **) &slack); free_and_null ((char **) &dj); free_and_null ((char **) &pi); /* Free up the problem as allocated by CPXcreateprob, if necessary */ if ( lp != NULL ) { status = CPXfreeprob (env, &lp); if ( status ) { fprintf (stderr, "CPXfreeprob failed, error code %d.\n", status); } } /* Free up the CPLEX environment, if necessary */ if ( env != NULL ) { status = CPXcloseCPLEX (&env); /* Note that CPXcloseCPLEX produces no output, so the only way to see the cause of the error is to use CPXgeterrorstring. For other CPLEX routines, the errors will be seen if the CPX_PARAM_SCRIND indicator is set to CPX_ON. */ if ( status ) { char errmsg[1024]; fprintf (stderr, "Could not close CPLEX environment.\n"); CPXgeterrorstring (env, status, errmsg); fprintf (stderr, "%s", errmsg); } } return (status); } /* END main */ /* This simple routine frees up the pointer *ptr, and sets *ptr to NULL */ static void free_and_null (ptr) char **ptr; { if ( *ptr != NULL ) { free (*ptr); *ptr = NULL; } } /* END free_and_null */ static void fillvec(int node, int *bitvector) { if(leftchild[node] == node) { bitvector[node] = 1; return; } fillvec(leftchild[node], bitvector); fillvec(rightchild[node], bitvector); return; } static void compute_margin(int l, int r, int *bitvector, double *m) { int *b1, *b2; int i, j; b1 = (int *) calloc(ntum, sizeof(int)); b2 = (int *) calloc(ntum, sizeof(int)); for(i = 0; i < ntum; i++) { bitvector[i] = 0; b1[i] = b2[i] = 0; } fillvec(l, b1); fillvec(r, b2); *m = margin[l]; *m = *m < margin[r] ? *m : margin[r]; for(i = 0; i < ntum; i++) if(b1[i] == 1) { bitvector[i] = 1; for(j = 0; j < ntum; j++) if(b2[j] == 1) { *m = *m < pair[i][j] ? *m : pair[i][j]; bitvector[j] = 1; } } free_and_null((char **) &b1); free_and_null((char **) &b2); } static void usage (progname) char *progname; { fprintf (stderr,"Usage: %s -X\n", progname); fprintf (stderr," where X is one of the following options: \n"); fprintf (stderr," r generate problem by row\n"); fprintf (stderr," c generate problem by column\n"); fprintf (stderr," n generate problem by nonzero\n"); fprintf (stderr," Exiting...\n"); } /* END usage */ /* These functions all populate the problem with data for the following linear program: Maximize obj: x1 + 2 x2 + 3 x3 Subject To c1: - x1 + x2 + x3 <= 20 c2: x1 - 3 x2 + x3 <= 30 Bounds 0 <= x1 <= 40 End */ /* To populate by row, we first create the columns, and then add the rows. */ static int populatebyrow (env, lp, norm, tum, nnorm, ntum, nge, active, numactive) CPXENVptr env; CPXLPptr lp; double **norm; double **tum; int nnorm; int ntum; int nge; int *active; int numactive; { int i,j; NUMROWS = nnorm + 1 + numactive; NUMNZ = NUMROWS * NUMCOLS; for(j = 0; j < nge; j++) for(i = 0; i < numactive; i++) { rmatval[(i+nnorm+1)*NUMCOLS+j] = tum[active[i]][j]; rmatind[(i+nnorm+1)*NUMCOLS+j] = j; } for(i = 0; i < NUMROWS; i++) { rhs[i] = 0.0; rowname[i] = ""; rmatbeg[i] = i*NUMCOLS; rmatind[i*NUMCOLS+NUMCOLS-2] = NUMCOLS-2; rmatind[i*NUMCOLS+NUMCOLS-1] = NUMCOLS-1; if(i > 0 && i <= nnorm) { sense[i] = 'L'; rmatval[i*NUMCOLS+NUMCOLS-2] = 0.0; rmatval[i*NUMCOLS+NUMCOLS-1] = -1.0; } else if(i > nnorm) { sense[i] = 'G'; rmatval[i*NUMCOLS+NUMCOLS-2] = -1.0; rmatval[i*NUMCOLS+NUMCOLS-1] = 0.0; } } for(j = 0; j < NUMCOLS; j++) { colname[j] = ""; rmatval[j] = 1.0; rmatind[j] = j; } rmatval[NUMCOLS-2] = 0.0; rmatval[NUMCOLS-1] = 0.0; sense[0] = 'E'; rhs[0] = 1.0; obj[NUMCOLS-2] = 1.0; obj[NUMCOLS-1] = -1.0; for(j = 0; j < NUMCOLS-2; j++) { obj[j] = 0.0; lb[j] = 0.0; ub[j] = 1.0; } lb[NUMCOLS-2] = -100.0; lb[NUMCOLS-1] = -100.0; ub[NUMCOLS-2] = 100.0; ub[NUMCOLS-1] = 100.0; CPXchgobjsen (env, lp, CPX_MAX); /* Problem is maximization */ status = CPXnewcols (env, lp, NUMCOLS, obj, lb, ub, NULL, colname); if ( status ) goto TERMINATE; status = CPXaddrows (env, lp, 0, NUMROWS, NUMNZ, rhs, sense, rmatbeg, rmatind, rmatval, NULL, rowname); if ( status ) goto TERMINATE; TERMINATE: return (status); } /* END populatebyrow */