package index;
//import ilog.opl_core.*;
//import ilog.opl_core.cppimpl.IloTupleSet;
import ilog.concert.IloException;
import ilog.concert.IloIntVarMap;
import ilog.concert.IloTuple;
import ilog.cplex.cppimpl.IloCplex;
import ilog.opl.IloOplDataSource;
import ilog.opl.IloOplErrorHandler;
import ilog.opl.IloOplException;
import ilog.opl.IloOplFactory;
import ilog.opl.IloOplModel;
import ilog.opl.IloOplModelDefinition;
import ilog.opl.IloOplModelSource;
import ilog.opl.IloOplSettings;

import java.io.FileOutputStream;
import java.io.IOException;
import java.io.OutputStreamWriter;
import java.io.Writer;
import java.util.*;

import org.jgrapht.generate.GraphGenerator;

import exp.RunMesasurements;

import graph.*;
import stats.*;
import misc.*;
import alg.*;

public class IntervalBounds {
	
	// --------------------UBs-----------------
	
	// UB without adjacency requirement
	public static double getUB(Graph g) {
		if(g.getn()==1) return 0;
		TreeMap<Double,Integer> pns= new TreeMap<Double, Integer>();
		TreeMap<Double,Integer> nes = new TreeMap<Double, Integer>();
		// get the distributions
		g.getPosNegDists(pns,nes);
		
		if(pns.isEmpty()) return 0.0;
		// connect the positive nodes using the best negative edges while the score improves
		double score = pns.lastKey();
		if (pns.get(pns.lastKey()) > 1) pns.put(pns.lastKey(), pns.get(pns.lastKey())-1);
		else pns.remove(pns.lastKey());
		
		if (pns.isEmpty() || nes.isEmpty()) return score;
		while (nes.lastKey() + pns.lastKey() > 0){
			score += nes.lastKey() + pns.lastKey();
			if (pns.get(pns.lastKey()) > 1) pns.put(pns.lastKey(), pns.get(pns.lastKey())-1);
			else pns.remove(pns.lastKey());
			if (nes.get(nes.lastKey()) > 1) nes.put(nes.lastKey(), nes.get(nes.lastKey())-1);
			else nes.remove(nes.lastKey());
			if(pns.size()==0 || nes.size()==0) break;
		}
		
		return score;
	}
	// UB with negative edges adjacent to clusters
	public static double getUBAdjNegative(Graph g) {
		if(g.getn()==1) return 0;
		ArrayList<ComparableCluster> cl = new ArrayList<ComparableCluster>();
		// get clusters
		g.getPosNegClusters(cl);
		
		// The UB is eithe one of the clusters on its own
		double ub = 0;
		for (ComparableCluster c: cl) { 
			if (ub < c.s) ub = c.s;
		}
		// Or the UB can be a combination of clusters
		// that may get positive profit when interconnected
		double score = 0.0;
		//double min_neg = 0.0;
		for (ComparableCluster c: cl) {
			if (c.s + c.n >= 0.0)
				score += (c.s + c.n);
		}
		if (ub < score) ub = score;
		// FIXED ERROR IN LOGIC: we should charge for all edges
		// otherwise we may obtain score lower than the actual
		//score -= min_neg;
		return ub;
	}
	
	// -------------------- NAIVE UBE + UB ------------------------------
	// Performs naive pruning of intervals of one and the same size
	private static double[] pruneSubIntervals(TimeGraph tg, int s, int e, 
											  double lb, 
											  HashMap<double[], Double> vcache,
											  boolean useMonotonicity) {
		assert(tg.gett() > e);
		double[] ubs = new double[e-s+1];
		Arrays.fill(ubs, Double.POSITIVE_INFINITY);
		
		double[] weights = new double[tg.getm()];
		Arrays.fill(weights, Double.POSITIVE_INFINITY);
		int[] pos = new int[tg.getm()];
		Arrays.fill(pos,-1);
		
		double domub,slack;
		int cntp;
		int cnt_ubruns = 0, cnt_all = 0;
		for (int sz=e-s+1;sz>0; sz--) {
//			// if this size is mostly pruned skip it
//			cntp = 0;
//			for (int j = s; j < e-sz+1; j++)
//				if (all_ubs[sz-1][j] > lb)
//					cntp++;
//			if (cntp < 2) continue;
			
			
			cnt_all++;
			tg.updateBestTimeLapsedEdges(s, e, sz, weights, pos);
			
			if (useMonotonicity) {
				domub = Double.POSITIVE_INFINITY;
				// try to reuse previous results
				for (double[] v:vcache.keySet()) {
					slack = Common.slack_dominates(v,weights);
					if ( slack +  vcache.get(v) <= lb) {
						domub = slack + vcache.get(v);
						break;
					}
				}
				if (domub > lb) {
					tg.we = weights;
					cnt_ubruns++;
					ubs[sz-1] = getUB(tg);
					if (ubs[sz-1] <= lb) 
						vcache.put(Arrays.copyOf(tg.we, tg.we.length), ubs[sz-1]);
				} else {
					ubs[sz-1] = domub;
				}
			} else {
				tg.we = weights;
				cnt_ubruns++;
				ubs[sz-1] = getUB(tg);
			}
		}
		//System.err.print("UBE runs:" + cnt_ubruns + " ( of " + cnt_all + " )\n");
		return ubs;
	}
	
