Java using Maven amateurishly OUTPUT
//*********************************************************************************************************************
//GraphFlow.java
package com.smkltd.app;
import java.util.BitSet; //to make a two bit two's complement array like object for totally unimodular matrix
public class GraphFlow
{
public int Rows;
public int Cols;
public int Iter;
public int Origin;
public int Destination;
public int smd;
public int MaxIter = 1000;
public int Objective = 0; // could be changed to floating point depending on cost and supply demand data types
public String minshipdescript;
public int[] ColumnHeading; // number the n edges 0,...,n-1
public int[] RowHeading; // number the m vertices n,...,n+m-1
public int[] SupplyDemand; // non negative cost vector, could be made floating point to be modified by simplex algo
public int[] B; //to be modified by simplex algo
public int[] TransportationCost; // non negative cost vector, could be made floating point to be modified by simplex algo
public int[] EdgeCapacity;
public int[] C; //to be modified by simplex algo
public int[] Transport; // units transported along each edge
public int[] SupplyChange; // column vector result of Incidence*Transport
public int[] InitSupplyState;
public int[] FinalSupplyState;
int[] walk;
public TwoBit Incidence; //two bit two's complement rep of incidence before simplex operation
private TwoBit Tableau;
public TwoBit MfTableau;
public GraphFlow(int nrows, int ncols)
{
this.Rows = nrows;
this.Cols = ncols;
this.ColumnHeading = new int[ncols];
this.RowHeading = new int[nrows];
for (int jj=0; jj<ncols; jj++)
{
this.ColumnHeading[jj] = jj;
}
for (int ii=0; ii<nrows; ii++)
{
this.RowHeading[ii] = ii+ncols;
}
this.TransportationCost = new int[ncols];
this.EdgeCapacity = new int[ncols];
this.Transport = new int[ncols];
this.SupplyDemand = new int[nrows];
this.SupplyChange = new int[nrows];
this.InitSupplyState = new int[nrows];
this.FinalSupplyState = new int[nrows];
this.Incidence = new TwoBit(nrows,ncols);
}
public void xchange(int pvtr, int pvtc)
{
// xchange operation preserves total unimodular property
int swap;
if (Tableau.Get(pvtr,pvtc)==0) {
throw new IllegalArgumentException("pivot element is 0 in xchange");
} else {
Objective = Objective-B[pvtr]*C[pvtc]/Tableau.Get(pvtr,pvtc);
for (int ii=0; ii<B.length; ii++)
{
if (ii!=pvtr) {
B[ii] = B[ii]-Tableau.Get(ii,pvtc)*B[pvtr]/Tableau.Get(pvtr,pvtc);
}
}
for(int jj=0; jj<C.length; jj++)
{
if (jj!=pvtc) {
C[jj] = C[jj]-Tableau.Get(pvtr,jj)*C[pvtc]/Tableau.Get(pvtr,pvtc);
}
}
C[pvtc] = C[pvtc]/Tableau.Get(pvtr,pvtc);
B[pvtr] = -B[pvtr]/Tableau.Get(pvtr,pvtc);
for(int ii=0; ii<Tableau.Rows; ii++)
{
for (int jj=0; jj<Tableau.Cols; jj++)
{
if ((ii!=pvtr) & (jj!=pvtc)) {
Tableau.Set(ii,jj,(Tableau.Get(ii,jj)-Tableau.Get(ii,pvtc)*Tableau.Get(pvtr,jj)/Tableau.Get(pvtr,pvtc)));
}
}
if (ii!=pvtr) {
Tableau.Set(ii,pvtc,(Tableau.Get(ii,pvtc)/Tableau.Get(pvtr,pvtc)));
}
}
for(int jj=0; jj<Tableau.Cols; jj++)
{
if (jj!=pvtc) {
Tableau.Set(pvtr,jj,(-Tableau.Get(pvtr,jj)/Tableau.Get(pvtr,pvtc)));
}
}
Tableau.Set(pvtr,pvtc,(1/Tableau.Get(pvtr,pvtc)));
swap = ColumnHeading[pvtc]; //relabel to indicate what row and column was xchanged
ColumnHeading[pvtc] = RowHeading[pvtr];
RowHeading[pvtr] = swap;
}
}
public void SimplexMin()
{
for (int kk=1; kk <= MaxIter; kk++)
{
int pr = -1; //arrays must be indexed from zero
int pc = -1; //ditto
int cr = 0;
int cr2 = 0;
for (int ii = 0; ii < Tableau.Rows; ii++)
{
if (B[ii] < 0) {
for (int jj = 0; jj < Tableau.Cols; jj++)
{
if (Tableau.Get(ii,jj) > 0) {
pr = ii;
pc = jj;
cr = C[pc]/Tableau.Get(pr,pc);
break;
}
}
}
if (pr != -1) {
break;
}
}
if ((pr != -1) & (pc != -1)) {
for (int jj = 0; jj < Tableau.Cols; jj++)
{
if (Tableau.Get(pr,jj) > 0) {
cr2 = C[jj]/Tableau.Get(pr,jj);
if (cr2 < cr) {
pc = jj;
cr = cr2;
}
}
}
xchange(pr,pc);
} else {
break;
}
Iter = kk;
}
}
public void MinShipSetUp()
{
int[] indeg = new int[Rows];
int[] outdeg = new int[Rows];
int[] deg = new int[Rows];
minshipdescript = "";
Tableau = new TwoBit(Rows,Cols);
B = new int[Tableau.Rows]; // B = -b
C = new int[Tableau.Cols];
ColumnHeading = new int[Tableau.Cols];
RowHeading = new int[Tableau.Rows];
for (int row = 0; row < Tableau.Rows; row++)
{
B[row] = SupplyDemand[row];
RowHeading[row] = row + Cols;
}
for (int col = 0; col < Tableau.Cols; col++)
{
C[col] = TransportationCost[col];
ColumnHeading[col] = col;
}
for (int ii = 0; ii < Rows; ii++)
{
indeg[ii] = 0;
outdeg[ii] = 0;
for (int jj = 0; jj < Cols; jj++)
{
if (Incidence.Get(ii,jj) > 0) {
indeg[ii]=indeg[ii]+1;
} else if (Incidence.Get(ii,jj) < 0) {
outdeg[ii]=outdeg[ii]+1;
}
}
deg[ii] = indeg[ii] + outdeg[ii];
if ((indeg[ii]==deg[ii]) & (B[ii]>0)) {
B[ii] = 0; //adjust for sink with supply
minshipdescript = "adjusted for sinks with supply: ";
}
if ((outdeg[ii]==deg[ii]) & (B[ii]<0)) {
B[ii] = 0; //adjust for source with demand
minshipdescript = "adjusted for sources with demand: ";
}
}
int supply = 0;
int demand = 0;
smd = 0; //supply minus demand
for (int ii = 0; ii < Rows; ii++)
{
if (B[ii] > 0) {
supply = supply + B[ii];
} else if (B[ii]<0) {
demand = demand - B[ii];
}
}
smd = supply - demand;
for (int row = 0; row < Tableau.Rows; row++)
{
if (smd<0) {
B[row] = -B[row];
}
for (int col = 0; col < Tableau.Cols; col++)
{
if (smd<0) {
Tableau.Set(row,col,-Incidence.Get(row,col));
} else {
Tableau.Set(row,col,Incidence.Get(row,col));
}
}
}
if (smd<0) {
minshipdescript = minshipdescript + "available supply less than demand";
} else if (smd>0) {
minshipdescript = minshipdescript + "available supply greater than demand";
} else {
minshipdescript = minshipdescript + "available supply equals demand";
}
Objective = 0;
Iter = 0;
Transport = new int[Cols];
SupplyChange = new int[Rows];
Change();
}//I hate curly braces!
public void MinShip()
{
SimplexMin();
//Transport = new int[Cols];
for (int ii=0; ii<Rows; ii++)
{
if (RowHeading[ii]<Tableau.Cols) {
Transport[RowHeading[ii]] = B[ii];
}
}
Change();
}
public void ShortestPath()
{
Tableau = new TwoBit(Rows,Cols);
B = new int[Tableau.Rows]; // B = -b
C = new int[Tableau.Cols];
ColumnHeading = new int[Tableau.Cols];
RowHeading = new int[Tableau.Rows];
//need to test if origin is a sink or destination is a source
for (int row = 0; row < Tableau.Rows; row++)
{
if (row==Origin) {
B[row] = 1;
SupplyDemand[row] = 1;
} else if (row==Destination) {
B[row] = -1;
SupplyDemand[row] = -1;
} else {
B[row] = 0;
SupplyDemand[row] = 0;
}
RowHeading[row] = row + Cols;
}
for (int col = 0; col < Tableau.Cols; col++)
{
C[col] = TransportationCost[col];
ColumnHeading[col] = col;
}
for (int row = 0; row < Tableau.Rows; row++)
{
for (int col = 0; col < Tableau.Cols; col++)
{
Tableau.Set(row,col,Incidence.Get(row,col));
}
}
MinShip();
//sequence the path
//howmany edges in path?
int lenth = 0;
for (int ii=0; ii<Cols; ii++)
{
if (Transport[ii]==1) {
lenth = lenth +1;
}
}
int lenth2 = 2*lenth+1;
walk = new int[lenth2];
int jj = 0;
walk[jj] = Origin;
while (jj<lenth2-1)
{
for (int ii=0; ii<Cols; ii++)
{
//System.out.println("1st "+jj+" "+walk[jj]+" "+ii+" "+Incidence.Get(walk[jj],ii));
if (Incidence.Get(walk[jj],ii)==-1) {
if (Transport[ii] == 1) {
jj=jj+1;
walk[jj] = ii;
break;
}
}
}
for (int ii=0; ii<Rows; ii++)
{
//System.out.println("2nd "+jj+" "+ii+" "+walk[jj]+" "+Incidence.Get(ii,walk[jj]));
if (Incidence.Get(ii,walk[jj])==1) {
jj=jj+1;
walk[jj] = ii;
break;
}
}
}
for (int ii=0; ii<walk.length; ii++) // add Cols to vertices which is how we're labeling
//of course we hven't made this adjustmenmt when input Origin and Destination
{
if (ii%2==0) {
walk[ii]=walk[ii]+Cols;
}
}
}
public void Change()
{
//just a matrix multiplication
for (int ii=0; ii<Rows; ii++)
{ SupplyChange[ii] = 0;
for (int jj = 0; jj<Cols; jj++)
{
SupplyChange[ii] = SupplyChange[ii] + Incidence.Get(ii,jj)*Transport[jj];
}
FinalSupplyState[ii] = InitSupplyState[ii] + SupplyChange[ii];
}
}
public void MaxFlowSetUp()
{
Tableau = new TwoBit(Cols, 2*(Rows-2)+Cols);
B = new int[Tableau.Rows];
C = new int[Tableau.Cols];
ColumnHeading = new int[Tableau.Cols];
RowHeading = new int[Tableau.Rows];
for (int ii=0; ii<Tableau.Rows; ii++)
{
RowHeading[ii] = ii;
}
int kk = -1;
for (int ii=0; ii<Rows; ii++)
{
if ((ii != Origin) & (ii != Destination))
{
kk=kk+1;
C[kk]=0;
ColumnHeading[kk]=ii+Tableau.Rows;
for (int jj = 0; jj<Tableau.Rows; jj++)
{
Tableau.Set(jj,kk,-Incidence.Get(ii,jj));
}
}
}
for (int ii=0; ii<Rows; ii++)
{
if ((ii != Origin) & (ii != Destination))
{
kk=kk+1;
C[kk]=0;
ColumnHeading[kk]=ii+Tableau.Rows+Rows;
for (int jj = 0; jj<Tableau.Rows; jj++)
{
Tableau.Set(jj,kk,Incidence.Get(ii,jj));
}
}
}
for (int ii=kk+1; ii<Tableau.Cols; ii++)
{
C[ii]=EdgeCapacity[ii-kk-1];
ColumnHeading[ii]=ii+kk+Rows;
for (int jj = 0; jj<Tableau.Rows; jj++)
{
if (jj==(ii-kk-1)) {
Tableau.Set(jj,ii,1);
} else {
Tableau.Set(jj,ii,0);
}
}
}
for (int jj = 0; jj <Tableau.Rows; jj++)
{
B[jj] = -Incidence.Get(Destination,jj);
}
Objective = 0;
Iter=0;
Transport = new int[Cols];
SupplyChange = new int[Rows];
Change();
}
public void MaxFlow()
{
SimplexMin();
for (int ii=0; ii<Tableau.Cols; ii++)
{
if (ColumnHeading[ii]<Tableau.Rows) {
Transport[ColumnHeading[ii]] = C[ii];
}
}
Change();
}
public void MaxFlowAtMinCostSetUp()
{
for (int ii=0; ii<Rows; ii++)
{
SupplyDemand[ii] = -SupplyChange[ii];
}
MinShipSetUp();
//MinShip();
}
public void PrintMinShipProgram()
{
System.out.println(" ");
System.out.println(" Transportation at minimum cost in a directed network linear program");
System.out.println(" ");
System.out.println(" x >= 0 Tableau");
System.out.println(" Ax - b >= 0 [A ][-b] pivot to make [-b] >= 0 and [c'] >= 0");
System.out.println(" c'x = min [c'][ 0]");
System.out.println(" ");
System.out.println(" "+minshipdescript);
System.out.println(" ");
System.out.println(" Simplex iterations = "+Iter);
System.out.printf(" V\\E ");
for (int col = 0; col < ColumnHeading.length; col++)
{
System.out.printf("%2d ", ColumnHeading[col]);
}
System.out.printf(" -b");
System.out.println("");
for (int row = 0; row < Rows; row++)
{
System.out.printf(" %3d ", RowHeading[row]);
for (int col = 0; col < Cols; col++)
{
System.out.printf("%2d ", Tableau.Get(row,col));
}
System.out.printf("%3d ", B[row]);
System.out.println("");
}
System.out.printf(" c' ");
for (int col = 0; col < C.length; col++) {
System.out.printf("%2d ", C[col]);
}
System.out.printf("%3d", Objective);
System.out.println(" ");
System.out.println(" ");
System.out.printf(" ship ");
for (int col = 0; col < Cols; col++)
{
System.out.printf("%2d ", Transport[col]);
}
System.out.println(" ");
System.out.println(" ");
System.out.printf(" VertexLabels ");
for(int ii=0; ii<Rows; ii++)
{
System.out.printf("%3d ", Cols+ii);
}
System.out.println("");
System.out.printf(" SupplyDemand ");
for(int ii=0; ii<Rows; ii++)
{
System.out.printf("%3d ", SupplyDemand[ii]);
}
System.out.println("");
System.out.printf(" SupplyChange ");
for(int ii=0; ii<Rows; ii++)
{
System.out.printf("%3d ", SupplyChange[ii]);
}
System.out.println("");
System.out.printf(" InitialSupply");
for(int ii=0; ii<Rows; ii++)
{
System.out.printf("%3d ", InitSupplyState[ii]);
}
System.out.println("");
System.out.printf(" FinalSupply ");
for(int ii=0; ii<Rows; ii++)
{
System.out.printf("%3d ", FinalSupplyState[ii]);
}
System.out.println(" ");
}
public void PrintMaxFlowProgram()
{
System.out.println(" ");
System.out.println(" Maximum flow from source to destination in a directed network linear program");
System.out.println(" ");
System.out.println(" z <= 0 Tableau");
System.out.println(" -A'z - c <= 0 [A ][-c] here c contains the edge capacities with remaining entries zero");
System.out.println(" b'z = min [b'][ 0] here b is destination row in Incidence matrix");
System.out.println(" ");
System.out.println(" setting up as dual program");
System.out.println(" ");
System.out.println(" x >= 0 Tableau");
System.out.println(" Ax - b >= 0 [A ][-b] pivot to make [-b] >= 0 and [c'] >= 0");
System.out.println(" c'x = min [c'][ 0]");
System.out.println(" ");
System.out.println(" Simplex iterations = "+Iter);
System.out.printf(" E\\V ");
for (int col = 0; col < Tableau.Cols; col++)
{
System.out.printf("%2d ", ColumnHeading[col]);
}
System.out.printf("%2d ",-(Destination+Cols));
System.out.println("");
for (int row = 0; row < Tableau.Rows; row++)
{
System.out.printf(" %3d ", RowHeading[row]);
for (int col = 0; col < Tableau.Cols; col++)
{
System.out.printf("%2d ", Tableau.Get(row,col));
}
System.out.printf("%3d ", B[row]);
System.out.println("");
}
System.out.printf(" c' ");
for (int col = 0; col < Tableau.Cols; col++) {
System.out.printf("%2d ", C[col]);
}
System.out.printf("%3d", Objective);
System.out.println(" ");
System.out.println(" ");
System.out.printf(" ship ");
for (int col = 0; col < Cols; col++)
{
System.out.printf("%2d ", Transport[col]);
}
System.out.println(" ");
System.out.println(" ");
System.out.printf(" VertexLabels ");
for(int ii=0; ii<Rows; ii++)
{
System.out.printf("%3d ", Cols+ii);
}
System.out.println("");
System.out.printf(" SupplyChange ");
for(int ii=0; ii<Rows; ii++)
{
System.out.printf("%3d ", SupplyChange[ii]);
}
System.out.println("");
System.out.printf(" InitialSupply");
for(int ii=0; ii<Rows; ii++)
{
System.out.printf("%3d ", InitSupplyState[ii]);
}
System.out.println("");
System.out.printf(" FinalSupply ");
for(int ii=0; ii<Rows; ii++)
{
System.out.printf("%3d ", FinalSupplyState[ii]);
}
System.out.println(" ");
}
public static void main( String[] args )
{
//incidence matrix of directed connected multigraph note: 2 bit two's complement integer matrix would save a lot of space
//the convention here is that +1 represents the arrow head and -1 represents the arrow tail
//so a non negative row represents a sink and a non positive row represents a source
//TwoBit.class allows two bit two's complement incidence representation
//the linear program we are solving to find transport through the network at minimum cost is:
//find c'x = min x >= 0 Tableau
//subject to Ax - b >= 0 Ax - b >= 0 [A ][-b] pivot to make [-b] >=0 and [c'] >= 0
//where x >= 0 c'x = min [c'][ 0]
//the dual program would be
//find b'y = max z <= 0 where z = -y
//subject to A'y - c <= 0 -A'z - c <= 0 [-A'][-c] pivot to make [-c] <= 0 and [b'] <= 0
//where y >= 0 b'z = min [ b'][ 0]
//note that Xr,c(-A') = -Xc,r(A)' for an xchange step, where r,c is the pivot
byte[][] MyIncidence = {
{1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{-1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, -1, 0, -1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, -1, 0, -1, 0, -1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, -1, 0, -1, 0, -1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, -1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, -1, 0, 0, 0, -1, 0, 1, 0, 0, 0},
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, -1, 0, 1, 1, 0},
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, -1, -1, 0, 1},
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, -1, -1}
};
int[] MyTransportationCost = {1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1}; // the cost of shipping one unit along edge
int[] MySupplyDemand = {-6,-14,-18,-6,-46,-14,24,-9,30,59};// an entry<0 represents a vertex with demand = abs(entry) and an entry>0 represents a vertex with supply = entry
int[] MyInitSupplyState = {15,2,23,0,7,9,36,1,55,72};
int[] MyEdgeCapacity = {13,17,28,7,23,17,9,19,10,18,8,5,21,6,27,16,5,29,19,7,18,19};
int[] MySupplyDemand2 = {6,-16,-20,-6,-46,-14,24,-9,30,51};
int[] MySupplyDemand3 = {-6,-14,-18,-6,-46,-14,40,-9,37,68};
GraphFlow MyGraphFlow = new GraphFlow(MyIncidence.length,MyIncidence[0].length);
for (int row = 0; row < MyGraphFlow.Rows; row++)
{
MyGraphFlow.SupplyDemand[row] = MySupplyDemand[row];
MyGraphFlow.InitSupplyState[row] = MyInitSupplyState[row];
for (int col = 0; col < MyGraphFlow.Cols; col++)
{
MyGraphFlow.Incidence.Set(row,col,MyIncidence[row][col]);
}
}
for (int col = 0; col < MyGraphFlow.Cols; col++)
{
MyGraphFlow.TransportationCost[col] = MyTransportationCost[col];
MyGraphFlow.EdgeCapacity[col] = MyEdgeCapacity[col];
}
MyGraphFlow.MinShipSetUp();
MyGraphFlow.PrintMinShipProgram();
MyGraphFlow.MinShip();
MyGraphFlow.PrintMinShipProgram();
MyGraphFlow.SupplyDemand = MySupplyDemand2;
MyGraphFlow.MinShipSetUp();
MyGraphFlow.PrintMinShipProgram();
MyGraphFlow.MinShip();
MyGraphFlow.PrintMinShipProgram();
MyGraphFlow.SupplyDemand = MySupplyDemand3;
MyGraphFlow.MinShipSetUp();
MyGraphFlow.PrintMinShipProgram();
MyGraphFlow.MinShip();
MyGraphFlow.PrintMinShipProgram();
MyGraphFlow.Origin=9;
MyGraphFlow.Destination=0;
MyGraphFlow.MaxFlowSetUp();
MyGraphFlow.PrintMaxFlowProgram();
MyGraphFlow.MaxFlow();
MyGraphFlow.PrintMaxFlowProgram();
MyGraphFlow.MaxFlowAtMinCostSetUp();
MyGraphFlow.PrintMinShipProgram();
MyGraphFlow.MinShip();
MyGraphFlow.PrintMinShipProgram();
MyGraphFlow.ShortestPath();
System.out.println(" ");
System.out.printf(" shortest path ");
for (int ii=0; ii<MyGraphFlow.walk.length; ii++)
{
System.out.printf("%2d",MyGraphFlow.walk[ii]);
System.out.printf(", ");
}
System.out.println(" ");
System.out.println(" ");
System.out.println(" ");
System.out.println(" Testing the TwoBit object rep of the incidence matrix:");
System.out.println(" ");
for (int row = 0; row < MyGraphFlow.Incidence.Rows; row++)
{
System.out.printf(" ");
for (int col = 0; col < MyGraphFlow.Incidence.Cols; col++)
{
System.out.printf("%2d ", MyGraphFlow.Incidence.Get(row,col));
}
System.out.println("");
}
System.out.println(" ");
System.out.println(" ");
MyGraphFlow.Incidence.LaplacianMatrix();
System.out.println(" Printing Laplacian Matrix for Network");
System.out.println(" ");
for (int row = 0; row < MyGraphFlow.Rows; row++)
{
System.out.printf(" ");
for (int col = 0; col < MyGraphFlow.Rows; col++)
{
System.out.printf("%2d ", MyGraphFlow.Incidence.LaplacianMatrix[row][col]);
}
System.out.println("");
}
}
}
//*******************************************************************************************************************************************
//TwoBit.java
package com.smkltd.app;
import java.util.BitSet; //to make a two bit two's complement array like object for totally unimodular matrix
public class TwoBit
{
public int Rows;
public int Cols;
public int[][] LaplacianMatrix;
public double[][] LaplacianMatrixR;
private BitSet[][] twobit;
public TwoBit(int nrows, int ncols)
{
this.Rows = nrows;
this.Cols = ncols;
this.twobit = new BitSet[2][nrows];
for (int ii=0;ii<nrows;ii++)
{
this.twobit[0][ii] = new BitSet(ncols);
this.twobit[1][ii] = new BitSet(ncols);
}
}
public int Get(int r, int c)
{
int ii;
if (r>=Rows) {
throw new IllegalArgumentException("row index out of bounds in TwoBit.Get");
} else if (c>=Cols) {
throw new IllegalArgumentException("col index out of bounds in TwoBit.Get");
} else {
if ((twobit[0][r].get(c))&(twobit[1][r].get(c))) {
ii = -1;
} else if (!(twobit[0][r].get(c)) & !(twobit[1][r].get(c))) {
ii=0;
} else if (!(twobit[0][r].get(c)) & (twobit[1][r].get(c))) {
ii=1;
} else {
ii=-2;
}
}
return ii;
}
public void Set(int r, int c, int pm1v0)
{
if ((pm1v0!=-1)&(pm1v0!=0)&(pm1v0!=1)) {
throw new IllegalArgumentException("trying to set value <> -1,0,1 in TwoBit.Set");
} else if (r>=Rows) {
throw new IllegalArgumentException("row index out of bounds in TwoBit.Set");
} else if (c>=Cols) {
throw new IllegalArgumentException("col index out of bounds in TwoBit.Set");
} else {
if (pm1v0 == -1) {
twobit[0][r].set(c);
twobit[1][r].set(c);
} else if (pm1v0 == 0) {
twobit[0][r].clear(c);
twobit[1][r].clear(c);
} else if (pm1v0 == 1) {
twobit[0][r].clear(c);
twobit[1][r].set(c);
}
}
}
public int Laplacian(int r, int c)
{
//also called the Kirchhoff matrix
int lp;
if (r>=Rows) {
throw new IllegalArgumentException("row index out of bounds in TwoBit.Laplacian");
} else if (c>=Rows) {
throw new IllegalArgumentException("col index out of bounds in TwoBit.Laplacian");
} else {
lp = 0;
for (int ii=0; ii<Cols; ii++)
{
lp = lp + Get(r,ii)*Get(c,ii);
}
}
return lp;
}
public int Laplacian(int r, int c, int[] K)
{
//also called the Kirchhoff matrix
int lp;
if (r>=Rows) {
throw new IllegalArgumentException("row index out of bounds in TwoBit.Laplacian");
} else if (c>=Rows) {
throw new IllegalArgumentException("col index out of bounds in TwoBit.Laplacian");
} else if (K.length!=Cols) {
throw new IllegalArgumentException("K[] has incorrect length in TwoBit.Laplacian");
} else {
lp = 0;
for (int ii=0; ii<Cols; ii++)
{
lp = lp + Get(r,ii)*Get(c,ii)*K[c];
}
}
return lp;
}
public double Laplacian(int r, int c, double[] R)
{
//also called the Kirchhoff matrix
double lp;
if (r>=Rows) {
throw new IllegalArgumentException("row index out of bounds in TwoBit.Laplacian");
} else if (c>=Rows) {
throw new IllegalArgumentException("col index out of bounds in TwoBit.Laplacian");
} else if (R.length!=Cols) {
throw new IllegalArgumentException("R[] has incorrect length in TwoBit.Laplacian");
} else {
lp = 0;
for (int ii=0; ii<Cols; ii++)
{
lp = lp + Get(r,ii)*Get(c,ii)*R[c];
}
}
return lp;
}
public void LaplacianMatrix()
{
LaplacianMatrix = new int[Rows][Rows];
for (int ii=0; ii<Rows; ii++)
{
for (int jj=0; jj<Rows; jj++)
{
LaplacianMatrix[ii][jj]=Laplacian(ii,jj);
}
}
}
public void LaplacianMatrixR(double[] CondR)
{
LaplacianMatrixR = new double[Rows][Rows];
for (int ii=0; ii<Rows; ii++)
{
for (int jj=0; jj<Rows; jj++)
{
LaplacianMatrixR[ii][jj]=Laplacian(ii,jj,CondR);
}
}
}
public void xchange(int pvtr, int pvtc)
{
if (pvtr>=Rows) {
throw new IllegalArgumentException("row index out of bounds in TwoBit.xchange()");
} else if (pvtc>=Cols) {
throw new IllegalArgumentException("col index out of bounds in TwoBit.xchange()");
} else if (Get(pvtr,pvtc)==0) {
throw new IllegalArgumentException("Illegal pivot element in TwoBit.xchange()");
} else {
for (int ii=0; ii<Rows; ii++) //need to catch a not totally unimodular condition though set() will bomb in this case
{
for (int jj=0; jj<Cols; jj++)
{
if ((ii!=pvtr) & (jj!=pvtc)) {
Set(ii,jj,(Get(ii,jj)-Get(ii,pvtc)*Get(pvtr,jj)/Get(pvtr,pvtc)));
}
}
if (ii!=pvtr) {
Set(ii,pvtc,(Get(ii,pvtc)/Get(pvtr,pvtc)));
}
}
for(int jj=0; jj<Cols; jj++)
{
if (jj!=pvtc) {
Set(pvtr,jj,(-Get(pvtr,jj)/Get(pvtr,pvtc)));
}
}
Set(pvtr,pvtc,(1/Get(pvtr,pvtc)));
}
}
public static void main( String[] args )
{
int v = 5;
int e;
e=(int)(v*(v-1)/2);
TwoBit mytwobit = new TwoBit(v,e);
//mytwobit.Create(v,e);
int kk=-1;
for (int ii=0; ii<v-1;ii++)
{
for (int jj=ii+1; jj<v;jj++)
{
kk=kk+1;
mytwobit.Set(ii,kk,1);
mytwobit.Set(jj,kk,-1);
}
}
System.out.println(" ");
for (int ii=0; ii<mytwobit.Rows;ii++)
{
for (int jj=0; jj<mytwobit.Cols;jj++)
{
System.out.printf("%2d ", mytwobit.Get(ii,jj));
}
System.out.println("");
}
}
}