public class StandardLDA {
	/**
	 * @param args
	 * @throws IOException 
	 * @throws MatrixException 
	 */
	public static void main(String[] args) throws MatrixException, IOException 
	{
		// Parse command-line args
		//
		String basefn = args[0];
		int numsamp = Integer.parseInt(args[1]);
		int topN = Integer.parseInt(args[2]); 
		int randseed = Integer.parseInt(args[3]); 

		// Load corpus and parameters (checking vocab dim agreement)
		//
		LDAParameters p = new LDAParameters(basefn, randseed);
		Corpus c = new Corpus(basefn);		
		assert(p.W == c.W);
		
		// Run standard LDA for numsamp 
		//
		DiscreteSample s = runStandardLDA(c, p, numsamp);
 		
		// Write out results
		//
		s.writePhiTheta(p, basefn);
		s.writeSample(basefn);
		c.writeTopics(basefn, s.getPhi(p), topN);
	}
	
	/**
	 * Standard (no logic) LDA, using collapsed Gibbs sampling
	 * 
	 * @param c
	 * @param p
	 * @param numsamp
	 * @return
	 */
	public static DiscreteSample runStandardLDA(Corpus c, LDAParameters p, int numsamp)
	{	
		// Do Collapsed Gibbs sampling and return final sample
		DiscreteSample s = CollapsedGibbs.doGibbs(c, p, numsamp);		
		return s;
	}
} public class CollapsedGibbs {

	/**
	 * Do a single Collapsed Gibbs sample, but with topic-words Phi fixed
	 * 
	 * @param c Contains words, documents
	 * @param p Contains hyperparameters
	 * @param s Sample object to be updated in place
	 * @param onlineInit If true, don't pre-subtract counts 
	 */
	public static void fixedPhiSample(Corpus c, LDAParameters p, DiscreteSample s,
			boolean onlineInit, double[][] phi)
	{
		int N = c.N;
		int T = p.T;
		
		double[] tmp = new double[T];
		for(int i = 0; i < N; i++)
		{
			// Reset sampling normalization sum to 0
			double normsum = 0;

			// Remove current assignment from counts
			// (unless we're doing 'online-style' init)
			if(!onlineInit)
			{
				s.updateCounts(c.w[i], s.z[i], c.d[i], -1);				
			}
		
			// Get un-normalized probabilities for each topic
			for(int j = 0; j < T; j++)
			{				
				double phiterm = phi[j][c.w[i]];
				double num2 = s.nd[c.d[i]][j] + p.alpha[j];				
				tmp[j] = phiterm * num2;
				normsum += tmp[j];
			}		
			
			// Sample the assignment
			s.z[i] = MiscUtil.multSample(p.rng, tmp, normsum);
			
			// Update the count matrices
			s.updateCounts(c.w[i], s.z[i], c.d[i], 1);			
		}		
		return;	
	}
	
	public static DiscreteSample runGroundGibbs(Corpus c, LDAParameters p, 
			DiscreteSample s, List<LogicRule> rules, int numsamp)
	{
		// Cast LogicRule to GroundableRule
		List<GroundableRule> grules = GroundRules.groundCast(rules);
		
		// Apply evidence
		for(GroundableRule gr : grules)
			gr.applyEvidence(c, p.T);
		// Ground the Logic Rules
		GroundRules gr = new GroundRules(grules, s.z, p.rng, p.T);
		
		// Best overall z-assignment found thus far
		//		
		int[] bestz = new int[s.z.length];
		System.arraycopy(s.z, 0, bestz, 0, s.z.length);
		MirrorDescent logicrules = new MirrorDescent(rules, p.rng);
		double bestobj = EvalLDA.ldaLoglike(s, p) + logicrules.satWeight(s.z);
		
		// Do the samples
		for(int si = 0; si < numsamp; si++)
		{
			System.out.println(String.format("Sample %d of %d", si, numsamp));
			CollapsedGibbs.groundGibbsSample(gr, c, p, s, false);
			
			// Do we have a new best sample?
			if((si+1) % 500 == 0)
			{
				double newobj = EvalLDA.ldaLoglike(s, p) + logicrules.satWeight(s.z);
				if(newobj > bestobj)
				{
					System.arraycopy(s.z, 0, bestz, 0, s.z.length);
					bestobj = newobj;
				}
			}
		}
		
		// Use the best sample
		s.repopZ(c, bestz);
		return s;
	}
	
	/**
	 * Externally called method for Ground LogicLDA sampling
	 * @param gr
	 * @param c
	 * @param p
	 * @param numsamp
	 * @return
	 */
	public static DiscreteSample doGroundGibbs(List<GroundableRule> grules,
			Corpus c, LDAParameters p, int numsamp)
	{
		// Get relevant dimensions
		int N = c.N;		
		int T = p.T;
		int W = p.W;
		int D = c.D;
		
		// Online initialization
		DiscreteSample s = new DiscreteSample(N, T, W, D);
		CollapsedGibbs.gibbsSample(c, p, s, true);
						
		// Ground the Logic Rules
		GroundRules gr = new GroundRules(grules, s.z, p.rng, p.T);		
		
		// Do the samples
		for(int si = 0; si < numsamp; si++)
		{
			System.out.println(String.format("Sample %d of %d", si, numsamp));
			CollapsedGibbs.groundGibbsSample(gr, c, p, s, false);
		}
		
		return s;
	}
	
	
	/**
	 * Do a single Ground LogicLDA sample
	 * 
	 * @param logicweights
	 * @param c
	 * @param p
	 * @param s
	 * @param onlineInit
	 */
	public static void groundGibbsSample(GroundRules gr,
			Corpus c, LDAParameters p, DiscreteSample s, 
			boolean onlineInit)            
	{				
		int N = c.N;
		int T = p.T;
		
		double[] tmp = new double[T];
		for(int i = 0; i < N; i++)
		{
			// Reset sampling normalization sum to 0
			double normsum = 0;

			// Remove current assignment from counts
			// (unless we're doing 'online-style' init)
			if(!onlineInit)
			{
				s.updateCounts(c.w[i], s.z[i], c.d[i], -1);				
			}
		
			// Get un-normalized probabilities for each topic
			for(int j = 0; j < T; j++)
			{
				double num1 = s.nw[c.w[i]][j] + p.beta[j][c.w[i]];
				double den1 = s.nwcolsums[j] + p.betasums[j];
				double num2 = s.nd[c.d[i]][j] + p.alpha[j];			
								
				tmp[j] = (num1 / den1) * num2;
				
				// If applicable, multiply the standard Collapsed Gibbs term 
				// by the exp of the logic contribution				
				if(gr.inLogic(i))
				{
					s.z[i] = j;
					tmp[j] *= Math.exp(gr.evalAssign(s.z, i));
				}					
								
				normsum += tmp[j];
			}		
			
			// Sample the assignment
			s.z[i] = MiscUtil.multSample(p.rng, tmp, normsum);
			
			// Update the count matrices
			s.updateCounts(c.w[i], s.z[i], c.d[i], 1);			
		}		
		return;
	}
	
	
	/**
	 * External method for doing online-init, then numsamp Logic Gibbs samples
	 * 
	 * @param logicweights N x T matrix of (independent) rule contributions
	 * @param c Contains words, documents
	 * @param p Contains hyperparameters
	 * @param numsamp How many samples to do
	 * @param randseed Seed for random number generator 
	 * @return The final sample from the Markov Chain
	 */
	public static DiscreteSample doLogicGibbs(double[][] logicweights, 
			Corpus c, LDAParameters p, int numsamp)
	{
		// Get relevant dimensions
		int N = c.N;		
		int T = p.T;
		int W = p.W;
		int D = c.D;
		
		// Online initialization
		DiscreteSample s = new DiscreteSample(N, T, W, D);
		CollapsedGibbs.logicGibbsSample(logicweights, c, p, s, true);
		
		// Do the samples
		for(int si = 0; si < numsamp; si++)
		{
			System.out.println(String.format("Sample %d of %d", si, numsamp));
			CollapsedGibbs.logicGibbsSample(logicweights, c, p, s, false);
		}
		
		return s;
	}
	
	/**
	 * Do a single Logic Collapsed Gibbs sample
	 * 
	 * @param logicweights N x T matrix of (independent) rule contributions
	 * @param c Contains words, documents
	 * @param p Contains hyperparameters
	 * @param s Sample object to be updated in place
	 * @param rng Random number generator for sampling
	 * @param onlineInit If true, don't pre-subtract counts 
	 */
	public static void logicGibbsSample(double[][] logicweights,
			Corpus c, LDAParameters p, DiscreteSample s, 
			boolean onlineInit)            
	{				
		int N = c.N;
		int T = p.T;
		
		double[] tmp = new double[T];
		for(int i = 0; i < N; i++)
		{
			// Reset sampling normalization sum to 0
			double normsum = 0;

			// Remove current assignment from counts
			// (unless we're doing 'online-style' init)
			if(!onlineInit)
			{
				s.updateCounts(c.w[i], s.z[i], c.d[i], -1);				
			}
		
			// Get un-normalized probabilities for each topic
			for(int j = 0; j < T; j++)
			{
				double num1 = s.nw[c.w[i]][j] + p.beta[j][c.w[i]];
				double den1 = s.nwcolsums[j] + p.betasums[j];
				double num2 = s.nd[c.d[i]][j] + p.alpha[j];			
				if(logicweights[i] != null)
					tmp[j] = (num1 / den1) * num2 * Math.exp(logicweights[i][j]);
				else
					tmp[j] = (num1 / den1) * num2;
				normsum += tmp[j];
			}		
			
			// Sample the assignment
			s.z[i] = MiscUtil.multSample(p.rng, tmp, normsum);
			
			// Update the count matrices
			s.updateCounts(c.w[i], s.z[i], c.d[i], 1);			
		}		
		return;
	}
	
	/**
	 * External method for doing online-init, then numsamp Gibbs samples
	 * @param c Contains words, documents
	 * @param p Contains hyperparameters
	 * @param numsamp How many samples to do
	 * @param randseed Seed for random number generator 
	 * @return The final sample from the Markov Chain
	 */
	public static DiscreteSample doGibbs(Corpus c, LDAParameters p, int numsamp)
	{
		// Get relevant dimensions
		int N = c.N;		
		int T = p.T;
		int W = p.W;
		int D = c.D;
		
		// Online initialization
		DiscreteSample s = new DiscreteSample(N, T, W, D);
		CollapsedGibbs.gibbsSample(c, p, s, true);
		
		// Do the samples
		for(int si = 0; si < numsamp; si++)
		{
			System.out.println(String.format("Sample %d of %d", si, numsamp));
			CollapsedGibbs.gibbsSample(c, p, s, false);
		}
		
		return s;
	}
	
	/**
	 * Do a single Collapsed Gibbs sample
	 * 
	 * @param c Contains words, documents
	 * @param p Contains hyperparameters
	 * @param s Sample object to be updated in place
	 * @param rng Random number generator for sampling
	 * @param onlineInit If true, don't pre-subtract counts 
	 */
	public static void gibbsSample(Corpus c, LDAParameters p, DiscreteSample s, 
			boolean onlineInit)            
	{				
		int N = c.N;
		int T = p.T;
		
		double[] tmp = new double[T];
		for(int i = 0; i < N; i++)
		{
			// Reset sampling normalization sum to 0
			double normsum = 0;

			// Remove current assignment from counts
			// (unless we're doing 'online-style' init)
			if(!onlineInit)
			{
				s.updateCounts(c.w[i], s.z[i], c.d[i], -1);				
			}
		
			// Get un-normalized probabilities for each topic
			for(int j = 0; j < T; j++)
			{
				double num1 = s.nw[c.w[i]][j] + p.beta[j][c.w[i]];
				double den1 = s.nwcolsums[j] + p.betasums[j];
				double num2 = s.nd[c.d[i]][j] + p.alpha[j];				
				tmp[j] = (num1 / den1) * num2;
				normsum += tmp[j];
			}		
			
			// Sample the assignment
			s.z[i] = MiscUtil.multSample(p.rng, tmp, normsum);
			
			// Update the count matrices
			s.updateCounts(c.w[i], s.z[i], c.d[i], 1);			
		}		
		return;
	}
} public class DiscreteSample implements Sample {
	
	// Minimum values for MAP estimates 
	public static final double MIN_PHI = 0.001;
	public static final double MIN_THETA = 0.001;
	
	public long[][] nw;
	public long[][] nd;
	public long[] nwcolsums;
	public int[] z;
	
	
	/**
	 * Init from a given *.init or *.sample file
	 * 
	 * @param N
	 * @param T
	 * @param W
	 * @param D
	 * @param filename
	 * @param c
	 */
	public DiscreteSample(int N, int T, int W, int D,
			String filename, Corpus c)
	{
		// Hidden topic assignments
		z = new int[N];		
		// Count matrices
		nw = new long[W][T];
		nd = new long[D][T];		
		// Column sums for word-topic count array
		nwcolsums = new long[T];
		for(int j = 0; j < T; j++)		
			nwcolsums[j] = 0;
		
		this.fromFile(filename, c);
	}
	
	/**
	 * Init to empty sample
	 * 
	 * @param N
	 * @param T
	 * @param W
	 * @param D
	 */
	public DiscreteSample(int N, int T, int W, int D)
	{
		// Hidden topic assignments
		z = new int[N];		
		// Count matrices
		nw = new long[W][T];
		nd = new long[D][T];		
		// Column sums for word-topic count array
		nwcolsums = new long[T];
		for(int j = 0; j < T; j++)		
			nwcolsums[j] = 0;				
	}
	
	/**
	 *  Init from a given z-assignment
	 * @param N
	 * @param T
	 * @param W
	 * @param D
	 * @param givenz
	 * @param c
	 */
	public DiscreteSample(int N, int T, int W, int D,
			int[] givenz, Corpus c)
	{
		assert(givenz.length == N && c.N == N);
		this.z = new int[N];		
						
		// Count matrices
		nw = new long[W][T];
		nd = new long[D][T];		
		// Column sums for word-topic count array
		nwcolsums = new long[T];
		for(int j = 0; j < T; j++)		
			nwcolsums[j] = 0;		
		
		// Set count matrices accordingly
		for(int i = 0; i < z.length; i++)
		{
			assert(givenz[i] > 0 && givenz[i] < T);
			this.z[i] = givenz[i];
			updateCounts(c.w[i], z[i], c.d[i], 1);
		}
	}
	
	/**
	 * Repopulate z-sample, making necessary data struct changes
	 * @param z
	 */
	public void repopZ(Corpus c, int[] newz)
	{
		assert(this.z.length == newz.length);
		for(int i = 0; i < this.z.length; i++)
			this.reassign(c, i, newz[i]);	
	}
	
	
	
	/** 
	 * Do bookkeeping assoc with reassigning z[idx] to newz
	 * @param c
	 * @param idx
	 * @param newz
	 */
	public void reassign(Corpus c, int idx, int newz)
	{
		this.updateCounts(c.w[idx], this.z[idx], c.d[idx], -1);
		this.z[idx] = newz;				
		this.updateCounts(c.w[idx], this.z[idx], c.d[idx], 1);	
	}
	
	
	/**
	 * Calculate MAP phi
	 * 
	 * @param p
	 * @param enw
	 */
	public double[][] mapPhi(LDAParameters p, double[][] phi)
	{				
		// Estimate entries
		for(int t = 0; t < p.T; t++)
		{
			double normsum = 0;
			for(int w = 0; w < p.W; w++)
			{
				// Cannot allow negative entries
				phi[t][w] = Math.max(MIN_PHI, this.nw[w][t] + p.beta[t][w] - 1);
				normsum += phi[t][w];
			}
			// Normalize
			for(int w = 0; w < p.W; w++)
				phi[t][w] /= normsum;		
		}
		return phi;
	}

	/**
	 * Calculate MAP theta
	 * 
	 * @param p
	 * @param end
	 */
	public double[][] mapTheta(LDAParameters p, double[][] theta)
	{
		// Estimate entries
		for(int d = 0; d < nd.length; d++)
		{
			double normsum = 0;
			for(int t = 0; t < p.T; t++)
			{
				// Cannot allow negative entries
				theta[d][t] = Math.max(MIN_THETA, nd[d][t] + p.alpha[t] - 1);
				normsum += theta[d][t];
			}
			// Normalize
			for(int t = 0; t < p.T; t++)
				theta[d][t] /= normsum;			
		}
		return theta;
	}
	
	
	/**
	 * Initialize from a *.sample file
	 * @param filename
	 * @param c
	 */
	public void fromFile(String filename, Corpus c)
	{
		// Read from file
		int[] zfile = MiscUtil.intListUnbox(FileUtil.readIntFile(filename));
		assert(zfile.length == z.length);
		z = zfile;
		
		// Set count matrices accordingly
		for(int i = 0; i < z.length; i++)
		{
			updateCounts(c.w[i], z[i], c.d[i], 1);
		}
	}

	/**
	 * Return the discrete z-vector
	 * 
	 * @return
	 */
	public int[] getZ()
	{
		return z;
	}
	
	/**
	 * Update count matrices used for Gibbs sampling
	 * @param w
	 * @param t
	 * @param d
	 * @param update
	 */
	public void updateCounts(int w, int t, int d, int update)
	{
		nw[w][t] += update;
		nwcolsums[t] += update;
		nd[d][t] += update;
	}
	
	/**
	 * Estimate mean of Phi posterior from current sample
	 * @param p
	 * @return
	 */
	public Matrix getPhi(LDAParameters p)
	{
		long W = nw.length;
		long T = nw[0].length;
		long D = nd.length;

		// Alloc phi matrix
		long[] msize2 = new long[2];
		msize2[0] = T;
		msize2[1] = W;		
		Matrix phi = MatrixFactory.zeros(ValueType.DOUBLE, msize2);

		// Estimate entries
		for(int ti = 0; ti < T; ti++)
		{
			double colsum = nwcolsums[ti];
			double bsum = p.betasums[ti];
			for(int wi = 0; wi < W; wi++)
			{
				double val = (nw[wi][ti] + p.beta[ti][wi]) / (colsum + bsum);
				phi.setAsDouble(val, ti, wi);
			}
		}				
		return phi;
	}

	/**
	 * Estimate mean of Theta posterior from current sample 
	 * @param p
	 * @return
	 */
	public Matrix getTheta(LDAParameters p)
	{
		long D =  nd.length;
		long T =  nd[0].length;
		long[] msize2 = new long[2];
		msize2[0] = D;
		msize2[1] = T;		
		Matrix theta = MatrixFactory.zeros(ValueType.DOUBLE, msize2);		

		// Pre-calculate alpha sum
		double asum = p.alphaSum();

		for(int di = 0; di < D; di++)
		{
			// Get sum of counts for this document			
			double doclen = 0;
			for(int ti = 0; ti < T; ti++)
				doclen += nd[di][ti];

			// Calculate entries
			for(int ti = 0; ti < T; ti++)
			{
				double val = (nd[di][ti] + p.alpha[ti]) / (doclen + asum);
				theta.setAsDouble(val, di, ti);
			}
		}				
		return theta;
	}
	
	public void writeSample(String basefn)
	{
		try
		{			
			ArrayList<Integer> zlist = new ArrayList<Integer>();
			for(int zi : z)
				zlist.add(zi);
			FileUtil.writeIntFile(String.format	("%s.sample", basefn), zlist);
		}
		catch(IOException ioe)
		{
			System.out.println("Problem writing sample out to file");
			System.out.println(ioe.toString());
		}
	}
	 	
	
	/**
	 * Estimate phi/theta and write out to %s.phi/theta   
	 * 
	 * @param p Hyperparameters alpha/beta are required to estimate phi/theta 
	 * @param basefn Base output filename
	 */
	public void writePhiTheta(LDAParameters p, String basefn)
	{		
		try
		{
			FileUtil.writeMatrix(String.format("%s.phi", basefn), 
					this.getPhi(p));
			FileUtil.writeMatrix(String.format("%s.theta", basefn), 
					this.getTheta(p));			
		}
		catch(IOException ioe)
		{
			System.out.println("Phi/theta writeout failed");
			System.out.println(ioe.toString());			
		}		
	}
} public class RelaxedSample implements Sample
{
	// MAP estimates should really have 
	public static final double MIN_PHI = 0.001;
	public static final double MIN_THETA = 0.001;
	
	public double[][] zrelax;
	public double[][] phi;
	public double[][] theta;
	
	public RelaxedSample(Corpus c, LDAParameters p, DiscreteSample s)
	{		
		// Init phi/theta as mean of posteriors from prev sample
		//
		zrelax = new double[c.N][p.T];		
		phi = (s.getPhi(p)).toDoubleArray();
		theta = (s.getTheta(p)).toDoubleArray();
		
		// Init z-relax as posterior given phi/theta
		for(int i = 0; i < s.z.length; i++)
		{			
			double normsum = 0;
			for(int t = 0; t < p.T; t++)
			{
				zrelax[i][t] = theta[c.d[i]][t] * phi[t][c.w[i]];
				assert(!Double.isNaN(zrelax[i][t]));
				normsum += zrelax[i][t]; 
			}
			assert(normsum > 0);
			for(int t = 0; t < p.T; t++)
				zrelax[i][t] /= normsum;						
		}		
	}
	
	/**
	 * Do argmax to convert relaxed z to a single "hard" z sample
	 * @return
	 */
	public int[] getZ()
	{
		int N = zrelax.length;
		int T = zrelax[0].length;
		int[] hardz = new int[N];
		for(int i = 0; i < N; i++)
		{
			int maxidx = 0;
			double maxval= -1;
			for(int t = 0; t < T; t++)
			{
				if(zrelax[i][t] > maxval)
				{
					maxidx = t;
					maxval = zrelax[i][t];					
				}				
			}
			hardz[i] = maxidx;			
		}
		return hardz;
	}
		
	/**
	 * Do an Entropic Mirror Descent step on a single relaxed z-entry
	 * 
	 * @param i Index to update
	 * @param gradSignMag Sign and magnitude of gradient
	 * @param tGradIdx Indices of entries w/ non-zero gradients
	 * @param stepSize Step size parameter
	 */
	public void emdaStep(Gradient stepGrad, double stepSize)
	{		
		// Each z-entry where we have gradient information
		for(int gi = 0; gi < stepGrad.gradients.length; gi++)
		{
			// Get true index and gradient
			int i = stepGrad.indices[gi];
			double[] curGrad = stepGrad.gradients[gi];
			// Multiply entries w/ non-zero gradients
			for(int t = 0; t < curGrad.length; t++)
				zrelax[i][t] *= Math.exp(stepSize * curGrad[t]);
			// Re-normalize
			double normsum = 0;
			for(double val : zrelax[i])
				normsum += val;
			assert(normsum > 0);
			for(int t = 0; t < curGrad.length; t++)
				zrelax[i][t] /= normsum;
		}
	}	
		
	public Matrix getPhi(LDAParameters p)
	{
		return MatrixFactory.importFromArray(phi);
	}
	
	public Matrix getTheta(LDAParameters p)
	{
		return MatrixFactory.importFromArray(theta);
	}
	
	/**
	 * Update MAP estimates of phi/theta
	 * @param c
	 * @param p
	 */
	public void updatePhiTheta(Corpus c, LDAParameters p)
	{
		// Calculate 'expected' NW / ND count matrices 
		double[][] enw = new double[p.W][p.T];
		double[][] end = new double[c.D][p.T];		
		for(int i = 0; i < c.N; i++)
		{
			for(int t = 0; t < p.T; t++)
			{
				enw[c.w[i]][t] += zrelax[i][t];
				end[c.d[i]][t] += zrelax[i][t];				
			}
		}
		
		// Update our MAP estimates
		updateMapPhi(p, enw);
		updateMapTheta(p, end);
	}

	/**
	 * Re-calculate MAP phi
	 * 
	 * @param p
	 * @param enw
	 */
	private void updateMapPhi(LDAParameters p, double[][] enw)
	{				
		// Estimate entries
		for(int t = 0; t < p.T; t++)
		{
			double normsum = 0;
			for(int w = 0; w < p.W; w++)
			{
				// Cannot allow negative entries
				phi[t][w] = Math.max(MIN_PHI, enw[w][t] + p.beta[t][w] - 1);
				normsum += phi[t][w];
			}
			// Normalize
			for(int w = 0; w < p.W; w++)
				phi[t][w] /= normsum;								
		}
	}

	/**
	 * Re-calculate MAP theta
	 * 
	 * @param p
	 * @param end
	 */
	private void updateMapTheta(LDAParameters p, double[][] end)
	{
		// Estimate entries
		for(int d = 0; d < end.length; d++)
		{
			double normsum = 0;
			for(int t = 0; t < p.T; t++)
			{
				// Cannot allow negative entries
				theta[d][t] = Math.max(MIN_THETA, end[d][t] + p.alpha[t] - 1);
				normsum += theta[d][t];
			}
			// Normalize
			for(int t = 0; t < p.T; t++)
				theta[d][t] /= normsum;			
		}
	}
	
	/**
	 * Convert to discrete z-vector and write out to file
	 * 
	 * @param basefn
	 */
	public void writeSample(String basefn)
	{
		try
		{			
			ArrayList<Integer> zlist = new ArrayList<Integer>();
			for(int zi : this.getZ())
				zlist.add(zi);
			FileUtil.writeIntFile(String.format	("%s.sample", basefn), zlist);
		}
		catch(IOException ioe)
		{
			System.out.println("Problem writing sample out to file");
			System.out.println(ioe.toString());
		}
	}
	
	/**
	 * Estimate phi/theta and write out to %s.phi/theta   
	 * 
	 * @param p Hyperparameters alpha/beta are required to estimate phi/theta 
	 * @param basefn Base output filename
	 */
	public void writePhiTheta(LDAParameters p, String basefn)
	{		
		try
		{
			Matrix matphi = MatrixFactory.importFromArray(phi);
			FileWriter phiout = new FileWriter(String.format("%s.phi",
					basefn));
			matphi.exportToWriter(FileFormat.TXT, phiout);
			phiout.close();
					
			Matrix mattheta = MatrixFactory.importFromArray(theta);
			FileWriter thetaout = new FileWriter(String.format("%s.theta",
					basefn));
			mattheta.exportToWriter(FileFormat.TXT, thetaout);
			thetaout.close();
		}
		catch(IOException ioe)
		{
			System.out.println("Phi/theta writeout failed");
			System.out.println(ioe.toString());			
		}		
	}
} public class MirrorDescent 
{

	private Random rng;
	private double[] ruleWeights;
	private double ruleWeightSum;
	public LogicRule[] rules;
	private double pullfactor;
	
	public MirrorDescent(List<LogicRule> lstRules, Random rng)
	{		
		// Rule weights will be used for sampling
		ruleWeightSum = 0;
		Vector<Double> vecRuleWeights = new Vector<Double>();
		for(LogicRule lr : lstRules)
		{
			double weight = lr.getTotalSamplingWeight();
			ruleWeightSum += weight;
			vecRuleWeights.add(weight);
		}
		ruleWeights = MiscUtil.doubleListUnbox(vecRuleWeights);
		
		// Save the actual rules and random number generator
		this.rng = rng;
		rules = new LogicRule[lstRules.size()];
		rules = lstRules.toArray(rules);
	}

	/**
	 * Calc total satisfied weight over all rules
	 * 
	 * @param z
	 * @return
	 */
	public double satWeight(int[] z)
	{
		double weight = 0;
		for(LogicRule lr : this.rules)
			weight += lr.getRuleWeight() * lr.numSat(z);
		return weight;
	}
	
	public double totalWeight()
	{
		double total = 0;
		for(LogicRule lr : this.rules)
			total += lr.getRuleWeight() * lr.numGroundings();
		return total;
	}
	
	/**
	 * Get the z-label style weights for all IndependentRule 
	 * (for use in Logic Collapsed Gibbs, etc)
	 *  
	 * @param N
	 * @param T
	 * @return
	 */
	public double[][] seedsToZL(int N, int T)
	{
		double[][] retval = new double[N][]; // default to null
		for(LogicRule lr : rules)	
		{
			if(lr instanceof IndependentRule)			
			{
				// Need to 'sparse-add' these double[][] 
				double[][] ruleweights = ((IndependentRule) lr).toZLabel(N, T); 
				for(int i = 0; i < N; i++)
				{
					if(ruleweights[i] != null && retval[i] != null)
					{
						// both non-null, need to add
						for(int t = 0; t < ruleweights[i].length; t++)
							retval[i][t] += ruleweights[i][t];
					}
					else if(ruleweights[i] != null)
					{
						// only ruleweights non-null, just copy over
						retval[i] = ruleweights[i];
					}
					// ruleweights null so just continue
				}
			}
		}
		return retval;
	}
	
	/**
	 * Write logic satisfaction report out to *.logic file
	 * 
	 * @param z
	 * @param basename
	 */
	public void satReport(int[] z, String basename)
	{
		try 
		{
			String filename = String.format("%s.logic", basename);
			FileUtil.fileSpit(filename,this.satReport(z));
		}
		catch (IOException ioe)
		{
			System.out.println("Problem writing sample out to file");
			System.out.println(ioe.toString());
		}
	}
	
	
	/**
	 * Return logic satisfaction report as String
	 * 
	 * @param z
	 * @return
	 */
	public String satReport(int[] z)
	{
		String retval = "";
		retval += String.format("%d rules\n-------\n\n", rules.length);
		for(LogicRule lr : rules)
		{
			retval += lr.toString();
			long numground = lr.numGroundings();
			retval += String.format("%d groundings (weight = %.1f)\n",
					numground, lr.getTotalSamplingWeight());
			int numsat = lr.numSat(z); 
			retval += String.format("%d groundings satisfied (%d unsat)\n\n",
					numsat, numground - numsat);
		}
		return retval;
	}
	
	/**
	 * Nice summary of rule set and nontrivial groundings
	 */
	public String toString()
	{
		int totalground = 0;
		String retval = "";
		retval += String.format("%d rules\n-------\n\n", rules.length);
		for(LogicRule lr : rules)
		{
			retval += lr.toString();
			retval += String.format("%d groundings (weight = %.1f)\n\n",
					lr.numGroundings(), lr.getTotalSamplingWeight());
			totalground += lr.numGroundings();
		}
		
		retval += String.format("Total number of nontrivial groundings = %d\n",
				totalground);
		
		return retval;
	}	
	
	public static RelaxedSample runSGD(Corpus c, LDAParameters p, 
			List<LogicRule> rules, int numouter, 
			int numinner, DiscreteSample s)
	{
		// Apply evidence
		for(LogicRule lr : rules)
			lr.applyEvidence(c, p.T);
		
		// Init rule set
		MirrorDescent rs = new MirrorDescent(rules, p.rng);
		
		// Init relaxed z-sample 
		RelaxedSample relax = new RelaxedSample(c, p, s);						
						
		// Do LogicLDA MAP inference via Stochastic Gradient Descent				
		double stepa = Math.sqrt(numinner);
		double stepb = (double) numinner;
		relax = rs.doSGD(c, p, relax, numouter, numinner, stepa, stepb);
		
		// Return final RelaxedSample
		return relax;
	}
	
	 /**
	  * Do stochastic gradient descent MAP inference
	  * 
	  * @param c
	  * @param p
	  * @param relax
	  * @param numouter
	  * @param numinner
	  * @param stepa
	  * @param stepb
	  * @return
	  */
	public RelaxedSample doSGD(Corpus c, LDAParameters p,
			RelaxedSample relax,
			int numouter, int numinner,
			double stepa, double stepb)							
	{
		// Each outer loop, re-estimate phi/theta
		for(int nout = 0; nout < numouter; nout++)
		{
			System.out.println(String.format("Outer loop %d of %d (%d inner)", 
					nout, numouter, numinner));
			relax.updatePhiTheta(c, p);
			// Each inner loop takes a single stochastic gradient step			
			for(int nin = 0; nin < numinner; nin++)
			{
				double stepSize = stepa / Math.sqrt(stepb + nin);
				// Randomly sample a rule
				int chosen = MiscUtil.multSample(rng, 
						ruleWeights, ruleWeightSum);
				
				// Randomly sample a grounding and gradient
				Gradient rGrad = rules[chosen].randomGradient(relax, rng);
				
				// Take an EMDA step
				relax.emdaStep(rGrad, stepSize);
			}			
		}				
		return relax;
	}
		
} public class LDARule implements LogicRule {

	private double ruleWeight;
	private int N;
	private Corpus c;
	
	public LDARule(Corpus c, double ruleWeight)
	{
		this.c = c;
		this.ruleWeight = ruleWeight;
		this.N = c.N;
	}
	
	public double getRuleWeight()
	{
		return this.ruleWeight;
	}
	
	public double getTotalSamplingWeight() 
	{
		return N * ruleWeight;
	}
	
	public long numGroundings()
	{
		return N;
	}

	public String toString() 
	{
		return String.format("LDA pseudorule (corpus length = %d)\n", N);
	}

	/**
	 * Calculate LDA gradient for a given index i
	 * @param i 
	 * @return
	 */
	public double[][] ldaGradient(Corpus c, double[][] phi, double[][] theta,
			int i)
	{
		int T = phi.length;
		double[][] gradient = new double[1][T];
		for(int t = 0; t < T; t++)
			gradient[0][t] = phi[t][c.w[i]] * theta[c.d[i]][t];
		return gradient;
	}
	
	public Gradient randomGradient(RelaxedSample relax, Random rng) 
	{
		int[] indices = new int[1];
		indices[0] = rng.nextInt(N);								
		return new Gradient(ldaGradient(c, relax.phi, relax.theta, indices[0]), 
				indices); 
	}

	public LogicRule[] docPartition(Corpus c, HashMap<Integer,Integer> docMap)
	{
		// TODO: Implement LDA rule splitting
		// (add another layer of indirection)
		return null;
	}
	
	/**
	 * Number of satisfied groundings is undefined for LDARule
	 */
	public int numSat(int[] z)
	{
		return 0;
	}
	
	public void applyEvidence(Corpus c, int T)
	{	
		return;
	}
} public interface LogicRule 
{
	
	/**
	 * Weight of this rule for random sampling purposes
	 * @return [total sampling weight] = [sampling weight] x [num groundings] 
	 */
	public double getTotalSamplingWeight();
	
	/**
	 * Actual weight of this logic rule
	 * @return [rule weight] = [sampling weight] * [step weight]
	 */
	public double getRuleWeight();	
		 	
	/**
	 * Randomly sample a grounding and take a gradient step wrt to it
	 * @param relax Relaxed z-array (will be updated in-place!)
	 */
	public Gradient randomGradient(RelaxedSample relax, Random rng);
	
	/**
	 * 
	 * @param c
	 * @param T
	 */
	public void applyEvidence(Corpus c, int T);
	
	/**
	 * How many (non-trivial) groundings  
	 * @return
	 */
	public long numGroundings();
	
	/**
	 * Given a sample z, how many groundings are satisfied?
	 * @param z
	 * @return
	 */
	public int numSat(int[] z);
} public enum RuleType {

	SEED("SEED"),
	SENTEXCL("SENTEXCL"),
	SENTINCL("SENTINCL"),	
	CL("CL"),
	ML("ML"),
	DOCSEED("DOCSEED");
	
	private final String ruleName;
	
	RuleType(String name)
	{
		this.ruleName = name;		
	}
		
	public boolean matchesName(String str)
	{
		return this.ruleName.equalsIgnoreCase(str);
	}
} public class SeedRule implements IndependentRule, GroundableRule {

	private Seeds seedWords;
	private HashSet<Integer> hashSeedTopics;
	private int[] groundings;
	
	private double sampWeight;
	private double stepWeight;
	
	private int[] indices;
	private double[][] gradient;
	
	// below members are only used if rule is grounded as GroundableRule
	private Map<Integer, Set<Grounding>> invIndex;
	private Set<Grounding> unsat;
	
	private Set<Grounding> topicDupe;
	
	/**
	 * If wi in seed words, then zi in seed topics
	 *   
	 * @param argToks
	 */
	public SeedRule(double sampWeight, double stepWeight,
			Vector<String> argToks) {
	
		// Error check input
		assert (argToks.size() == 2);
	
		// Save weights 		
		this.sampWeight = sampWeight; // used for rule sampling 
		this.stepWeight = stepWeight; // used for gradient step-size
		
		// Get HashSet of seed words (or n-grams)
		//
		// Format:
		// 3434, 456 732, 563
		// means that our seed terms are:
		// 3434
		// 456 732 (only when appearing as bigram)
		// 563
		// 
		//
		seedWords = new Seeds(argToks.get(0));
		
		// Get HashSet of seed topics
		//
		Vector<Integer> vecSeedTopics = new Vector<Integer>();
		StringTokenizer seedtok = new StringTokenizer(argToks.get(1), ",");
		while (seedtok.hasMoreTokens())
			vecSeedTopics.add(Integer.parseInt(seedtok.nextToken()));
		hashSeedTopics = new HashSet<Integer>();
		for(Integer zi : vecSeedTopics)		
			hashSeedTopics.add(zi);		
		
		// Init groundings to null (will be set later by applyEvidence)
		this.indices = new int[1];
		this.groundings = null;
		
		// Init yet other members
		// (will be set by groundRule)
		// below members are only used if rule is grounded as GroundableRule
		this.invIndex = null;
		this.unsat = null;
	}

	public double getRuleWeight()
	{
		return this.sampWeight * this.stepWeight;
	}
	
	public double[][] toZLabel(int N, int T)
	{
		evidenceCheck("toZLabel()");

		 // entries *not* associated with weights will be null
		double[][] retval = new double[N][];
		
		for(int i : groundings)	
		{
			retval[i] = new double[T];		
			for(int t = 0; t < T; t++)			
				retval[i][t] = gradient[0][t] * sampWeight;
		}
		return retval;
	}
	
	public String toString()
	{
		return String.format("Seed Rule (samp weight=%.1f, step weight=%.1f)" +
				"\nWords = %s\nTopics = %s\n",
				sampWeight, stepWeight,
				seedWords.toString(),
				Arrays.toString(hashSeedTopics.toArray()));
	}
	
	/**
	 * Find non-trivial groundings of this rule
	 * @param c
	 * @param p
	 */
	public void applyEvidence(Corpus c, int T)
	{
		// Get all non-trival ground indices
		//
		HashSet<Integer> hashGround = new HashSet<Integer>();
		for(Integer wi = 0; wi < c.N; wi++)
		{
			int hit = seedWords.hit(c.w, wi);
			if(hit > 0)
				for(int hiti = 0; hiti < hit; hiti++)
					hashGround.add(wi + hiti);			
		}		
		 
		// Convert from HashSet to int[]
		//
		Integer[] boxground = new Integer[hashGround.size()];
		boxground = hashGround.toArray(boxground);	
		groundings = new int[hashGround.size()];
		for(int i = 0; i < groundings.length; i++)
			groundings[i] = boxground[i];
		
		// Pre-calc gradients calculations
		// (not dependent on zrelax so will never change)
		//
		gradient = new double[1][T];
		for(Integer t = 0; t < T; t++)		
			if(hashSeedTopics.contains(t))		
				gradient[0][t] = stepWeight;
			else
				gradient[0][t] = 0;		
	}
	
	//
	// Methods below should only be called after applyEvidence()
	//
	
	/**
	 * Ensure that evidence has been applied 
	 */	
	private void evidenceCheck(String methodname)
	{
		if(groundings == null)
		{
			String errmsg = String.format("ERROR: %s called without applying evidence", 
					methodname);
			System.out.println(errmsg);					
			System.exit(1);		
		}		
	}
	
	public double getTotalSamplingWeight() 
	{
		evidenceCheck("getWeight()");
		
		return (this.groundings.length) * this.sampWeight;		
	}

	public long numGroundings()
	{
		evidenceCheck("numGroundings()");
		
		return this.groundings.length;
	}
	
	public Gradient randomGradient(RelaxedSample relax, Random rng)
	{
		evidenceCheck("randomGradient()");
		
		// Sample a random grounding
		indices[0] = this.groundings[rng.nextInt(this.groundings.length)];		
		
		// Return random gradient
		return new Gradient(gradient, indices);		
	}
	
	/**
	 * Count satisfied groundings for this rule
	 * (easy to calc)
	 */
	public int numSat(int[] z)
	{
		evidenceCheck("numSat()");
				
		int sat = 0;
		for(int gi : this.groundings)
			if(hashSeedTopics.contains(z[gi]))
				sat += 1;
		return sat;
	}
	
	//
	// GroundableRule methods
	//

	/**
	 * Ensure that evidence has been applied 
	 */	
	private void groundCheck(String methodname)
	{
		if(this.invIndex == null)
		{
			String errmsg = String.format(
					"ERROR: %s called before grounding rule", 
					methodname);
			System.out.println(errmsg);					
			System.exit(1);		
		}		
	}
	
	/**
	 * For this rule, is a particular grounding satisfied?
	 * 
	 * @param z
	 * @param g
	 * @return
	 */
	private boolean groundingSat(int[] z, Grounding g)
	{
		return(this.hashSeedTopics.contains(z[g.get(0)]));		
	}
	
	public void groundRule(int[] z)
	{
		// Must have applied evidence before we can 
		// generate non-trivial groundings
		evidenceCheck("groundRule()");
		
		// Init data struct
		this.invIndex = new HashMap<Integer, Set<Grounding>>();
		this.unsat = new HashSet<Grounding>();
		
		for(int idxa : this.groundings)
		{
			// Ensure we have inverted index entries for idxa			
			if(!this.invIndex.containsKey(idxa))
				this.invIndex.put(idxa, new HashSet<Grounding>());				
			
			// Add this grounding to each entry
			Grounding newg = new Grounding(idxa);
			this.invIndex.get(idxa).add(newg);
				
			// Initialize unsat
			if(!groundingSat(z, newg))
				this.unsat.add(newg);				
		}
	}
	
	public void groundPenalty(int T)
	{
		groundCheck("groundPenalty()");
		
		// This will contain T copies of each grounding
		// (one per topic)
		topicDupe = new HashSet<Grounding>();
		
		// Get the first values element
		Iterator<Set<Grounding>> iterGround = 
			this.invIndex.values().iterator();
		Set<Grounding> allGround = iterGround.next();
		
		for(Grounding g : allGround)		
			for(int ti = 0; ti < T; ti++)			
				topicDupe.add(new Grounding(g.values[0], ti));
	}
	
	public double evalAssign(int[] z, int idx)
	{
		groundCheck("evalAssign()");
		
		if(!this.invIndex.containsKey(idx))
			return 0;
		
		double satweight = 0;
		for(Grounding g : this.invIndex.get(idx))
		{
			if(groundingSat(z, g))
				satweight += this.sampWeight * this.stepWeight;
		}
		return satweight;
	}

	public Map<Integer, Set<Grounding>> getInvIndex()
	{
		groundCheck("getInvIndex()");
		return this.invIndex;
	}
		
	public Set<Grounding> getUnSat()
	{
		groundCheck("getUnSat()");
		return this.unsat;
	}
	
	public void updateUnSat(int[] z, int idx)
	{			
		groundCheck("updateUnSat()");
		
		if(!this.invIndex.containsKey(idx))
			return;
		
		for(Grounding g : this.invIndex.get(idx))
		{
			if(!groundingSat(z,g))			
				this.unsat.add(g);
			else
				this.unsat.remove(g);
		}
	}

} public class SentExclRule implements LogicRule {

	private int excluder;
	private int excludee;
	
	private double sampWeight;
	private double stepWeight;
	
	// Maps sentence idx --> list of corpus idxs in that sentence
	private HashMap<Integer, ArrayList<Integer>> sentences;
	private int[] sentWeightIdx;
	private double[] sentWeights; 
	private long numGround;
	
	private int[] indices;
	private int T;
	
	public SentExclRule(double sampWeight, double stepWeight,
			Vector<String> argToks) 
	{	
		// Error check input
		assert (argToks.size() == 2);
	
		// Save weights 		
		this.sampWeight = sampWeight; // used for rule sampling 
		this.stepWeight = stepWeight; // used for gradient step-size
		
		// Get "excluder" topic
		//
		this.excluder = new Integer(argToks.get(0));
		
		// Get "excludee" topics
		//
		this.excludee = new Integer(argToks.get(1));
		
		// Init sentences to null (will be set later by applyEvidence)
		//
		this.indices = new int[2];		
		this.sentences = null;
	}
	
	public double getRuleWeight()
	{
		return this.sampWeight * this.stepWeight;
	}
	
	public String toString()
	{
		return String.format("Sentence Exclusion Rule " +
				"(samp weight=%.1f, step weight=%.1f)" +
				"\nExcluder = %d\nExcludee = %d\n",
				sampWeight, stepWeight, excluder, excludee);				
	}
	
	public void applyEvidence(Corpus c, int T) {
		String sentStr = SideInfoType.SENTENCE.infoName;
		assert(c.sideInfo.containsKey(sentStr));
		assert(excluder < T && excludee < T);
		
		// Not sure how to make this cast type-safe?
		//
		// Idea: put (int[]) type in SideInfoType SENTENCE enum,
		// (but I don't know how to do that in Java...)
		int[] sentidx = (int[]) c.sideInfo.get(sentStr); 
		
		// Build hash: [sentence idx] --> [associated corpus idx]
		//
		sentences = new HashMap<Integer, ArrayList<Integer>>();
		for(int i = 0; i < sentidx.length; i++)
		{
			int si = sentidx[i];
			
			if(!sentences.containsKey(si))
				sentences.put(si, new ArrayList<Integer>());			
			sentences.get(si).add(i);			
		}
		
		// We will sample groundings uniformly from sentence length squared
		//
		// sentWeightIdx is a new layer of indirection to allow non-sequential
		// sentence indices (eg, 3 3 3 4 4 7 7 0 0 0), for example in a 
		// cross-fold validation context
		//  
		sentWeightIdx = new int[sentences.keySet().size()];
		sentWeights = new double[sentences.keySet().size()];
		numGround = 0;
		
		int sidx = 0;
		for(Map.Entry<Integer, ArrayList<Integer>> entry : 
			sentences.entrySet())
		{
			double current = Math.pow(entry.getValue().size(), 2);
			
			sentWeightIdx[sidx] = entry.getKey();
			sentWeights[sidx] = current;
			
			numGround += current;
			sidx += 1;
		}				
		
		this.T = T;		
	}

	//
	// Methods below should only be called after applyEvidence()
	//
	
	/**
	 * Ensure that evidence has been applied 
	 */	
	private void evidenceCheck(String methodname)
	{
		if(sentences == null)
		{
			String errmsg = String.format("ERROR: %s called without applying evidence", 
					methodname);
			System.out.println(errmsg);					
			System.exit(1);		
		}		
	}
	
	public double getTotalSamplingWeight() 
	{
		evidenceCheck("getWeight()");
		return sampWeight * numGround;		
	}

	public Gradient randomGradient(RelaxedSample relax, Random rng) 
	{	
		evidenceCheck("randomGradient()");
		
		// Sample a random sentence
		//
		int si = this.sentWeightIdx[MiscUtil.multSample(rng, sentWeights, 
				(double) numGround)];
		
		// Sample a pair of corpus indices within this sentence
		//
		ArrayList<Integer> sentIdx = sentences.get(si);
		int excluderIdx = sentIdx.get(rng.nextInt(sentIdx.size()));
		int excludeeIdx = sentIdx.get(rng.nextInt(sentIdx.size()));
		indices[0] = excluderIdx;
		indices[1] = excludeeIdx;
		
		// Construct gradient
		//
		double[][] gradient = new double[2][this.T]; // default value is zero 
		gradient[0][excluder] = -1 * stepWeight * relax.zrelax[excludeeIdx][excludee];
		gradient[1][excludee] = -1 * stepWeight * relax.zrelax[excluderIdx][excluder];
		
		return new Gradient(gradient, indices);
	}
	
	public long numGroundings() 	
	{		
		evidenceCheck("numGroundings()");		
		return numGround;
	}

	public int numSat(int[] z) 
	{
		evidenceCheck("numSat()");		
		
		// Count number of unsatisifed groundings
		long numUnsat = 0;
		for(Map.Entry<Integer, ArrayList<Integer>> entry : 
			sentences.entrySet())
		{	
			int nExcluder = 0;
			int nExcludee = 0;
			for(int idx : entry.getValue())
				if(z[idx] == excluder)
					nExcluder += 1;
				else if(z[idx] == excludee)
					nExcludee += 1;
			numUnsat += nExcluder * nExcludee;			
		}
		
		return (int) (numGroundings() - numUnsat); 				
	}

} public class SentInclRule implements LogicRule 
{
	private int includer;
	private int includee;
	
	private double sampWeight;
	private double stepWeight;
	
	// Maps sentence idx --> list of corpus idxs in that sentence
	private HashMap<Integer, ArrayList<Integer>> sentences;
	private int[] sentidx;
	private long numGround;

	public SentInclRule(double sampWeight, double stepWeight,
			Vector<String> argToks) 
	{	
		// Error check input
		assert (argToks.size() == 2);
	
		// Save weights 		
		this.sampWeight = sampWeight; // used for rule sampling 
		this.stepWeight = stepWeight; // used for gradient step-size
		
		// Get "includer" topic
		//
		this.includer = new Integer(argToks.get(0));
		
		// Get "includee" topics
		//
		this.includee = new Integer(argToks.get(1));
		
		// Init sentences and sentidx to null 
		// (will be set later by applyEvidence)
		//		
		this.sentences = null;
	}
	
	public double getRuleWeight()
	{
		return this.sampWeight * this.stepWeight;
	}
	
	public String toString()
	{
		return String.format("Sentence Inclusion Rule " +
				"(samp weight=%.1f, step weight=%.1f)" +
				"\nIncluder = %d\nIncludee = %d\n",
				sampWeight, stepWeight, includer, includee);				
	}
	
	public void applyEvidence(Corpus c, int T) 
	{
		// Extract sentence information from Corpus side info
		String sentStr = SideInfoType.SENTENCE.infoName;
		assert(c.sideInfo.containsKey(sentStr));
		assert(includer < T && includee < T);
		
		// Not sure how to make this cast type-safe?
		//
		// Idea: put (int[]) type in SideInfoType SENTENCE enum,
		// (but I don't know how to do that in Java...)
		sentidx = (int[]) c.sideInfo.get(sentStr); 
		
		// Build hash: [sentence idx] --> [associated corpus idx]
		//
		sentences = new HashMap<Integer, ArrayList<Integer>>();
		for(int i = 0; i < sentidx.length; i++)
		{
			int si = sentidx[i];
			
			if(!sentences.containsKey(si))
				sentences.put(si, new ArrayList<Integer>());			
			sentences.get(si).add(i);			
		}
		
		// The number of groundings is simply equal to the corpus length
		numGround = sentidx.length;		
	}	
	
	/**
	 * Ensure that evidence has been applied 
	 */	
	private void evidenceCheck(String methodname)
	{
		if(sentences == null)
		{
			String errmsg = 
				String.format("ERROR: %s called without applying evidence",
						methodname);
			System.out.println(errmsg);					
			System.exit(1);		
		}		
	}
	
	
	public double getTotalSamplingWeight() 
	{
		evidenceCheck("getWeight()");				
		return numGround * sampWeight;
	}

	public long numGroundings() 
	{
		evidenceCheck("numGroundings()");
		return numGround;
	}

	public Gradient randomGradient(RelaxedSample relax, Random rng) 
	{
		evidenceCheck("randomGradient()");
		
		// Randomly sample any index (this will be the includee)
		//
		int incleeIdx = rng.nextInt(sentidx.length);
		
		// Get the indices of all other words in that sentence
		//
		ArrayList<Integer> coSentence = sentences.get(sentidx[incleeIdx]);
		
		// Init gradients and indices
		//
		int S = coSentence.size();		
		int T = relax.zrelax[0].length;
		int[] indices = MiscUtil.intListUnbox(coSentence);
		double[][] gradients = new double[S][T];
		
		// Populate gradient
		//
		// Poly: 1 - (1 - zjTr) (1 - zkTr) ... (ziTe)
		// 
		// (i is incleeIdx, Te is inclee topic, Tr is incler topic)
		//
		for(int i = 0; i < indices.length; i++)
		{
			// Which topic has a non-zero gradient for this index?
			int t = -1;
			if(indices[i] == incleeIdx)
			{
				t = this.includee;
				gradients[i][t] = -1 * stepWeight;
			}
			else
			{
				t = this.includer;
				gradients[i][t] = 1 * stepWeight;
			}
						
			// To calc gradient, sum entries for all entries OTHER THAN i			
			for(int j = 0; j < indices.length; j++)
			{
				if(j == i)
					continue;
				
				// Need to handle the includee index differently 
				if(indices[j] == incleeIdx)
				{
					gradients[i][t] *= relax.zrelax[indices[j]][includee];
				}
				else
				{
					gradients[i][t] *= (1 - relax.zrelax[indices[j]][includer]);
				}				
			}
		}
		
		return new Gradient(gradients, indices);
	}
	
	public int numSat(int[] z) 
	{
		evidenceCheck("numSat()");
		
		long numUnsat = 0;

		// Look for unsatisifed occurrences in each sentence
		for(Map.Entry<Integer, ArrayList<Integer>> entry : 
			sentences.entrySet())
		{	
			// Have we seen the includER in this sentence?
			boolean includerPresent = false;
			// How many includEEs have we seen in this sentence?
			int nIncludee = 0;
			
			for(int idx : entry.getValue())
				if(z[idx] == includer)
					includerPresent = true;
				else if(z[idx] == includee)
					nIncludee += 1;
			
			// If we did not see the includER, then
			// all includEEs represent rule violations
			if(!includerPresent)
				numUnsat += nIncludee;
		}
		
		return (int) (numGroundings() - numUnsat); 				
	}

	//
	// TODO: Implement GroundableRule so we can try all runs
	//
	
	
	
} public class CLRule implements LogicRule, GroundableRule 
{
	// The two word types to be Cannot-linked
	private int wordA;
	private int wordB;
	
	// Gradient sampling and stepsize weights 
	private double sampWeight;
	private double stepWeight;
	
	// indices of occurrences of word A and word B
	private ArrayList<Integer> idxA;
	private ArrayList<Integer> idxB;
	
	private int[] indices;
	private double[][] gradient;
	private int T;

	// below members are only used if rule is grounded as GroundableRule
	private Map<Integer, Set<Grounding>> invIndex;
	private Set<Grounding> unsat;

	private Set<Grounding> topicDupe;
	
	public CLRule(double sampWeight, double stepWeight,
			Vector<String> argToks)
	{
		// Error check input
		assert (argToks.size() == 2);
	
		// Save weights 		
		this.sampWeight = sampWeight; // used for rule sampling 
		this.stepWeight = stepWeight; // used for gradient step-size
		
		// Get the two Cannot-Linked words
		//
		this.wordA = new Integer(argToks.get(0));
		this.wordB = new Integer(argToks.get(1));
		
		// Init other members
		// (will be set by applyEvidence)
		//
		this.indices = null;
		this.gradient = null;
		this.idxA = null;
		this.idxB = null;
		
		// Init yet other members
		// (will be set by groundRule)
		// below members are only used if rule is grounded as GroundableRule
		this.invIndex = null;
		this.unsat = null;		
	}
	
	public double getRuleWeight()
	{
		return this.sampWeight * this.stepWeight;
	}
	
	/**
	 * Ensure that evidence has been applied 
	 */	
	private void evidenceCheck(String methodname)
	{
		if(this.idxA == null)
		{
			String errmsg = String.format("ERROR: %s called without applying evidence", 
					methodname);
			System.out.println(errmsg);					
			System.exit(1);		
		}		
	}
	
	/**
	 * Total weight of rule (for sampling purposes) 
	 */
	public double getTotalSamplingWeight() 
	{
		evidenceCheck("getWeight()");
		return this.sampWeight * this.numGroundings();
	}

	/**
	 * Sample a random rule grounding
	 */
	public Gradient randomGradient(RelaxedSample relax, Random rng) 
	{
		evidenceCheck("randomGradient()");
		// Sample a word A idx and a word B idx
		int ai = this.idxA.get(rng.nextInt(this.idxA.size()));
		int bi = this.idxB.get(rng.nextInt(this.idxB.size()));		
		indices[0] = ai;
		indices[1] = bi;
		
		// CL logic clause polynomial is 
		// 1 - (zi0*zj0 + zi1*zj1 + ... + zi(T-1)*zj(T-1))
		//
		for(int ti = 0; ti < this.T; ti++)
		{		
			// Construct the gradient
			gradient[0][ti] = -1 * stepWeight * relax.zrelax[bi][ti];
			gradient[1][ti] = -1 * stepWeight * relax.zrelax[ai][ti];
		}
		
		return new Gradient(gradient, indices);		
	}

	/**
	 * Given the corpus, find occurrences of word A and word B
	 */
	public void applyEvidence(Corpus c, int T) 
	{
		this.idxA = new ArrayList<Integer>();
		this.idxB = new ArrayList<Integer>();
		
		for(int idx = 0; idx < c.w.length; idx++)
		{
			if(c.w[idx] == this.wordA)
				this.idxA.add(idx);
			else if(c.w[idx] == this.wordB)
				this.idxB.add(idx);
		}
		
		// Init other fields
		this.T = T;
		this.indices = new int[2];
		this.gradient = new double[2][T];		
	}

	/**
	 * Nice summary of rule for logic report
	 */
	public String toString()
	{
		return String.format("Cannot-Link Rule " +
				"(samp weight=%.1f, step weight=%.1f)" +
				"\nWord A = %d\nWord B = %d\n",
				sampWeight, stepWeight, wordA, wordB);
	}
	
	/**
	 * How many (non-trivial) groundings does this rule have?
	 */
	public long numGroundings() 
	{
		evidenceCheck("numGroundings()");
		return this.idxA.size() * this.idxB.size();
	}

	/**
	 * For this z-assignment, how many groundings are satisfied?
	 */
	public int numSat(int[] z) 
	{

		// Get topic assignment counts for all word A idx 
		int[] za = new int[this.T]; // default value is zero
		for(int idx : this.idxA)
			za[z[idx]] += 1;
		
		// Get topic assignment counts for all word B idx 
		int[] zb = new int[this.T]; // default value is zero
		for(int idx : this.idxB)
			zb[z[idx]] += 1;
		
		// For each topic, all PAIRS of idxA-idxB entries assigned 
		// to that topic correspond to unsat groundings
		long numUnsat = 0;
		for(int zi = 0; zi < this.T; zi++)
			numUnsat += za[zi] * zb[zi];
		
		return (int) (numGroundings() - numUnsat);		
	}

	//
	// GroundableRule methods
	//

	/**
	 * Ensure that evidence has been applied 
	 */	
	private void groundCheck(String methodname)
	{
		if(this.invIndex == null)
		{
			String errmsg = String.format(
					"ERROR: %s called before grounding rule", 
					methodname);
			System.out.println(errmsg);					
			System.exit(1);		
		}		
	}
	
	/**
	 * For this rule, is a particular grounding satisfied?
	 * 
	 * @param z
	 * @param g
	 * @return
	 */
	private boolean groundingSat(int[] z, Grounding g)
	{
		int idxa = g.get(0);
		int idxb = g.get(1);
		return (z[idxa] != z[idxb]);
	}
	
	public void groundRule(int[] z)
	{
		// Must have applied evidence before we can 
		// generate non-trivial groundings
		evidenceCheck("groundRule()");
		
		// Init data struct
		this.invIndex = new HashMap<Integer, Set<Grounding>>();
		this.unsat = new HashSet<Grounding>();
		
		for(int idxa : this.idxA)
			for(int idxb : this.idxB)
			{
				// Ensure we have inverted index entries for idxa and idxb
				if(!this.invIndex.containsKey(idxa))
					this.invIndex.put(idxa, new HashSet<Grounding>());
				if(!this.invIndex.containsKey(idxb))
					this.invIndex.put(idxb, new HashSet<Grounding>());
				
				// Add this grounding to each entry
				Grounding newg = new Grounding(idxa, idxb);								
				this.invIndex.get(idxa).add(newg);
				this.invIndex.get(idxb).add(newg);
					
				// Initialize unsat
				if(!groundingSat(z, newg))
					this.unsat.add(newg);				
			}
	}

	public void groundPenalty(int T)
	{
		groundCheck("groundPenalty()");
		
		// This will contain T copies of each grounding
		// (one per topic)
		topicDupe = new HashSet<Grounding>();
		
		for(Set<Grounding> gset : this.invIndex.values())
			for(Grounding g : gset)
				for(int ti = 0; ti < T; ti++)			
					topicDupe.add(new Grounding(g.values[0], g.values[1], ti));						
	}
	
	public double evalAssign(int[] z, int idx)
	{
		groundCheck("evalAssign()");
		
		if(!this.invIndex.containsKey(idx))
			return 0;
		
		double satweight = 0;
		for(Grounding g : this.invIndex.get(idx))
		{
			if(groundingSat(z, g))
				satweight += this.sampWeight * this.stepWeight;
		}
		return satweight;
	}

	public Map<Integer, Set<Grounding>> getInvIndex()
	{
		groundCheck("getInvIndex()");
		return this.invIndex;
	}
		
	public Set<Grounding> getUnSat()
	{
		groundCheck("getUnSat()");
		return this.unsat;
	}
	
	public void updateUnSat(int[] z, int idx)
	{	
		groundCheck("updateUnSat()");
		
		System.out.println(String.format("dupesize=%d", 
				this.topicDupe.size()));
		
		if(!this.invIndex.containsKey(idx))
			return;
		
		for(Grounding g : this.invIndex.get(idx))
		{
			if(!groundingSat(z,g))			
				this.unsat.add(g);
			else
				this.unsat.remove(g);
		}
	}
} public class MLRule implements GroundableRule
{
	// The two word types to be Must-Linked
	private int wordA;
	private int wordB;

	// Gradient sampling and stepsize weights 
	private double sampWeight;
	private double stepWeight;

	// indices of occurrences of word A and word B
	private ArrayList<Integer> idxA;
	private ArrayList<Integer> idxB;

	private int[] indices;
	private double[][] gradient;
	private int T;

	// below members are only used if rule is grounded as GroundableRule
	private Map<Integer, Set<Grounding>> invIndex;
	private Set<Grounding> unsat;
	
	private Set<Grounding> topicDupe;
	
	public MLRule(double sampWeight, double stepWeight,
			Vector<String> argToks)
	{
		// Error check input
		assert (argToks.size() == 2);

		// Save weights 		
		this.sampWeight = sampWeight; // used for rule sampling 
		this.stepWeight = stepWeight; // used for gradient step-size

		// Get the two Must-Linked words
		//
		this.wordA = new Integer(argToks.get(0));
		this.wordB = new Integer(argToks.get(1));

		// Init other members
		// (will be set by applyEvidence)
		//
		this.indices = null;
		this.gradient = null;
		this.idxA = null;
		this.idxB = null;
		
		// Init yet other members
		// (will be set by groundRule)
		// below members are only used if rule is grounded as GroundableRule
		this.invIndex = null;
		this.unsat = null;		
	}

	
	public double getRuleWeight()
	{
		return this.sampWeight * this.stepWeight;
	}
	
	
	/**
	 * Ensure that evidence has been applied 
	 */	
	private void evidenceCheck(String methodname)
	{
		if(this.idxA == null)
		{
			String errmsg = String.format("ERROR: %s called without applying evidence", 
					methodname);
			System.out.println(errmsg);					
			System.exit(1);		
		}		
	}

	/**
	 * Total weight of rule (for sampling purposes) 
	 */
	public double getTotalSamplingWeight() 
	{
		evidenceCheck("getWeight()");
		return this.sampWeight * this.numGroundings();
	}

	/**
	 * Sample a random rule grounding
	 */
	public Gradient randomGradient(RelaxedSample relax, Random rng) 
	{
		evidenceCheck("randomGradient()");
		// Sample a word A idx and a word B idx
		int ai = this.idxA.get(rng.nextInt(this.idxA.size()));
		int bi = this.idxB.get(rng.nextInt(this.idxB.size()));		
		indices[0] = ai;
		indices[1] = bi;

		// ML logic clause to corresponding polynomial is
		// (zi0 => zj0) ^ (zi1 => zj1) 
		// 1 - ( (1-zi0)*zj0 + (1-zi1)*zj1 + ...)
		//
		// We need to choose which 'side' of our bidirectional implication
		// to take the gradient of (shouldn't really matter...)		
		//
		int bidir = rng.nextInt(2);
		if(bidir == 0)
		{
			for(int ti = 0; ti < this.T; ti++)
			{		
				// Construct the gradient zit => zjt
				this.gradient[0][ti] = 
					stepWeight * (relax.zrelax[bi][ti] - 1);
				this.gradient[1][ti] = 
					stepWeight * relax.zrelax[ai][ti];
			}
		}
		else
		{
			for(int ti = 0; ti < this.T; ti++)
			{		
				// Construct the gradient for zit <= zjt
				this.gradient[0][ti] = 
					stepWeight * relax.zrelax[bi][ti];
				this.gradient[1][ti] = 
					stepWeight * (relax.zrelax[ai][ti] - 1);
			}
		}	
				
		return new Gradient(gradient, indices);		
	}

	/**
	 * Given the corpus, find occurrences of word A and word B
	 */
	public void applyEvidence(Corpus c, int T) 
	{
		this.idxA = new ArrayList<Integer>();
		this.idxB = new ArrayList<Integer>();

		for(int idx = 0; idx < c.w.length; idx++)
		{
			if(c.w[idx] == this.wordA)
				this.idxA.add(idx);
			else if(c.w[idx] == this.wordB)
				this.idxB.add(idx);
		}

		// Init other fields
		this.T = T;
		this.indices = new int[2];
		this.gradient = new double[2][T];		
	}

	/**
	 * Nice summary of rule for logic report
	 */
	public String toString()
	{
		return String.format("Must-Link Rule " +
				"(samp weight=%.1f, step weight=%.1f)" +
				"\nWord A = %d\nWord B = %d\n",
				sampWeight, stepWeight, wordA, wordB);
	}
	
	/**
	 * How many (non-trivial) groundings does this rule have?
	 * 
	 * (could multiply by 2 to account for each side of bidirectional,
	 * but do not, in order to be consistent with Cannot-Link rule)
	 */
	public long numGroundings() 
	{
		evidenceCheck("numGroundings()");
		return (this.idxA.size() * this.idxB.size());
	}

	/**
	 * For this z-assignment, how many groundings are satisfied?
	 */
	public int numSat(int[] z) 
	{		
		// Get topic assignment counts for all word A idx 
		int[] za = new int[this.T]; // default value is zero
		for(int idx : this.idxA)
			za[z[idx]] += 1;

		// Get topic assignment counts for all word B idx 
		int[] zb = new int[this.T]; // default value is zero
		for(int idx : this.idxB)
			zb[z[idx]] += 1;

		// For each topic, all PAIRS of idxA-idxB entries assigned 
		// to that topic correspond to satisfied groundings
		long numSat = 0;
		for(int zi = 0; zi < this.T; zi++)
			numSat += za[zi] * zb[zi];

		return (int) numSat;
	}

	//
	// GroundableRule methods
	//

	/**
	 * Ensure that evidence has been applied 
	 */	
	private void groundCheck(String methodname)
	{
		if(this.invIndex == null)
		{
			String errmsg = String.format(
					"ERROR: %s called before grounding rule", 
					methodname);
			System.out.println(errmsg);					
			System.exit(1);		
		}		
	}
	
	/**
	 * For this rule, is a particular grounding satisfied?
	 * 
	 * @param z
	 * @param g
	 * @return
	 */
	private boolean groundingSat(int[] z, Grounding g)
	{
		int idxa = g.get(0);
		int idxb = g.get(1);
		return (z[idxa] == z[idxb]);
	}
	
	public void groundRule(int[] z)
	{
		// Must have applied evidence before we can 
		// generate non-trivial groundings
		evidenceCheck("groundRule()");
		
		// Init data struct
		this.invIndex = new HashMap<Integer, Set<Grounding>>();
		this.unsat = new HashSet<Grounding>();
		
		for(int idxa : this.idxA)
			for(int idxb : this.idxB)
			{
				// Ensure we have inverted index entries for idxa and idxb
				if(!this.invIndex.containsKey(idxa))
					this.invIndex.put(idxa, new HashSet<Grounding>());
				if(!this.invIndex.containsKey(idxb))
					this.invIndex.put(idxb, new HashSet<Grounding>());
				
				// Add this grounding to each entry
				Grounding newg = new Grounding(idxa, idxb);								
				this.invIndex.get(idxa).add(newg);
				this.invIndex.get(idxb).add(newg);
					
				// Initialize unsat
				if(!groundingSat(z, newg))
					this.unsat.add(newg);				
			}
	}
	
	public void groundPenalty(int T)
	{
		groundCheck("groundPenalty()");
		
		// This will contain T copies of each grounding
		// (one per topic)
		topicDupe = new HashSet<Grounding>();
		
		// Get the first values element
		Iterator<Set<Grounding>> iterGround = 
			this.invIndex.values().iterator();
		Set<Grounding> allGround = iterGround.next();
		
		for(Grounding g : allGround)		
			for(int ti = 0; ti < T; ti++)			
				topicDupe.add(new Grounding(g.values[0], g.values[1], ti));
	}
	
	public double evalAssign(int[] z, int idx)
	{
		groundCheck("evalAssign()");
		
		if(!this.invIndex.containsKey(idx))
			return 0;
		
		double satweight = 0;
		for(Grounding g : this.invIndex.get(idx))
		{
			if(groundingSat(z, g))
				satweight += this.sampWeight * this.stepWeight;
		}
		return satweight;
	}

	public Map<Integer, Set<Grounding>> getInvIndex()
	{
		groundCheck("getInvIndex()");
		return this.invIndex;
	}
		
	public Set<Grounding> getUnSat()
	{
		groundCheck("getUnSat()");
		return this.unsat;
	}
	
	public void updateUnSat(int[] z, int idx)
	{	
		groundCheck("updateUnSat()");
		
		if(!this.invIndex.containsKey(idx))
			return;
		
		for(Grounding g : this.invIndex.get(idx))
		{
			if(!groundingSat(z,g))			
				this.unsat.add(g);
			else
				this.unsat.remove(g);
		}
	}
} public class DocRule implements IndependentRule {

	private int docLabel;
	private HashSet<Integer> hashSeedTopics;
	private int[] groundings;
	
	private double sampWeight;
	private double stepWeight;
	
	private int[] indices;
	private double[][] gradient;
	
	// below members are only used if rule is grounded as GroundableRule
	private Map<Integer, Set<Grounding>> invIndex;
	private Set<Grounding> unsat;
	
	private Set<Grounding> topicDupe;
	
	/**
	 * If doc has label, then z-label words in doc
	 * to provided topics
	 */
	public DocRule(double sampWeight, double stepWeight,
			Vector<String> argToks)
	{
		// Error check input
		assert (argToks.size() == 2);

		// Save weights 		
		this.sampWeight = sampWeight; // used for rule sampling 
		this.stepWeight = stepWeight; // used for gradient step-size

		// Get docLabel of interest	
		docLabel = Integer.parseInt(argToks.get(0));

		// Get HashSet of seed topics
		Vector<Integer> vecSeedTopics = new Vector<Integer>();
		StringTokenizer seedtok = new StringTokenizer(argToks.get(1), ",");
		while (seedtok.hasMoreTokens())
			vecSeedTopics.add(Integer.parseInt(seedtok.nextToken()));
		hashSeedTopics = new HashSet<Integer>();
		for(Integer zi : vecSeedTopics)		
			hashSeedTopics.add(zi);		

		// Init groundings to null (will be set later by applyEvidence)
		this.indices = new int[1];
		this.groundings = null;			
		
		// Init yet other members
		// (will be set by groundRule)
		// below members are only used if rule is grounded as GroundableRule
		this.invIndex = null;
		this.unsat = null;
	}
	
	@Override
	public void applyEvidence(Corpus c, int T) 
	{
		// error-check 
		String doclabelStr = SideInfoType.DOCLABEL.infoName;
		assert(c.sideInfo.containsKey(doclabelStr));
		for(int t : hashSeedTopics)
			assert(t >= 0 && t < T);
		
		// get document labels 
		int[] doclabels = (int[]) c.sideInfo.get(doclabelStr); 

		// find all idx where corresponding doc has label doclabel
		List<Integer >groundList = new ArrayList<Integer>();
		for(int i = 0; i < c.d.length; i++)
			if(doclabels[c.d[i]] == docLabel)
				groundList.add(i);
		this.groundings = new int[groundList.size()];
		for(int i = 0; i < groundList.size(); i++)
			this.groundings[i] = groundList.get(i);
		
		// Pre-calc gradients calculations
		// (not dependent on zrelax so will never change)
		//
		this.gradient = new double[1][T];
		for(Integer t = 0; t < T; t++)		
			if(hashSeedTopics.contains(t))		
				gradient[0][t] = stepWeight;
			else
				gradient[0][t] = 0;	
	}
	
	/**
	 * Ensure that evidence has been applied 
	 */	
	private void evidenceCheck(String methodname)
	{
		if(this.groundings == null)
		{
			String errmsg = String.format("ERROR: %s called without applying evidence", 
					methodname);
			System.out.println(errmsg);					
			System.exit(1);		
		}		
	}	
	
	//
	// Methods below should only be called after applyEvidence()
	//	
	
	public double[][] toZLabel(int N, int T)
	{
		evidenceCheck("toZLabel()");

		// entries *not* associated with weights will be null
		double[][] retval = new double[N][];
		
		// Populate z-label weights
		for(int i : groundings)		
		{			
			retval[i] = new double[T];
			for(int t = 0; t < T; t++)			
				retval[i][t] = gradient[0][t] * sampWeight;
		}
		return retval;
	}
	
	@Override
	public double getTotalSamplingWeight() 
	{		
		evidenceCheck("getTotalSamplingWeight()");

		return this.sampWeight * this.groundings.length;
	}

	@Override
	public double getRuleWeight() 
	{
		evidenceCheck("getRuleWeight()");

		return this.sampWeight * this.stepWeight;
	}

	@Override
	public Gradient randomGradient(RelaxedSample relax, Random rng) 
	{
		evidenceCheck("randomGradient");

		// Sample a random grounding
		indices[0] = this.groundings[rng.nextInt(this.groundings.length)];		
		
		// Return random gradient
		return new Gradient(gradient, indices);
	}	

	@Override
	public long numGroundings() 
	{
		evidenceCheck("numGroundings()");

		return this.groundings.length;
	}

	@Override
	public int numSat(int[] z) 
	{
		evidenceCheck("numSat()");
				
		int sat = 0;
		for(int gi : this.groundings)
			if(hashSeedTopics.contains(z[gi]))
				sat += 1;
		return sat;		
	}

	//
	// GroundableRule methods
	//

	/**
	 * Ensure that evidence has been applied 
	 */	
	private void groundCheck(String methodname)
	{
		if(this.invIndex == null)
		{
			String errmsg = String.format(
					"ERROR: %s called before grounding rule", 
					methodname);
			System.out.println(errmsg);					
			System.exit(1);		
		}		
	}
	
	/**
	 * For this rule, is a particular grounding satisfied?
	 * 
	 * @param z
	 * @param g
	 * @return
	 */
	private boolean groundingSat(int[] z, Grounding g)
	{
		return(this.hashSeedTopics.contains(z[g.get(0)]));		
	}
	
	public void groundRule(int[] z)
	{
		// Must have applied evidence before we can 
		// generate non-trivial groundings
		evidenceCheck("groundRule()");
		
		// Init data struct
		this.invIndex = new HashMap<Integer, Set<Grounding>>();
		this.unsat = new HashSet<Grounding>();
		
		for(int idxa : this.groundings)
		{
			// Ensure we have inverted index entries for idxa			
			if(!this.invIndex.containsKey(idxa))
				this.invIndex.put(idxa, new HashSet<Grounding>());				
			
			// Add this grounding to each entry
			Grounding newg = new Grounding(idxa);
			this.invIndex.get(idxa).add(newg);
				
			// Initialize unsat
			if(!groundingSat(z, newg))
				this.unsat.add(newg);				
		}
	}
	
	public void groundPenalty(int T)
	{
		groundCheck("groundPenalty()");
		
		// This will contain T copies of each grounding
		// (one per topic)
		topicDupe = new HashSet<Grounding>();
		
		// Get the first values element
		Iterator<Set<Grounding>> iterGround = 
			this.invIndex.values().iterator();
		Set<Grounding> allGround = iterGround.next();
		
		for(Grounding g : allGround)		
			for(int ti = 0; ti < T; ti++)			
				topicDupe.add(new Grounding(g.values[0], ti));
	}
	
	public double evalAssign(int[] z, int idx)
	{
		groundCheck("evalAssign()");
		
		if(!this.invIndex.containsKey(idx))
			return 0;
		
		double satweight = 0;
		for(Grounding g : this.invIndex.get(idx))
		{
			if(groundingSat(z, g))
				satweight += this.sampWeight * this.stepWeight;
		}
		return satweight;
	}

	public Map<Integer, Set<Grounding>> getInvIndex()
	{
		groundCheck("getInvIndex()");
		return this.invIndex;
	}
		
	public Set<Grounding> getUnSat()
	{
		groundCheck("getUnSat()");
		return this.unsat;
	}
	
	public void updateUnSat(int[] z, int idx)
	{			
		groundCheck("updateUnSat()");
		
		if(!this.invIndex.containsKey(idx))
			return;
		
		for(Grounding g : this.invIndex.get(idx))
		{
			if(!groundingSat(z,g))			
				this.unsat.add(g);
			else
				this.unsat.remove(g);
		}
	}

	
}
import logiclda.rules.SentExclRule;
import logiclda.rules.SentInclRule;
import logiclda.rules.CLRule;
import logiclda.rules.MLRule;
import logiclda.rules.DocRule;
import org.ujmp.core.exceptions.MatrixException;
import java.util.*;
import java.io.*;
import logiclda.StandardLDA;
import logiclda.infer.CollapsedGibbs;
import logiclda.infer.DiscreteSample;
import logiclda.infer.RelaxedSample;
import logiclda.infer.MirrorDescent;
import logiclda.rules.LDARule;
import logiclda.rules.LogicRule;
import logiclda.rules.RuleType;
import logiclda.rules.SeedRule;
/**
    LogicLDA - Topic modeling with First-Order Logic (FOL) domain knowledge        
    Copyright (C) 2010, David Andrzejewski (andrzeje@cs.wisc.edu)

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

package logiclda;




public class LogicLDA {

	/**
	 * @param args
	 * @throws IOException 
	 * @throws MatrixException 
	 */
	public static void main(String[] args) throws MatrixException, IOException 
	{
		// Parse command-line args
		//
		String basefn = args[0];
		int numsamp = Integer.parseInt(args[1]);
		int numouter= Integer.parseInt(args[2]); 
		int numinner = Integer.parseInt(args[3]); 
		int topN = Integer.parseInt(args[4]); 
		int randseed = Integer.parseInt(args[5]); 

		// Load corpus and parameters (checking vocab dim agreement)
		//
		LDAParameters p = new LDAParameters(basefn, randseed);
		Corpus c = new Corpus(basefn);		
		assert(p.W == c.W);
				
		// Load logic rules 
		// (ldaRule controls whether to include LDA phi/theta rule)
		//
		boolean ldaRule = true;
		MirrorDescent rs = constructRuleSet(basefn, c, p.T, randseed, ldaRule);		
		System.out.println(rs.toString());
		
		// Initialize Logic LDA
		//
		File init = new File(String.format("%s.%s", basefn, "init"));
		DiscreteSample s;
		if(init.exists())
		{
			// from *.init file
			s = new DiscreteSample(c.N, p.T, c.W, c.D, init.getCanonicalPath(), c);
		}
		else
		{
			// If we are going to run Logic SGD inference later, 
			// run Logic Collapsed Gibbs for numsamp
			if(numouter > 0)
			{
				double[][] logicweights = rs.seedsToZL(c.N, p.T);
				s = CollapsedGibbs.doLogicGibbs(logicweights, c, p, numsamp);
			}
			else
			{
				// else assume we were meant to use no logic whatsoever...
				s = StandardLDA.runStandardLDA(c, p, numsamp);
			}
		}
		
		// Run LogicLDA MAP inference
 		//
 		RelaxedSample relax = runLogicLDA(c, p, rs, s, numouter, numinner);

		// Write out results
		//
		relax.writePhiTheta(p, basefn);
		relax.writeSample(basefn);
		rs.satReport(relax.getZ(), basefn);
		c.writeTopics(basefn, relax.getPhi(p), topN);
	}
	
	/**
	 * Run LogicLDA MAP inference, starting from init DiscreteSample
	 * 
	 * @param c
	 * @param p
	 * @param rs
	 * @param s
	 * @param numouter
	 * @param numinner
	 * @return
	 */
	public static RelaxedSample runLogicLDA(Corpus c, LDAParameters p, MirrorDescent rs, 
			DiscreteSample s, int numouter, int numinner)
	{				
		// Init relaxed z-sample and RuleSet
		RelaxedSample relax = new RelaxedSample(c, p, s);						
		
		// Do LogicLDA MAP inference via Stochastic Gradient Descent				
		double stepa = Math.sqrt(numinner);
		double stepb = (double) numinner;
		relax = rs.doSGD(c, p, relax, numouter, numinner, stepa, stepb);
		
		// Return final RelaxedSample
		return relax;
	}

	/**
	 * Build rule set from *.rules file
	 * 
	 * @param basefn
	 * @param c
	 * @param T
	 * @param randseed
	 * @param ldaRule
	 * @return
	 */
	public static MirrorDescent constructRuleSet(String basefn, Corpus c, 
			int T, int randseed, boolean ldaRule)
	{		
		// Read logic rules from *.rules file
		Vector<LogicRule> rules = readRules(String.format("%s.rules",basefn));
		// Add LDA rule
		if(ldaRule)
			rules.add(new LDARule(c, 1));			
		// Apply the evidence
 		for(LogicRule lr : rules)
			lr.applyEvidence(c, T);
		// Create RuleSet
 		return new MirrorDescent(rules, new Random(randseed + 1));	
	}
	
	/**
	 * Static factory for constructing rules from file
	 */
	public static Vector<LogicRule> readRules(String rulefn)
	{
		Vector<LogicRule> rules = new Vector<LogicRule>();
				
		try
	    {
			BufferedReader in = 
				new BufferedReader(new FileReader(rulefn));
			String curLine = in.readLine();
			while(curLine != null)
		    {				
				Vector<String> ruleToks= new Vector<String>();
				
				StringTokenizer stok = new StringTokenizer(curLine,"-");
		    	while(stok.hasMoreTokens())
		    		ruleToks.add(stok.nextToken());
		    	
		    	// Pop sampling and step weights 
		    	double sampWeight = Double.parseDouble(ruleToks.remove(0));
		    	double stepWeight = Double.parseDouble(ruleToks.remove(0));
		    	// Rest of string contains rule parameters
		    	String rtStr = ruleToks.remove(0);
		    	
		    	for(RuleType rt : RuleType.values())
		    	{
		    		if(rt.matchesName(rtStr))
		    		{
		    			switch(rt) 
		    			{
		    				case SEED:
		    					rules.add(new SeedRule(sampWeight, stepWeight,