	// OLD combination of NAIVE UBE + UB on the rest
	public static void getPruningUBAlign(TimeGraph tg, int s, int e, 
										int num_ube_levels, 
										double lb, 
										boolean useMonotonicity, 
										double[] pruned,
										double[][] ubs) {
		int totsz = e-s+1;
		int cnt_ube = 0;
		int sz = totsz, szsub, fnp, npruned;
		int from, to;
		pruned[0] = 0; pruned[1] = 0;
		
		// initialize the cache of dominanT vectors 
		double perc = 0;
		HashMap<double[], Double> vcache = new HashMap<double[],Double>();
		
		// Prune using exhaustive ub
		if (num_ube_levels > 0) {
			int inc = 0;
			while (sz > 0 && cnt_ube < num_ube_levels) {
				//System.err.print("Size ------------------->" + sz + "\n");
				int cnt = 0;
				do {
					from = s + cnt*inc;
					to = from + sz - 1;
					if (to > e) { // handle the last interval
						to = e;
						from = to-sz+1;
					}
	
					double[] p = pruneSubIntervals(tg, from, to, lb, vcache, false);
					//update the ub estimates of all subintervals
					for(szsub=0; szsub < p.length; szsub++) {
						for (int pos=from-s; pos <= to-s-szsub; pos++) {
							ubs[szsub][pos] = Math.min(p[szsub],ubs[szsub][pos]);
						}
					}
					//System.err.print("Size: " + sz + "\tPos: " + (from -s )+ "\tPruned: " + percPruned(ubs, lb) + "\n");
					//printPruningPerc(lb, ubs);
					cnt++;
				} while (to < e);
				
				
				// Find the first non-pruned size
				fnp = sz;
				while (fnp > 1) {
					npruned = 0; 
					for(int i=0;i<=totsz - fnp; i++) 
						if(ubs[fnp-1][i] <= lb) 
							npruned++;
					if ( npruned == totsz - fnp + 1 ) {
						fnp--;
					}else { 
						break;
					}
				}
				if (1 == fnp) {
					break;
				} else if (sz == fnp) {
					sz = (int)(sz/2.0);
					inc = sz /4;
				} else {
					if ((sz-fnp)/2 == 0 ) inc = 1;
					else inc = (sz-fnp)/2;
					sz = fnp + inc;
					//sz = first_sz_nonpruned;
				}
				cnt_ube++;
			}
			//printPruningPerc(99.0, lbs);
			//printPruningPerc(lb, ubs);
			perc = percPruned(ubs, lb);
			pruned[0] = perc;
			//ExtToolsDriver.plotMatrix(ubs,lb,1.0,ubs.length);
			printPruningPerc(lb, ubs);
			System.err.print("UBe" + (cnt_ube) + " pruned: \t" + ((int)(99.0*perc) +1) + "\n");
		}
		
		int vcache_size_ube = vcache.size();
		
		// use the UB
		int cnt_ub = 0, cnt_all = 0;
		for(sz = 0; sz < ubs.length; sz++) {
			for (int j = 0; j < ubs[sz].length-sz; j++) {
				if (ubs[sz][j] > lb) {
					cnt_all++;
					tg.aggregateByTime(j+s, j+s+sz);
					
					if (useMonotonicity) {
						double domub = Double.POSITIVE_INFINITY;
						double slack;
						
						// try to reuse previous results
						for (double[] v:vcache.keySet()) {
							slack = Common.slack_dominates(v,tg.we);
							if ( slack +  vcache.get(v) <= lb) {
								domub = slack + vcache.get(v);
								break;
							}
						}
						if (domub > lb) {
							ubs[sz][j] = getUB(tg); 
							cnt_ub++;
							if (ubs[sz][j] <= lb) 
								vcache.put(Arrays.copyOf(tg.we, tg.we.length), ubs[sz][j]);
						} else {
							ubs[sz][j] = domub;
						}
					} else { // no monotonicity
						ubs[sz][j] = getUB(tg); 
						cnt_ub++;
					}
				}
			}
		}
		//System.err.print("#UBs: " + cnt_ub + "\t dominated:" + (cnt_all - cnt_ub) + "\t vector cache:" +vcache.size() +  "( " + vcache_size_ube + "from UBE) \n");
		//printPruningPerc(lb, ubs);
		perc = percPruned(ubs, lb);
		pruned[1] = perc;
		System.err.print("UB pruning: " + ((int)(99.0*perc) +1) + "\n");
//		System.err.print("pruned ub asub:" + pruned_ub_al_ub.size() + "\n");
//		ExtToolsDriver.plotIntersvals(pruned_ub_al_ub, e, "ub");
	}
	
	
	// UB Pruning with/without an index
	public static int pruneUB(ITimeGraph tg, int s, int e, double lb,  
										double[][] ubs,
										IntervalIndexExactCeiling ind,
										boolean useMonotonicity,
										double alpha,
										TreeSet<ComparableInterval> notpruned, 
										RunMesasurements rm) {
		int sz = e-s+1;
		notpruned.clear();
		// initialize the cache of dominanT vectors 
		HashMap<double[], Double> vcache = null;
		if(useMonotonicity) vcache = new HashMap<double[],Double>();
		boolean grouped = false;
		int y=-1, xfrom=-1, xto = -1; 
		
		// use the UB
		int cnt_ub = 0, cnt_all = 0, cnt_lb = 0;
		for(sz = 0; sz < ubs.length; sz++) {
			for (int j = 0; j < ubs[sz].length-sz; j++) {
				 // check for both previous and curr ub since it might have been refined by groupoing
				if (ubs[sz][j] > lb) {
					cnt_all++;
					grouped = false; y=-1; xfrom=-1; xto = -1;
					if (ind != null) {
						// grouping lines
						int ssz = sz;
						for ( ssz = sz; ssz <= sz/alpha && ssz < ubs.length - j; ssz++) {
							if(ubs[ssz][j] <= lb) break;
						}
						if (ssz-1 > sz) {
							// group intervals
							grouped = true; y=j; xfrom=sz; xto=ssz-1;
							ind.getCeiling(y,xfrom+y, y+xto, tg.getGraph().we);
						} else {
							ind.getAggregate(j,j+sz,tg.getGraph().we);
						}
					} else {
						// no index so just plain old aggregate
						tg.aggregateByTime(j+s, j+s+sz);
					}
					if (useMonotonicity) {
						double domub = Double.POSITIVE_INFINITY;
						double slack;
						
						// try to reuse previous results
						for (double[] v:vcache.keySet()) {
							slack = Common.slack_dominates(v,tg.getGraph().we);
							if ( slack +  vcache.get(v) <= lb) {
								domub = slack + vcache.get(v);
								break;
							}
						}
						if (domub > lb) { // could not hit the dominant set
							if (! grouped) {
								//ubs[sz][j] = getUB(tg); 
								if (ind != null) {
									ubs[sz][j] = Common.sumPos(tg.getGraph().we) / 2.0;
								}
								if (ubs[sz][j] > lb) {
									ubs[sz][j] = getUBAdjNegative(tg.getGraph());
									cnt_ub++;
									if (ubs[sz][j] > lb) {
										cnt_lb++;
										notpruned.add(new ComparableInterval(j, sz, ubs[sz][j], tg.getGraph().we));
										rm.veri_p++;
									}else { // pruned
										rm.ub_p++;
										vcache.put(Arrays.copyOf(tg.getGraph().we, tg.getGraph().we.length), ubs[sz][j]);
									}
								} else {
									rm.sop_p++;
								}
							} else { // grouped
								double v = getUBAdjNegative(tg.getGraph()); cnt_ub++;
								if (v <= lb) { // ceiling pruned
									vcache.put(Arrays.copyOf(tg.getGraph().we, tg.getGraph().we.length), v);
									// update ubs
									for (int x = xfrom; x <= xto; x++) {
										ubs[x][y] = v;
										rm.ubcoarse_p++;
									}
								} else { // did not prune the ceiling
									// break group 
									for (int x = xfrom; x <= xto; x++) {
										ind.getAggregate(y,y+x,tg.getGraph().we);
										ubs[x][y] = Common.sumPos(tg.getGraph().we)  / 2.0;
										if ( ubs[x][y] > lb) {	
											ubs[x][y] = getUBAdjNegative(tg.getGraph());
											cnt_ub++;
											if (ubs[x][y] > lb) {
												cnt_lb++;
												notpruned.add(new ComparableInterval(y, x, ubs[x][y], tg.getGraph().we));
												rm.veri_p++;
											} else {
												rm.ub_p++;
											}
										} else {
											rm.sop_p++;
										}
										if (ubs[x][y] <= lb) {
											vcache.put(Arrays.copyOf(tg.getGraph().we, tg.getGraph().we.length), ubs[x][y]);
										}
									}
								}
							}
						} else { // hit the dominant set of pruned
							if (! grouped) {
								ubs[sz][j] = domub;
							} else {
								for (int x = xfrom; x <= xto; x++) ubs[x][y] = domub;
							}
						}
					} else { // no monotonicity
						if (! grouped) {
							//ubs[sz][j] = getUB(tg); 
							if (ind != null) {
								ubs[sz][j] = Common.sumPos(tg.getGraph().we)  / 2.0;
							}
							if (ubs[sz][j] > lb) {
								ubs[sz][j] = getUBAdjNegative(tg.getGraph()); cnt_ub++;
								if (ubs[sz][j] > lb) {
									cnt_lb++;
									notpruned.add(new ComparableInterval(j, sz, ubs[sz][j], tg.getGraph().we));
									rm.veri_p++;
								} else {
									rm.ub_p++;
								}
							} else {
								rm.sop_p++;
							}
						} else { // grouped
							// TODO(remove)
							assert(tg.getGraph().isUndirected());
							double v = getUBAdjNegative(tg.getGraph()); cnt_ub++;
							if (v <= lb) { // ceiling pruned
								// update ubs
								for (int x = xfrom; x <= xto; x++) { 
									ubs[x][y] = v;
								    rm.ubcoarse_p++;
								}
							} else { // did not prune the ceiling
								// break group 
								for (int x = xfrom; x <= xto; x++) {
									ind.getAggregate(y,y+x,tg.getGraph().we);
									// try sum of positives
									ubs[x][y] = Common.sumPos(tg.getGraph().we) / 2.0;
									if ( ubs[x][y] > lb) {
										ubs[x][y] = getUBAdjNegative(tg.getGraph());
										if (ubs[x][y] > lb) {
											cnt_lb++;
											notpruned.add(new ComparableInterval(y, x, ubs[x][y], tg.getGraph().we));
											rm.veri_p++;
										} else {
											rm.ub_p++;
										}
										cnt_ub++;
									} else {
										rm.sop_p++;
									}
								}
							}
						}
					}
				}
			}
		}
		return cnt_ub;
	}
	
