From: adam Date: Sun, 4 Apr 2004 03:22:57 +0000 (+0000) Subject: new linear constraint solver for layout X-Git-Url: http://git.megacz.com/?p=org.ibex.core.git;a=commitdiff_plain;h=678a052b8a88a204220b57071360775ffa7b8c2b new linear constraint solver for layout darcs-hash:20040404032257-5007d-bca57114e1842dcdf4e79b01ce67cf04f5208351.gz --- diff --git a/src/org/ibex/Box.java b/src/org/ibex/Box.java index 524cb7c..3892ce9 100644 --- a/src/org/ibex/Box.java +++ b/src/org/ibex/Box.java @@ -214,11 +214,6 @@ public final class Box extends JSScope implements Scheduler.Task { // static stuff so we don't have to keep reallocating private static int[] numRowsInCol = new int[65535]; - private static LENGTH[] colWidth = new LENGTH[65535]; - private static LENGTH[] colMaxWidth = new LENGTH[65535]; - private static LENGTH[] rowHeight = new LENGTH[65535]; - private static LENGTH[] rowMaxHeight = new LENGTH[65535]; - static { for(int i=0; icoeff.length) coeff = new float[nc+1]; + LinearProgramming.Simplex lp_h = new LinearProgramming.Simplex(); + LinearProgramming.Problem lpr_h = new LinearProgramming.Problem(nc, nc); + LinearProgramming.Simplex lp_v = new LinearProgramming.Simplex(); + LinearProgramming.Problem lpr_v = new LinearProgramming.Problem(nc, nc); + + //#repeat col/row colspan/rowspan contentwidth/contentheight width/height colMaxWidth/rowMaxHeight colWidth/rowHeight \ + // HSHRINK/VSHRINK maxwidth/maxheight cols/rows minwidth/minheight colWidth/rowHeight x_slack/y_slack lp_h/lp_v lpr_h/lpr_v + do { + // objective function + coeff[cols*2+numkids] = coeff[cols*2+numkids+1] = (float)-10000.0; // attempt to make sum of columns equal to parent width + for(int i=cols*2; i=child.col && i=child.col && i 0 && cols > 0 && startslack != x_slack;) { - int increment = max(1, x_slack / cols); - startslack = x_slack; - for(short col=0; col < cols; col++) { - // FIXME: double check this - int diff = min(min(colMaxWidth[col], colWidth[col] + increment) - colWidth[col], x_slack); - x_slack -= diff; - colWidth[col] += diff; + + // third priority: try to make columns similar size + for(int i=0 ; i unbounded) child_x -= (child_width - unbounded) / 2; + // colWidth/rowHeight child_width/child_height ALIGN_RIGHT/ALIGN_BOTTOM ALIGN_LEFT/ALIGN_TOP lpr_h/lpr_v + child_width = 0; + child_x = 0; + for(int i=0; i < child.col; i++) child_x += Math.round(lpr_h.solution[lpr_h.rows+i+1]); + for(int i = child.col; i < child.col + child.colspan; i++) child_width += Math.round(lpr_h.solution[lpr_h.rows+i+1]); + child_x += (child_width - min(child_width, child.test(HSHRINK) ? child.contentwidth : child.maxwidth)) / 2; + child_width = min(child_width, child.test(HSHRINK) ? child.contentwidth : child.maxwidth); //#end } child.resize(child_x, child_y, child_width, child_height); } - // cleanup - for(int i=0; i 0) + child.place_children(); } @@ -685,7 +702,6 @@ public final class Box extends JSScope implements Scheduler.Task { fillcolor = newfillcolor; } else if(value instanceof JS) { texture = Picture.load((JS)value, this); - if (texture != null) perform(); } else { throw new JSExn("fill must be null, a String, or a stream, not a " + value.getClass()); } diff --git a/src/org/ibex/Surface.java b/src/org/ibex/Surface.java index 394f299..b553056 100644 --- a/src/org/ibex/Surface.java +++ b/src/org/ibex/Surface.java @@ -248,7 +248,7 @@ public abstract class Surface extends PixelBuffer implements Scheduler.Task { dirty(0, root.maxheight - Main.scarImage.height, Main.scarImage.width, Main.scarImage.height); } root.resize(root.x, root.y, root.maxwidth, root.maxheight); - root.resize_children(); + root.place_children(); setSize(root.width, root.height); String oldcursor = cursor; cursor = "default"; diff --git a/src/org/ibex/util/LinearProgramming.java b/src/org/ibex/util/LinearProgramming.java new file mode 100644 index 0000000..869737f --- /dev/null +++ b/src/org/ibex/util/LinearProgramming.java @@ -0,0 +1,3128 @@ +package org.ibex.util; +import java.io.* ; +import java.util.* ; + +public class LinearProgramming { + + public final static short FAIL = -1; + + public final static short NULL = 0; + public final static short FALSE = 0; + public final static short TRUE = 1; + + public final static short DEFNUMINV = 50; + + /* solve status values */ + public final static short OPTIMAL = 0; + public final static short MILP_FAIL = 1; + public final static short INFEASIBLE = 2; + public final static short UNBOUNDED = 3; + public final static short FAILURE = 4; + public final static short RUNNING = 5; + + /* lag_solve extra status values */ + public final static short FEAS_FOUND = 6; + public final static short NO_FEAS_FOUND = 7; + public final static short BREAK_BB = 8; + + public final static short FIRST_NI = 0; + public final static short RAND_NI = 1; + + public final static short LE = 0; + public final static short EQ = 1; + public final static short GE = 2; + public final static short OF = 3; + + public final static short MAX_WARN_COUNT = 20; + + public final static float DEF_INFINITE = (float)1e24; /* limit for dynamic range */ + public final static float DEF_EPSB = (float)5.01e-7; /* for rounding RHS values to 0 determine + infeasibility basis */ + public final static float DEF_EPSEL = (float)1e-8; /* for rounding other values (vectors) to 0 */ + public final static float DEF_EPSD = (float)1e-6; /* for rounding reduced costs to zero */ + public final static float DEF_EPSILON = (float)1e-3; /* to determine if a float value is integer */ + + public final static float PREJ = (float)1e-3; /* pivot reject (try others first) */ + + public final static int HASHSIZE = 10007; /* prime number is better, MB */ + public final static int ETA_START_SIZE = 10000; /* start size of array Eta. Realloced if needed */ + public final static String STD_ROW_NAME_PREFIX = "r_"; + + static class Ref { + float value; + public Ref(float v) { value = v; } + } + + public static class Simplex { + /* Globals */ + Problem Lp; /* pointer to active problem */ + int Rows; + int columns; + int Sum; + int Non_zeros; + int Level; + Matrix[] Mat; + int[] Col_no; + int[] Col_end; + int[] Row_end; + float[] Orig_rh; + float[] Rh; + float[] Rhs; + float[] Orig_upbo; + float[] Orig_lowbo; + float[] Upbo; + float[] Lowbo; + int[] Bas; + short[] Basis; + short[] Lower; + int Eta_alloc; + int Eta_size; + float[] Eta_value; + int[] Eta_row_nr; + int[] Eta_col_end; + int Num_inv; + float[] Solution; + float[] Best_solution; + float Infinite; + float Epsilon; + float Epsb; + float Epsd; + float Epsel; + + float TREJ; + float TINV; + + short Maximise; + short Floor_first; + float Extrad; + + int Warn_count; /* used in CHECK version of rounding macro */ + + void inc_mat_space(Problem lp, int maxextra) + { + if(lp.non_zeros + maxextra > lp.mat_alloc) + { + int i; + lp.mat_alloc = lp.non_zeros + maxextra; + Matrix[] old_mat = lp.mat; + lp.mat = new Matrix[lp.mat_alloc]; + for (i = old_mat.length; i < lp.mat.length; i++) + lp.mat[i] = new Matrix(0,0); + System.arraycopy(old_mat, 0, lp.mat, 0, old_mat.length); + int[] old_col_no = lp.col_no; + lp.col_no = new int[lp.mat_alloc]; + System.arraycopy(old_col_no, 0, lp.col_no, 0, old_col_no.length); + if (lp.active != FALSE) + { + Mat=lp.mat; + Col_no=lp.col_no; + } + } + } // end of inc_mat_space + + void inc_row_space(Problem lp) + { + if(lp.rows > lp.rows_alloc) + { + lp.rows_alloc=lp.rows+10; + lp.sum_alloc=lp.rows_alloc+lp.columns_alloc; + + float[] db_ptr = lp.orig_rh; + lp.orig_rh = new float[lp.rows_alloc + 1]; + System.arraycopy(db_ptr, 0, lp.orig_rh, 0, db_ptr.length); + + db_ptr = lp.rh; + lp.rh = new float[lp.rows_alloc + 1]; + System.arraycopy(db_ptr, 0, lp.rh, 0, db_ptr.length); + + db_ptr = lp.rhs; + lp.rhs = new float[lp.rows_alloc + 1]; + System.arraycopy(db_ptr, 0, lp.rhs, 0, db_ptr.length); + + db_ptr = lp.orig_upbo; + lp.orig_upbo = new float[lp.sum_alloc + 1]; + System.arraycopy(db_ptr, 0, lp.orig_upbo, 0, db_ptr.length); + + db_ptr = lp.upbo; + lp.upbo = new float[lp.sum_alloc + 1]; + System.arraycopy(db_ptr, 0, lp.upbo, 0, db_ptr.length); + + db_ptr = lp.orig_lowbo; + lp.orig_lowbo = new float[lp.sum_alloc + 1]; + System.arraycopy(db_ptr, 0, lp.orig_lowbo, 0, db_ptr.length); + + db_ptr = lp.lowbo; + lp.lowbo = new float[lp.sum_alloc + 1]; + System.arraycopy(db_ptr, 0, lp.lowbo, 0, db_ptr.length); + + db_ptr = lp.solution; + lp.solution = new float[lp.sum_alloc + 1]; + System.arraycopy(db_ptr, 0, lp.solution, 0, db_ptr.length); + + db_ptr = lp.best_solution; + lp.best_solution = new float[lp.sum_alloc + 1]; + System.arraycopy(db_ptr, 0, lp.best_solution, 0, db_ptr.length); + + int[] int_ptr = lp.row_end; + lp.row_end = new int[lp.rows_alloc + 1]; + System.arraycopy(int_ptr, 0, lp.row_end, 0, int_ptr.length); + + short[] short_ptr = lp.basis; + lp.basis = new short[lp.sum_alloc + 1]; + System.arraycopy(short_ptr, 0, lp.basis, 0, short_ptr.length); + + short_ptr = lp.lower; + lp.lower = new short[lp.sum_alloc + 1]; + System.arraycopy(short_ptr, 0, lp.lower, 0, short_ptr.length); + + int_ptr = lp.bas; + lp.bas = new int[lp.rows_alloc + 1]; + System.arraycopy(int_ptr, 0, lp.bas, 0, int_ptr.length); + + db_ptr = lp.duals; + lp.duals = new float[lp.rows_alloc + 1]; + System.arraycopy(db_ptr, 0, lp.duals, 0, db_ptr.length); + + short_ptr = lp.ch_sign; + lp.ch_sign = new short[lp.rows_alloc + 1]; + System.arraycopy(short_ptr, 0, lp.ch_sign, 0, short_ptr.length); + + int_ptr = lp.eta_col_end; + lp.eta_col_end = new int[lp.rows_alloc + lp.max_num_inv]; + System.arraycopy(int_ptr, 0, lp.eta_col_end, 0, int_ptr.length); + + if(lp.names_used != FALSE) { + String[] str_ptr = lp.row_name; + lp.row_name = new String[lp.rows_alloc + 1]; + System.arraycopy(str_ptr, 0, lp.row_name, 0, str_ptr.length); + } + + if(lp.scaling_used != FALSE) { + db_ptr = lp.scale; + lp.scale = new float[lp.sum_alloc + 1]; + System.arraycopy(db_ptr, 0, lp.scale, 0, db_ptr.length); + } + + if(lp.active != FALSE) + set_globals(lp); + } + } + + void inc_col_space(Problem lp) + { + if(lp.columns >= lp.columns_alloc) + { + int[] int_ptr; + lp.columns_alloc=lp.columns+10; + lp.sum_alloc=lp.rows_alloc+lp.columns_alloc; + + float[] float_ptr = lp.orig_upbo; + lp.orig_upbo = new float[lp.sum_alloc + 1]; + System.arraycopy(float_ptr, 0, lp.orig_upbo, 0, float_ptr.length); + + float_ptr = lp.upbo; + lp.upbo = new float[lp.sum_alloc + 1]; + System.arraycopy(float_ptr, 0, lp.upbo, 0, float_ptr.length); + + float_ptr = lp.orig_lowbo; + lp.orig_lowbo = new float[lp.sum_alloc + 1]; + System.arraycopy(float_ptr, 0, lp.orig_lowbo, 0, float_ptr.length); + + float_ptr = lp.lowbo; + lp.lowbo = new float[lp.sum_alloc + 1]; + System.arraycopy(float_ptr, 0, lp.lowbo, 0, float_ptr.length); + + float_ptr = lp.solution; + lp.solution = new float[lp.sum_alloc + 1]; + System.arraycopy(float_ptr, 0, lp.solution, 0, float_ptr.length); + + float_ptr = lp.best_solution; + lp.best_solution = new float[lp.sum_alloc + 1]; + System.arraycopy(float_ptr, 0, lp.best_solution, 0, float_ptr.length); + + short[] short_ptr = lp.basis; + lp.basis = new short[lp.sum_alloc + 1]; + System.arraycopy(short_ptr, 0, lp.basis, 0, short_ptr.length); + + short_ptr = lp.lower; + lp.lower = new short[lp.sum_alloc + 1]; + System.arraycopy(short_ptr, 0, lp.lower, 0, short_ptr.length); + + if(lp.names_used != FALSE) { + int[] str_ptr = lp.col_name; + lp.col_name = new int[lp.columns_alloc + 1]; + System.arraycopy(str_ptr, 0, lp.col_name, 0, str_ptr.length); + } + + if(lp.scaling_used != FALSE) { + float_ptr = lp.scale; + lp.scale = new float[lp.sum_alloc + 1]; + System.arraycopy(float_ptr, 0, lp.scale, 0, float_ptr.length); + } + + int_ptr = lp.col_end; + lp.col_end = new int[lp.columns_alloc + 1]; + System.arraycopy(int_ptr, 0, lp.col_end, 0, int_ptr.length); + + if(lp.active != FALSE) + set_globals(lp); + } + } // end of inc_col_space + + public void set_mat(Problem lp, int Row, int column, float Value) + { + int elmnr, lastelm, i; + // System.err.println("lp.mat.length = " + lp.mat.length); + + if(Row > lp.rows || Row < 0) + System.err.print("Row out of range"); + if(column > lp.columns || column < 1) + System.err.print("column out of range"); + if(lp.scaling_used != FALSE) + Value *= lp.scale[Row] * lp.scale[lp.rows + column]; + + if(true /*abs(Value) > lp.epsilon*/) + { + if (lp.basis[column] == TRUE && Row > 0) + lp.basis_valid = FALSE; + lp.eta_valid = FALSE; + elmnr = lp.col_end[column-1]; + while((elmnr < lp.col_end[column]) ? + (lp.mat[elmnr].row_nr != Row) : false) + elmnr++; + + if((elmnr != lp.col_end[column]) ? + (lp.mat[elmnr].row_nr == Row) : false ) + if (lp.scaling_used != FALSE) + { + if (lp.ch_sign[Row] != FALSE) + lp.mat[elmnr].value = + -Value * lp.scale[Row] * lp.scale[column]; + else + lp.mat[elmnr].value = + Value * lp.scale[Row] * lp.scale[column]; + } + else + { + if (lp.ch_sign[Row] != FALSE) + lp.mat[elmnr].value = -Value; + else + lp.mat[elmnr].value = Value; + } + else + { + /* check if more space is needed for matrix */ + inc_mat_space(lp,1); + + /* Shift the matrix */ + lastelm=lp.non_zeros; + for(i = lastelm; i > elmnr ; i--) + lp.mat[i]=lp.mat[i-1]; + for(i = column; i <= lp.columns; i++) + lp.col_end[i]++; + + /* Set new element */ + lp.mat[elmnr].row_nr=Row; + + if (lp.scaling_used != FALSE) + { + if (lp.ch_sign[Row] != FALSE) + lp.mat[elmnr].value=-Value*lp.scale[Row]*lp.scale[column]; + else + lp.mat[elmnr].value=Value*lp.scale[Row]*lp.scale[column]; + } + else + { + if (lp.ch_sign[Row] != FALSE) + lp.mat[elmnr].value=-Value; + else + lp.mat[elmnr].value=Value; + } + + lp.row_end_valid=FALSE; + + lp.non_zeros++; + if (lp.active != FALSE) + Non_zeros=lp.non_zeros; + } + } + } + + public void set_obj_fn(Problem lp, float[] row) + { + for(int i=row.length-1; i>0; i--) row[i] = row[i-1]; + row[0] = (float)0.0; + + int i; + for(i = 1; i <= lp.columns; i++) + set_mat(lp, 0, i, row[i]); + } + + + public void add_constraint(Problem lp, float[] row, short constr_type, float rh) + { + for(int i=row.length-1; i>0; i--) row[i] = row[i-1]; + row[0] = (float)0.0; + + Matrix[] newmat; + int i, j; + int elmnr; + int stcol; + short[] addtoo; + + addtoo = new short[lp.columns + 1]; + + for(i = 1; i <= lp.columns; i++) + if(row[i]!=0) + { + addtoo[i]=TRUE; + lp.non_zeros++; + } + else + addtoo[i]=FALSE; + + newmat = new Matrix[lp.non_zeros]; + for (i = 0; i < newmat.length; i++) + newmat[i] = new Matrix(0, (float)0); + + inc_mat_space(lp, 0); + lp.rows++; + lp.sum++; + inc_row_space(lp); + + if (lp.scaling_used != FALSE) + { + /* shift scale */ + for(i=lp.sum; i > lp.rows; i--) + lp.scale[i]=lp.scale[i-1]; + lp.scale[lp.rows]=1; + } + + if (lp.names_used != FALSE) + lp.row_name[lp.rows] = new String("r_" + lp.rows); + + if(lp.scaling_used != FALSE && lp.columns_scaled != FALSE) + for(i = 1; i <= lp.columns; i++) + row[i] *= lp.scale[lp.rows+i]; + + if(constr_type==GE) + lp.ch_sign[lp.rows] = TRUE; + else + lp.ch_sign[lp.rows] = FALSE; + + elmnr = 0; + stcol = 0; + for(i = 1; i <= lp.columns; i++) + { + for(j = stcol; j < lp.col_end[i]; j++) + { + newmat[elmnr].row_nr=lp.mat[j].row_nr; + newmat[elmnr].value=lp.mat[j].value; + elmnr++; + } + if(addtoo[i] != FALSE) + { + if(lp.ch_sign[lp.rows] != FALSE) + newmat[elmnr].value = -row[i]; + else + newmat[elmnr].value = row[i]; + newmat[elmnr].row_nr = lp.rows; + elmnr++; + } + stcol=lp.col_end[i]; + lp.col_end[i]=elmnr; + } + + lp.mat = newmat; + + for(i=lp.sum ; i > lp.rows; i--) + { + lp.orig_upbo[i]=lp.orig_upbo[i-1]; + lp.orig_lowbo[i]=lp.orig_lowbo[i-1]; + lp.basis[i]=lp.basis[i-1]; + lp.lower[i]=lp.lower[i-1]; + } + + for(i= 1 ; i <= lp.rows; i++) + if(lp.bas[i] >= lp.rows) + lp.bas[i]++; + + if(constr_type==LE || constr_type==GE) + { + lp.orig_upbo[lp.rows]=lp.infinite; + } + else if(constr_type==EQ) + { + lp.orig_upbo[lp.rows]=0; + } + else + { + System.err.print("Wrong constraint type\n"); + System.exit(FAIL); + } + + lp.orig_lowbo[lp.rows]=0; + + if(constr_type==GE && rh != 0) + lp.orig_rh[lp.rows]=-rh; + else + lp.orig_rh[lp.rows]=rh; + + lp.row_end_valid=FALSE; + + lp.bas[lp.rows]=lp.rows; + lp.basis[lp.rows]=TRUE; + lp.lower[lp.rows]=TRUE; + + if (lp.active != FALSE) + set_globals(lp); + lp.eta_valid=FALSE; + } + + public void del_constraint(Problem lp, int del_row) + { + int i, j; + int elmnr; + int startcol; + + if(del_row < 1 || del_row > lp.rows) + { + System.err.println("There is no constraint nr. " + del_row); + System.exit(FAIL); + } + + elmnr=0; + startcol=0; + + for(i = 1; i <= lp.columns; i++) + { + for(j=startcol; j < lp.col_end[i]; j++) + { + if(lp.mat[j].row_nr!=del_row) + { + lp.mat[elmnr]=lp.mat[j]; + if(lp.mat[elmnr].row_nr > del_row) + lp.mat[elmnr].row_nr--; + elmnr++; + } + else + lp.non_zeros--; + } + startcol=lp.col_end[i]; + lp.col_end[i]=elmnr; + } + for(i = del_row; i < lp.rows; i++) + { + lp.orig_rh[i] = lp.orig_rh[i+1]; + lp.ch_sign[i] = lp.ch_sign[i+1]; + lp.bas[i] = lp.bas[i+1]; + if (lp.names_used != FALSE) + lp.row_name[i] = lp.row_name[i+1]; + } + for(i = 1; i < lp.rows; i++) + if(lp.bas[i] > del_row) + lp.bas[i]--; + + for(i=del_row; i < lp.sum; i++) + { + lp.lower[i]=lp.lower[i+1]; + lp.basis[i]=lp.basis[i+1]; + lp.orig_upbo[i]=lp.orig_upbo[i+1]; + lp.orig_lowbo[i]=lp.orig_lowbo[i+1]; + if (lp.scaling_used != FALSE) + lp.scale[i]=lp.scale[i+1]; + } + + lp.rows--; + lp.sum--; + + lp.row_end_valid=FALSE; + + if (lp.active != FALSE) + set_globals(lp); + lp.eta_valid=FALSE; + lp.basis_valid=FALSE; + } + + + public void add_column(Problem lp, float[] column) + { + int i, elmnr; + + lp.columns++; + lp.sum++; + inc_col_space(lp); + inc_mat_space(lp, lp.rows+1); + + if (lp.scaling_used != FALSE) + { + for(i = 0; i <= lp.rows; i++) + column[i]*=lp.scale[i]; + lp.scale[lp.sum]=1; + } + + elmnr=lp.col_end[lp.columns-1]; + for(i = 0 ; i <= lp.rows ; i++) + if(column[i] != 0) + { + lp.mat[elmnr].row_nr=i; + if(lp.ch_sign[i] != FALSE) + lp.mat[elmnr].value=-column[i]; + else + lp.mat[elmnr].value=column[i]; + lp.non_zeros++; + elmnr++; + } + lp.col_end[lp.columns]=elmnr; + lp.orig_lowbo[lp.sum]=0; + lp.orig_upbo[lp.sum]=lp.infinite; + lp.lower[lp.sum]=TRUE; + lp.basis[lp.sum]=FALSE; + if (lp.names_used != FALSE) + lp.col_name[lp.columns] = 0; + + + lp.row_end_valid=FALSE; + + if (lp.active != FALSE) + { + Sum=lp.sum; + columns=lp.columns; + Non_zeros=lp.non_zeros; + } + } + + public void del_column(Problem lp, int column) + { + int i, j, from_elm, to_elm, elm_in_col; + if(column > lp.columns || column < 1) + System.err.print("column out of range in del_column"); + for(i = 1; i <= lp.rows; i++) + { + if(lp.bas[i]==lp.rows+column) + lp.basis_valid=FALSE; + else if(lp.bas[i] > lp.rows+column) + lp.bas[i]--; + } + for(i = lp.rows+column; i < lp.sum; i++) + { + if (lp.names_used != FALSE) + lp.col_name[i-lp.rows] = lp.col_name[i-lp.rows+1]; + lp.orig_upbo[i]=lp.orig_upbo[i+1]; + lp.orig_lowbo[i]=lp.orig_lowbo[i+1]; + lp.upbo[i]=lp.upbo[i+1]; + lp.lowbo[i]=lp.lowbo[i+1]; + lp.basis[i]=lp.basis[i+1]; + lp.lower[i]=lp.lower[i+1]; + if (lp.scaling_used != FALSE) + lp.scale[i]=lp.scale[i+1]; + } + for(i = 0; i < lp.nr_lagrange; i++) + for(j = column; j <= lp.columns; j++) + lp.lag_row[i][j]=lp.lag_row[i][j+1]; + to_elm=lp.col_end[column-1]; + from_elm=lp.col_end[column]; + elm_in_col=from_elm-to_elm; + for(i = from_elm; i < lp.non_zeros; i++) + { + lp.mat[to_elm]=lp.mat[i]; + to_elm++; + } + for(i = column; i < lp.columns; i++) + lp.col_end[i]=lp.col_end[i+1]-elm_in_col; + lp.non_zeros -= elm_in_col; + lp.row_end_valid=FALSE; + lp.eta_valid=FALSE; + + lp.sum--; + lp.columns--; + if (lp.active != FALSE) + set_globals(lp); + } + + public void bound_sum(Problem lp, int column1, int column2, float bound, short type, float[] scratch) { + for(int i=0; i lp.columns || column < 1) + System.err.print("column out of range"); + if (lp.scaling_used != FALSE) + value /= lp.scale[lp.rows + column]; + if(value < lp.orig_lowbo[lp.rows + column]) + System.err.print("UpperBound must be >= lowerBound"); + lp.eta_valid = FALSE; + lp.orig_upbo[lp.rows+column] = value; + } + + public void set_lowbo(Problem lp, int column, float value) + { + if(column > lp.columns || column < 1) + System.err.print("column out of range"); + if (lp.scaling_used != FALSE) + value /= lp.scale[lp.rows + column]; + if(value > lp.orig_upbo[lp.rows + column]) + System.err.print("UpperBound must be >= lowerBound"); + lp.eta_valid = FALSE; + lp.orig_lowbo[lp.rows+column] = value; + } + + public void set_rh(Problem lp, int row, float value) + { + if(row > lp.rows || row < 0) + System.err.print("Row out of Range"); + if(row == 0) /* setting of RHS of OF not meaningful */ + { + System.err.println("Warning: attempt to set RHS of objective function, ignored"); + return; + } + if (lp.scaling_used != FALSE) + if (lp.ch_sign[row] != FALSE) + lp.orig_rh[row] = -value * lp.scale[row]; + else + lp.orig_rh[row] = value * lp.scale[row]; + else + if (lp.ch_sign[row] != FALSE) + lp.orig_rh[row] = -value; + else + lp.orig_rh[row] = value; + lp.eta_valid = FALSE; + } + + public void set_rh_vec(Problem lp, float[] rh) + { + int i; + if (lp.scaling_used != FALSE) + for(i = 1; i <= lp.rows; i++) + if(lp.ch_sign[i] != FALSE) + lp.orig_rh[i]=-rh[i]*lp.scale[i]; + else + lp.orig_rh[i]=rh[i]*lp.scale[i]; + else + for(i=1; i <= lp.rows; i++) + if(lp.ch_sign[i] != FALSE) + lp.orig_rh[i]=-rh[i]; + else + lp.orig_rh[i]=rh[i]; + lp.eta_valid=FALSE; + } + + public void set_maxim(Problem lp) + { + int i; + if(lp.maximise==FALSE) + { + for(i = 0; i < lp.non_zeros; i++) + if(lp.mat[i].row_nr==0) + lp.mat[i].value*=-1; + lp.eta_valid=FALSE; + } + lp.maximise=TRUE; + lp.ch_sign[0]=TRUE; + if (lp.active != FALSE) + Maximise=TRUE; + } + + public void set_minim(Problem lp) + { + int i; + if(lp.maximise==TRUE) + { + for(i = 0; i < lp.non_zeros; i++) + if(lp.mat[i].row_nr==0) + lp.mat[i].value = -lp.mat[i].value; + lp.eta_valid=FALSE; + } + lp.maximise=FALSE; + lp.ch_sign[0]=FALSE; + if (lp.active != FALSE) + Maximise=FALSE; + } + + public void set_constr_type(Problem lp, int row, short con_type) + { + int i; + if(row > lp.rows || row < 1) + System.err.print("Row out of Range"); + if(con_type==EQ) + { + lp.orig_upbo[row]=0; + lp.basis_valid=FALSE; + if (lp.ch_sign[row] != FALSE) + { + for(i = 0; i < lp.non_zeros; i++) + if(lp.mat[i].row_nr==row) + lp.mat[i].value*=-1; + lp.eta_valid=FALSE; + lp.ch_sign[row]=FALSE; + if(lp.orig_rh[row]!=0) + lp.orig_rh[row]*=-1; + } + } + else if(con_type==LE) + { + lp.orig_upbo[row]=lp.infinite; + lp.basis_valid=FALSE; + if (lp.ch_sign[row] != FALSE) + { + for(i = 0; i < lp.non_zeros; i++) + if(lp.mat[i].row_nr==row) + lp.mat[i].value*=-1; + lp.eta_valid=FALSE; + lp.ch_sign[row]=FALSE; + if(lp.orig_rh[row]!=0) + lp.orig_rh[row]*=-1; + } + } + else if(con_type==GE) + { + lp.orig_upbo[row]=lp.infinite; + lp.basis_valid=FALSE; + if(lp.ch_sign[row] == FALSE) + { + for(i = 0; i < lp.non_zeros; i++) + if(lp.mat[i].row_nr==row) + lp.mat[i].value*=-1; + lp.eta_valid=FALSE; + lp.ch_sign[row]=TRUE; + if(lp.orig_rh[row]!=0) + lp.orig_rh[row]*=-1; + } + } + else + System.err.print("Constraint type not (yet) implemented"); + } + + public float mat_elm(Problem lp, int row, int column) + { + float value; + int elmnr; + if(row < 0 || row > lp.rows) + System.err.print("Row out of range in mat_elm"); + if(column < 1 || column > lp.columns) + System.err.print("column out of range in mat_elm"); + value=0; + elmnr=lp.col_end[column-1]; + while(lp.mat[elmnr].row_nr != row && elmnr < lp.col_end[column]) + elmnr++; + if(elmnr != lp.col_end[column]) + { + value = lp.mat[elmnr].value; + if (lp.ch_sign[row] != FALSE) + value = -value; + if (lp.scaling_used != FALSE) + value /= lp.scale[row] * lp.scale[lp.rows + column]; + } + return(value); + } + + + public void get_row(Problem lp, int row_nr, float[] row) + { + int i, j; + + if(row_nr < 0 || row_nr > lp.rows) + System.err.print("Row nr. out of range in get_row"); + for(i = 1; i <= lp.columns; i++) + { + row[i]=0; + for(j=lp.col_end[i-1]; j < lp.col_end[i]; j++) + if(lp.mat[j].row_nr==row_nr) + row[i]=lp.mat[j].value; + if (lp.scaling_used != FALSE) + row[i] /= lp.scale[lp.rows+i] * lp.scale[row_nr]; + } + if(lp.ch_sign[row_nr] != FALSE) + for(i=0; i <= lp.columns; i++) + if(row[i]!=0) + row[i] = -row[i]; + } + + public void get_column(Problem lp, int col_nr, float[] column) + { + int i; + + if(col_nr < 1 || col_nr > lp.columns) + System.err.print("Col. nr. out of range in get_column"); + for(i = 0; i <= lp.rows; i++) + column[i]=0; + for(i = lp.col_end[col_nr-1]; i < lp.col_end[col_nr]; i++) + column[lp.mat[i].row_nr]=lp.mat[i].value; + for(i = 0; i <= lp.rows; i++) + if(column[i] !=0) + { + if(lp.ch_sign[i] != FALSE) + column[i]*=-1; + if (lp.scaling_used != FALSE) + column[i]/=(lp.scale[i] * lp.scale[lp.rows+col_nr]); + } + } + + public void get_reduced_costs(Problem lp, float[] rc) + { + int varnr, i, j; + float f; + + if(lp.basis_valid == FALSE) + System.err.print("Not a valid basis in get_reduced_costs"); + set_globals(lp); + if(lp.eta_valid == FALSE) + invert(); + for(i = 1; i <= lp.sum; i++) + rc[i] = 0; + rc[0] = 1; + btran(rc); + for(i = 1; i <= lp.columns; i++) + { + varnr = lp.rows + i; + if(Basis[varnr] == FALSE) + if(Upbo[varnr] > 0) + { + f = 0; + for(j = Col_end[i - 1]; j < Col_end[i]; j++) + f += rc[Mat[j].row_nr] * Mat[j].value; + rc[varnr] = f; + } + } + for(i = 1; i <= Sum; i++) + rc[i] = round(rc[i], Epsd); + } + + short is_feasible(Problem lp, float[] values) + { + int i, elmnr; + float[] this_rhs; + float dist; + + if (lp.scaling_used != FALSE) + { + for(i = lp.rows+1; i <= lp.sum; i++) + if( values[i - lp.rows] < lp.orig_lowbo[i]*lp.scale[i] + || values[i - lp.rows] > lp.orig_upbo[i]*lp.scale[i]) + return(FALSE); + } + else + { + for(i = lp.rows+1; i <= lp.sum; i++) + if( values[i - lp.rows] < lp.orig_lowbo[i] + || values[i - lp.rows] > lp.orig_upbo[i]) + return(FALSE); + } + this_rhs = new float[lp.rows+1]; + for (i = 0; i <= lp.rows; i++) + this_rhs[i] = 0; + + for(i = 1; i <= lp.columns; i++) + for(elmnr = lp.col_end[i - 1]; elmnr < lp.col_end[i]; elmnr++) + this_rhs[lp.mat[elmnr].row_nr] += lp.mat[elmnr].value * values[i]; + for(i = 1; i <= lp.rows; i++) + { + dist = lp.orig_rh[i] - this_rhs[i]; + dist = round(dist, (float)0.001); + if((lp.orig_upbo[i] == 0 && dist != 0) || dist < 0) + { + return(FALSE); + } + } + return(TRUE); + } + + short column_in_lp(Problem lp, float[] testcolumn) + { + int i, j; + short ident; + float value; + + if (lp.scaling_used != FALSE) + for(i = 1; i <= lp.columns; i++) + { + ident = TRUE; + j = lp.col_end[i-1]; + while(ident != FALSE && (j < lp.col_end[i])) + { + value = lp.mat[j].value; + if(lp.ch_sign[lp.mat[j].row_nr] != FALSE) + value = -value; + value /= lp.scale[lp.rows+i]; + value /= lp.scale[lp.mat[j].row_nr]; + value -= testcolumn[lp.mat[j].row_nr]; + if(Math.abs(value) > (float)0.001) /* should be some epsilon? MB */ + ident=FALSE; + j++; + } + if(ident != FALSE) + return(TRUE); + } + else + for(i = 1; i <= lp.columns; i++) + { + ident = TRUE; + j = lp.col_end[i-1]; + while(ident != FALSE && j < lp.col_end[i]) + { + value = lp.mat[j].value; + if(lp.ch_sign[lp.mat[j].row_nr] != FALSE) + value *= -1; + value -= testcolumn[lp.mat[j].row_nr]; + if( Math.abs(value) > (float)0.001 ) + ident=FALSE; + j++; + } + if(ident != FALSE) + return(TRUE); + } + return(FALSE); + } + + private float minmax_to_scale(float min, float max) + { + float scale; + + /* should do something sensible when min or max is 0, MB */ + if((min == 0) || (max == 0)) + return((float)1); + + scale = (float)(1 / Math.pow(Math.E, (Math.log(min) + Math.log(max)) / 2)); + return(scale); + } + + public void unscale_columns(Problem lp) + { + int i, j; + + /* unscale mat */ + for(j = 1; j <= lp.columns; j++) + for(i = lp.col_end[j - 1]; i < lp.col_end[j]; i++) + lp.mat[i].value /= lp.scale[lp.rows + j]; + + /* unscale Bounds as well */ + for(i = lp.rows + 1; i < lp.sum; i++) + { + if(lp.orig_lowbo[i] != 0) + lp.orig_lowbo[i] *= lp.scale[i]; + if(lp.orig_upbo[i] != lp.infinite) + lp.orig_upbo[i] *= lp.scale[i]; + } + + for(i=lp.rows+1; i<= lp.sum; i++) + lp.scale[i]=1; + lp.columns_scaled=FALSE; + lp.eta_valid=FALSE; + } + + public void unscale(Problem lp) + { + int i, j; + + if (lp.scaling_used != FALSE) + { + + /* unscale mat */ + for(j = 1; j <= lp.columns; j++) + for(i = lp.col_end[j - 1]; i < lp.col_end[j]; i++) + lp.mat[i].value /= lp.scale[lp.rows + j]; + + /* unscale Bounds */ + for(i = lp.rows + 1; i < lp.sum; i++) + { + if(lp.orig_lowbo[i] != 0) + lp.orig_lowbo[i] *= lp.scale[i]; + if(lp.orig_upbo[i] != lp.infinite) + lp.orig_upbo[i] *= lp.scale[i]; + } + + /* unscale the matrix */ + for(j = 1; j <= lp.columns; j++) + for(i = lp.col_end[j-1]; i < lp.col_end[j]; i++) + lp.mat[i].value /= lp.scale[lp.mat[i].row_nr]; + + /* unscale the rhs! */ + for(i = 0; i <= lp.rows; i++) + lp.orig_rh[i] /= lp.scale[i]; + + lp.scaling_used=FALSE; + lp.eta_valid=FALSE; + } + } + + + public void auto_scale(Problem lp) + { + int i, j, row_nr, IntUsed; + float row_max[], row_min[], scalechange[], absval; + + if(lp.scaling_used == 0) + { + lp.scale = new float[lp.sum_alloc + 1]; + for(i = 0; i <= lp.sum; i++) + lp.scale[i]=1; + } + + row_max = new float[lp.rows + 1]; + row_min = new float[lp.rows + 1]; + scalechange = new float[lp.sum + 1]; + + /* initialise min and max values */ + for(i = 0; i <= lp.rows; i++) + { + row_max[i]=0; + row_min[i]=lp.infinite; + } + + /* calculate min and max absolute values of rows */ + for(j = 1; j <= lp.columns; j++) + for(i = lp.col_end[j - 1]; i < lp.col_end[j]; i++) + { + row_nr = lp.mat[i].row_nr; + absval = Math.abs(lp.mat[i].value); + if(absval!=0) + { + row_max[row_nr] = Math.max(row_max[row_nr], absval); + row_min[row_nr] = Math.min(row_min[row_nr], absval); + } + } + /* calculate scale factors for rows */ + for(i = 0; i <= lp.rows; i++) + { + scalechange[i] = minmax_to_scale(row_min[i], row_max[i]); + lp.scale[i] *= scalechange[i]; + } + + /* now actually scale the matrix */ + for(j = 1; j <= lp.columns; j++) + for(i = lp.col_end[j - 1]; i < lp.col_end[j]; i++) + lp.mat[i].value *= scalechange[lp.mat[i].row_nr]; + + /* and scale the rhs! */ + for(i = 0; i <= lp.rows; i++) + lp.orig_rh[i] *= scalechange[i]; + + IntUsed=FALSE; + i=lp.rows+1; + while(IntUsed == FALSE && i <= lp.sum) + { + IntUsed=FALSE; + i++; + } + if(IntUsed == FALSE) + { + float col_max, col_min; + + /* calculate scales */ + for(j = 1; j <= lp.columns; j++) + { + col_max = 0; + col_min = lp.infinite; + for(i = lp.col_end[j - 1]; i < lp.col_end[j]; i++) + { + if(lp.mat[i].value!=0) + { + col_max = Math.max(col_max, Math.abs(lp.mat[i].value)); + col_min = Math.min(col_min, Math.abs(lp.mat[i].value)); + } + } + scalechange[lp.rows + j] = minmax_to_scale(col_min, col_max); + lp.scale[lp.rows + j] *= scalechange[lp.rows + j]; + } + + /* scale mat */ + for(j = 1; j <= lp.columns; j++) + for(i = lp.col_end[j - 1]; i < lp.col_end[j]; i++) + lp.mat[i].value *= scalechange[lp.rows + j]; + + /* scale Bounds as well */ + for(i = lp.rows + 1; i < lp.sum; i++) + { + if(lp.orig_lowbo[i] != 0) + lp.orig_lowbo[i] /= scalechange[i]; + if(lp.orig_upbo[i] != lp.infinite) + lp.orig_upbo[i] /= scalechange[i]; + } + lp.columns_scaled=TRUE; + } + lp.scaling_used=TRUE; + lp.eta_valid=FALSE; + } + + public void reset_basis(Problem lp) { lp.basis_valid=FALSE; } + + /* Globals used by solver */ + short JustInverted; + short Status; + short Doiter; + short DoInvert; + short Break_bb; + + void set_globals(Problem lp) + { + if(Lp != null) + Lp.active = FALSE; + Lp = lp; + Rows = lp.rows; + columns = lp.columns; + Sum = Rows + columns; + Non_zeros = lp.non_zeros; + Mat = lp.mat; + Col_no = lp.col_no; + Col_end = lp.col_end; + Row_end = lp.row_end; + Rh = lp.rh; + Rhs = lp.rhs; + Orig_rh = lp.orig_rh; + Orig_upbo = lp.orig_upbo; + Orig_lowbo = lp.orig_lowbo; + Upbo = lp.upbo; + Lowbo = lp.lowbo; + Bas = lp.bas; + Basis = lp.basis; + Lower = lp.lower; + Eta_alloc = lp.eta_alloc; + Eta_size = lp.eta_size; + Num_inv = lp.num_inv; + Eta_value = lp.eta_value; + Eta_row_nr = lp.eta_row_nr; + Eta_col_end = lp.eta_col_end; + Solution = lp.solution; + Best_solution = lp.best_solution; + Infinite = lp.infinite; + Epsilon = lp.epsilon; + Epsb = lp.epsb; + Epsd = lp.epsd; + Epsel = lp.epsel; + + /* ?? MB */ + TREJ = TREJ; + TINV = TINV; + + Maximise = lp.maximise; + Floor_first = lp.floor_first; + lp.active = TRUE; + + // System.out.println("Sum = " + Sum); + } // end of set_global + + private void ftran(int start, + int end, + float[] pcol) + { + int i, j, k, r; + float theta; + + for(i = start; i <= end; i++) + { + k = Eta_col_end[i] - 1; + r = Eta_row_nr[k]; + theta = pcol[r]; + if(theta != 0) + for(j = Eta_col_end[i - 1]; j < k; j++) + pcol[Eta_row_nr[j]] += theta * Eta_value[j]; /* cpu expensive line */ + pcol[r] *= Eta_value[k]; + } + /* round small values to zero */ + for(i = 0; i <= Rows; i++) + round(pcol[i], Epsel); + } /* ftran */ + + private void btran(float[] row) + { + int i, j, k; + float f; + + for(i = Eta_size; i >= 1; i--) + { + f = 0; + k = Eta_col_end[i] - 1; + for(j = Eta_col_end[i - 1]; j <= k; j++) + f += row[Eta_row_nr[j]] * Eta_value[j]; + f = round(f, Epsel); + row[Eta_row_nr[k]] = f; + } + } /* btran */ + + + static short Isvalid(Problem lp) + { + int i, j, rownum[], colnum[]; + int num[], row_nr; + + if(lp.row_end_valid == FALSE) + { + num = new int[lp.rows + 1]; + rownum = new int[lp.rows + 1]; + for(i = 0; i <= lp.rows; i++) + { + num[i] = 0; + rownum[i] = 0; + } + for(i = 0; i < lp.non_zeros; i++) + rownum[lp.mat[i].row_nr]++; + lp.row_end[0] = 0; + for(i = 1; i <= lp.rows; i++) + lp.row_end[i] = lp.row_end[i - 1] + rownum[i]; + for(i = 1; i <= lp.columns; i++) + for(j = lp.col_end[i - 1]; j < lp.col_end[i]; j++) + { + row_nr = lp.mat[j].row_nr; + if(row_nr != 0) + { + num[row_nr]++; + lp.col_no[lp.row_end[row_nr - 1] + num[row_nr]] = i; + } + } + lp.row_end_valid = TRUE; + } + if(lp.valid != FALSE) + return(TRUE); + rownum = new int[lp.rows + 1]; + for (i = 0; i <= lp.rows; i++) + rownum[i] = 0; + colnum = new int[lp.columns + 1]; + for (i = 0; i <= lp.columns; i++) + colnum[i] = 0; + for(i = 1 ; i <= lp.columns; i++) + for(j = lp.col_end[i - 1]; j < lp.col_end[i]; j++) + { + colnum[i]++; + rownum[lp.mat[j].row_nr]++; + } + for(i = 1; i <= lp.columns; i++) + if(colnum[i] == 0) + if (lp.names_used != FALSE) + System.err.print("Warning: Variable " + lp.col_name[i] + + " not used in any constraints\n"); + else + System.err.print("Warning: Variable " + i + " not used in any constaints\n"); + lp.valid = TRUE; + return(TRUE); + } + + private void resize_eta() + { + Eta_alloc *= 2; + float[] db_ptr = Eta_value; + Eta_value = new float[Eta_alloc]; + System.arraycopy(db_ptr, 0, Eta_value, 0, db_ptr.length); + Lp.eta_value = Eta_value; + + int[] int_ptr = Eta_row_nr; + Eta_row_nr = new int[Eta_alloc]; + System.arraycopy(int_ptr, 0, Eta_row_nr, 0, int_ptr.length); + Lp.eta_row_nr = Eta_row_nr; + } /* resize_eta */ + + + private void condensecol(int row_nr, + float[] pcol) + { + int i, elnr; + + elnr = Eta_col_end[Eta_size]; + + if(elnr + Rows + 2 > Eta_alloc) /* maximum local growth of Eta */ + resize_eta(); + + for(i = 0; i <= Rows; i++) + if(i != row_nr && pcol[i] != 0) + { + Eta_row_nr[elnr] = i; + Eta_value[elnr] = pcol[i]; + elnr++; + } + Eta_row_nr[elnr] = row_nr; + Eta_value[elnr] = pcol[row_nr]; + elnr++; + Eta_col_end[Eta_size + 1] = elnr; + } /* condensecol */ + + + private void addetacol() + { + int i, j, k; + float theta; + + j = Eta_col_end[Eta_size]; + Eta_size++; + k = Eta_col_end[Eta_size]; + theta = 1 / (float) Eta_value[k - 1]; + Eta_value[k - 1] = theta; + for(i = j; i < k - 1; i++) + Eta_value[i] *= -theta; + JustInverted = FALSE; + } /* addetacol */ + + private void setpivcol(short lower, + int varin, + float[] pcol) + { + int i, colnr; + float f; + + if(lower != FALSE) + f = 1; + else + f = -1; + for(i = 0; i <= Rows; i++) + pcol[i] = 0; + if(varin > Rows) + { + colnr = varin - Rows; + for(i = Col_end[colnr - 1]; i < Col_end[colnr]; i++) + pcol[Mat[i].row_nr] = Mat[i].value * f; + pcol[0] -= Extrad * f; + } + else + if(lower != FALSE) + pcol[varin] = 1; + else + pcol[varin] = -1; + ftran(1, Eta_size, pcol); + } /* setpivcol */ + + + private void minoriteration(int colnr, + int row_nr) + { + int i, j, k, wk, varin, varout, elnr; + float piv = 0, theta; + + varin = colnr + Rows; + elnr = Eta_col_end[Eta_size]; + wk = elnr; + Eta_size++; + if(Extrad != 0) + { + Eta_row_nr[elnr] = 0; + Eta_value[elnr] = -Extrad; + elnr++; + } + for(j = Col_end[colnr - 1] ; j < Col_end[colnr]; j++) + { + k = Mat[j].row_nr; + if(k == 0 && Extrad != 0) + Eta_value[Eta_col_end[Eta_size -1]] += Mat[j].value; + else if(k != row_nr) + { + Eta_row_nr[elnr] = k; + Eta_value[elnr] = Mat[j].value; + elnr++; + } + else + piv = Mat[j].value; + } + Eta_row_nr[elnr] = row_nr; + Eta_value[elnr] = 1 / (float) piv; + elnr++; + theta = Rhs[row_nr] / (float) piv; + Rhs[row_nr] = theta; + for(i = wk; i < elnr - 1; i++) + Rhs[Eta_row_nr[i]] -= theta * Eta_value[i]; + varout = Bas[row_nr]; + Bas[row_nr] = varin; + Basis[varout] = FALSE; + Basis[varin] = TRUE; + for(i = wk; i < elnr - 1; i++) + Eta_value[i] /= - (float) piv; + Eta_col_end[Eta_size] = elnr; + } /* minoriteration */ + + private void rhsmincol(float theta, + int row_nr, + int varin) + { + int i, j, k, varout; + float f; + + if(row_nr > Rows + 1) + { + System.err.print("Error: rhsmincol called with row_nr: " + + row_nr + ", rows: " + Rows + "\n"); + System.err.print("This indicates numerical instability\n"); + System.exit(FAIL); + } + j = Eta_col_end[Eta_size]; + k = Eta_col_end[Eta_size + 1]; + for(i = j; i < k; i++) + { + f = Rhs[Eta_row_nr[i]] - theta * Eta_value[i]; + f = round(f, Epsb); + Rhs[Eta_row_nr[i]] = f; + } + Rhs[row_nr] = theta; + varout = Bas[row_nr]; + Bas[row_nr] = varin; + Basis[varout] = FALSE; + Basis[varin] = TRUE; + } /* rhsmincol */ + + + void invert() + { + int i, j, v, wk, numit, varnr, row_nr, colnr, varin; + float theta; + float[] pcol; + short[] frow; + short[] fcol; + int[] rownum, col, row; + int[] colnum; + + rownum = new int[Rows + 1]; + for (i = 0; i <= Rows; i++) + rownum[i] = 0; + col = new int[Rows + 1]; + for (i = 0; i <= Rows; i++) + col[i] = 0; + row = new int[Rows + 1]; + for (i = 0; i <= Rows; i++) + row[i] = 0; + pcol = new float[Rows + 1]; + for (i = 0; i <= Rows; i++) + pcol[i] = 0; + frow = new short[Rows + 1]; + for(i = 0; i <= Rows; i++) + frow[i] = TRUE; + fcol = new short[columns + 1]; + for(i = 0; i < columns; i++) + fcol[i] = FALSE; + colnum = new int[columns + 1]; + for(i = 0; i <= columns; i++) + colnum[i] = 0; + + for(i = 0; i <= Rows; i++) + if(Bas[i] > Rows) + fcol[Bas[i] - Rows - 1] = TRUE; + else + frow[Bas[i]] = FALSE; + + for(i = 1; i <= Rows; i++) + if(frow[i] != FALSE) + for(j = Row_end[i - 1] + 1; j <= Row_end[i]; j++) + { + wk = Col_no[j]; + if(fcol[wk - 1] != FALSE) + { + colnum[wk]++; + rownum[i - 1]++; + } + } + for(i = 1; i <= Rows; i++) + Bas[i] = i; + for(i = 1; i <= Rows; i++) + Basis[i] = TRUE; + for(i = 1; i <= columns; i++) + Basis[i + Rows] = FALSE; + + for(i = 0; i <= Rows; i++) + Rhs[i] = Rh[i]; + + for(i = 1; i <= columns; i++) + { + varnr = Rows + i; + if(Lower[varnr] == FALSE) + { + theta = Upbo[varnr]; + for(j = Col_end[i - 1]; j < Col_end[i]; j++) + Rhs[Mat[j].row_nr] -= theta * Mat[j].value; + } + } + for(i = 1; i <= Rows; i++) + if(Lower[i] == FALSE) + Rhs[i] -= Upbo[i]; + Eta_size = 0; + v = 0; + row_nr = 0; + Num_inv = 0; + numit = 0; + while(v < Rows) + { + row_nr++; + if(row_nr > Rows) + row_nr = 1; + v++; + if(rownum[row_nr - 1] == 1) + if(frow[row_nr] != FALSE) + { + v = 0; + j = Row_end[row_nr - 1] + 1; + while(fcol[Col_no[j] - 1] == FALSE) + j++; + colnr = Col_no[j]; + fcol[colnr - 1] = FALSE; + colnum[colnr] = 0; + for(j = Col_end[colnr - 1]; j < Col_end[colnr]; j++) + if(frow[Mat[j].row_nr] != FALSE) + rownum[Mat[j].row_nr - 1]--; + frow[row_nr] = FALSE; + minoriteration(colnr, row_nr); + } + } + v = 0; + colnr = 0; + while(v columns) + colnr = 1; + v++; + if(colnum[colnr] == 1) + if(fcol[colnr - 1] != FALSE) + { + v = 0; + j = Col_end[colnr - 1] + 1; + while(frow[Mat[j - 1].row_nr] == FALSE) + j++; + row_nr = Mat[j - 1].row_nr; + frow[row_nr] = FALSE; + rownum[row_nr - 1] = 0; + for(j = Row_end[row_nr - 1] + 1; j <= Row_end[row_nr]; j++) + if(fcol[Col_no[j] - 1] != FALSE) + colnum[Col_no[j]]--; + fcol[colnr - 1] = FALSE; + numit++; + col[numit - 1] = colnr; + row[numit - 1] = row_nr; + } + } + for(j = 1; j <= columns; j++) + if(fcol[j - 1] != FALSE) + { + fcol[j - 1] = FALSE; + setpivcol(Lower[Rows + j], j + Rows, pcol); + row_nr = 1; + while((frow[row_nr] == FALSE || pcol[row_nr] == FALSE) + && row_nr <= Rows) + row_nr++; /* this sometimes sets row_nr to Rows + 1 and makes + rhsmincol crash. Solved in 2.0? MB */ + if(row_nr == Rows + 1) + System.err.print("Inverting failed"); + frow[row_nr] = FALSE; + condensecol(row_nr, pcol); + theta = Rhs[row_nr] / (float) pcol[row_nr]; + rhsmincol(theta, row_nr, Rows + j); + addetacol(); + } + for(i = numit - 1; i >= 0; i--) + { + colnr = col[i]; + row_nr = row[i]; + varin = colnr + Rows; + for(j = 0; j <= Rows; j++) + pcol[j] = 0; + for(j = Col_end[colnr - 1]; j < Col_end[colnr]; j++) + pcol[Mat[j].row_nr] = Mat[j].value; + pcol[0] -= Extrad; + condensecol(row_nr, pcol); + theta = Rhs[row_nr] / (float) pcol[row_nr]; + rhsmincol(theta, row_nr, varin); + addetacol(); + } + for(i = 1; i <= Rows; i++) + Rhs[i] = round(Rhs[i], Epsb); + JustInverted = TRUE; + DoInvert = FALSE; + } /* invert */ + + private short colprim(Ref colnr, + short minit, + float[] drow) + { + int varnr, i, j; + float f, dpiv; + + dpiv = -Epsd; + colnr.value = 0; + if(minit == FALSE) + { + for(i = 1; i <= Sum; i++) + drow[i] = 0; + drow[0] = 1; + btran(drow); + for(i = 1; i <= columns; i++) + { + varnr = Rows + i; + if(Basis[varnr] == FALSE) + if(Upbo[varnr] > 0) + { + f = 0; + for(j = Col_end[i - 1]; j < Col_end[i]; j++) + f += drow[Mat[j].row_nr] * Mat[j].value; + drow[varnr] = f; + } + } + for(i = 1; i <= Sum; i++) + drow[i] = round(drow[i], Epsd); + } + /* + System.out.println("dpiv = " + dpiv); + System.out.println("drow[] at colprim"); + for(i = 1; i <= Sum; i++) // DEBUG + { + System.out.print(drow[i] + " "); + } + System.out.println(); + */ + for(i = 1; i <= Sum; i++) + if(Basis[i] == FALSE) + if(Upbo[i] > 0) + { + if(Lower[i] != FALSE) + f = drow[i]; + else + f = -drow[i]; + if(f < dpiv) + { + dpiv = f; + colnr.value = i; + } + } + if(Lp.trace != FALSE) + if(colnr.value > 0) + System.err.print("col_prim:" + colnr.value + ", reduced cost: " + dpiv + "\n"); + else + System.err.print("col_prim: no negative reduced costs found, optimality!\n"); + if(colnr.value == 0) + { + Doiter = FALSE; + DoInvert = FALSE; + Status = OPTIMAL; + } + return(colnr.value > 0 ? (short)1 : (short)0); + } /* colprim */ + + private short rowprim(int colnr, + Ref row_nr, + Ref theta, + float[] pcol) + { + int i; + float f = 0, quot; + + row_nr.value = 0; + theta.value = Infinite; + for(i = 1; i <= Rows; i++) + { + f = pcol[i]; + if(Math.abs(f) < TREJ) + f = 0; + if(f != 0) + { + quot = 2 * Infinite; + if(f > 0) + quot = Rhs[i] / (float) f; + else + if(Upbo[Bas[i]] < Infinite) + quot = (Rhs[i] - Upbo[Bas[i]]) / (float) f; + round(quot, Epsel); + if(quot < theta.value) + { + theta.value = quot; + row_nr.value = i; + } + } + } + if(row_nr.value == 0) + for(i = 1; i <= Rows; i++) + { + f = pcol[i]; + if(f != 0) + { + quot = 2 * Infinite; + if(f > 0) + quot = Rhs[i] / (float) f; + else + if(Upbo[Bas[i]] < Infinite) + quot = (Rhs[i] - Upbo[Bas[i]]) / (float) f; + quot = round(quot, Epsel); + if(quot < theta.value) + { + theta.value = quot; + row_nr.value = i; + } + } + } + + if(theta.value < 0) + { + System.err.println("Warning: Numerical instability, qout = " + theta.value); + System.err.println("pcol[" + row_nr.value + "] = " + f + ", rhs[" + + row_nr.value + "] = " + Rhs[(int)row_nr.value] + + " , upbo = " + Upbo[Bas[(int)row_nr.value]]); + } + if(row_nr.value == 0) + { + if(Upbo[colnr] == Infinite) + { + Doiter = FALSE; + DoInvert = FALSE; + Status = UNBOUNDED; + } + else + { + i = 1; + while(pcol[i] >= 0 && i <= Rows) + i++; + if(i > Rows) /* empty column with upperBound! */ + { + Lower[colnr] = FALSE; + Rhs[0] += Upbo[colnr]*pcol[0]; + Doiter = FALSE; + DoInvert = FALSE; + } + else if(pcol[i]<0) + { + row_nr.value = i; + } + } + } + if(row_nr.value > 0) + Doiter = TRUE; + if(Lp.trace != FALSE) + System.out.println("row_prim:" + row_nr.value + ", pivot element:" + + pcol[(int)row_nr.value]); + + return((row_nr.value > 0) ? (short)1 : (short)0); + } /* rowprim */ + + private short rowdual(Ref row_nr) + { + int i; + float f, g, minrhs; + short artifs; + + row_nr.value = 0; + minrhs = -Epsb; + i = 0; + artifs = FALSE; + while(i < Rows && artifs == FALSE) + { + i++; + f = Upbo[Bas[i]]; + if(f == 0 && (Rhs[i] != 0)) + { + artifs = TRUE; + row_nr.value = i; + } + else + { + if(Rhs[i] < f - Rhs[i]) + g = Rhs[i]; + else + g = f - Rhs[i]; + if(g < minrhs) + { + minrhs = g; + row_nr.value = i; + } + } + } + + if(Lp.trace != FALSE) + { + if(row_nr.value > 0) + { + System.out.println("row_dual:" + row_nr.value + + ", rhs of selected row: " + + Rhs[(int)row_nr.value]); + if(Upbo[Bas[(int)row_nr.value]] < Infinite) + System.out.println(" upper Bound of basis variable: " + + Upbo[Bas[(int)row_nr.value]]); + } + else + System.out.print("row_dual: no infeasibilities found\n"); + } + + return(row_nr.value > 0 ? (short)1 : (short)0); + } /* rowdual */ + + private short coldual(int row_nr, + Ref colnr, + short minit, + float[] prow, + float[] drow) + { + int i, j, r, varnr; + float theta, quot, pivot, d, f, g; + + Doiter = FALSE; + if(minit == FALSE) + { + for(i = 0; i <= Rows; i++) + { + prow[i] = 0; + drow[i] = 0; + } + drow[0] = 1; + prow[row_nr] = 1; + for(i = Eta_size; i >= 1; i--) + { + d = 0; + f = 0; + r = Eta_row_nr[Eta_col_end[i] - 1]; + for(j = Eta_col_end[i - 1]; j < Eta_col_end[i]; j++) + { + /* this is where the program consumes most cpu time */ + f += prow[Eta_row_nr[j]] * Eta_value[j]; + d += drow[Eta_row_nr[j]] * Eta_value[j]; + } + f = round(f, Epsel); + prow[r] = f; + d = round(d, Epsel); + drow[r] = d; + } + for(i = 1; i <= columns; i++) + { + varnr = Rows + i; + if(Basis[varnr] == FALSE) + { + d = - Extrad * drow[0]; + f = 0; + for(j = Col_end[i - 1]; j < Col_end[i]; j++) + { + d = d + drow[Mat[j].row_nr] * Mat[j].value; + f = f + prow[Mat[j].row_nr] * Mat[j].value; + } + drow[varnr] = d; + prow[varnr] = f; + } + } + for(i = 0; i <= Sum; i++) + { + prow[i] = round(prow[i], Epsel); + drow[i] = round(drow[i], Epsd); + } + } + if(Rhs[row_nr] > Upbo[Bas[row_nr]]) + g = -1; + else + g = 1; + pivot = 0; + colnr.value = 0; + theta = Infinite; + for(i = 1; i <= Sum; i++) + { + if(Lower[i] != FALSE) + d = prow[i] * g; + else + d = -prow[i] * g; + if((d < 0) && (Basis[i] == FALSE) && (Upbo[i] > 0)) + { + if(Lower[i] == FALSE) + quot = -drow[i] / (float) d; + else + quot = drow[i] / (float) d; + if(quot < theta) + { + theta = quot; + pivot = d; + colnr.value = i; + } + else if((quot == theta) && (Math.abs(d) > Math.abs(pivot))) + { + pivot = d; + colnr.value = i; + } + } + } + if(Lp.trace != FALSE) + System.out.println("col_dual:" + colnr.value + ", pivot element: " + + prow[(int)colnr.value]); + + if(colnr.value > 0) + Doiter = TRUE; + + return(colnr.value > 0 ? (short)1 : (short)0); + } /* coldual */ + + private void iteration(int row_nr, + int varin, + Ref theta, + float up, + Ref minit, + Ref low, + short primal, + float[] pcol) + { + int i, k, varout; + float f; + float pivot; + + Lp.iter++; + minit.value = theta.value > (up + Epsb) ? 1 : 0; + if(minit.value != 0) + { + theta.value = up; + low.value = low.value == 0 ? 1 : 0; + } + k = Eta_col_end[Eta_size + 1]; + pivot = Eta_value[k - 1]; + for(i = Eta_col_end[Eta_size]; i < k; i++) + { + f = Rhs[Eta_row_nr[i]] - theta.value * Eta_value[i]; + f = round(f, Epsb); + Rhs[Eta_row_nr[i]] = f; + } + if(minit.value == 0) + { + Rhs[row_nr] = theta.value; + varout = Bas[row_nr]; + Bas[row_nr] = varin; + Basis[varout] = FALSE; + Basis[varin] = TRUE; + if(primal != FALSE && pivot < 0) + Lower[varout] = FALSE; + if(low.value == 0 && up < Infinite) + { + low.value = TRUE; + Rhs[row_nr] = up - Rhs[row_nr]; + for(i = Eta_col_end[Eta_size]; i < k; i++) + Eta_value[i] = -Eta_value[i]; + } + addetacol(); + Num_inv++; + } + if(Lp.trace != FALSE) + { + System.out.print("Theta = " + theta.value + " "); + if(minit.value != 0) + { + if(Lower[varin] == FALSE) + System.out.print("Iteration:" + Lp.iter + + ", variable" + varin + " changed from 0 to its upper Bound of " + + Upbo[varin] + "\n"); + else + System.out.print("Iteration:" + Lp.iter + ", variable" + varin + + " changed its upper Bound of " + Upbo[varin] + " to 0\n"); + } + else + System.out.print("Iteration:" + Lp.iter + ", variable" + varin + + " entered basis at:" + Rhs[row_nr] + "\n"); + if(primal == 0) + { + f = 0; + for(i = 1; i <= Rows; i++) + if(Rhs[i] < 0) + f -= Rhs[i]; + else + if(Rhs[i] > Upbo[Bas[i]]) + f += Rhs[i] - Upbo[Bas[i]]; + System.out.println("feasibility gap of this basis:" + (float)f); + } + else + System.out.println("objective function value of this feasible basis: " + Rhs[0]); + } + } /* iteration */ + + + private int solvelp() + { + int i, j, varnr; + float f = 0, theta = 0; + short primal; + float[] drow, prow, Pcol; + short minit; + int colnr, row_nr; + colnr = 0; + row_nr = 0; + short flag; + Ref ref1, ref2, ref3; + ref1 = new Ref(0); + ref2 = new Ref(0); + ref3 = new Ref(0); + + drow = new float[Sum + 1]; + prow = new float[Sum + 1]; + Pcol = new float[Rows + 1]; + for (i = 0; i <= Sum; i++) { + drow[i] = 0; + prow[i] = 0; + } + for (i = 0; i <= Rows; i++) + Pcol[i] = 0; + + + Lp.iter = 0; + minit = FALSE; + Status = RUNNING; + DoInvert = FALSE; + Doiter = FALSE; + i = 0; + primal = TRUE; + while(i != Rows && primal != FALSE) + { + i++; + primal = (Rhs[i] >= 0 && Rhs[i] <= Upbo[Bas[i]]) ? + (short)1: (short)0; + } + if(Lp.trace != FALSE) + { + if(primal != FALSE) + System.out.print("Start at feasible basis\n"); + else + System.out.print("Start at infeasible basis\n"); + } + if(primal == FALSE) + { + drow[0] = 1; + for(i = 1; i <= Rows; i++) + drow[i] = 0; + Extrad = 0; + for(i = 1; i <= columns; i++) + { + varnr = Rows + i; + drow[varnr] = 0; + for(j = Col_end[i - 1]; j < Col_end[i]; j++) + if(drow[Mat[j].row_nr] != 0) + drow[varnr] += drow[Mat[j].row_nr] * Mat[j].value; + if(drow[varnr] < Extrad) + Extrad = drow[varnr]; + } + } + else + Extrad = 0; + if(Lp.trace != FALSE) + System.out.println("Extrad = " + Extrad); + minit = FALSE; + + while(Status == RUNNING) + { + Doiter = FALSE; + DoInvert = FALSE; + + construct_solution(Solution); + if(primal != FALSE) + { + ref1.value = colnr; + flag = colprim(ref1, minit, drow); + colnr = (int)ref1.value; + if (flag != FALSE) + { + setpivcol(Lower[colnr], colnr, Pcol); + ref1.value = row_nr; + ref2.value = theta; + flag = rowprim(colnr, ref1, ref2, Pcol); + row_nr = (int)ref1.value; + theta = ref2.value; + if (flag != FALSE) + condensecol(row_nr, Pcol); + + } + } + else /* not primal */ + { + if(minit == FALSE) { + ref1.value = row_nr; + flag = rowdual(ref1); + row_nr = (int)ref1.value; + } + if(row_nr > 0 ) + { + ref1.value = colnr; + flag = coldual(row_nr, ref1, minit, prow, drow); + colnr = (int)ref1.value; + if (flag != FALSE) + { + setpivcol(Lower[colnr], colnr, Pcol); + /* getting div by zero here ... MB */ + if(Pcol[row_nr] == 0) + { + System.err.println("An attempt was made to divide by zero (Pcol[" + + row_nr + "])"); + System.err.println("This indicates numerical instability"); + Doiter = FALSE; + if(JustInverted == FALSE) + { + System.out.println("Reinverting Eta"); + DoInvert = TRUE; + } + else + { + System.out.println("Can't reinvert, failure"); + Status = FAILURE; + } + } + else + { + condensecol(row_nr, Pcol); + f = Rhs[row_nr] - Upbo[Bas[row_nr]]; + if(f > 0) + { + theta = f / (float) Pcol[row_nr]; + if(theta <= Upbo[colnr]) + Lower[Bas[row_nr]] = + (Lower[Bas[row_nr]] == FALSE)? + (short)1:(short)0; + } + else /* f <= 0 */ + theta = Rhs[row_nr] / (float) Pcol[row_nr]; + } + } + else + Status = INFEASIBLE; + } + else + { + primal = TRUE; + Doiter = FALSE; + Extrad = 0; + DoInvert = TRUE; + } + } + if(Doiter != FALSE) { + ref1.value = theta; + ref2.value = minit; + ref3.value = Lower[colnr]; + iteration(row_nr, colnr, ref1, Upbo[colnr], ref2, ref3, primal, Pcol); + theta = ref1.value; + minit = (short)ref2.value; + Lower[colnr] = (short)ref3.value; + } + if(Num_inv >= Lp.max_num_inv) + DoInvert = TRUE; + if(DoInvert != FALSE) + { + invert(); + } + } + + Lp.total_iter += Lp.iter; + + return(Status); + } /* solvelp */ + + + private short is_int(float value) + { + float tmp; + + tmp = (float)(value - Math.floor(value)); + if(tmp < Epsilon) + return(TRUE); + if(tmp > (1 - Epsilon)) + return(TRUE); + return(FALSE); + } /* is_int */ + + private void construct_solution(float[] sol) + { + int i, j, basi; + float f; + + /* zero all results of rows */ + for (i = 0; i <= Rows; i++) + sol[i] = 0; + + if (Lp.scaling_used != FALSE) + { + for(i = Rows + 1; i <= Sum; i++) + sol[i] = Lowbo[i] * Lp.scale[i]; + for(i = 1; i <= Rows; i++) + { + basi = Bas[i]; + if(basi > Rows) + sol[basi] += Rhs[i] * Lp.scale[basi]; + } + for(i = Rows + 1; i <= Sum; i++) + if(Basis[i] == FALSE && Lower[i] == FALSE) + sol[i] += Upbo[i] * Lp.scale[i]; + + for(j = 1; j <= columns; j++) + { + f = sol[Rows + j]; + if(f != 0) + for(i = Col_end[j - 1]; i < Col_end[j]; i++) + sol[Mat[i].row_nr] += (f / Lp.scale[Rows+j]) + * (Mat[i].value / Lp.scale[Mat[i].row_nr]); + } + + for(i = 0; i <= Rows; i++) + { + if(Math.abs(sol[i]) < Epsb) + sol[i] = 0; + else + if(Lp.ch_sign[i] != FALSE) + sol[i] = -sol[i]; + } + } + else + { + for(i = Rows + 1; i <= Sum; i++) + sol[i] = Lowbo[i]; + for(i = 1; i <= Rows; i++) + { + basi = Bas[i]; + if(basi > Rows) + sol[basi] += Rhs[i]; + } + for(i = Rows + 1; i <= Sum; i++) + if(Basis[i] == FALSE && Lower[i] == FALSE) + sol[i] += Upbo[i]; + for(j = 1; j <= columns; j++) + { + f = sol[Rows + j]; + if(f != 0) + for(i = Col_end[j - 1]; i < Col_end[j]; i++) + sol[Mat[i].row_nr] += f * Mat[i].value; + } + + for(i = 0; i <= Rows; i++) + { + if(Math.abs(sol[i]) < Epsb) + sol[i] = 0; + else + if(Lp.ch_sign[i] != FALSE) + sol[i] = -sol[i]; + } + } + } /* construct_solution */ + + private void calculate_duals() + { + int i; + + /* initialise */ + for(i = 1; i <= Rows; i++) + Lp.duals[i] = 0; + Lp.duals[0] = 1; + btran(Lp.duals); + if (Lp.scaling_used != FALSE) + for(i = 1; i <= Rows; i++) + Lp.duals[i] *= Lp.scale[i]/Lp.scale[0]; + + /* the dual values are the reduced costs of the slacks */ + /* When the slack is at its upper Bound, change the sign. Can this happen? */ + for(i = 1; i <= Rows; i++) + { + if(Lp.basis[i] != FALSE) + Lp.duals[i] = 0; + else if( Lp.ch_sign[0] == Lp.ch_sign[i]) + Lp.duals[i] = -Lp.duals[i]; + } + } + + private int milpsolve(float[] upbo, + float[] lowbo, + short[] sbasis, + short[] slower, + int[] sbas) + { + int i, j, failure, notint, is_worse; + float theta, tmpfloat; + Random rdm = new Random(); + notint = 0; + + if(Break_bb != FALSE) + return(BREAK_BB); + Level++; + Lp.total_nodes++; + if(Level > Lp.max_level) + Lp.max_level = Level; + /* make fresh copies of upbo, lowbo, rh as solving changes them */ + System.arraycopy(upbo, 0, Upbo, 0, Sum + 1); + System.arraycopy(lowbo, 0, Lowbo, 0, Sum + 1); + System.arraycopy(sbasis, 0, Basis, 0, Sum + 1); + System.arraycopy(slower, 0, Lower, 0, Sum + 1); + System.arraycopy(sbas, 0, Bas, 0, Rows + 1); + System.arraycopy(Orig_rh, 0, Rh, 0, Rows + 1); + + if(Lp.anti_degen != FALSE) + { + for(i = 1; i <= columns; i++) + { + tmpfloat = rdm.nextFloat()*(float)0.001; + if(tmpfloat > Epsb) + Lowbo[i + Rows] -= tmpfloat; + tmpfloat = rdm.nextFloat()*(float)0.001; + if(tmpfloat > Epsb) + Upbo[i + Rows] += tmpfloat; + } + Lp.eta_valid = FALSE; + } + + if(Lp.eta_valid == FALSE) + { + for(i = 1; i <= columns; i++) + if(Lowbo[Rows + i] != 0) + { + theta = Lowbo[ Rows + i]; + if(Upbo[Rows + i]= Best_solution[0]) ? 1:0; + + if(is_worse != FALSE) + { + if(Lp.verbose != FALSE) + System.out.println("level" + Level + " OPT NOB value " + Solution[0] + + " Bound " + Best_solution[0]); + Level--; + return(MILP_FAIL); + } + + /* check if solution contains enough ints */ + if(Lp.bb_rule == FIRST_NI) + { + notint = 0; + i = Rows + 1; + while(i <= Sum && notint == 0) + { + i++; + } + } + if(Lp.bb_rule == RAND_NI) + { + int nr_not_int, select_not_int; + nr_not_int = 0; + for(i = Rows + 1; i <= Sum; i++) + if(nr_not_int == 0) + notint = 0; + else + { + select_not_int=(rdm.nextInt() % nr_not_int) + 1; + i = Rows + 1; + while(select_not_int > 0) + { + i++; + } + notint = i - 1; + } + } + + if(Lp.verbose == TRUE) + if(notint != FALSE) + System.out.println("level " + Level + " OPT value " + Solution[0]); + else + System.out.println("level " + Level + " OPT INT value " + Solution[0]); + + if(notint != FALSE) /* there is at least one value not yet int */ + { + /* set up two new problems */ + float[] new_upbo, new_lowbo; + float new_bound; + short[] new_lower, new_basis; + int[] new_bas; + int resone, restwo; + + /* allocate room for them */ + new_upbo = new float[Sum + 1]; + new_lowbo = new float[Sum + 1]; + new_lower = new short[Sum + 1]; + new_basis = new short[Sum + 1]; + new_bas = new int[Rows + 1]; + System.arraycopy(upbo, 0, new_upbo, 0, Sum + 1); + System.arraycopy(lowbo, 0, new_lowbo, 0, Sum + 1); + System.arraycopy(Lower, 0, new_lower, 0, Sum + 1); + System.arraycopy(Basis, 0, new_basis, 0, Sum + 1); + System.arraycopy(Bas, 0, new_bas, 0, Rows +1); + + if(Floor_first != FALSE) + { + new_bound = (float)(Math.ceil(Solution[notint]) - 1); + /* this Bound might conflict */ + if(new_bound < lowbo[notint]) + { + resone = MILP_FAIL; + } + else /* Bound feasible */ + { + new_upbo[notint] = new_bound; + Lp.eta_valid = FALSE; + resone = milpsolve(new_upbo, lowbo, new_basis, new_lower, + new_bas); + Lp.eta_valid = FALSE; + } + new_bound += 1; + if(new_bound > upbo[notint]) + { + restwo = MILP_FAIL; + } + else /* Bound feasible */ + { + new_lowbo[notint] = new_bound; + Lp.eta_valid = FALSE; + restwo = milpsolve(upbo, new_lowbo, new_basis, new_lower, + new_bas); + Lp.eta_valid = FALSE; + } + } + else /* take ceiling first */ + { + new_bound = (float)Math.ceil(Solution[notint]); + /* this Bound might conflict */ + if(new_bound > upbo[notint]) + { + resone = MILP_FAIL; + } + else /* Bound feasible */ + { + new_lowbo[notint] = new_bound; + Lp.eta_valid = FALSE; + resone = milpsolve(upbo, new_lowbo, new_basis, new_lower, + new_bas); + Lp.eta_valid = FALSE; + } + new_bound -= 1; + if(new_bound < lowbo[notint]) + { + restwo = MILP_FAIL; + } + else /* Bound feasible */ + { + new_upbo[notint] = new_bound; + Lp.eta_valid = FALSE; + restwo = milpsolve(new_upbo, lowbo, new_basis, new_lower, + new_bas); + Lp.eta_valid = FALSE; + } + } + if(resone != FALSE && restwo != FALSE) + /* both failed and must have been infeasible */ + failure = INFEASIBLE; + else + failure = OPTIMAL; + + } + else /* all required values are int */ + { + if(Maximise != FALSE) + is_worse = (Solution[0] < Best_solution[0]) ? 1:0; + else + is_worse = (Solution[0] > Best_solution[0]) ? 1:0; + + if(is_worse == FALSE) /* Current solution better */ + { + System.arraycopy(Solution, 0, Best_solution, 0, Sum + 1); + calculate_duals(); + if(Lp.break_at_int != FALSE) + { + if(Maximise != FALSE && (Best_solution[0] > Lp.break_value)) + Break_bb = TRUE; + if(Maximise == FALSE && (Best_solution[0] < Lp.break_value)) + Break_bb = TRUE; + } + } + } + } + + /* failure can have the values OPTIMAL, UNBOUNDED and INFEASIBLE. */ + Level--; + return(failure); + } /* milpsolve */ + + public int solve(Problem lp) + { + int result, i; + + if(lp.active == FALSE) + set_globals(lp); + + lp.total_iter = 0; + lp.max_level = 1; + lp.total_nodes = 0; + + if(Isvalid(lp) != FALSE) + { + if(Maximise != FALSE && lp.obj_bound == Infinite) + Best_solution[0]=-Infinite; + else if(Maximise == FALSE && lp.obj_bound==-Infinite) + Best_solution[0] = Infinite; + else + Best_solution[0] = lp.obj_bound; + + Level = 0; + + if(lp.basis_valid == FALSE) + { + for(i = 0; i <= lp.rows; i++) + { + lp.basis[i] = TRUE; + lp.bas[i] = i; + } + for(i = lp.rows+1; i <= lp.sum; i++) + lp.basis[i] = FALSE; + for(i = 0; i <= lp.sum; i++) + lp.lower[i] = TRUE; + lp.basis_valid = TRUE; + } + + lp.eta_valid = FALSE; + Break_bb = FALSE; + result = milpsolve(Orig_upbo, Orig_lowbo, Basis, Lower, Bas); + lp.eta_size = Eta_size; + lp.eta_alloc = Eta_alloc; + lp.num_inv = Num_inv; + + return(result); + } + + /* if we get here, Isvalid(lp) failed. I suggest we return FAILURE */ + return(FAILURE); + } + + public int lag_solve(Problem lp, float start_bound, int num_iter, short verbose) + { + int i, j, result, citer; + short status, OrigFeas, AnyFeas, same_basis; + float[] OrigObj, ModObj, SubGrad, BestFeasSol; + float Zub, Zlb, Ztmp, pie; + float rhsmod, Step, SqrsumSubGrad; + int[] old_bas; + short[] old_lower; + + /* allocate mem */ + OrigObj = new float[lp.columns + 1]; + ModObj = new float[lp.columns + 1]; + for (i = 0; i <= lp.columns; i++) + ModObj[i] = 0; + SubGrad = new float[lp.nr_lagrange]; + for (i = 0; i < lp.nr_lagrange; i++) + SubGrad[i] = 0; + BestFeasSol = new float[lp.sum + 1]; + for (i = 0; i <= lp.sum; i++) + BestFeasSol[i] = 0; + old_bas = new int[lp.rows + 1]; + System.arraycopy(lp.bas, 0, old_bas, 0, lp.rows + 1); + old_lower = new short[lp.sum + 1]; + System.arraycopy(lp.lower, 0, old_lower, 0, lp.sum + 1); + + get_row(lp, 0, OrigObj); + + pie = 2; + + if(lp.maximise != FALSE) + { + Zub = DEF_INFINITE; + Zlb = start_bound; + } + else + { + Zlb = -DEF_INFINITE; + Zub = start_bound; + } + status = RUNNING; + Step = 1; + OrigFeas = FALSE; + AnyFeas = FALSE; + citer = 0; + + for(i = 0 ; i < lp.nr_lagrange; i++) + lp.lambda[i] = 0; + + while(status == RUNNING) + { + citer++; + + for(i = 1; i <= lp.columns; i++) + { + ModObj[i] = OrigObj[i]; + for(j = 0; j < lp.nr_lagrange; j++) + { + if(lp.maximise != FALSE) + ModObj[i] -= lp.lambda[j] * lp.lag_row[j][i]; + else + ModObj[i] += lp.lambda[j] * lp.lag_row[j][i]; + } + } + for(i = 1; i <= lp.columns; i++) + { + set_mat(lp, 0, i, ModObj[i]); + /* fSystem.out.print(stderr, "%f ", ModObj[i]); */ + } + rhsmod = 0; + for(i = 0; i < lp.nr_lagrange; i++) + if(lp.maximise != FALSE) + rhsmod += lp.lambda[i] * lp.lag_rhs[i]; + else + rhsmod -= lp.lambda[i] * lp.lag_rhs[i]; + + if(verbose != FALSE) + { + System.out.println("Zub: " + Zub + " Zlb: " + Zlb + " Step: " + Step + + " pie: " + pie + " Feas " + OrigFeas); + for(i = 0; i < lp.nr_lagrange; i++) + System.out.println(i + " SubGrad " + SubGrad[i] + " lambda " + lp.lambda[i]); + } + + + result = solve(lp); + + same_basis = TRUE; + i = 1; + while(same_basis != FALSE && i < lp.rows) + { + same_basis = (old_bas[i] == lp.bas[i])? (short)1: (short)0; + i++; + } + i = 1; + while(same_basis != FALSE && i < lp.sum) + { + same_basis=(old_lower[i] == lp.lower[i])? (short)1:(short)0; + i++; + } + if(same_basis == FALSE) + { + System.arraycopy(lp.lower, 0, old_lower, 0, lp.sum+1); + System.arraycopy(lp.bas, 0, old_bas, 0, lp.rows+1); + pie *= 0.95; + } + + if(verbose != FALSE) + System.out.println("result: " + result + " same basis: " + same_basis); + + if(result == UNBOUNDED) + { + for(i = 1; i <= lp.columns; i++) + System.out.print(ModObj[i] + " "); + System.exit(FAIL); + } + + if(result == FAILURE) + status = FAILURE; + + if(result == INFEASIBLE) + status = INFEASIBLE; + + SqrsumSubGrad = 0; + for(i = 0; i < lp.nr_lagrange; i++) + { + SubGrad[i]= -lp.lag_rhs[i]; + for(j = 1; j <= lp.columns; j++) + SubGrad[i] += lp.best_solution[lp.rows + j] * lp.lag_row[i][j]; + SqrsumSubGrad += SubGrad[i] * SubGrad[i]; + } + + OrigFeas = TRUE; + for(i = 0; i < lp.nr_lagrange; i++) + if(lp.lag_con_type[i] != FALSE) + { + if(Math.abs(SubGrad[i]) > lp.epsb) + OrigFeas = FALSE; + } + else if(SubGrad[i] > lp.epsb) + OrigFeas = FALSE; + + if(OrigFeas != FALSE) + { + AnyFeas = TRUE; + Ztmp = 0; + for(i = 1; i <= lp.columns; i++) + Ztmp += lp.best_solution[lp.rows + i] * OrigObj[i]; + if((lp.maximise != FALSE) && (Ztmp > Zlb)) + { + Zlb = Ztmp; + for(i = 1; i <= lp.sum; i++) + BestFeasSol[i] = lp.best_solution[i]; + BestFeasSol[0] = Zlb; + if(verbose != FALSE) + System.out.println("Best feasible solution: " + Zlb); + } + else if(Ztmp < Zub) + { + Zub = Ztmp; + for(i = 1; i <= lp.sum; i++) + BestFeasSol[i] = lp.best_solution[i]; + BestFeasSol[0] = Zub; + if(verbose != FALSE) + System.out.println("Best feasible solution: " + Zub); + } + } + + if(lp.maximise != FALSE) + Zub = Math.min(Zub, rhsmod + lp.best_solution[0]); + else + Zlb = Math.max(Zlb, rhsmod + lp.best_solution[0]); + + if(Math.abs(Zub-Zlb) < (float)0.001) + { + status = OPTIMAL; + } + Step = (float)(pie * ((1.05*Zub) - Zlb) / SqrsumSubGrad); + + for(i = 0; i < lp.nr_lagrange; i++) + { + lp.lambda[i] += Step * SubGrad[i]; + if(lp.lag_con_type[i] == FALSE && lp.lambda[i] < 0) + lp.lambda[i] = 0; + } + + if(citer == num_iter && status==RUNNING) + if(AnyFeas != FALSE) + status = FEAS_FOUND; + else + status = NO_FEAS_FOUND; + } + + for(i = 0; i <= lp.sum; i++) + lp.best_solution[i] = BestFeasSol[i]; + + for(i = 1; i <= lp.columns; i++) + set_mat(lp, 0, i, OrigObj[i]); + + if(lp.maximise != FALSE) + lp.lag_bound = Zub; + else + lp.lag_bound = Zlb; + return(status); + } + + } // end of class solve + + + public static class Matrix + { + int row_nr; + float value; + public Matrix(int r, float v) { + row_nr = r; + value = v; + } + } + + + public static class Problem { + String lp_name; /* the name of the lp */ + + public short active; /*TRUE if the globals point to this structure*/ + public short verbose; /* ## Verbose flag */ + public short debug; /* ## Print B&B information */ + public short trace; /* ## Print information on pivot selection */ + public short anti_degen; /* ## Do perturbations */ + + public int rows; /* Nr of constraint rows in the problem */ + int rows_alloc; /* The allocated memory for Rows sized data */ + int columns; /* The number of columns (= variables) */ + int columns_alloc; + int sum; /* The size of the variables + the slacks */ + int sum_alloc; + + short names_used; /* Flag to indecate if names for rows and + columns are used */ + String[] row_name; /* rows_alloc+1 */ + int[] col_name; /* columns_alloc+1 */ + + /* Row[0] of the sparce matrix is the objective function */ + + int non_zeros; /* The number of elements in the sparce matrix*/ + int mat_alloc; /* The allocated size for matrix sized + structures */ + Matrix[] mat; /* mat_alloc :The sparse matrix */ + int[] col_end; /* columns_alloc+1 :Cend[i] is the index of the + first element after column i. + column[i] is stored in elements + col_end[i-1] to col_end[i]-1 */ + int[] col_no; /* mat_alloc :From Row 1 on, col_no contains the + column nr. of the + nonzero elements, row by row */ + short row_end_valid; /* true if row_end & col_no are valid */ + int[] row_end; /* rows_alloc+1 :row_end[i] is the index of the + first element in Colno after row i */ + float[] orig_rh; /* rows_alloc+1 :The RHS after scaling & sign + changing, but before `Bound transformation' */ + float[] rh; /* rows_alloc+1 :As orig_rh, but after Bound + transformation */ + float[] rhs; /* rows_alloc+1 :The RHS of the curent simplex + tableau */ + float[] orig_upbo; /* sum_alloc+1 :Bound before transformations */ + float[] orig_lowbo; /* " " */ + float[] upbo; /* " " :Upper bound after transformation + & B&B work*/ + float[] lowbo; /* " " :Lower bound after transformation + & B&B work */ + + short basis_valid; /* TRUE is the basis is still valid */ + int[] bas; /* rows_alloc+1 :The basis column list */ + short[] basis; /* sum_alloc+1 : basis[i] is TRUE if the column + is in the basis */ + short[] lower; /* " " :TRUE is the variable is at its + lower bound (or in the basis), it is FALSE + if the variable is at its upper bound */ + + short eta_valid; /* TRUE if current Eta structures are valid */ + int eta_alloc; /* The allocated memory for Eta */ + int eta_size; /* The number of Eta columns */ + int num_inv; /* The number of float pivots */ + int max_num_inv; /* ## The number of float pivots between + reinvertions */ + float[] eta_value; /* eta_alloc :The Structure containing the + values of Eta */ + int[] eta_row_nr; /* " " :The Structure containing the Row + indexes of Eta */ + int[] eta_col_end; /* rows_alloc + MaxNumInv : eta_col_end[i] is + the start index of the next Eta column */ + + short bb_rule; /* what rule for selecting B&B variables */ + + short break_at_int; /* TRUE if stop at first integer better than + break_value */ + float break_value; + + float obj_bound; /* ## Objective function bound for speedup of + B&B */ + int iter; /* The number of iterations in the simplex + solver (LP) */ + int total_iter; /* The total number of iterations (B&B) (ILP)*/ + int max_level; /* The Deepest B&B level of the last solution */ + int total_nodes; /* total number of nodes processed in b&b */ + public float[] solution; /* sum_alloc+1 :The Solution of the last LP, + 0 = The Optimal Value, + 1..rows The Slacks, + rows+1..sum The Variables */ + public float[] best_solution; /* " " :The Best 'Integer' Solution */ + float[] duals; /* rows_alloc+1 :The dual variables of the + last LP */ + + short maximise; /* TRUE if the goal is to maximise the + objective function */ + short floor_first; /* TRUE if B&B does floor bound first */ + short[] ch_sign; /* rows_alloc+1 :TRUE if the Row in the matrix + has changed sign + (a`x > b, x>=0) is translated to + s + -a`x = -b with x>=0, s>=0) */ + + short scaling_used; /* TRUE if scaling is used */ + short columns_scaled; /* TRUE is the columns are scaled too, Only use + if all variables are non-integer */ + float[] scale; /* sum_alloc+1 :0..Rows the scaling of the Rows, + Rows+1..Sum the scaling of the columns */ + + int nr_lagrange; /* Nr. of Langrangian relaxation constraints */ + float[][]lag_row; /* NumLagrange, columns+1:Pointer to pointer of + rows */ + float[] lag_rhs; /* NumLagrange :Pointer to pointer of Rhs */ + float[] lambda; /* NumLagrange :Lambda Values */ + short[] lag_con_type; /* NumLagrange :TRUE if constraint type EQ */ + float lag_bound; /* the lagrangian lower bound */ + + short valid; /* Has this lp pased the 'test' */ + float infinite; /* ## numercal stuff */ + float epsilon; /* ## */ + float epsb; /* ## */ + float epsd; /* ## */ + float epsel; /* ## */ + + + public Problem (int nrows, int ncolumns) + { + int i, nsum; + + nsum=nrows+ncolumns; + if(rows < 0 || columns < 0) + System.err.print("rows < 0 or columns < 0"); + + lp_name = new String("unnamed"); + active=FALSE; + verbose=FALSE; + debug=FALSE; + trace=FALSE; + + rows=nrows; + columns=ncolumns; + sum=nsum; + rows_alloc=rows; + columns_alloc=columns; + sum_alloc=sum; + names_used=FALSE; + + obj_bound=DEF_INFINITE; + infinite=DEF_INFINITE; + epsilon=DEF_EPSILON; + epsb=DEF_EPSB; + epsd=DEF_EPSD; + epsel=DEF_EPSEL; + non_zeros=0; + mat_alloc=1; + + mat = new Matrix[mat_alloc]; + for (i = 0; i < mat_alloc; i++) + mat[i] = new Matrix(0, 0); + + col_no = new int[mat_alloc]; + for (i = 0; i < mat_alloc; i++) + col_no[i] = 0; + + col_end = new int[columns + 1]; + for (i = 0; i < columns + 1; i++) + col_end[i] = 0; + + row_end = new int[rows + 1]; + for (i = 0; i < rows + 1; i++) + row_end[i] = 0; + + row_end_valid=FALSE; + + orig_rh = new float[rows + 1]; + for (i = 0; i < rows + 1; i++) + orig_rh[i] = 0; + + rh = new float[rows + 1]; + for (i = 0; i < rows + 1; i++) + rh[i] = 0; + + rhs = new float[rows + 1]; + for (i = 0; i < rows + 1; i++) + rhs[i] = 0; + + orig_upbo = new float[sum + 1]; + for(i = 0; i <= sum; i++) + orig_upbo[i]=infinite; + + upbo = new float[sum + 1]; + for (i = 0; i < sum + 1; i++) + upbo[i] = 0; + + orig_lowbo = new float[sum + 1]; + for (i = 0; i < sum + 1; i++) + orig_lowbo[i] = 0; + + lowbo = new float[sum + 1]; + for (i = 0; i < sum + 1; i++) + lowbo[i] = 0; + + basis_valid=TRUE; + + bas = new int[rows+1]; + for (i = 0; i <= rows; i++) + bas[i] = 0; + + basis = new short[sum + 1]; + for (i = 0; i <= sum; i++) + basis[i] = 0; + + lower = new short[sum + 1]; + for(i = 0; i <= rows; i++) + { + bas[i]=i; + basis[i]=TRUE; + } + for(i = rows + 1; i <= sum; i++) + basis[i]=FALSE; + for(i = 0 ; i <= sum; i++) + lower[i]=TRUE; + + eta_valid=TRUE; + eta_size=0; + eta_alloc=10000; + max_num_inv=DEFNUMINV; + + nr_lagrange=0; + + eta_value = new float[eta_alloc]; + for (i = 0; i < eta_alloc; i++) + eta_value[i] = 0; + + eta_row_nr = new int[eta_alloc]; + for (i = 0; i < eta_alloc; i++) + eta_row_nr[i] = 0; + + eta_col_end = new int[rows_alloc + max_num_inv]; + for (i = 0; i < rows_alloc + max_num_inv; i++) + eta_col_end[i] = 0; + + bb_rule=FIRST_NI; + break_at_int=FALSE; + break_value=0; + + iter=0; + total_iter=0; + + solution = new float[sum + 1]; + for (i = 0; i <= sum; i++) + solution[i] = 0; + + best_solution = new float[sum + 1]; + for (i = 0; i <= sum; i++) + best_solution[i] = 0; + + duals = new float[rows + 1]; + for (i = 0; i <= rows; i++) + duals[i] = 0; + + maximise = FALSE; + floor_first = TRUE; + + scaling_used = FALSE; + columns_scaled = FALSE; + + ch_sign = new short[rows + 1]; + for(i = 0; i <= rows; i++) + ch_sign[i] = FALSE; + + valid = FALSE; + } // end of constructor from row and column + + //*************************************** + // return the ith member of the best_solution[] + // + public float getBestSolution(int i) { + return best_solution[i]; + } + + //*************************************** + // get the number of rows + // + public int getRows() { + return rows; + } + + //*************************************** + // get the number of columns + // + public int getcolumns() { + return columns; + } + + } // end of class Problem + + private final static float round( float val, float eps) { return (Math.abs(val) < eps) ? 0 : val; } + +}