	// -------------------- Greedy Cluster formation --------------------
	
	// performs a neighborhood base expansion of interval groups
	public static void growUBEPrunable(TimeGraph tg, int s, int e, double lb,
			double[][] posit, double[][] ubs, double[][] partit) {
		
		// coefficient of the STDEV
		double ss = 0.1 ;
		// positive mass that is not included in UB because it is not combinable
		double P = 0.3; 
		// the fraction of considered negative edges in UB
		double N = 0;
		
		double[] curr = new double[tg.getm()];
		int cid = 1;
		for(int sz = partit.length-1; sz >= 0; sz--) {
			for (int j = 0; j < partit.length-sz; j++) {
				if (ubs[sz][j] > lb) partit[sz][j] = -30;
				else if (posit[sz][j] > lb && partit[sz][j] == Double.POSITIVE_INFINITY) {
					if (sz < j) {
						System.currentTimeMillis();
					}
					Arrays.fill(curr,Double.NEGATIVE_INFINITY);
					ArrayList<String> q = new ArrayList<String>();
					ArrayList<String> added = new ArrayList<String>();
					q.add(sz + "_" + j);
					Cluster c = new Cluster(0, 0, 0);
					int cnt=0;
					while (!q.isEmpty()) {
						int ssz = Integer.parseInt(q.get(0).split("_")[0]);
						int jj = Integer.parseInt(q.get(0).split("_")[1]);
						tg.aggregateByTime(s+jj, jj+s+ssz);
						for (int i =0 ; i < tg.getm(); i++) {
							tg.we[i] = Math.max(tg.we[i], curr[i]);
						}
						
						if (getUB(tg) <= lb) {
							added.add(ssz+"_"+jj);
							cnt++;
							curr = Arrays.copyOf(tg.we, tg.we.length);
							addNeighbors(ssz,jj,lb,posit,ubs,partit,q,added);
						} 
						q.remove(0);
					}
					if (cnt > 0) {
						for (String key: added) 
							partit[Integer.parseInt(key.split("_")[0])]
							      [Integer.parseInt(key.split("_")[1])] = cid;
						cid++;
					}
					else {
						for (String key: added) 
							partit[Integer.parseInt(key.split("_")[0])]
							      [Integer.parseInt(key.split("_")[1])] = -20;
					}
				}
			}
		}
		permuteClusterIds(partit, cid);
		
		ExtToolsDriver.plotMatrix(partit,1.0,-30,partit.length);
	}
	
	// performs a neighborhood base expansion of interval groups
	// returns the number of computed upper bounds
	public static int prunePartitinons(TimeGraph tg, int s, int e, double lb,
			double[][] posit,  IntervalIndexExactCeiling ind,
			double[][] ubs, double[][] partit,
			double ss, double P) {
		
		// coefficient of the STDEV
		//double ss = 0.2;
		// positive mass that is not included in UB because it is not combinable
		//double P = 0.98; 
		// the fraction of considered negative edges in UB
		double N = 0;
		
		double[] curr = new double[tg.getm()];
		int num_ub = 0;
		int cid = 1;
		
		for(int sz = partit.length-1; sz >= 0; sz--) {
			for (int j = 0; j < partit.length-sz; j++) {
				//if(posit[sz][j] > lb && ubs[sz][j] <= lb) {
				if(posit[sz][j] > lb ) {
					if (partit[sz][j] == Double.POSITIVE_INFINITY) {
						Arrays.fill(curr,Double.NEGATIVE_INFINITY);
						ArrayList<String> q = new ArrayList<String>();
						ArrayList<String> added = new ArrayList<String>();
						q.add(sz + "_" + j);
						Cluster c = new Cluster(0, 0, 0);
						int cnt=0;
						while (!q.isEmpty()) {
							int ssz = Integer.parseInt(q.get(0).split("_")[0]);
							int jj = Integer.parseInt(q.get(0).split("_")[1]);
							
							ind.getAggregate(jj,jj+ssz, tg.we);
							//tg.aggregateByTime(s+jj, jj+s+ssz);
							double sumpos=0, sumneg=0;
							for (int i =0 ; i < tg.getm(); i++) {
								tg.we[i] = Math.max(tg.we[i], curr[i]);
								if (tg.we[i] >0) sumpos+= tg.we[i];
								if (tg.we[i] <0) sumneg+= tg.we[i];
							}
							
							Cluster cn = new Cluster(c.m_size + 1, sumpos, sumneg);
							if (1.0 == cn.m_size ||
								cn.cost(lb, P, N, ss) < (c.cost(lb, P, N, ss) + 1)) {
								c = cn;
								added.add(ssz+"_"+jj);
								cnt++;
								curr = Arrays.copyOf(tg.we, tg.we.length);
								addNeighbors(ssz,jj,lb,posit,ubs,partit,q,added);
							} 
							q.remove(0);
						}
						if (cnt > 0) {
							num_ub++;
							for(int i = 0; i < tg.getm(); i++)tg.we[i] = curr[i];
							double ub = getUBAdjNegative(tg);
							if (ub > lb) {
								num_ub += added.size();
								for (String key: added) {
									int ssz = Integer.parseInt(key.split("_")[0]);
									int jj = Integer.parseInt(key.split("_")[1]);
									ind.getAggregate(jj, jj+ssz, tg.we);
									ubs[ssz][jj] = getUBAdjNegative(tg);
								}
							}
							if (true) {
								// book keeping - remove in production
								for (String key: added) { 
									int ssz = Integer.parseInt(key.split("_")[0]);
									int jj = Integer.parseInt(key.split("_")[1]);
									partit[ssz][jj] = cid;
								    if (ub <= lb) ubs[ssz][jj] = ub;  
								}
								
							}
							cid++;
						}
						else {
							for (String key: added) 
								partit[Integer.parseInt(key.split("_")[0])]
								      [Integer.parseInt(key.split("_")[1])] = -20;
						}
					}
				}
//				else if (ubs[sz][j] > lb) 
//					partit[sz][j] = -30;
			}
		}
		
		//permuteClusterIds(partit, cid);
		return num_ub;
	}

	private static void addNeighbors(int sz, int j, double lb,
			double[][] posit, double[][] ubs, double[][] partit,
			ArrayList<String> q, ArrayList<String> added) {
		for (int pos = j-1; pos <= j+1; pos++) {
			if(pos>=0 && pos < partit.length - sz && 
			  (!q.contains(sz + "_" + pos)) && (!added.contains(sz + "_" + pos))){
				if(posit[sz][pos] > lb && 
				   //ubs[sz][pos] <= lb &&
				   partit[sz][pos] == Double.POSITIVE_INFINITY)
					q.add(sz + "_" + pos);
			}
		}
		if (sz-1>= 0){ 
			for (int pos = j-1; pos <= j+3; pos++) {
				if(pos>=0 && (pos < partit.length - sz+1) && 
				  (!q.contains((sz-1) + "_" + pos)) && (!added.contains((sz -1) + "_" + pos))){
					if(posit[sz-1][pos] > lb && 
					   //ubs[sz-1][pos] <= lb && 
					   partit[sz-1][pos] == Double.POSITIVE_INFINITY)
						q.add((sz-1) + "_" + pos);
				}
			}
		}
	}

	// --------------------------- STATS and UTILS -------------------------
	public static double percPruned(double[][] ubs, double lb) {
		int pr=0, total=0;
		for(int i=0; i < ubs.length; i++) {
			for (int j = 0; j < ubs[i].length - i; j++) {
				if (ubs[i][j] <= lb) pr++;
				total++;
			}
		}
		return 1.0*pr/total;
	}
	
	public static int numPruned(double[][] ubs, double lb) {
		int pr=0;
		for(int i=0; i < ubs.length; i++) {
			for (int j = 0; j < ubs[i].length - i; j++) {
				if (ubs[i][j] <= lb) pr++;
			}
		}
		return pr;
	}
	
	public static int biggestUnprunedSize(double[][] ubs, double lb) {
		int pr=0, total=0;
		for(int i=ubs.length-1; i > -1; i--) {
			for (int j = 0; j < ubs[i].length - i; j++) {
				if (ubs[i][j] > lb) 
					return i;
			}
		}
		return -1;
	}
	
	private static boolean expandable(int sz, int pos, double[][] lbs) {
		double lb = 0;
		if (sz>0 ) {
			// sz-1, pos-1
			if (pos-1 >= 0) {
				if (lbs[sz-1][pos-1] == Double.NEGATIVE_INFINITY){
					return true;
				}
			}
			// sz-1, pos
			if (lbs[sz-1][pos] == Double.NEGATIVE_INFINITY){
				return true;
			}
			// sz-1, pos+1
			if (pos+1 < lbs.length-(sz-1)) { 
				if (lbs[sz-1][pos+1] == Double.NEGATIVE_INFINITY){
					return true;
				}
			}
		}
		// same size
		// sz, pos-1
		if (pos-1 >= 0) {
			if (lbs[sz][pos-1] == Double.NEGATIVE_INFINITY){
				return true;
			}
		}
		// sz, pos+1
		if (pos+1 < lbs.length-sz){ 
			if (lbs[sz][pos+1] == Double.NEGATIVE_INFINITY){
				return true;
			}
		}
		// size + 1
		if (sz+1 < lbs.length) {
			// sz+1, pos-1
			if (pos-1 >= 0) {
				if (lbs[sz+1][pos-1] == Double.NEGATIVE_INFINITY){
					return true;
				}
			}
			// sz+1, pos
			if (pos < lbs.length-(sz+1) ) {
				if (lbs[sz+1][pos] == Double.NEGATIVE_INFINITY){
					return true;
				}
			}
			// sz+1, pos+1
			if (pos+1 < lbs.length-(sz+1) ) { 
				if (lbs[sz+1][pos+1] == Double.NEGATIVE_INFINITY){
					return true;
				}
			}
		}
		return false;
	}
	
	private static void printPruningPerc(double norm,
			double[][] ubs) {
		int val;
		for(int sz=0; sz < ubs.length; sz++) {
			System.err.print((sz + 1)+ ":\t");
			for (int j=0; j < ubs.length - sz; j++) {
					val = (int)(99.0*ubs[sz][j]/norm);
					if (val < 0) val = 0;
					System.err.print( ((val<10)?("  " + val):((val<100)?(" " + val):(val))) + " ");
			}
			System.err.print("\n");
		}
	}
	
	public static TreeMap<Double,Integer> getMaxSTClusterScores(Graph g) {
		if(g.getn()==1) return null;
		TreeMap<Double,Integer> pns= new TreeMap<Double, Integer>();
		TreeMap<Double,Integer> nes = new TreeMap<Double, Integer>();
		// get the distributions
		g.getPosNegDists(pns,nes);
		return pns;
	}
	
	public static void permuteClusterIds(double[][] partitioning,
			int num_clusters) {
		Random r = new Random();
		int[] m = new int[num_clusters];
		for (int i=0; i < m.length; i++) m[i] = i;
		for (int i=0; i < 100*m.length; i++) {
			int i1 = r.nextInt(m.length-1)+1;
			int i2 = r.nextInt(m.length-1)+1;
			int tmp = m[i1];
			m[i1] = m[i2];
			m[i2] = tmp;
		}
		for (int i=0; i < partitioning.length; i++)
			for(int j=0; j < partitioning.length; j++)
				if (partitioning[i][j]>0 && partitioning[i][j] < num_clusters)
					partitioning[i][j] = m[(int)partitioning[i][j]];
	}
	
	// -------------------------- LB --------------------------------------
	
	// use this for edge-weighted
	public static Tree topDownLBold(Graph g, BitSet edges) {
		if(g.getn()==1) return null;
		
		long start = System.currentTimeMillis();
		Tree t = ST.computeMaxST(g);
		long mstt = System.currentTimeMillis() - start;
		
		start = System.currentTimeMillis();
		// distribute positive edge mass to adjacent nodes
		// the doubled edge score here is skipped
		for (int i=0; i<g.getn();i++) {
			double sum=0;
			for (int j=g.ind[i]; j < g.ind[i+1]; j++) {
				if (g.we[j] > 0 && g.endv[j] > i) // ONLY IN ONE DIRECTION
					sum+=g.we[j];
			}
			t.wn[t.getNodeIndex(g.names[i])] = sum;
		}
		// keep only negative edges in the tree, 
		// take their absolute value
		for (int i=0; i < t.getm(); i++) {
			t.we[i] = Math.abs(Math.min(0.0, t.we[i]));
		}
		long transformt = System.currentTimeMillis() - start;
		
		start = System.currentTimeMillis();
		// Tree PCST
		Tree ot = PCST.computeOptTreeNW(t);
		long pcstt = System.currentTimeMillis() - start;
		
		start = System.currentTimeMillis();
		if (null != edges) {
			// first add connector edges from the spanning tree
			ArrayList<Integer> el = g.getSubgraphEdgesStructureOnly(ot);
			for (int e:el) edges.set(e);
			// now add all positive edges adjacent to included nodes
			HashSet<String> nds = new HashSet<String>();
			for(String n:ot.names) nds.add(n);
			for (int i=0; i < g.getn(); i++){
				if(!nds.contains(g.names[i])) continue;
				for (int j=g.ind[i]; j < g.ind[i+1]; j++) {
					if (g.we[j] >=0) 
						edges.set(j);
				}
			}
		}
		long reconstructt = System.currentTimeMillis() - start;
		
		
		//ExtToolsDriver.showInGoogleEarth(ot.names);
		start = System.currentTimeMillis();
		double s =  ot.getScore();
		long scoret = System.currentTimeMillis() - start;
		
		
//		System.err.print("" + s + "\tmst: " + mstt +  "ms\ttrans: " + transformt + 
//						 "ms\tpcst: " + pcstt + "ms\treconstr: " + reconstructt + 
//						 "ms\tscore: " + scoret + "ms\n");
		return ot;
	}
	
	
    public static void MaximumScoreSubgraph(Graph g, BitSet edgs) throws IOException {
        String datFileName = "graph.dat";
        Tree subgraph = null;
        IO.writeILPInput(g,datFileName);
        
        // Run ILP
        int status = 127;
        try {
        	IloOplFactory.setDebugMode(false);
        	IloOplFactory oplF = new IloOplFactory();
        	IloOplErrorHandler errHandler = oplF.createOplErrorHandler();
        	IloOplModelSource modelSource = oplF.createOplModelSource("MS.mod");
        	IloOplSettings settings = oplF.createOplSettings(errHandler);
        	IloOplModelDefinition def = oplF.createOplModelDefinition(modelSource,settings);
        	ilog.cplex.IloCplex cplex = oplF.createCplex();
        	cplex.setOut(null);
        	IloOplModel opl = oplF.createOplModel(def, cplex);
        	IloOplDataSource dataSource = oplF.createOplDataSource(datFileName);
        	opl.addDataSource(dataSource);
        	opl.generate();
        	
        	if (cplex.solve())
        	{
        		ilog.concert.IloTupleSet edges = opl.getElement("edges").asTupleSet();
        		IloIntVarMap y = opl.getElement("y").asIntVarMap();
        		for (Iterator it = edges.iterator(); it.hasNext(); ) {
        			IloTuple tuple = (IloTuple) it.next();
        			int u = tuple.getIntValue("u");
        			int v = tuple.getIntValue("v");
        			if ((opl.getCplex().getValue(y.get(tuple)) == 1) && (u <= v))
        			{
        				edgs.set(g.getEdgeIndex(u-1, v-1));
        				edgs.set(g.getEdgeIndex(v-1, u-1));
        			}
        		}
        	}
        	else
        	{
        		System.out.println("No solution!");
        	}
        	oplF.end();
        	status = 0;
        }
        catch (IloOplException e) {
        	e.printStackTrace();
        	status = 2;
        }
        catch (IloException e) {
        	e.printStackTrace();
        	status = 3;
        }
        catch (Exception e) {
        	e.printStackTrace();
        	status = 4;
        }
    }


	
	// use this for edge-weighted fast version
	public static void topDownLBFast(Graph g, BitSet edges) {
		if(g.getn()==1) return;
		Arrays.fill(g.wn, 0);
		
		long start = System.currentTimeMillis();
		TreeSet<ComparableEdge> edg = new TreeSet<ComparableEdge>();
		// Order edges by increasing cost
		for (int i = 0; i < g.getn(); i++) {
			for( int j = g.ind[i]; j < g.ind[i+1]; j++ ) {
				if (i < g.endv[j]) {
					edg.add(new ComparableEdge(g.we[j], j, g.getReverseEdgeIndex(j), i, g.endv[j]));
				}
			}
		}
		long ordert = System.currentTimeMillis() - start;
		
		start = System.currentTimeMillis();
		Tree t = ST.computeMaxSTFast(g, edg);t.buildN2I();
		long mstt = System.currentTimeMillis() - start;
		
		
		// distribute positive edge mass to adjacent nodes
		// the doubled edge score here is skipped
		for (ComparableEdge ce : edg.descendingSet()) {
			if (ce.w > 0) {
				if(t.n2i.containsKey(g.names[ce.v1])){
					t.wn[t.n2i.get(g.names[ce.v1])] += ce.w;
				}else if (t.n2i.containsKey(g.names[ce.v2])){
					t.wn[t.n2i.get(g.names[ce.v2])] += ce.w;
				}
			} else {
				break;
			}
		}
	
		// keep only negative edges in the tree, 
		// take their absolute value
		for (int i=0; i < t.getm(); i++) {
			if(t.we[i] >= 0) t.we[i] = 0;
			else t.we[i] = (-1)*t.we[i];
		}

		start = System.currentTimeMillis();
		// Tree PCST
		Tree ot = PCST.computeOptTreeNW(t);
		
		if (null != edges) {
			// first add connector edges from the spanning tree
			ArrayList<Integer> el = g.getSubgraphEdgesStructureOnly(ot);
			for (int e:el) edges.set(e);
			// now add all positive edges adjacent to included nodes
			HashSet<String> nds = new HashSet<String>();
			for(String n:ot.names) nds.add(n);
			for (int i=0; i < g.getn(); i++){
				if(!nds.contains(g.names[i])) continue;
				for (int j=g.ind[i]; j < g.ind[i+1]; j++) {
					if (g.we[j] >=0) { 
						edges.set(j);
						edges.set(g.getReverseEdgeIndex(j));
					}
				}
			}
		}
		// TODO(remove)
		for (int j = 0; j < g.getm(); j++) {
			assert(edges.get(j) == edges.get(g.getReverseEdgeIndex(j)));
		}
		
		//ExtToolsDriver.showInGoogleEarth(ot.names);
		return;
	}
	
	// topK LB in interpolated matrix
	public static double getLBtopKInMatrix(ITimeGraph tg, int s, int e, 
			   SplineInterpolation si, int size, int k,
			   IntervalIndexExactCeiling ind, double[][] lbs,
			   boolean exact) throws IOException{
		
		long start = System.currentTimeMillis();
		TreeSet<ComparableInterval> topk = getTopKIntervals(si.m, k);
		long topkt = System.currentTimeMillis() - start;
		
		start = System.currentTimeMillis();
		updateTopKIntervals(si, topk);
		long updatet = System.currentTimeMillis() - start;
		
		start = System.currentTimeMillis();
		double score = evaluateTopKLB(tg, s, ind, topk, lbs, exact);
		long evalt = System.currentTimeMillis() - start;
		//System.err.print("LB topk:" + topkt  + "\teval:" + evalt + "\n" );
		
		return score;
	}
	
	// topK LB in real matrix
	public static double getLBtopKInMatrix(TimeGraph tg, int s, int e, 
			   double[][] mat, int k, 
			   IntervalIndexExactCeiling ind, double[][] lbs,
			   boolean exact) throws IOException{
		long start = System.currentTimeMillis();
		TreeSet<ComparableInterval> topk = getTopKIntervals(mat, k);
		long topkt = System.currentTimeMillis() - start;
		
		start = System.currentTimeMillis();
		double score = evaluateTopKLB(tg, s, ind, topk, lbs, exact);
		long evalt = System.currentTimeMillis() - start;
		System.err.print("LB topk:" + topkt  + "\teval:" + evalt + "\n" );
		return score;
	}
	
	// evaluates a LB for each of the topk intervals
	public static double evaluateTopKLB(ITimeGraph tg, int s,
			IntervalIndexExactCeiling ind, TreeSet<ComparableInterval> topk,
			double[][] lbs, boolean exact) throws IOException {
		Double neginf = Double.NEGATIVE_INFINITY;
		double score = 0;
		for (ComparableInterval ci: topk.descendingSet()) {
			if (neginf.equals(lbs[ci.size][ci.pos])){
				if (ind==null) {
					tg.aggregateByTime(ci.pos+s, ci.pos + s + ci.size);
				} else {
					ind.getAggregate(ci.pos, ci.pos + ci.size, tg.getGraph().we);
				}
				BitSet edges = new BitSet(tg.getm());
				if(exact)MaximumScoreSubgraph(tg.getGraph(), edges);
				else topDownLBFast(tg.getGraph(), edges);
				lbs[ci.size][ci.pos] = tg.getGraph().getScore(edges);
				
//				System.err.print("LB:" + "[" + (ci.pos +s) + "," + 
//						         (ci.pos+ci.size-1+s) + "]" + 
//						         lbs[ci.size][ci.pos] + "\n");
			}
			score = Math.max(score, lbs[ci.size][ci.pos]);
		}
		return score;
	}
	
	// top k points from an interpolated matrix
	public static TreeSet<ComparableInterval> getTopKIntervals(SplineInterpolation si, int size, int k) {
		TreeSet<ComparableInterval> topk = new TreeSet<ComparableInterval>();
		double v = 0;
		for (int sz = 0; sz < size; sz++) {
			for (int pos= 0; pos < size - sz; pos++) {
				v = si.interpolate(sz, pos);
				if(topk.size() < k || v > topk.first().val) {
					ComparableInterval intvl = new ComparableInterval(pos,sz,v);
					topk.add(intvl);
					if (topk.size()>k) topk.remove(topk.first());
				}
			}
		}
		return topk;
	}
	
	// updates the topK intervals, using the interpolation by sampling "near" the
	// top  intervals obtained from the coarse matrix
	private static void updateTopKIntervals(SplineInterpolation si,	TreeSet<ComparableInterval> topk) {
		int k = topk.size();
		TreeSet<ComparableInterval> newtopk = new TreeSet<ComparableInterval>();
		for (ComparableInterval ci: topk) {
			si.getInterpolatedNeighborsGreaterThan(ci, topk, newtopk);
		}
		while (newtopk.size() > 0) {
			// add the new intervals
			newtopk.addAll(topk);
			topk.clear();
			int cnt=0;
			// keep only the top in k
			for(ComparableInterval ci: newtopk.descendingSet()) {
				topk.add(ci); cnt++; 
				if(cnt==k) break;
			}
			newtopk.clear();
			// add more neighbors based on new topk
			for (ComparableInterval ci: topk) {
				si.getInterpolatedNeighborsGreaterThan(ci, topk, newtopk);
			};
		}
	}
	
	// top k points from actual matrix
	public static TreeSet<ComparableInterval> getTopKIntervals(double[][] mat, int k) {
		TreeSet<ComparableInterval> topk = new TreeSet<ComparableInterval>();
		for (int sz = 0; sz < mat.length; sz++) {
			for (int pos= 0; pos < mat.length - sz; pos++) {
				if(topk.size() < k || mat[sz][pos] > topk.first().val) {
					ComparableInterval intvl = new ComparableInterval(pos,sz,mat[sz][pos]);
					topk.add(intvl);
					if (topk.size()>k) topk.remove(topk.first());
				}
			}
		}
		return topk;
	}
	
	// methods to compare spline with
	// LB greedy: runs topK in descenig SOP in aggregated intervals
	public static double getLBGreedy(TimeGraph tg, int s, int e, int size, int k,
			IntervalIndexExactCeiling ind, double[][] coarse, long limit,
			boolean exact ) throws IOException {
		
		long start = System.currentTimeMillis();
		TreeSet<ComparableInterval> topk = getTopKIntervals(coarse, k*2);
		
		double mscore = 0;
		for (ComparableInterval ci: topk.descendingSet()) {
			ind.getAggregate(ci.pos, ci.pos + ci.size, tg.we);
			BitSet edges = new BitSet(tg.getm());
			if (exact) MaximumScoreSubgraph(tg, edges);
			else topDownLBFast(tg, edges);
			mscore = Math.max(tg.getScore(edges), mscore);
			if (System.currentTimeMillis() - start > limit) break;
		}
		return mscore;
	}
	
	public static double getLBRandomSampling(TimeGraph tg, int s, int e, int size, 
											 IntervalIndexExactCeiling ind, 
											 long limit, boolean exact ) throws IOException {		
		long start = System.currentTimeMillis();
		double mscore = 0;
		while (limit > System.currentTimeMillis() - start) {
			int sz = Common.rand.nextInt(size);
			int pos = Common.rand.nextInt(size-sz);
			ind.getAggregate(pos, pos + sz, tg.we);
			BitSet edges = new BitSet(tg.getm());
			if (exact) MaximumScoreSubgraph(tg, edges);
			else topDownLBFast(tg, edges);
			mscore = Math.max(tg.getScore(edges), mscore);
		}
		return mscore;
	}
	
	
	
	public static double getLBParticleFiltering(TimeGraph tg, int s, int e, int size, int k,
			IntervalIndexExactCeiling ind, double[][] coarse, long limit) {
		long start = System.currentTimeMillis();
		TreeSet<ComparableInterval> topk = getTopKIntervals(coarse, k);
		TreeSet<ComparableInterval> newtop = new TreeSet<ComparableInterval>();
		
		double mscore = topk.last().val;
		while (limit > System.currentTimeMillis() - start) {
			double sum = 0;
			for (ComparableInterval ci: topk) sum += ci.val;
			for (ComparableInterval ci: topk) ci.Probability = ci.val/sum;
			//select k particles 
			
			  // sample around each
			
			// update top k and mscore
			
			mscore = topk.last().val;
		}
		return mscore;
	}
	
	
	
	
	// runs the TopDown on the remaining in a descending order of their UBs. if lb improves prunes additionally
	public static int handleNonPruned(ITimeGraph tg, int s, int e, 
			RunMesasurements rm, double[][] lbs, 
			TreeSet<ComparableInterval> toverify, 
			IntervalIndexExactCeiling ind, boolean exact ) throws IOException {
		int cnt_lb=0;
		double updated_lb = rm.best_lb;
		for (ComparableInterval ci: toverify.descendingSet()){
			if ( ci.val > updated_lb && lbs[ci.size][ci.pos] < 0) {
				//for (int i=0; i < tg.getm(); i++) tg.we[i] = ci.w[i];
				tg.aggregateByTime(ci.pos+s, s+ci.pos+ci.size);
				BitSet edges = new BitSet(tg.getm());
				if (exact) MaximumScoreSubgraph(tg.getGraph(), edges);
				else topDownLBFast(tg.getGraph(), edges);
				lbs[ci.size][ci.pos] = tg.getGraph().getScore(edges); cnt_lb++;
				//assert(lbs[ci.size][ci.pos] == getSolutionScoreInOriginalGraph(tg,ci,t,false));
				if (lbs[ci.size][ci.pos] > updated_lb) 
					updated_lb = lbs[ci.size][ci.pos];
			}
			ci.act = lbs[ci.size][ci.pos];
		}
		rm.best_lb = updated_lb;
		return cnt_lb;
	}
	
	public static ComparableInterval handleNonPruned(TimeGraph tg, int s, 
			int e, boolean exact) throws IOException {
		int size = (e-s+1);
		Tree t = null;
		double score = 0;
		ComparableInterval maxi = new ComparableInterval(-1, -1, -1);
		for (int sz = 0; sz < size; sz++) {
			for (int pos = 0; pos < size-sz; pos++) {
				tg.aggregateByTime(s+pos, s + pos+sz );
				BitSet edges = new BitSet(tg.getm());
				if (exact) MaximumScoreSubgraph(tg, edges);
				else topDownLBFast(tg, edges);
				score = tg.getScore(edges);
				if (score > maxi.val) {
					maxi.pos = pos;
					maxi.size = sz;
					maxi.val = score;
				}
			}
		}
		return maxi;
	}
	
	
	public static double getSolutionScoreInOriginalGraph(TimeGraph tg,
			ComparableInterval ci, Graph fast, boolean print) {
		for (int i=0; i < tg.getm(); i++) tg.we[i] = ci.w[i];
		double sc = 0;
		for (String name:fast.names) {
			int i = tg.getNodeIndex(name);
			for (int j=tg.ind[i]; j<tg.ind[i+1];j++) {
				if(i < tg.endv[j] && tg.we[j] > 0) {
					sc += tg.we[j]; 
				}
			}
		}
		if(fast.getn() > 1) {
			// add tree edges
			for (int i = 0; i < fast.getn(); i++) {
				String n1 = fast.names[i];
				int ind1 = tg.getNodeIndex(n1);
				for (int j = fast.ind[i]; j < fast.ind[i+1]; j++) {
					String n2 = fast.names[fast.endv[j]];
					int ind2 = tg.getNodeIndex(n2);
					if(ind1 < ind2) continue;
					if(tg.we[tg.getEdgeIndex(ind1, ind2)] < 0) {
						sc += tg.we[tg.getEdgeIndex(ind1, ind2)];
					}
				}
			}
		}
		if(print) {
			for (String name:fast.names) {			
				System.err.print(tg.neighToString(name) + "\n");
			}
		}
		
		return sc;
	}
	
	// --------OLD AND UNUSED ----------
	
	// computes the heuristic in the Highest LB position
	private static double refineLBBUp(TimeGraph tg, int s, double[][] lbs) {
		double lb = 0;
		int sz = -1, pos = -1;
		// find the position of the best lb that has not been extended
		for(int i = 0; i < lbs.length; i++) {
			for (int j = 0; j < lbs[i].length-i; j++) {
				if (lbs[i][j] > lb) {
					lb = lbs[i][j];
					sz = i;
					pos = j;
				}
			}
		}
		
		// Bottom up heuristic refinement
		tg.aggregateByTime(s+pos,s+pos+sz);
		BUp b = new BUp(tg.getAggNodeWeighted());
		b.resetClusters = false;
		b.useTM = false;
		b.computeOptimalLMSMTM();
		lbs[s+pos][s+pos+sz] = b.getBestClusterScore(true);
		return Math.max(lb, lbs[s+pos][s+pos+sz]); 
	}
	
	// DEPRECATED computes lb around the max lb found so far OLD
	private static double refineLBGreedyLB(TimeGraph tg, int s, double[][] lbs) {
		double lb = 0;
		int sz = -1, pos = -1;
		// find the position of the best lb that has not been extended
		for(int i = 0; i < lbs.length; i++) {
			for (int j = 0; j < lbs[i].length-i; j++) {
				if (lbs[i][j] > lb && expandable(i,j,lbs)) {
					sz = i; pos = j; lb = lbs[i][j]; 
				}
			}
		}
		
		// try all adjacent positions to the max
		int offset_pos = s + pos;
		// size-1
		if (sz>0 ) {
			// sz-1, pos-1
			if (pos-1 >= 0) {
				if (lbs[sz-1][pos-1] == Double.NEGATIVE_INFINITY){
					tg.aggregateByTime(offset_pos-1,offset_pos-1+sz-1);
					Tree t = topDownLBold(tg, null);
					lbs[sz-1][pos-1] = t.getScore();
				}
			}
			// sz-1, pos
			if (lbs[sz-1][pos] == Double.NEGATIVE_INFINITY){
				tg.aggregateByTime(offset_pos,offset_pos+sz-1);
				Tree t =  topDownLBold(tg, null);
				lbs[sz-1][pos] = t.getScore();
			}
			// sz-1, pos+1
			if (pos+1 < lbs.length-(sz-1)) { 
				if (lbs[sz-1][pos+1] == Double.NEGATIVE_INFINITY){
					tg.aggregateByTime(offset_pos+1,(offset_pos+1)+(sz-1));
					Tree t =  topDownLBold(tg, null);
					lbs[sz-1][pos+1] = t.getScore();
				}
			}
		}
		
		// same size
		// sz, pos-1
		if (pos-1 >= 0) {
			if (lbs[sz][pos-1] == Double.NEGATIVE_INFINITY){
				tg.aggregateByTime(offset_pos-1,offset_pos-1+sz);
				Tree t =  topDownLBold(tg, null);
				lbs[sz][pos-1] = t.getScore();
			}
		}
		// sz, pos+1
		if (pos+1 < lbs.length-sz ) { 
			if (lbs[sz][pos+1] == Double.NEGATIVE_INFINITY){
				tg.aggregateByTime(offset_pos+1,offset_pos+sz);
				Tree t =  topDownLBold(tg, null);
				lbs[sz][pos+1] = t.getScore();
			}
		}
		
		// size + 1
		if (sz+1 < lbs.length) {
			// sz+1, pos-1
			if (pos-1 >= 0) {
				if (lbs[sz+1][pos-1] == Double.NEGATIVE_INFINITY){
					tg.aggregateByTime(offset_pos-1,(offset_pos-1)+(sz+1));
					Tree t =  topDownLBold(tg, null);
					lbs[sz+1][pos-1] = t.getScore();
				}
			}
			// sz+1, pos
			if (pos < lbs.length-(sz+1) ) {
				if (lbs[sz+1][pos] == Double.NEGATIVE_INFINITY){
					tg.aggregateByTime(offset_pos,offset_pos+sz+1);
					Tree t =  topDownLBold(tg, null);
					lbs[sz+1][pos] = t.getScore();
				}
			}
			// sz+1, pos+1
			if (pos+1 < lbs.length-(sz+1) ) { 
				if (lbs[sz+1][pos+1] == Double.NEGATIVE_INFINITY){
					tg.aggregateByTime(offset_pos+1,(offset_pos+1)+(sz+1));
					Tree t =  topDownLBold(tg, null);
					lbs[sz+1][pos+1] = t.getScore();
				}
			}
		}
		// find the position of the best lb
		for(int i = 0; i < lbs.length; i++) {
			for (int j = 0; j < lbs[i].length-i; j++) {
				if (lbs[i][j] > lb) {
					lb = lbs[i][j]; 
				}
			}
		}
		return lb;
	}

	// DEPRECATED get a good upper bound within the interval
	// initially just compute all
	public static double getAGoodLB(TimeGraph tg, int s, int e, double[][] lbs, boolean useActual) {
		double lb = 0.0;
		int pos = -1;
		for (int i=s; i<=e; i++) {
			tg.aggregateByTime(i, i);
			Tree t =  topDownLBold(tg, null);
			lbs[0][i-s] = t.getScore();
			if (lb<lbs[0][i-s]) {
				lb = lbs[0][i-s];
				pos = i;
			}
			//System.err.print(i + "\t" + lbs[i-s] + "\n");
		}
		System.err.print("Best LB in a single interval:" + lb + "\n");
		
		// Refine the lower bound in a greedy manner
		double newlb = 0;
		while((newlb = refineLBGreedyLB(tg, s, lbs)) > lb) {
			lb = newlb;
			//printPruningPerc(99.0, lbs);
		} while (newlb > lb);
		//printPruningPerc(99.0, lbs);
		System.err.print("Greedy extend around biggest: " + lb + "\n");	
		
		if (useActual) {
			// Compute the heuristic UB 
			double act = refineLBBUp(tg, s, lbs);
			System.err.print("Actual at biggest: " + act + "\n");
			lb = Math.max(act, lb);
		}

		return lb;
	}
	
	// DEPRECATED enumerates all subintervals
	public static double getBestLB(TimeGraph tg, int s, int e,  double[][] lbs) {
		int l = -1, r = -1; 
		double stmp, score=0;
		for (int i=s; i<=e; i++) {
			for (int j=i; j <=e; j++) {
				tg.aggregateByTime(i, j);
				Tree t =  topDownLBold(tg, null);
				stmp = t.getScore();
				lbs[j-i][i-s] = stmp;
				if(score<stmp) {l = i; r = j; score = stmp;}
			}
		}
		ExtToolsDriver.plotMatrix(lbs, 1.0,0.0,lbs.length);
		System.err.print("Actual Interval [" + l + "-" + r+ "]: " + 
				score);
		return score;
	}
	
	
	
	// ------------------- Sum of positive based solutions -----------------------------------
	
	// Computes and approximation to the sum of positive, 
	// grouping intervals with the same starting position that overlap by at least alpha fraction 
	public static double coarse(ITimeGraph tg, int s, int e, double alpha, 
									IntervalIndexExactCeiling ind,
									double[][] ceiling,
									double[][] actual) {
		int size =  e-s+1;
		// ceiling lines
		double[] ceil = new double[tg.getm()];
		double[] act = new double[tg.getm()];
		for (int sz=0; sz < size; sz++) {
			
			for (int pos=0; pos + sz < size; pos++) {
				if (ceiling[sz][pos] >= 0) continue;
				Arrays.fill(ceil, 0);
				ind.getCeiling(pos, pos+sz, pos+(int)Math.min( sz/alpha,size-pos-1), ceil);
				ind.getAggregate(pos, pos+sz, act);
				
				//int a = tg.slices[pos].cardinality();
				double sp  = 0, sp_act = 0;
				for (int i =0; i < tg.getn(); i++) {
					for (int j= tg.ind(i); j < tg.ind(i+1); j++) {
						if (i < tg.endv(j)) {
							if(ceil[j] > 0)  
								sp += ceil[j];
							if(act[j]  > 0) 
								sp_act += act[j];
						}
					}
				}
				
				ceiling[sz][pos] = sp_act;
				actual[sz][pos] = sp_act;
				for (int ssz = sz; ssz <= sz/alpha && ssz < size-pos; ssz++) {
					if(ceiling[ssz][pos] < 0) 
						ceiling[ssz][pos] = sp; 
				}
			}
		}
		
		// Compute the whole interval(needed for the interpolation)
		tg.aggregateByTime(s, e);
		actual[actual.length-1][0] = 0;
		for (int i =0; i < tg.getn(); i++) {
			for (int j= tg.ind(i); j < tg.ind(i+1); j++) {
				if (i < tg.endv(j)) {
					if(tg.we(j)  > 0) 
						actual[actual.length-1][0] += tg.we(j);
				}
			}
		}
		
		//permuteClusterIds(partit, cid);
		return 0;
	}
	
	// Computes the sum of positives in all intervals
	public static void sumPositive(TimeGraph tg, int s, int e, double[][] sumpos) {
		double spos=0;
		// hold the single interval score vectors
		double[][] singles = new double[e-s+1][tg.getm()];
		double[] curr = new double[tg.getm()];
		for (int pos=s; pos<=e; pos++) {
			tg.aggregateByTime(pos, pos);
			spos = 0;
			for (int i = 0; i < tg.getn(); i++) {
				for (int j=tg.ind[i]; j < tg.ind[i+1]; j++) {
					singles[pos-s][j] = tg.we[j];
					if (i < tg.endv[j]) { // Count each edge only once
						spos += Math.max(0, tg.we[j]);
					}
				}
			}
			sumpos[0][pos-s] = spos;
		}
		for (int pos = s; pos <= e; pos++) {
			curr = Arrays.copyOf(singles[pos-s], curr.length);
			for (int sz = 1; sz <= e-pos; sz++) {
				spos = 0;
				for (int i = 0; i < tg.getn(); i++) {
					for (int j=tg.ind[i]; j < tg.ind[i+1]; j++) {
						curr[j] += singles[pos-s+sz][j];
						if (i < tg.endv[j]) { // Count each edge only once
							spos += Math.max(0, curr[j]);
						}
					}
				}
				sumpos[sz][pos-s] = spos;
			}
		}
	}
	
	// experimental only;Computes the sum of positives in all intervals; obtains a LB based on top-k sum-of-pos intervals
	public static void probePositive(TimeGraph tg, int s, int e, 
									   double[][] sumpos,
									   double[][] sumneg,
									   double[][] stdev) {
		
		double spos=0, sneg=0, score=0, std = 0, avg = 0;
		
		// hold the single interval score vectors
		double[][] singles = new double[e-s+1][tg.getm()];
		double[] curr = new double[tg.getm()];
		for (int pos=s; pos<=e; pos++) {
			tg.aggregateByTime(pos, pos);
			spos = 0; sneg = 0;
			for (int i = 0; i < tg.getn(); i++) {
				for (int j=tg.ind[i]; j < tg.ind[i+1]; j++) {
					singles[pos-s][j] = tg.we[j];
					if (i < tg.endv[j]) { // Count each edge only once
						spos += Math.max(0, tg.we[j]);
						sneg += Math.min(0, tg.we[j]);
					}
				}
			}
//			std = 0;
//			avg = (spos + sneg)/(tg.getm()/2.0);
//			for (int i = 0; i < tg.getn(); i++) {
//				for (int j=tg.ind[i]; j < tg.ind[i+1]; j++) {
//					if (i < tg.endv[j]) { // Count each edge only once
//						std +=(avg - tg.we[j])*(avg - tg.we[j]);
//					}
//				}
//			}
//	
//			stdev[0][pos-s] = Math.sqrt(std);
			sumpos[0][pos-s] = spos;
//			sumneg[0][pos-s] = sneg;
		}
		for (int pos = s; pos <= e; pos++) {
			curr = Arrays.copyOf(singles[pos-s], curr.length);
			for (int sz = 1; sz <= e-pos; sz++) {
				spos = 0; sneg = 0;
				for (int i = 0; i < tg.getn(); i++) {
					for (int j=tg.ind[i]; j < tg.ind[i+1]; j++) {
						curr[j] += singles[pos-s+sz][j];
						if (i < tg.endv[j]) { // Count each edge only once
							spos += Math.max(0, curr[j]);
							sneg += Math.min(0, curr[j]);
						}
					}
				}
//				std = 0;
//				avg = (spos + sneg)/(tg.getm()/2.0);
//				for (int i = 0; i < tg.getn(); i++) {
//					for (int j=tg.ind[i]; j < tg.ind[i+1]; j++) {
//						if (i < tg.endv[j]) { // Count each edge only once
//							std +=(avg - tg.we[j])*(avg - tg.we[j]);
//						}
//					}
//				}
//				stdev[sz][pos-s] = Math.sqrt(std);
				sumpos[sz][pos-s] = spos;
//				sumneg[sz][pos-s] = sneg;
			}
		}
		//ExtToolsDriver.plotMatrix(sumpos,1.0,0.0,sumpos.length);
	}
	
	// use this for node-weighted edge-weighted
	public static double getMinSTLB(Graph g) {
		if(g.getn()==1) return g.wn[0];
		Tree t = ST.computeMST(g);
		// get score returns the proper score (no double counting)
		return PCST.computeOptTreeNW(t).getScore();
	}
	
	public static void testUBStructure() {
		Graph g = generator.GraphGenerator.getTestGraphUBTest();
		double ub = getUBAdjNegative(g);
		assert(5.0 == ub);
	}
		
	public static void main(String[] args) {
		testUBStructure();
	}
	
	
}
