diff --git a/hippyflow/collectives/collective.py b/hippyflow/collectives/collective.py
index da857f7..acc5015 100644
--- a/hippyflow/collectives/collective.py
+++ b/hippyflow/collectives/collective.py
@@ -87,7 +87,7 @@ def allReduce(self, v, op):
             self._allReduce_array(v_array, op)
             return v_array[0]
                 
-        elif type(v) in [int, np.int, np.int32]:
+        elif type(v) in [int, np.int32]:
             v_array = np.array([v], dtype=np.int32)
             self._allReduce_array(v_array, op)
             return v_array[0]
@@ -126,7 +126,7 @@ def bcast(self, v, root = 0):
         broadcasted lives.
         """
         
-        if type(v) in [float, np.float64,int, np.int, np.int32]:
+        if type(v) in [float, np.float64, int, np.int32]:
             v_array = np.array([v])
             self.comm.Bcast(v_array,root = root)
             return v_array[0]
diff --git a/hippyflow/modeling/KLEProjector.py b/hippyflow/modeling/KLEProjector.py
index 7ce8fed..c98280d 100644
--- a/hippyflow/modeling/KLEProjector.py
+++ b/hippyflow/modeling/KLEProjector.py
@@ -14,10 +14,12 @@
 import dolfin as dl
 import numpy as np
 import os
-from hippylib import *
 from mpi4py import MPI 
 import time
 
+import hippylib as hp
+
+from ..collectives.collective import NullCollective
 from ..collectives.collectiveOperator import CollectiveOperator
 from ..collectives.comm_utils import checkMeshConsistentPartitioning
 from .jacobian import *
@@ -30,16 +32,14 @@ def KLEParameterList():
 	This function implements the parameter list for the KLE object
 	"""
 	parameters = {}
-	parameters['sample_per_process'] 	= [100, 'Number of samples per process']
 	parameters['error_test_samples'] 		= [50, 'Number of samples for error test']
-	parameters['rank'] 				 	= [128, 'Rank of subspace']
-	parameters['oversampling'] 		 	= [10, 'Oversampling parameter for randomized algorithms']
-	parameters['verbose']				= [True, 'Boolean for printing']
-
+	parameters['rank'] 				 		= [128, 'Rank of subspace']
+	parameters['oversampling'] 		 		= [10, 'Oversampling parameter for randomized algorithms']
+	parameters['verbose']					= [True, 'Boolean for printing']
 	parameters['output_directory']			= [None,'output directory for saving arrays and plots']
 	parameters['plot_label_suffix']			= ['', 'suffix for plot label']
 
-	return ParameterList(parameters)
+	return hp.ParameterList(parameters)
 
 class MRinvM:
 	"""
@@ -106,10 +106,10 @@ def random_input_projector(self):
 		"""
 		m_KLE = dl.Vector(self.mesh_constructor_comm)
 		self.prior.M.init_vector(m_KLE,0)
-		Omega = MultiVector(m_KLE,self.parameters['rank'] + self.parameters['oversampling'])
+		Omega = hp.MultiVector(m_KLE,self.parameters['rank'] + self.parameters['oversampling'])
 
 		if self.collective.rank() == 0:
-			parRandom.normal(1.,Omega)
+			hp.parRandom.normal(1.,Omega)
 			Omega.orthogonalize()
 		else:
 			Omega.zero()
@@ -133,22 +133,22 @@ def construct_input_subspace(self,M_orthogonal = True):
 
 		m_KLE = dl.Vector(self.mesh_constructor_comm)
 		KLE_Operator.init_vector(m_KLE,0)
-		Omega = MultiVector(m_KLE,self.parameters['rank'] + self.parameters['oversampling'])
+		Omega = hp.MultiVector(m_KLE,self.parameters['rank'] + self.parameters['oversampling'])
 
 		if self.collective.rank() == 0:
-			parRandom.normal(1.,Omega)
+			hp.parRandom.normal(1.,Omega)
 		else:
 			Omega.zero()
 
 		self.collective.bcast(Omega,root = 0)
 
 		if M_orthogonal:
-			self.d_KLE, self.V_KLE = doublePassG(KLE_Operator,\
+			self.d_KLE, self.V_KLE = hp.doublePassG(KLE_Operator,\
 				self.prior.M, self.prior.Msolver, Omega,self.parameters['rank'],s=1)
 			self.M_orthogonal = True
 		else:
-			RsolverOperator = Solver2Operator(self.prior.Rsolver)
-			self.d_KLE, self.V_KLE = doublePass(RsolverOperator, Omega,self.parameters['rank'],s=1)
+			RsolverOperator = hp.Solver2Operator(self.prior.Rsolver)
+			self.d_KLE, self.V_KLE = hp.doublePass(RsolverOperator, Omega,self.parameters['rank'],s=1)
 			self.M_orthogonal = False
 
 		self._subspace_construction_time = time.time() - t0
@@ -202,15 +202,15 @@ def test_errors(self, ranks = [None],cut_off = 1e-12):
 		projection_vector = dl.Vector(self.mesh_constructor_comm)
 		self.prior.init_vector(projection_vector,0)
 
-		LocalParameters = MultiVector(projection_vector,self.parameters['error_test_samples'])
+		LocalParameters = hp.MultiVector(projection_vector,self.parameters['error_test_samples'])
 		LocalParameters.zero()
 		# Generate samples
 		for i in range(self.parameters['error_test_samples']):
 			t0 = time.time()
-			parRandom.normal(1,self.noise)
+			hp.parRandom.normal(1,self.noise)
 			self.prior.sample(self.noise,LocalParameters[i])
 
-		LocalErrors = MultiVector(projection_vector,self.parameters['error_test_samples'])
+		LocalErrors = hp.MultiVector(projection_vector,self.parameters['error_test_samples'])
 		
 
 		for rank_index,rank in enumerate(ranks):
@@ -227,7 +227,7 @@ def test_errors(self, ranks = [None],cut_off = 1e-12):
 			if self.M_orthogonal:
 				InputProjectorOperator = PriorPreconditionedProjector(V_KLE,self.prior.M, input_init_vector_lambda)
 			else:
-				InputProjectorOperator = LowRankOperator(np.ones_like(d_KLE),V_KLE, input_init_vector_lambda)
+				InputProjectorOperator = hp.LowRankOperator(np.ones_like(d_KLE),V_KLE, input_init_vector_lambda)
 		
 			rel_errors = np.zeros(LocalErrors.nvec())
 			for i in range(LocalErrors.nvec()):
diff --git a/hippyflow/modeling/PODProjector.py b/hippyflow/modeling/PODProjector.py
index 0ad4bfa..c429908 100644
--- a/hippyflow/modeling/PODProjector.py
+++ b/hippyflow/modeling/PODProjector.py
@@ -14,7 +14,7 @@
 import dolfin as dl
 import numpy as np
 import time
-from hippylib import *
+import hippylib as hp
 from mpi4py import MPI 
 import os
 
@@ -24,6 +24,8 @@
 from ..utilities.plotting import *
 from .priorPreconditionedProjector import PriorPreconditionedProjector
 
+CONTROL = 3
+
 
 def PODParameterList():
 	"""
@@ -39,14 +41,14 @@ def PODParameterList():
 	parameters['output_directory']			= [None,'output directory for saving arrays and plots']
 	parameters['plot_label_suffix']			= ['', 'suffix for plot label']
 
-	return ParameterList(parameters)
+	return hp.ParameterList(parameters)
 
 
 class PODProjector:
 	"""
 	Projector class based on proper orthogonal decomposition
 	"""
-	def __init__(self,observable, prior, mesh_constructor_comm = None ,collective = None, parameters = PODParameterList()):
+	def __init__(self,observable, prior, control_distribution = None, mesh_constructor_comm = None ,collective = None, parameters = PODParameterList()):
 		"""
 		Constructor
 			- :code:`observable` - object that implements the observable mapping :math:`m -> q(m)`
@@ -57,6 +59,8 @@ def __init__(self,observable, prior, mesh_constructor_comm = None ,collective =
 		"""
 		self.observable = observable
 		self.prior = prior
+		self.control_distribution = control_distribution
+
 		if mesh_constructor_comm is not None:
 			self.mesh_constructor_comm = mesh_constructor_comm
 		else:
@@ -77,8 +81,12 @@ def __init__(self,observable, prior, mesh_constructor_comm = None ,collective =
 		self.noise = dl.Vector(self.mesh_constructor_comm)
 		self.prior.init_vector(self.noise,"noise")
 
-		self.u = self.observable.generate_vector(STATE)
-		self.m = self.observable.generate_vector(PARAMETER)
+		self.u = self.observable.generate_vector(hp.STATE)
+		self.m = self.observable.generate_vector(hp.PARAMETER)
+		if self.control_distribution is None:
+			self.z = None
+		else:
+			self.z = self.observable.generate_vector(CONTROL)
 
 
 		self.d = None
@@ -90,15 +98,20 @@ def solve_at_mean(self):
 		"""
 		Solve the PDE at the mean
 		"""
-		try:
-			m_mean = self.prior.mean
-			self.u_at_mean = self.observable.problem.generate_state()
+		m_mean = self.prior.mean
+		self.u_at_mean = self.observable.problem.generate_state()
+		if self.control_distribution is not None:
+			if hasattr(self.control_distribution,'mean'):
+				z_mean = self.control_distribution.mean
+			else:
+				z_mean = self.observable.generate_vector(CONTROL)
+				self.control_distribution.sample(z_mean)
+			self.observable.problem.solveFwd(self.u_at_mean,[self.u_at_mean,m_mean,None,z_mean])
+		else:
 			self.observable.problem.solveFwd(self.u_at_mean,[self.u_at_mean,m_mean,None])
-		except:
-			self.u_at_mean = None
 
-	def generate_training_data(self,output_directory = 'data/',check_for_data = True,sequential = True,\
-																			compress_files = True):
+	def generate_training_data(self,check_for_data = True,sequential = True,\
+										compress_files = True):
 		"""
 		This method generates training data
 			- :code:`output_directory` - a string specifying the path to the directory where data
@@ -108,107 +121,207 @@ def generate_training_data(self,output_directory = 'data/',check_for_data = True
 		"""
 		self.solve_at_mean()
 		my_rank = int(self.collective.rank())
-		try:
-			os.makedirs(output_directory)
-		except:
-			pass
+		output_directory = self.parameters['output_directory']
+
+		os.makedirs(output_directory,exist_ok = True)
+
 		observable_vector = dl.Vector(self.mesh_constructor_comm)
 		self.observable.init_vector(observable_vector,dim = 0)
 		last_datum_generated = 0
-		m_shape = self.m.get_local().shape[0]
-		q_shape = observable_vector.get_local().shape[0]
-		print('m_shape = ',m_shape)
-		print('q_shape = ',q_shape)
+		parameter_dimension = self.m.get_local().shape[0]
+		output_dimension = observable_vector.get_local().shape[0]
+		print('dM = ',parameter_dimension)
+		print('dQ = ',output_dimension)
+		if self.control_distribution is not None:
+			control_dimension = self.z.get_local().shape[0]
+			print('dZ = ',control_dimension)
+		
 		if sequential:
 			rank_specific_directory = output_directory+'data_on_rank_'+str(my_rank)+'/'
 			os.makedirs(rank_specific_directory,exist_ok = True)
 			if check_for_data:
-				if os.path.isfile(rank_specific_directory+'ms_sample_'+str(my_rank)+'.npy') and \
-					os.path.isfile(rank_specific_directory+'qs_sample_'+str(my_rank)+'.npy'):
-					ms_generated = os.listdir(rank_specific_directory)
-					qs_generated = os.listdir(rank_specific_directory)
-					m_indices = [int(m_.split('m_sample_')[-1].split('.npy')[0]) for m_ in ms_generated]
-					last_m = max(m_indices)
-					q_indices = [int(m_.split('q_sample_')[-1].split('.npy')[0]) for q_ in qs_generated]
-					last_q = max(q_indices)
-					last_datum_generated = min(last_m,last_q)
+				if self.control_distribution is not None:
+					if os.path.isfile(rank_specific_directory+'m_sample_0.npy') and \
+						os.path.isfile(rank_specific_directory+'q_sample_0.npy') and \
+						os.path.isfile(rank_specific_directory+'z_sample_0.npy'):
+						all_files = os.listdir(rank_specific_directory)
+						ms_generated = []
+						qs_generated = []
+						zs_generated = []
+						for file in all_files:
+							if 'm_' in file:
+								ms_generated.append(file)
+							elif 'q_' in file:
+								qs_generated.append(file)
+							elif 'z_' in file:
+								zs_generated.append(file)
+						m_indices = [int(m_.split('m_sample_')[-1].split('.npy')[0]) for m_ in ms_generated]
+						last_m = max(m_indices)
+						q_indices = [int(q_.split('q_sample_')[-1].split('.npy')[0]) for q_ in qs_generated]
+						last_q = max(q_indices)
+						z_indices = [int(z_.split('z_sample_')[-1].split('.npy')[0]) for z_ in zs_generated]
+						last_z = max(z_indices)
+						last_datum_generated = min(last_m,last_q,last_z)
+				else:
+					if os.path.isfile(rank_specific_directory+'m_sample_0.npy') and \
+						os.path.isfile(rank_specific_directory+'q_sample_0.npy'):
+						all_files = os.listdir(rank_specific_directory)
+						ms_generated = []
+						qs_generated = []
+						for file in all_files:
+							if 'm_' in file:
+								ms_generated.append(file)
+							elif 'q_' in file:
+								qs_generated.append(file)
+						m_indices = [int(m_.split('m_sample_')[-1].split('.npy')[0]) for m_ in ms_generated]
+						last_m = max(m_indices)
+						q_indices = [int(q_.split('q_sample_')[-1].split('.npy')[0]) for q_ in qs_generated]
+						last_q = max(q_indices)
+						last_datum_generated = min(last_m,last_q)
 
 			t0 = time.time()
 			for i in range(last_datum_generated,self.parameters['data_per_process']):
 				print('Generating data number '+str(i))
-				converged = False
-				while not converged:
-					parRandom.normal(1,self.noise)
-					self.prior.sample(self.noise,self.m)
-					self.u.zero()
-					if self.u_at_mean is not None:
-						self.u.axpy(1.,self.u_at_mean)
+				solved = False
+				while not solved:
 					try:
-						this_q = self.observable.eval(self.m).get_local()
+						self.m.zero()
+						self.noise.zero()
+						hp.parRandom.normal(1,self.noise)
+						self.prior.sample(self.noise,self.m)
+						if self.control_distribution is None:
+							linearization_x = [self.u,self.m,None]
+						else:
+							self.z.zero()
+							self.control_distribution.sample(self.z)
+							linearization_x = [self.u,self.m,None,self.z]
+
+						self.observable.solveFwd(self.u,linearization_x)
 						this_m = self.m.get_local()
+						this_q = self.observable.evalu(self.u).get_local()
+						
 						np.save(rank_specific_directory+'m_sample_'+str(i)+'.npy',this_m)
 						np.save(rank_specific_directory+'q_sample_'+str(i)+'.npy',this_q)
-						converged = True
-						# If there is an issue with the solve move on
+
+						if self.control_distribution is not None:
+							this_z = self.z.get_local()
+							np.save(rank_specific_directory+'z_sample_'+str(i)+'.npy',this_z)
+
+						solved = True
 					except:
-						print('Issue with the nonlinear solve, moving on')
-						pass
-					
-				if self.parameters['verbose']:
-					print('On datum generated every ',(time.time() -t0)/(i - last_datum_generated+1),' s, on average.')
+						print('Issue with the forward solution, moving on.')
+
+					if self.parameters['verbose']:
+						print('On datum generated every ',(time.time() -t0)/(i - last_datum_generated+1),' s, on average.')
+
 			self._data_generation_time = time.time() - t0
 			if compress_files:
-				local_ms = np.zeros((self.parameters['data_per_process'],m_shape))
-				local_qs = np.zeros((self.parameters['data_per_process'],q_shape))
+				local_ms = np.zeros((self.parameters['data_per_process'],parameter_dimension))
+				local_qs = np.zeros((self.parameters['data_per_process'],output_dimension))
+				if self.control_distribution is not None:
+					local_zs = np.zeros((self.parameters['data_per_process'],control_dimension))
 				for i in range(0,self.parameters['data_per_process']):
 					local_ms[i] = np.load(rank_specific_directory+'m_sample_'+str(i)+'.npy')
 					local_qs[i] = np.load(rank_specific_directory+'q_sample_'+str(i)+'.npy')
-				np.savez_compressed(output_directory+'mq_on_rank'+str(my_rank)+'.npz',m_data = local_ms,q_data = local_qs)
+					if self.control_distribution is not None:
+						local_zs[i] = np.load(rank_specific_directory+'z_sample_'+str(i)+'.npy')
+				if self.control_distribution is not None:
+					np.savez_compressed(output_directory+'mqz_on_rank'+str(my_rank)+'.npz',m_data = local_ms,\
+										q_data = local_qs, z_data = local_zs)
+				else:
+					np.savez_compressed(output_directory+'mq_on_rank'+str(my_rank)+'.npz',m_data = local_ms,q_data = local_qs)
 
 		else:
 
-			local_ms = np.zeros((0,m_shape))
-			local_qs = np.zeros((0,q_shape))
+			local_ms = np.zeros((0,parameter_dimension))
+			local_qs = np.zeros((0,output_dimension))
+			if self.control_distribution is not None:
+				local_zs = np.zeros((0, self.control_distribution))
+
 			if check_for_data:
 				if os.path.isfile(output_directory+'ms_on_rank_'+str(my_rank)+'.npy') and \
 					os.path.isfile(output_directory+'qs_on_rank_'+str(my_rank)+'.npy'):
 					local_ms = np.load(output_directory+'ms_on_rank_'+str(my_rank)+'.npy')
 					local_qs = np.load(output_directory+'qs_on_rank_'+str(my_rank)+'.npy')
-					last_datum_generated = min(local_ms.shape[0],local_qs.shape[0])
+					if self.control_distribution is not None:
+						local_zs = np.load(output_directory+'zs_on_rank_'+str(my_rank)+'.npy')
+						last_datum_generated = min(local_ms.shape[0],local_qs.shape[0],local_zs.shape[0])
+					else:
+						last_datum_generated = min(local_ms.shape[0],local_qs.shape[0])
 
 			t0 = time.time()
 			# I think this is all hard coded for a single serial mesh, check if 
 			# the arrays need to be communicated to mesh rank 0 before being saved
 			for i in range(last_datum_generated,self.parameters['data_per_process']):
-				print('Generating data number '+str(i))
-				parRandom.normal(1,self.noise)
-				self.prior.sample(self.noise,self.m)
-				
-				self.u.zero()
-				self.u.axpy(1.,self.u_at_mean)
-				x = [self.u,self.m,None]
-				self.observable.setLinearizationPoint(x)
-				solution = self.observable.eval(self.m).get_local()
-				# If there is an issue with the solve move on
-				# local_qs = np.concatenate((local_qs,np.expand_dims(solution,0)))
-				# local_ms = np.concatenate((local_ms,np.expand_dims(self.m.get_local(),0)))
-				# np.save(output_directory+'ms_on_rank_'+str(my_rank)+'.npy',np.array(local_ms))
-				# np.save(output_directory+'qs_on_rank_'+str(my_rank)+'.npy',np.array(local_qs))
-				try:
-					solution = self.observable.eval(self.m).get_local()
-					# If there is an issue with the solve move on
-					local_qs = np.concatenate((local_qs,np.expand_dims(solution,0)))
-					local_ms = np.concatenate((local_ms,np.expand_dims(self.m.get_local(),0)))
-					np.save(output_directory+'ms_on_rank_'+str(my_rank)+'.npy',np.array(local_ms))
-					np.save(output_directory+'qs_on_rank_'+str(my_rank)+'.npy',np.array(local_qs))
-				except:
-					print('Issue with the nonlinear solve, moving on')
-					pass
-				
+				solved = False
+				while not solved:
+					try:
+						print('Generating data number '+str(i))
+						hp.parRandom.normal(1,self.noise)
+						self.prior.sample(self.noise,self.m)
+						self.u.zero()
+						self.u.axpy(1.,self.u_at_mean)
+
+
+						if self.control_distribution is not None:
+							self.z.zero()
+							self.control_distribution.sample(self.z)
+							linearization_x = [self.u,self.m,None,self.z]
+						else:
+							linearization_x = [self.u,self.m,None]
+
+
+						self.observable.solveFwd(self.u,linearization_x)
+						this_m = self.m.get_local()
+						this_q = self.observable.evalu(self.u).get_local()
+						local_qs = np.concatenate((local_qs,np.expand_dims(this_q,0)))
+						local_ms = np.concatenate((local_ms,np.expand_dims(this_m,0)))
+						np.save(output_directory+'ms_on_rank_'+str(my_rank)+'.npy',np.array(local_ms))
+						np.save(output_directory+'qs_on_rank_'+str(my_rank)+'.npy',np.array(local_qs))
+						if self.control_distribution is not None:
+							this_z = self.z.get_local()
+							local_zs = np.concatenate((local_zs,np.expand_dims(this_z,0)))
+							np.save(output_directory+'zs_on_rank_'+str(my_rank)+'.npy',np.array(local_zs))
+						solved = True
+
+					except:
+						print('Issue with nonlinear solve, moving on')
 				if self.parameters['verbose']:
 					print('On datum generated every ',(time.time() -t0)/(i - last_datum_generated+1),' s, on average.')
 			self._data_generation_time = time.time() - t0
 
+	def save_mass_and_stiffness_matrices(self):
+		"""
+		This method saves mass and stiffness matrices 
+		"""
+		if MPI.COMM_WORLD.rank == 0:
+			# Save mass matrix
+			output_directory = self.parameters['output_directory']
+			os.makedirs(output_directory,exist_ok = True)
+			import scipy.sparse as sp
+
+			u_trial = dl.TrialFunction(self.observable.problem.Vh[hp.STATE])
+			u_test = dl.TestFunction(self.observable.problem.Vh[hp.STATE])
+			M = dl.PETScMatrix(self.mesh_constructor_comm)
+			dl.assemble(dl.inner(u_trial,u_test)*dl.dx, tensor=M)
+			
+			# from scipy.sparse import csc_matrix, csr_matrix, save_npz
+			# from scipy.sparse import linalg as spla
+
+			M_mat = dl.as_backend_type(M).mat()
+			row,col,val = M_mat.getValuesCSR()
+			M_csr = sp.csr_matrix((val,col,row)) 
+			sp.save_npz(output_directory+'mass_csr',M_csr)
+
+			# Save stiffness matrix
+			K = dl.PETScMatrix(self.mesh_constructor_comm)
+			dl.assemble(dl.inner(dl.grad(u_trial),dl.grad(u_test))*dl.dx, tensor=K)
+			K_mat = dl.as_backend_type(K).mat()
+			row,col,val = K_mat.getValuesCSR()
+			K_csr = sp.csr_matrix((val,col,row)) 
+			sp.save_npz(output_directory+'stiffness_csr',K_csr)
+
+
 
 	def construct_subspace(self):
 		"""
@@ -218,7 +331,7 @@ def construct_subspace(self):
 		self.solve_at_mean()
 		observable_vector = dl.Vector(self.mesh_constructor_comm)
 		self.observable.init_vector(observable_vector,dim = 0)
-		LocalObservables = MultiVector(observable_vector,self.parameters['sample_per_process'])
+		LocalObservables = hp.MultiVector(observable_vector,self.parameters['sample_per_process'])
 		LocalObservables.zero()
 
 		#Read data from file and build subspace option
@@ -226,32 +339,36 @@ def construct_subspace(self):
 			# if self.parameters['verbose']:
 				# print('Starting observable generation for draw ',i)
 			observable_vector.zero()
-			parRandom.normal(1,self.noise)
+			hp.parRandom.normal(1,self.noise)
 			self.prior.sample(self.noise,self.m)
-			x = [self.u,self.m,None]
-			self.observable.setLinearizationPoint(x)
-			observable_vector.axpy(1.,self.observable.eval(self.m))
+			if self.control_distribution is not None:
+				self.z.zero()
+				self.control_distribution.sample(self.z)
+				x = [self.u,self.m,None,self.z]
+			else:
+				x = [self.u,self.m,None]
+			self.observable.solveFwd(self.u,x)
+			observable_vector.axpy(1.,self.observable.evalu(self.u))
 			LocalObservables[i].axpy(1.,observable_vector)
 
 		init_vector_lambda = lambda x, dim: self.observable.init_vector(x,dim = 0)
-
-		LocalPODOperator = LowRankOperator(np.ones(LocalObservables.nvec())/self.parameters['sample_per_process']\
+		LocalPODOperator = hp.LowRankOperator(np.ones(LocalObservables.nvec())/self.parameters['sample_per_process']\
 																			,LocalObservables,init_vector_lambda)
 
 		GlobalPODOperator = CollectiveOperator(LocalPODOperator, self.collective,mpi_op = 'avg')
 
 		x_POD = dl.Vector(self.mesh_constructor_comm)
 		LocalPODOperator.init_vector(x_POD,dim = 0)
-		Omega_POD = MultiVector(x_POD,self.parameters['rank'] + self.parameters['oversampling'])
+		Omega_POD = hp.MultiVector(x_POD,self.parameters['rank'] + self.parameters['oversampling'])
 
 		if self.collective.rank() == 0:
-			parRandom.normal(1.,Omega_POD)
+			hp.parRandom.normal(1.,Omega_POD)
 		else:
 			Omega_POD.zero()
 
 		self.collective.bcast(Omega_POD,root = 0)
 
-		self.d, self.U_MV = doublePass(GlobalPODOperator,Omega_POD,self.parameters['rank'],s = 1)
+		self.d, self.U_MV = hp.doublePass(GlobalPODOperator,Omega_POD,self.parameters['rank'],s = 1)
 
 		self._subspace_construction_time = time.time() - t0
 		if self.parameters['verbose'] and (self.mesh_constructor_comm.rank ==0):
@@ -273,6 +390,7 @@ def test_output_errors(self,ranks = [None],cut_off = 1e-10):
 			- :code:`ranks` - a python list of integers specifying ranks for projection tests
 			- :code:`cut_off` - where to truncate the ranks based on the spectral decay of the POD operator
 		"""
+		assert self.control_distribution is None, 'Not worked out yet for control problems'
 		ranks.sort()
 		if (self.d is None) or (self.U_MV is None):
 			if self.mesh_constructor_comm.rank == 0:
@@ -290,19 +408,24 @@ def test_output_errors(self,ranks = [None],cut_off = 1e-10):
 
 		observable_vector = dl.Vector(self.mesh_constructor_comm)
 		self.observable.init_vector(observable_vector,dim = 0)
-		LocalObservables = MultiVector(observable_vector,self.parameters['sample_per_process'])
+		LocalObservables = hp.MultiVector(observable_vector,self.parameters['sample_per_process'])
 		LocalObservables.zero()
 		for i in range(LocalObservables.nvec()):
 			t0 = time.time()
-			parRandom.normal(1,self.noise)
+			hp.parRandom.normal(1,self.noise)
 			self.prior.sample(self.noise,self.m)
-			x = [self.u,self.m,None]
-			self.observable.setLinearizationPoint(x)
-			LocalObservables[i].axpy(1.,self.observable.eval(self.m))
+			if self.control_distribution is not None:
+				self.z.zero()
+				self.control_distribution.sample(self.z)
+				x = [self.u,self.m,None,self.z]
+			else:
+				x = [self.u,self.m,None]
+			self.observable.solveFwd(self.u,x)
+			LocalObservables[i].axpy(1.,self.observable.evalu(self.u))
 			# if self.parameters['verbose'] and (self.mesh_constructor_comm.rank == 0):
 			# 	print('Generating local observable ',i,' for POD error test took',time.time() -t0, 's')
 
-		LocalErrors = MultiVector(observable_vector,self.parameters['sample_per_process'])
+		LocalErrors = hp.MultiVector(observable_vector,self.parameters['sample_per_process'])
 
 		projection_vector = dl.Vector(self.mesh_constructor_comm)
 		self.observable.init_vector(projection_vector,dim = 0)
@@ -314,7 +437,7 @@ def test_output_errors(self,ranks = [None],cut_off = 1e-10):
 				U_MV = self.U_MV
 				d = self.d
 			else:
-				U_MV = MultiVector(self.U_MV[0],rank)
+				U_MV = hp.MultiVector(self.U_MV[0],rank)
 				d = self.d[0:rank]
 				for i in range(rank):
 					U_MV[i].axpy(1.,self.U_MV[i])
@@ -322,7 +445,7 @@ def test_output_errors(self,ranks = [None],cut_off = 1e-10):
 					# U_MV[i].apply('')
 
 			init_vector_lambda = lambda x, dim: self.observable.init_vector(x,dim = 0)
-			PODOperator = LowRankOperator(np.ones_like(d),U_MV,init_vector_lambda)
+			PODOperator = hp.LowRankOperator(np.ones_like(d),U_MV,init_vector_lambda)
 
 			rel_errors = np.zeros(LocalErrors.nvec())
 			for i in range(LocalErrors.nvec()):
@@ -360,33 +483,53 @@ def two_state_solution(self):
 		if self.parameters['verbose']:
 			print('||m_mean|| = ',m_mean.norm('l2'))
 		m_mean_pvd = dl.File(save_states_dir+'m_mean.pvd')
-		m_mean_pvd << vector2Function(m_mean,self.observable.problem.Vh[PARAMETER])
-
+		m_mean_pvd << hp.vector2Function(m_mean,self.observable.problem.Vh[hp.PARAMETER])
 		u_at_mean = self.observable.problem.generate_state()
-		self.observable.problem.solveFwd(u_at_mean,[u_at_mean,m_mean,None])
+
+		if self.control_distribution is not None:
+			if hasattr(self.control_distribution,'mean'):
+				z_mean = self.control_distribution.mean
+			else:
+				print('Substituting a sample of z for z_mean, since ')
+				print('control_distribution did not have a mean attribute.')
+				z_mean = self.observable.generate_vector(CONTROL)
+				self.control_distribution.sample(z_mean)
+
+			linearization_x = [u_at_mean,m_mean,None,z_mean]
+		else:
+			linearization_x = [u_at_mean,m_mean,None]
+
+		self.observable.problem.solveFwd(u_at_mean,linearization_x)
 
 		if self.parameters['verbose']:
 			print('||v_at_mean|| = ',u_at_mean.norm('l2'))
 		v_at_mean_pvd = dl.File(save_states_dir+'v_at_mean.pvd')
-		v_at_mean_pvd << vector2Function(u_at_mean,self.observable.problem.Vh[STATE])
+		v_at_mean_pvd << hp.vector2Function(u_at_mean,self.observable.problem.Vh[hp.STATE])
 
 		# Sample from prior:
-		parRandom.normal(1,self.noise)
-		m_sample = self.observable.generate_vector(PARAMETER)
+		hp.parRandom.normal(1,self.noise)
+		m_sample = self.observable.generate_vector(hp.PARAMETER)
 		self.prior.sample(self.noise,m_sample)
 
 		if self.parameters['verbose']:
 			print('||m_sample|| = ',m_sample.norm('l2'))
 		m_sample_pvd = dl.File(save_states_dir+'m_sample.pvd')
-		m_sample_pvd << vector2Function(m_sample,self.observable.problem.Vh[PARAMETER])
+		m_sample_pvd << hp.vector2Function(m_sample,self.observable.problem.Vh[hp.PARAMETER])
 
 		u_at_sample = self.observable.problem.generate_state()
-		self.observable.problem.solveFwd(u_at_sample,[u_at_sample,m_sample,None])
+		if self.control_distribution is not None:
+			z_sample = self.observable.generate_vector(CONTROL)
+			self.control_distribution.sample(z_sample)
+			linearization_x = [u_at_mean,m_mean,None,z_sample]
+		else:
+			linearization_x = [u_at_mean,m_mean,None]
+
+		self.observable.problem.solveFwd(u_at_sample,linearization_x)
 
 		if self.parameters['verbose']:
 			print('||v_at_sample|| = ',u_at_sample.norm('l2'))
 		v_at_sample_pvd = dl.File(save_states_dir+'v_at_sample.pvd')
-		v_at_sample_pvd << vector2Function(u_at_sample,self.observable.problem.Vh[STATE])
+		v_at_sample_pvd << hp.vector2Function(u_at_sample,self.observable.problem.Vh[hp.STATE])
 
 
 
@@ -399,6 +542,7 @@ def input_output_error_test(self,V_MV,Cinv = None,rank_pairs = [None]):
 			- :code:`rank_pairs` - a python list of 2-tuples of ints specifying the input and output ranks 
 				to be used in the projection error test.
 		"""
+		assert self.control_distribution is None, 'Not worked out yet for control problems'
 		for (rank_in,rank_out) in rank_pairs:
 			assert rank_in <=V_MV.nvec()
 			assert rank_out <= self.U_MV.nvec()
@@ -408,31 +552,36 @@ def input_output_error_test(self,V_MV,Cinv = None,rank_pairs = [None]):
 		# truncate eigenvalues for numerical stability
 
 		# Instantiate parameter vector for requisite data structures
-		LocalParameters = MultiVector(self.m,self.parameters['sample_per_process'])
+		LocalParameters = hp.MultiVector(self.m,self.parameters['sample_per_process'])
 		LocalParameters.zero()
 
-		input_projection_vector = self.observable.generate_vector(PARAMETER)
+		input_projection_vector = self.observable.generate_vector(hp.PARAMETER)
 
 		# Instantiate an observable vector for requisite data structures
 		observable_vector = dl.Vector(self.mesh_constructor_comm)
 		self.observable.init_vector(observable_vector,dim = 0)
 
-		LocalObservables = MultiVector(observable_vector,self.parameters['sample_per_process'])
+		LocalObservables = hp.MultiVector(observable_vector,self.parameters['sample_per_process'])
 		LocalObservables.zero()
 		t0 = time.time()
 		for i in range(LocalObservables.nvec()):
-			parRandom.normal(1,self.noise)
+			hp.parRandom.normal(1,self.noise)
 			self.prior.sample(self.noise,LocalParameters[i])
-			x = [self.u,LocalParameters[i],None]
-			self.observable.setLinearizationPoint(x)
-			LocalObservables[i].axpy(1.,self.observable.eval(LocalParameters[i]))
+			if self.control_distribution is not None:
+				self.z.zero()
+				self.control_distribution.sample(self.z)
+				linearization_x = [self.u,LocalParameters[i],None,self.z]
+			else:
+				linearization_x = [self.u,LocalParameters[i],None]
+			self.observable.solveFwd(self.u,linearization_x)
+			LocalObservables[i].axpy(1.,self.observable.evalu(self.u))
 		if self.parameters['verbose'] and (self.mesh_constructor_comm.rank == 0):
 			print('Generating ',LocalObservables.nvec(),' local parameters and observables for POD error test took',time.time() -t0, 's')
 
 		# LocalProjectedObservables = MultiVector(observable_vector,self.parameters['sample_per_process'])
 		# LocalProjectedObservables.zero()
 
-		LocalErrors = MultiVector(observable_vector,self.parameters['sample_per_process'])
+		LocalErrors = hp.MultiVector(observable_vector,self.parameters['sample_per_process'])
 
 
 		output_projection_vector = dl.Vector(self.mesh_constructor_comm)
@@ -444,22 +593,22 @@ def input_output_error_test(self,V_MV,Cinv = None,rank_pairs = [None]):
 		global_std_rel_errors = []
 		for (rank_in,rank_out) in rank_pairs:
 			# Define input projector operator for rank_in
-			V_r = MultiVector(V_MV[0],rank_in)
+			V_r = hp.MultiVector(V_MV[0],rank_in)
 			for i in range(rank_in):
 				V_r[i].axpy(1.,V_MV[i])
 			input_init_vector_lambda = lambda x, dim: self.observable.init_vector(x,dim = 1)
 			if Cinv is not None:
 				InputProjectorOperator = PriorPreconditionedProjector(V_r,Cinv, input_init_vector_lambda)
 			else:
-				InputProjectorOperator = LowRankOperator(np.ones(rank_in),V_r, input_init_vector_lambda)
+				InputProjectorOperator = hp.LowRankOperator(np.ones(rank_in),V_r, input_init_vector_lambda)
 
 			# Define output projector operator for rank_out
-			U_MV = MultiVector(self.U_MV[0],rank_out)
+			U_MV = hp.MultiVector(self.U_MV[0],rank_out)
 			for i in range(rank_out):
 				U_MV[i].axpy(1.,self.U_MV[i])
 
 			init_vector_lambda = lambda x, dim: self.observable.init_vector(x,dim = 0)
-			PODOperator = LowRankOperator(np.ones(rank_out),U_MV,init_vector_lambda)
+			PODOperator = hp.LowRankOperator(np.ones(rank_out),U_MV,init_vector_lambda)
 
 
 			LocalErrors.zero()
@@ -472,9 +621,15 @@ def input_output_error_test(self,V_MV,Cinv = None,rank_pairs = [None]):
 				LocalErrors[i].axpy(1.,LocalObservables[i])
 				denominator = LocalErrors[i].norm('l2')
 				InputProjectorOperator.mult(LocalParameters[i],input_projection_vector)
-				x = [self.u,input_projection_vector,None]
-				self.observable.setLinearizationPoint(x)
-				reduced_q_vector.axpy(1.,self.observable.eval(input_projection_vector))
+				if self.control_distribution is not None:
+					self.z.zero()
+					self.control_distribution.sample(self.z)
+					linearization_x = [self.u,input_projection_vector,None, self.z]
+				else:
+					linearization_x = [self.u,input_projection_vector,None]
+
+				self.observable.solveFwd(self.u,linearization_x)
+				reduced_q_vector.axpy(1.,self.observable.evalu(self.u))
 
 				PODOperator.mult(reduced_q_vector,output_projection_vector)
 				LocalErrors[i].axpy(-1.,output_projection_vector)
diff --git a/hippyflow/modeling/__init__.py b/hippyflow/modeling/__init__.py
index 481d0cc..ba1e332 100644
--- a/hippyflow/modeling/__init__.py
+++ b/hippyflow/modeling/__init__.py
@@ -17,6 +17,8 @@
 
 from .cMinimization import ConstrainedNSolver
 
+from .controlJacobian import ObservableControlJacobian
+
 from .fullStateObservable import StateSpaceIdentityOperator
 
 from .hippylibModelWrapper import hippylibModelWrapper, hippylibModelWrapperSettings
@@ -33,7 +35,7 @@
 
 from .observable import LinearStateObservable, DomainRestrictedOperator, hippylibModelLinearStateObservable
 
-from .operatorWrappers import npToDolfinOperator
+from .operatorWrappers import npToDolfinOperator, MeanJTJfromDataOperator
 
 from .PODProjector import PODProjector, PODParameterList
 
diff --git a/hippyflow/modeling/activeSubspaceProjector.py b/hippyflow/modeling/activeSubspaceProjector.py
index a154965..7325b9a 100644
--- a/hippyflow/modeling/activeSubspaceProjector.py
+++ b/hippyflow/modeling/activeSubspaceProjector.py
@@ -11,24 +11,25 @@
 # terms of the GNU General Public License (as published by the Free
 # Software Foundation) version 2.0 dated June 1991.
 
-
 import dolfin as dl
 import numpy as np
 import os
-from hippylib import *
+import hippylib as hp
 from mpi4py import MPI 
 import time
 
-
 from ..collectives.collectiveOperator import CollectiveOperator, MatrixMultCollectiveOperator
 from ..collectives.collective import NullCollective
 from ..collectives.comm_utils import checkMeshConsistentPartitioning
 from .jacobian import *
+from .controlJacobian import ObservableControlJacobian
 from ..utilities.mv_utilities import mv_to_dense
 from ..utilities.plotting import *
 
 from .priorPreconditionedProjector import PriorPreconditionedProjector
 
+CONTROL = 3
+
 def ActiveSubspaceParameterList():
 	"""
 	This function implements a parameter list for the ActiveSubspaceProjector
@@ -38,6 +39,8 @@ def ActiveSubspaceParameterList():
 	parameters['jacobian_data_per_process'] = [512, 'Number of samples per process']
 	parameters['error_test_samples'] 		= [50, 'Number of samples for error test']
 	parameters['rank'] 				 		= [128, 'Rank of subspace']
+	parameters['jacobian_rank']				= [128, 'Rank of Jacobians generated']
+	parameters['control_jacobian_rank'] 	= [None, 'Rank of control Jacobians generated']
 	parameters['oversampling'] 		 		= [10, 'Oversampling parameter for randomized algorithms']
 	parameters['double_loop_samples']		= [20, 'Number of samples used in double loop MC approximation']
 	parameters['verbose']					= [True, 'Boolean for printing']
@@ -56,7 +59,7 @@ def ActiveSubspaceParameterList():
 	parameters['save_and_plot']				= [True,'Boolean for saving data and plots (only False for unit testing)']
 	parameters['store_Omega']				= [False,'Boolean for storing Gaussian random matrix (only True for unit testing)']
 	parameters['ms_given']					= [False,'Boolean for passing ms into serialized AS construction (only True for unit testing)']
-	return ParameterList(parameters)
+	return hp.ParameterList(parameters)
 
 
 class SummedListOperator:
@@ -91,7 +94,8 @@ class SeriallySampledJacobianOperator:
 	'''
 	Alterantive to SummedListOperator when memory is an issue for active subspace
 	'''
-	def __init__(self,observable,noise,prior,operation = 'JTJ',nsamples = None,ms = None,communicator=None,average=True):
+	def __init__(self,observable,noise,prior,control_distribution = None,operation = 'JTJ',\
+							nsamples = None,ms = None,zs = None,communicator=None,average=True):
 		'''
 		'''
 		assert operation in ['JTJ','JJT']
@@ -99,18 +103,32 @@ def __init__(self,observable,noise,prior,operation = 'JTJ',nsamples = None,ms =
 		self.observable = observable
 		self.noise = noise
 		self.prior = prior
+		self.control_distribution = control_distribution
 		self.operation = operation
 		self.nsamples = nsamples
 		self.average = average
 		self.ms = ms
+		if zs is not None:
+			self.zs = zs
+		else:
+			if type(self.ms) is list:
+				# This is a unit-test case
+				self.zs = len(self.ms)*[None]
+			else: 
+				self.zs = zs
+
 		
 		if communicator is None:
 			self.temp = None
 		else:
 			self.temp = dl.Vector(communicator)
 
-		self.u = self.observable.generate_vector(STATE)
-		self.m = self.observable.generate_vector(PARAMETER)
+		self.u = self.observable.generate_vector(hp.STATE)
+		self.m = self.observable.generate_vector(hp.PARAMETER)
+		if self.control_distribution is None:
+			self.z = None
+		else:
+			self.z = self.observable.generate_vector(CONTROL)
 
 	def init_vector(self,x):
 		"""
@@ -151,10 +169,17 @@ def matMvMult(self,x,y):
 						# Sample from the prior
 						self.m.zero()
 						self.noise.zero()
-						parRandom.normal(1,self.noise)
+						hp.parRandom.normal(1,self.noise)
 						self.prior.sample(self.noise,self.m)
+						if self.control_distribution is None:
+							linearization_x = [self.u,self.m,None]
+						else:
+							self.z.zero()
+							self.control_distribution.sample(self.z)
+							linearization_x = [self.u,self.m,None,self.z]
+
+
 						# Solve the PDE
-						linearization_x = [self.u,self.m,None]
 						print('Attempting to solve')
 						t0 = time.time()
 						self.observable.solveFwd(self.u,linearization_x)
@@ -165,7 +190,7 @@ def matMvMult(self,x,y):
 					except:
 						print('Issue with the solution, moving on')
 						pass
-				# Define action on matrix (as represented by MultiVector)
+				# Define action on matrix (as represented by hp.MultiVector)
 				for j in range(x.nvec()):
 					temp.zero()
 					t0 = time.time()
@@ -182,16 +207,18 @@ def matMvMult(self,x,y):
 		else:
 			# Each m should already not run into solver issues
 			nsamples = len(self.ms)
-			for m in self.ms:
+			for m,z in zip(self.ms,self.zs):
 				# Solve the PDE
 				print('Attempting to solve')
 				linearization_x = [self.u,m,None]
+				if z is not None:
+					linearization_x.append(z)
 				self.observable.solveFwd(self.u,linearization_x)
 				print('Solution succesful')
 				# set linearization point
 				self.observable.setLinearizationPoint(linearization_x)
 				solved = True
-				# Define action on matrix (as represented by MultiVector)
+				# Define action on matrix (as represented by hp.MultiVector)
 				for j in range(x.nvec()):
 					temp.zero()
 					operator_i.mult(x[j],temp)
@@ -211,7 +238,7 @@ class ActiveSubspaceProjector:
 	Output active subspace: :math:`J^*J = US^2U^*`
 	Input active subspace: :math:`JJ^* = VS^2V^*`
 	"""
-	def __init__(self,observable, prior, mesh_constructor_comm = None ,collective = NullCollective(),\
+	def __init__(self,observable, prior, control_distribution = None, mesh_constructor_comm = None ,collective = NullCollective(),\
 								  parameters = ActiveSubspaceParameterList()):
 		"""
 		Constructor
@@ -227,6 +254,7 @@ def __init__(self,observable, prior, mesh_constructor_comm = None ,collective =
 			print('Active Subspace object being created'.center(80))
 		self.observable = observable
 		self.prior = prior
+		self.control_distribution = control_distribution
 
 
 		if mesh_constructor_comm is not None:
@@ -265,10 +293,17 @@ def __init__(self,observable, prior, mesh_constructor_comm = None ,collective =
 			self.u = None
 			self.m = None
 			self.J = None
+			self.z = None
+			self.Jz = None
+			
+
 		else:		
 			self.us = None
 			self.ms = None
 			self.Js = None
+			self.zs = None
+			self.Jzs = None
+			
 
 		# Draw a new sample and set linearization point.
 		if self.parameters['initialize_samples']:
@@ -294,22 +329,31 @@ def _initialize_batched_samples(self):
 		This method initializes the samples from the prior used in sampling
 		"""
 		t0 = time.time()
-		self.us = [self.observable.generate_vector(STATE) for i in range(self.parameters['samples_per_process'])]
-		self.ms = [self.observable.generate_vector(PARAMETER) for i in range(self.parameters['samples_per_process'])]
-		for u,m,observable in zip(self.us,self.ms,self.observables):
+		self.us = [self.observable.generate_vector(hp.STATE) for i in range(self.parameters['samples_per_process'])]
+		self.ms = [self.observable.generate_vector(hp.PARAMETER) for i in range(self.parameters['samples_per_process'])]
+		if self.control_distribution is not None:
+			self.zs = [self.observable.generate_vector(CONTROL) for i in range(self.parameters['samples_per_process'])]
+		else:
+			self.zs = self.parameters['samples_per_process']*[None]
+		for u,m,z,observable in zip(self.us,self.ms,self.zs,self.observables):
 			solved = False
 			while not solved:
 				try:
 					self.noise.zero()
-					parRandom.normal(1,self.noise)
+					hp.parRandom.normal(1,self.noise)
 					# set linearization point
 					self.prior.sample(self.noise,m)
-					x = [u,m,None]
+					if self.control_distribution is not None:
+						self.control_distribution.sample(z)
+						x = [u,m,None,z]
+					else:
+						x = [u,m,None]
 					observable.solveFwd(u,x)
 					observable.setLinearizationPoint(x)
 					solved = True
 				except:
 					m.zero()
+					z.zero()
 					print('Issue with the solution, moving on')
 					pass
 		if self.parameters['verbose']:
@@ -319,6 +363,10 @@ def _initialize_batched_samples(self):
 				print(('Size of '+str(self.parameters['samples_per_process'])+' observable is '+str(asizeof.asizeof(self.observables)/1e6)+' MB').center(80))
 			except:
 				print('Install pympler and run again: pip install pympler'.center(80))
+		# This should be the same independent of whether their is a control variable
+		# dependence or not, since this Jacobian is just for the parameter.
+		# When adding mixed derivatives the Jacobian operator additionally needs
+		# to handle control adjoints.
 		self.Js = [ObservableJacobian(observable) for observable in self.observables]
 		total_init_time = time.time() - t0
 		for i in range(100):
@@ -326,18 +374,18 @@ def _initialize_batched_samples(self):
 			print('Initializing all batched samples took ',total_init_time, 's ')
 
 
-	def construct_input_subspace(self,prior_preconditioned = True):
+	def construct_input_subspace(self,prior_preconditioned = True,name_suffix = None):
 		if self.parameters['serialized_sampling']:
 			print('Construction via serialized AS construction'.center(80))
 			self._construct_serialized_jacobian_subspace(prior_preconditioned = prior_preconditioned,operation = 'JTJ')
 		else:
 			print('Construction via batched AS construction'.center(80))
-			self._construct_input_subspace_batched(prior_preconditioned = prior_preconditioned)
+			self._construct_input_subspace_batched(prior_preconditioned = prior_preconditioned,name_suffix = name_suffix)
 
 
 
 
-	def _construct_input_subspace_batched(self,prior_preconditioned = True):
+	def _construct_input_subspace_batched(self,prior_preconditioned = True,name_suffix = None):
 		"""
 		This method implements the input subspace constructor 
 			-:code:`prior_preconditioned` - a Boolean to decide whether to include the prior covariance in the decomposition
@@ -346,6 +394,7 @@ def _construct_input_subspace_batched(self,prior_preconditioned = True):
 		if self.Js is None:
 			self._initialize_batched_samples()
 
+
 		if self.parameters['verbose']:
 			print(80*'#')
 			print('Building derivative informed input subspace'.center(80))
@@ -360,10 +409,10 @@ def _construct_input_subspace_batched(self,prior_preconditioned = True):
 
 		x_GN = dl.Vector(self.mesh_constructor_comm)
 		GN_Hessians[0].init_vector(x_GN)
-		Omega = MultiVector(x_GN,self.parameters['rank'] + self.parameters['oversampling'])
+		Omega = hp.MultiVector(x_GN,self.parameters['rank'] + self.parameters['oversampling'])
 
 		if self.collective.rank() == 0:
-			parRandom.normal(1.,Omega)
+			hp.parRandom.normal(1.,Omega)
 			if self.parameters['store_Omega']:
 				self.Omega_GN = Omega
 		else:
@@ -374,13 +423,13 @@ def _construct_input_subspace_batched(self,prior_preconditioned = True):
 
 		if prior_preconditioned:
 			if hasattr(self.prior, "R"):
-				self.d_GN, self.V_GN = doublePassG(Average_GN_Hessian,\
+				self.d_GN, self.V_GN = hp.doublePassG(Average_GN_Hessian,\
 					self.prior.R, self.prior.Rsolver, Omega,self.parameters['rank'],s=1)
 			else:
-				self.d_GN, self.V_GN = doublePassG(Average_GN_Hessian,\
+				self.d_GN, self.V_GN = hp.doublePassG(Average_GN_Hessian,\
 					self.prior.Hlr, self.prior.Hlr, Omega,self.parameters['rank'],s=1)
 		else:
-			self.d_GN, self.V_GN = doublePass(Average_GN_Hessian,Omega,self.parameters['rank'],s=1)
+			self.d_GN, self.V_GN = hp.doublePass(Average_GN_Hessian,Omega,self.parameters['rank'],s=1)
 		total_init_time = time.time() - t0
 		for i in range(100):
 			print(80*'#')
@@ -391,15 +440,19 @@ def _construct_input_subspace_batched(self,prior_preconditioned = True):
 		if self.parameters['verbose'] and (self.mesh_constructor_comm.rank == 0):	
 			print(('Input subspace construction took '+str(self._input_subspace_construction_time)[:5]+' s').center(80))
 		if self.parameters['save_and_plot'] and MPI.COMM_WORLD.rank == 0:
-			np.save(self.parameters['output_directory']+'AS_input_projector',mv_to_dense(self.V_GN))
-			np.save(self.parameters['output_directory']+'AS_d_GN',self.d_GN)
-
-			out_name = self.parameters['output_directory']+'AS_input_eigenvalues_'+str(self.parameters['rank'])+'.pdf'
+			name = 'AS_'+str(int(self.parameters['samples_per_process']*self.collective.size()))
+			if name_suffix is not None:
+				assert type(name_suffix) is str
+				name += name_suffix
+			np.save(self.parameters['output_directory']+name+'_input_projector',mv_to_dense(self.V_GN))
+			np.save(self.parameters['output_directory']+name+'_d_GN',self.d_GN)
+
+			plot_out_name = self.parameters['output_directory']+name+'_input_eigenvalues_'+str(self.parameters['rank'])+'.pdf'
 			_ = spectrum_plot(self.d_GN,\
 				axis_label = ['i',r'$\lambda_i$',\
-				r'Eigenvalues of $\mathbb{E}_{\nu}[C{\nabla} q^T {\nabla} q]$'+self.parameters['plot_label_suffix']], out_name = out_name)
+				r'Eigenvalues of $\mathbb{E}_{\nu}[C{\nabla} q^T {\nabla} q]$'+self.parameters['plot_label_suffix']], out_name = plot_out_name)
 
-	def _construct_serialized_jacobian_subspace(self,prior_preconditioned = True, operation = 'JTJ'):
+	def _construct_serialized_jacobian_subspace(self,prior_preconditioned = True, operation = 'JTJ',name_suffix = None):
 		"""
 		This method implements the input subspace constructor 
 			-:code:`prior_preconditioned` - a Boolean to decide whether to include the prior covariance in the decomposition
@@ -410,31 +463,48 @@ def _construct_serialized_jacobian_subspace(self,prior_preconditioned = True, op
 		# ms_given is a unit testing case
 		if self.parameters['ms_given']:
 			assert self.ms is not None
+			if self.control_distribution is not None:
+				assert self.zs is not None
+				assert self.zs[0] is not None
 			Local_Average_Jacobian_Operator = SeriallySampledJacobianOperator(self.observable,self.noise,self.prior,operation = operation,\
-																				ms = self.ms)
+																				ms = self.ms,zs = self.zs)
 		else:
-			Local_Average_Jacobian_Operator = SeriallySampledJacobianOperator(self.observable,self.noise,self.prior,operation = operation,\
+			Local_Average_Jacobian_Operator = SeriallySampledJacobianOperator(self.observable,self.noise,self.prior,\
+														 control_distribution = self.control_distribution,operation = operation,\
 																				nsamples = self.parameters['samples_per_process'])
 		# This averaging assumes every process has an equal number of samples
 		# Otherwise it will bias towards a process with the fewest samples
 		Average_Jacobian_Operator = MatrixMultCollectiveOperator(Local_Average_Jacobian_Operator, self.collective, mpi_op = 'avg')
-		# Instantiate Gaussian random matrix
+		
 		if self.observable.problem.C is None:
+			# If this is the case, then the KKT blocks have not been built yet
+			# we overcome this by solving somewhere and setting the linearization pt.
 			m_mean = self.prior.mean
 			u_at_mean = self.observable.problem.generate_state()
-			self.observable.problem.solveFwd(u_at_mean,[u_at_mean,m_mean,None])
-			self.observable.setLinearizationPoint([u_at_mean,m_mean,None])
+			if self.control_distribution is not None:
+				if hasattr(self.control_distribution,'mean'):
+					z_mean = self.control_distribution.mean
+				else:
+					# Sample it somewhere
+					z_mean = self.observable.generate_vector(CONTROL)
+					self.control_distribution.sample(z_mean)
+				x_lin = [u_at_mean,m_mean,None,z_mean]
+			else:
+				x_lin = [u_at_mean,m_mean,None]
+			
+			self.observable.problem.solveFwd(u_at_mean,x_lin)
+			self.observable.setLinearizationPoint(x_lin)
 
+		# Instantiate Gaussian random matrix
 		x_Omega_construction = dl.Vector(self.mesh_constructor_comm)
-
 		Local_Average_Jacobian_Operator.init_vector(x_Omega_construction)
-		Omega = MultiVector(x_Omega_construction,self.parameters['rank'] + self.parameters['oversampling'])
+		Omega = hp.MultiVector(x_Omega_construction,self.parameters['rank'] + self.parameters['oversampling'])
 		Omega.zero()
 
 		if self.collective.rank() == 0:
 
 			if (operation == 'JTJ' and self.Omega_GN is None) or (operation == 'JJT' and self.Omega_NG is None):
-				parRandom.normal(1.,Omega)
+				hp.parRandom.normal(1.,Omega)
 			else:
 				if operation == 'JTJ':
 					Omega = self.Omega_GN
@@ -446,48 +516,52 @@ def _construct_serialized_jacobian_subspace(self,prior_preconditioned = True, op
 		if operation == 'JTJ':
 			if prior_preconditioned:
 				if hasattr(self.prior, "R"):
-					self.d_GN, self.V_GN = doublePassG(Average_Jacobian_Operator,\
+					self.d_GN, self.V_GN = hp.doublePassG(Average_Jacobian_Operator,\
 						self.prior.R, self.prior.Rsolver, Omega,self.parameters['rank'],s=1)
 				else:
-					self.d_GN, self.V_GN = doublePassG(Average_Jacobian_Operator,\
+					self.d_GN, self.V_GN = hp.doublePassG(Average_Jacobian_Operator,\
 						self.prior.Hlr, self.prior.Hlr, Omega,self.parameters['rank'],s=1)
 			else:
-				self.d_GN, self.V_GN = doublePass(Average_Jacobian_Operator,Omega,self.parameters['rank'],s=1)
+				self.d_GN, self.V_GN = hp.doublePass(Average_Jacobian_Operator,Omega,self.parameters['rank'],s=1)
 			self.prior_preconditioned = prior_preconditioned
 			self._input_subspace_construction_time = time.time() - t0
 			if self.parameters['verbose'] and (self.mesh_constructor_comm.rank == 0):	
 				print(('Input subspace construction took '+str(self._input_subspace_construction_time)[:5]+' s').center(80))
 		elif operation == 'JJT':
-			self.d_NG, self.U_NG = doublePass(Average_Jacobian_Operator,Omega,self.parameters['rank'],s=1)
+			self.d_NG, self.U_NG = hp.doublePass(Average_Jacobian_Operator,Omega,self.parameters['rank'],s=1)
 			self._output_subspace_construction_time = time.time() - t0
 			if self.parameters['verbose'] and (self.mesh_constructor_comm.rank ==0):	
 				print(('Output subspace construction took '+str(self._output_subspace_construction_time)[:5]+' s').center(80))
 		
 		if self.parameters['save_and_plot'] and MPI.COMM_WORLD.rank == 0:
+			name = 'AS_'+str(int(self.parameters['samples_per_process']*self.collective.size()))
+			if name_suffix is not None:
+				assert type(name_suffix) is str
+				name += name_suffix
 			if operation == 'JTJ':
-				np.save(self.parameters['output_directory']+'AS_input_projector',mv_to_dense(self.V_GN))
-				np.save(self.parameters['output_directory']+'AS_d_GN',self.d_GN)
-				out_name = self.parameters['output_directory']+'AS_input_eigenvalues_'+str(self.parameters['rank'])+'.pdf'
+				np.save(self.parameters['output_directory']+name+'_input_projector',mv_to_dense(self.V_GN))
+				np.save(self.parameters['output_directory']+name+'_d_GN',self.d_GN)
+				plot_out_name = self.parameters['output_directory']+name+'_input_eigenvalues_'+str(self.parameters['rank'])+'.pdf'
 				_ = spectrum_plot(self.d_GN,\
 					axis_label = ['i',r'$\lambda_i$',\
-					r'Eigenvalues of $\mathbb{E}_{\nu}[C{\nabla} q^T {\nabla} q]$'+self.parameters['plot_label_suffix']], out_name = out_name)
+					r'Eigenvalues of $\mathbb{E}_{\nu}[C{\nabla} q^T {\nabla} q]$'+self.parameters['plot_label_suffix']], out_name = plot_out_name)
 			if operation == 'JJT':
-				np.save(self.parameters['output_directory']+'AS_output_projector',mv_to_dense(self.U_NG))
-				np.save(self.parameters['output_directory']+'AS_d_NG',self.d_NG)
+				np.save(self.parameters['output_directory']+name+'_output_projector',mv_to_dense(self.U_NG))
+				np.save(self.parameters['output_directory']+name+'_d_NG',self.d_NG)
 
-				out_name = self.parameters['output_directory']+'AS_output_eigenvalues_'+str(self.parameters['rank'])+'.pdf'
+				plot_out_name = self.parameters['output_directory']+name+'_output_eigenvalues_'+str(self.parameters['rank'])+'.pdf'
 				_ = spectrum_plot(self.d_NG,\
 					axis_label = ['i',r'$\lambda_i$',\
-					r'Eigenvalues of $\mathbb{E}_{\nu}[{\nabla} q {\nabla} q^T]$'+self.parameters['plot_label_suffix']], out_name = out_name)
+					r'Eigenvalues of $\mathbb{E}_{\nu}[{\nabla} q {\nabla} q^T]$'+self.parameters['plot_label_suffix']], out_name = plot_out_name)
 
-	def construct_output_subspace(self):
+	def construct_output_subspace(self,name_suffix = None):
 		if self.parameters['serialized_sampling']:
 			print('Construction via serialized construction'.center(80))
-			self._construct_serialized_jacobian_subspace(operation = 'JJT')
+			self._construct_serialized_jacobian_subspace(operation = 'JJT',name_suffix = name_suffix)
 		else:
-			self._construct_output_subspace_batched()
+			self._construct_output_subspace_batched(name_suffix = name_suffix)
 
-	def _construct_output_subspace_batched(self):
+	def _construct_output_subspace_batched(self,name_suffix = None):
 		"""
 		This method implements the output subspace constructor 
 		"""
@@ -506,28 +580,32 @@ def _construct_output_subspace_batched(self):
 
 		x_NG = dl.Vector(self.mesh_constructor_comm)
 		NG_Hessians[0].init_vector(x_NG)
-		Omega = MultiVector(x_NG,self.parameters['rank'] + self.parameters['oversampling'])
+		Omega = hp.MultiVector(x_NG,self.parameters['rank'] + self.parameters['oversampling'])
 
 		if self.collective.rank() == 0:
-			parRandom.normal(1.,Omega)
+			hp.parRandom.normal(1.,Omega)
 			if self.parameters['store_Omega']:
 				self.Omega_NG = Omega
 		else:
 			Omega.zero()
 
 		self.collective.bcast(Omega,root = 0)
-		self.d_NG, self.U_NG = doublePass(Average_NG_Hessian,Omega,self.parameters['rank'],s=1)
+		self.d_NG, self.U_NG = hp.doublePass(Average_NG_Hessian,Omega,self.parameters['rank'],s=1)
 		self._output_subspace_construction_time = time.time() - t0
 		if self.parameters['verbose'] and (self.mesh_constructor_comm.rank ==0):	
 			print(('Output subspace construction took '+str(self._output_subspace_construction_time)[:5]+' s').center(80))
 		if self.parameters['save_and_plot'] and MPI.COMM_WORLD.rank == 0:
-			np.save(self.parameters['output_directory']+'AS_output_projector',mv_to_dense(self.U_NG))
-			np.save(self.parameters['output_directory']+'AS_d_NG',self.d_NG)
-
-			out_name = self.parameters['output_directory']+'AS_output_eigenvalues_'+str(self.parameters['rank'])+'.pdf'
+			name = 'AS_'+str(int(self.parameters['samples_per_process']*self.collective.size()))
+			if name_suffix is not None:
+				assert type(name_suffix) is str
+				name += name_suffix
+			np.save(self.parameters['output_directory']+name+'_output_projector',mv_to_dense(self.U_NG))
+			np.save(self.parameters['output_directory']+name+'_d_NG',self.d_NG)
+
+			plot_out_name = self.parameters['output_directory']+name+'_output_eigenvalues_'+str(self.parameters['rank'])+'.pdf'
 			_ = spectrum_plot(self.d_NG,\
 				axis_label = ['i',r'$\lambda_i$',\
-				r'Eigenvalues of $\mathbb{E}_{\nu}[{\nabla} q {\nabla} q^T]$'+self.parameters['plot_label_suffix']], out_name = out_name)
+				r'Eigenvalues of $\mathbb{E}_{\nu}[{\nabla} q {\nabla} q^T]$'+self.parameters['plot_label_suffix']], out_name = plot_out_name)
 
 	def construct_low_rank_Jacobians(self,check_for_data = True,compress_files = True):
 		if self.parameters['serialized_sampling']:
@@ -535,7 +613,16 @@ def construct_low_rank_Jacobians(self,check_for_data = True,compress_files = Tru
 		else:
 			self._construct_low_rank_Jacobians_batched(check_for_data = check_for_data)
 
-	def _construct_low_rank_Jacobians_serial(self,check_for_data = True,compress_files = True):
+	def construct_low_rank_control_Jacobians(self,check_for_data = True,compress_files = True):
+		if self.parameters['serialized_sampling']:
+			self._construct_low_rank_Jacobians_serial(check_for_data = check_for_data,compress_files = compress_files,\
+					parameter_jacobian = False,control_jacobian = True)
+		else:
+			raise
+			self._construct_low_rank_Jacobians_batched(check_for_data = check_for_data, control_jacobian = True)
+
+	def _construct_low_rank_Jacobians_serial(self,check_for_data = True,compress_files = True,\
+						parameter_jacobian = True, control_jacobian = False):
 		"""
 		This method generates low rank Jacobians for training (and also saves input output data in tandem)
 			- :code:`check_for_data` - a boolean to decide whether to check to see if the training
@@ -544,34 +631,83 @@ def _construct_low_rank_Jacobians_serial(self,check_for_data = True,compress_fil
 		is by default a jacobian_data/ directory
 		This allows for separate sampling for l2 loss and h1 seminorm loss
 		"""
-		my_rank = int(self.collective.rank())
-		jacobian_rank_specific_directory = self.parameters['output_directory']+'jacobian_data/proc_'+str(my_rank)+'/'
-		os.makedirs(jacobian_rank_specific_directory,exist_ok = True)
-		rank_specific_directory = self.parameters['output_directory']+'data_on_rank_'+str(my_rank)+'/'
-		os.makedirs(rank_specific_directory,exist_ok = True)
+		if control_jacobian:
+			assert (self.control_distribution is not None)
+
+		parameter_dimension = None
+
+		proc_id = int(self.collective.rank())
+		jacobian_process_specific_directory = self.parameters['output_directory']+'jacobian_data/proc_'+str(proc_id)+'/'
+		os.makedirs(jacobian_process_specific_directory,exist_ok = True)
+		process_specific_directory = self.parameters['output_directory']+'data_on_proc_'+str(proc_id)+'/'
+		os.makedirs(process_specific_directory,exist_ok = True)
 		if self.u is None:
-			self.u = self.observable.generate_vector(STATE)
+			self.u = self.observable.generate_vector(hp.STATE)
 		if self.m is None:
-			self.m = self.observable.generate_vector(PARAMETER)
+			self.m = self.observable.generate_vector(hp.PARAMETER)
+		if self.control_distribution is not None and self.z is None:
+			self.z = self.observable.generate_vector(CONTROL)
+
+		if self.control_distribution is not None:
+			assert self.observable.problem.Vh[hp.STATE].mesh().mpi_comm().size == 1, print('Only worked out for serial codes')
+			control_dimension = self.z.get_local().shape[0]
+
+		if parameter_jacobian:
+			self.J = ObservableJacobian(self.observable)
+			output_dimension,parameter_dimension = self.J.shape
+			parameter_rank = min(self.parameters['jacobian_rank'],output_dimension,parameter_dimension)
+
+		if parameter_dimension is None:
+			parameter_dimension = self.m.size()
+
+		if control_jacobian:
+			self.Jz = ObservableControlJacobian(self.observable)
+			output_dimension,control_dimension = self.Jz.shape
+			control_rank = min(self.parameters['control_jacobian_rank'],output_dimension,control_dimension)
 
-		self.J = ObservableJacobian(self.observable)
 
 		last_datum_generated = 0
-		output_dimension,input_dimension = self.J.shape
-		rank = min(self.parameters['rank'],output_dimension,input_dimension)
-		# Initialize randomized Omega
+		
 		if self.observable.problem.C is None:
+			# If this is the case, then the KKT blocks have not been built yet
+			# we overcome this by solving somewhere and setting the linearization pt.
 			m_mean = self.prior.mean
 			u_at_mean = self.observable.problem.generate_state()
-			self.observable.problem.solveFwd(u_at_mean,[u_at_mean,m_mean,None])
-			self.observable.setLinearizationPoint([u_at_mean,m_mean,None])
-		input_vector = dl.Vector(self.mesh_constructor_comm)
-		self.J.init_vector(input_vector,1)
-		Omega = MultiVector(input_vector,rank + self.parameters['oversampling'])
-		# Omega does not need to be communicated across processes in this case
-		# like with the global reduction collectives
-		# Omega can be the same for all samples
-		parRandom.normal(1.,Omega)
+			if self.control_distribution is not None:
+				if hasattr(self.control_distribution,'mean'):
+					z_mean = self.control_distribution.mean
+				else:
+					# Sample it somewhere
+					z_mean = self.observable.generate_vector(CONTROL)
+					self.control_distribution.sample(z_mean)
+				x_lin = [u_at_mean,m_mean,None,z_mean]
+			else:
+				x_lin = [u_at_mean,m_mean,None]
+			
+			self.observable.problem.solveFwd(u_at_mean,x_lin)
+			self.observable.setLinearizationPoint(x_lin)
+			if self.control_distribution is not None:
+				assert self.observable.problem.Cz is not None
+		
+		# Initialize randomized Omega
+		if parameter_jacobian:
+			parameter_vector = dl.Vector(self.mesh_constructor_comm)
+			self.J.init_vector(parameter_vector,1)
+			nvec_Omega_m = min(parameter_rank + self.parameters['oversampling'],output_dimension,parameter_dimension)
+			Omega_m = hp.MultiVector(parameter_vector,nvec_Omega_m)
+			# Omega does not need to be communicated across processes in this case
+			# like with the global reduction collectives
+			# Omega can be the same for all samples
+			hp.parRandom.normal(1.,Omega_m)
+		if control_jacobian:
+			control_vector = dl.Vector(self.mesh_constructor_comm)
+			self.Jz.init_vector(control_vector,1)
+			nvec_Omega_z = min(control_rank + self.parameters['oversampling'],output_dimension,control_dimension)
+			Omega_z = hp.MultiVector(control_vector,nvec_Omega_z)
+			# Omega does not need to be communicated across processes in this case
+			# like with the global reduction collectives
+			# Omega can be the same for all samples
+			hp.parRandom.normal(1.,Omega_z)
 
 		if check_for_data:
 			# Find largest mq pair generated
@@ -582,62 +718,119 @@ def _construct_low_rank_Jacobians_serial(self,check_for_data = True,compress_fil
 		t0 = time.time()
 		for i in range(last_datum_generated,self.parameters['jacobian_data_per_process']):
 			print(('Generating data number '+str(i)).center(80))
-			parRandom.normal(1,self.noise)
+			hp.parRandom.normal(1,self.noise)
 			self.m.zero() # This is probably redundant
 			self.prior.sample(self.noise,self.m)
-			x = [self.u,self.m,None]
+			if self.control_distribution is not None:
+				self.control_distribution.sample(self.z)
+				x = [self.u,self.m,None,self.z]
+			else:
+				x = [self.u,self.m,None]
 			self.observable.solveFwd(self.u,x)
 			self.observable.setLinearizationPoint(x)
 			this_m = self.m.get_local()
 			this_q = self.observable.evalu(self.u).get_local()
-
-			np.save(rank_specific_directory+'m_sample_'+str(i)+'.npy',this_m)
-			np.save(rank_specific_directory+'q_sample_'+str(i)+'.npy',this_q)
-
-			Omega.zero() # probably unecessary
-			parRandom.normal(1.,Omega)
-			U, sigma, V = accuracyEnhancedSVD(self.J,Omega,rank,s=1)
-			Unp = mv_to_dense(U)
-			Vnp = mv_to_dense(V)
-
-			np.save(jacobian_rank_specific_directory+'U_sample_'+str(i)+'.npy',Unp)
-			np.save(jacobian_rank_specific_directory+'sigma_sample_'+str(i)+'.npy',sigma)
-			np.save(jacobian_rank_specific_directory+'V_sample_'+str(i)+'.npy',Vnp)
+			np.save(process_specific_directory+'m_sample_'+str(i)+'.npy',this_m)
+			np.save(process_specific_directory+'q_sample_'+str(i)+'.npy',this_q)
+			if self.control_distribution is not None:
+				this_z = self.z.get_local()
+				np.save(process_specific_directory+'z_sample_'+str(i)+'.npy',this_z)
+
+			if parameter_jacobian:
+				Omega_m.zero() # probably unecessary
+				hp.parRandom.normal(1.,Omega_m)
+				U, sigma, V = hp.accuracyEnhancedSVD(self.J,Omega_m,parameter_rank,s=1)
+				Unp = mv_to_dense(U)
+				Vnp = mv_to_dense(V)
+
+				np.save(jacobian_process_specific_directory+'U_sample_'+str(i)+'.npy',Unp)
+				np.save(jacobian_process_specific_directory+'sigma_sample_'+str(i)+'.npy',sigma)
+				np.save(jacobian_process_specific_directory+'V_sample_'+str(i)+'.npy',Vnp)
+
+			if control_jacobian:
+				Omega_z.zero() # probably unecessary
+				hp.parRandom.normal(1.,Omega_z)
+
+				Uz, sigmaz, Vz = hp.accuracyEnhancedSVD(self.Jz,Omega_z,control_rank,s=1)
+				Uznp = mv_to_dense(Uz)
+				Vznp = mv_to_dense(Vz)
+
+				np.save(jacobian_process_specific_directory+'Uz_sample_'+str(i)+'.npy',Uznp)
+				np.save(jacobian_process_specific_directory+'sigmaz_sample_'+str(i)+'.npy',sigmaz)
+				np.save(jacobian_process_specific_directory+'Vz_sample_'+str(i)+'.npy',Vznp)
 
 			if self.parameters['verbose']:
-				print('One (m,q,J) generated every ',(time.time() -t0)/(i - last_datum_generated+1),' s, on average.')
+				if self.control_distribution is None:
+					print('One (m,q(m),J(m)) generated every ',(time.time() -t0)/(i - last_datum_generated+1),' s, on average.')
+				else:
+					pref = 'One (m,z,q(m,z),'
+					if parameter_jacobian:
+						pref += 'J(m,z)'
+					if control_jacobian:
+						pref += 'J_z(m,z)'
+					print(pref+') generated every ',(time.time() -t0)/(i - last_datum_generated+1),' s, on average.')
 
 		if compress_files:
 			print('Compressing mq data'.center(80))
 			t_start_mq = time.time()
-			local_ms = np.zeros((self.parameters['jacobian_data_per_process'],input_dimension))
+			local_ms = np.zeros((self.parameters['jacobian_data_per_process'],parameter_dimension))
 			local_qs = np.zeros((self.parameters['jacobian_data_per_process'],output_dimension))
+			if self.control_distribution is not None:
+				local_zs = np.zeros((self.parameters['jacobian_data_per_process'],control_dimension))
 			for i in range(0,self.parameters['jacobian_data_per_process']):
-				local_ms[i] = np.load(rank_specific_directory+'m_sample_'+str(i)+'.npy')
-				local_qs[i] = np.load(rank_specific_directory+'q_sample_'+str(i)+'.npy')
-			np.savez_compressed(self.parameters['output_directory']+'mq_on_rank'+str(my_rank)+'.npz',m_data = local_ms,q_data = local_qs)
-			print(('mq compression took '+str(time.time()-t_start_mq)+' s '))
-			t_start_J = time.time()
-			print('Compressing Jacobian data'.center(80))
-			local_Us = np.zeros((self.parameters['jacobian_data_per_process'],output_dimension,rank))
-			local_sigmas = np.zeros((self.parameters['jacobian_data_per_process'],rank))
-			local_Vs = np.zeros((self.parameters['jacobian_data_per_process'],input_dimension,rank))
-			for i in range(0,self.parameters['jacobian_data_per_process']):
-				local_Us[i] = np.load(jacobian_rank_specific_directory+'U_sample_'+str(i)+'.npy')
-				local_sigmas[i] = np.load(jacobian_rank_specific_directory+'sigma_sample_'+str(i)+'.npy')
-				local_Vs[i] = np.load(jacobian_rank_specific_directory+'V_sample_'+str(i)+'.npy')
-			np.savez_compressed(self.parameters['output_directory']+'J_on_rank'+str(my_rank)+'.npz',\
-						U_data = local_Us,sigma_data = local_sigmas,V_data = local_Vs)
-			print(('Jacobian compression took '+str(time.time()-t_start_J)+' s '))
-
-			if True:
-				out_name = self.parameters['output_directory']+'jacobian_singular_values_'+str(rank)+'.pdf'
-				plot_singular_values_with_std(np.mean(local_sigmas,axis=0),np.std(local_sigmas,axis=0),outname= out_name)
+				local_ms[i] = np.load(process_specific_directory+'m_sample_'+str(i)+'.npy')
+				local_qs[i] = np.load(process_specific_directory+'q_sample_'+str(i)+'.npy')
+				if self.control_distribution is not None:
+					local_zs[i] = np.load(process_specific_directory+'z_sample_'+str(i)+'.npy')
+			if self.control_distribution is None:
+				np.savez_compressed(self.parameters['output_directory']+'mq_on_proc'+str(proc_id)+'.npz',m_data = local_ms,q_data = local_qs)
+				print(('mq compression took '+str(time.time()-t_start_mq)+' s '))
+			else:
+				np.savez_compressed(self.parameters['output_directory']+'mzq_on_proc'+str(proc_id)+'.npz',m_data = local_ms,z_data = local_zs,\
+																															q_data = local_qs)
+				print(('mzq compression took '+str(time.time()-t_start_mq)+' s '))
+
+			if parameter_jacobian:
+				t_start_J = time.time()
+				print('Compressing Jacobian data'.center(80))
+				local_Us = np.zeros((self.parameters['jacobian_data_per_process'],output_dimension,parameter_rank))
+				local_sigmas = np.zeros((self.parameters['jacobian_data_per_process'],parameter_rank))
+				local_Vs = np.zeros((self.parameters['jacobian_data_per_process'],parameter_dimension,parameter_rank))
+				for i in range(0,self.parameters['jacobian_data_per_process']):
+					local_Us[i] = np.load(jacobian_process_specific_directory+'U_sample_'+str(i)+'.npy')
+					local_sigmas[i] = np.load(jacobian_process_specific_directory+'sigma_sample_'+str(i)+'.npy')
+					local_Vs[i] = np.load(jacobian_process_specific_directory+'V_sample_'+str(i)+'.npy')
+				np.savez_compressed(self.parameters['output_directory']+'J_on_proc'+str(proc_id)+'.npz',\
+							U_data = local_Us,sigma_data = local_sigmas,V_data = local_Vs)
+				print(('Jacobian compression took '+str(time.time()-t_start_J)+' s '))
+
+				if True:
+					out_name = self.parameters['output_directory']+'jacobian_singular_values_'+str(parameter_rank)+'.pdf'
+					plot_singular_values_with_std(np.mean(local_sigmas,axis=0),np.std(local_sigmas,axis=0),outname= out_name)
+
+			if control_jacobian:
+				t_start_Jz = time.time()
+				print('Compressing Jacobian data'.center(80))
+				local_Uzs = np.zeros((self.parameters['jacobian_data_per_process'],output_dimension,control_rank))
+				local_sigmazs = np.zeros((self.parameters['jacobian_data_per_process'],control_rank))
+				local_Vzs = np.zeros((self.parameters['jacobian_data_per_process'],control_dimension,control_rank))
+				for i in range(0,self.parameters['jacobian_data_per_process']):
+					local_Uzs[i] = np.load(jacobian_process_specific_directory+'Uz_sample_'+str(i)+'.npy')
+					local_sigmazs[i] = np.load(jacobian_process_specific_directory+'sigmaz_sample_'+str(i)+'.npy')
+					local_Vzs[i] = np.load(jacobian_process_specific_directory+'Vz_sample_'+str(i)+'.npy')
+				np.savez_compressed(self.parameters['output_directory']+'Jz_on_proc'+str(proc_id)+'.npz',\
+							Uz_data = local_Uzs,sigmaz_data = local_sigmazs,Vz_data = local_Vzs)
+				print(('Control Jacobian compression took '+str(time.time()-t_start_Jz)+' s '))
+
+				if True:
+					out_name = self.parameters['output_directory']+'control_jacobian_singular_values_'+str(control_rank)+'.pdf'
+					plot_singular_values_with_std(np.mean(local_sigmazs,axis=0),np.std(local_sigmazs,axis=0),outname= out_name)
 
 		self._jacobian_data_generation_time = time.time() - t0
 
 
-	def _construct_low_rank_Jacobians_batched(self,check_for_data = True):
+	def _construct_low_rank_Jacobians_batched(self,check_for_data = True,\
+						parameter_jacobian = True, control_jacobian = False):
 		"""
 		This method generates low rank Jacobians for training (and also saves input output data in tandem)
 			- :code:`check_for_data` - a boolean to decide whether to check to see if the training
@@ -646,37 +839,42 @@ def _construct_low_rank_Jacobians_batched(self,check_for_data = True):
 		is by default a jacobian_data/ directory
 		This allows for separate sampling for l2 loss and h1 seminorm loss
 		"""
+		assert not control_jacobian
 
 		if self.Js is None:
 			self.initialize_samples()
 
-		my_rank = int(self.collective.rank())
+		proc_id = int(self.collective.rank())
 		output_directory = self.parameters['output_directory']+'jacobian_data/'
 		os.makedirs(output_directory,exist_ok = True)
 
 		last_datum_generated = 0
-		output_dimension,input_dimension = self.Js[0].shape
+		output_dimension,parameter_dimension = self.Js[0].shape
+		if self.control_distribution is not None:
+			assert self.observable.problem.Vh[hp.STATE].mesh().mpi_comm().size == 1, print('Only worked out for serial codes')
+			control_dimension = self.z.get_local().shape[0]
 
-		rank = min(self.parameters['rank'],output_dimension,input_dimension)
+		rank = min(self.parameters['rank'],output_dimension,parameter_dimension)
 		local_Us = np.zeros((0,output_dimension, rank))	
 		local_sigmas = np.zeros((0,rank))
-		local_Vs = np.zeros((0,input_dimension, rank))
-		local_ms = np.zeros((0,input_dimension))
+		local_Vs = np.zeros((0,parameter_dimension, rank))
+		local_ms = np.zeros((0,parameter_dimension))
+		local_zs = np.zeros((0,control_dimension))
 		local_qs = np.zeros((0,output_dimension))
 		# Initialize arrays
 		if check_for_data:
 			# Save all five or restart sampling and saving
-			if os.path.isfile(output_directory+'Us_on_rank_'+str(my_rank)+'.npy') and \
-				os.path.isfile(output_directory+'sigmas_on_rank_'+str(my_rank)+'.npy') and \
-				os.path.isfile(output_directory+'Vs_on_rank_'+str(my_rank)+'.npy') and \
-				os.path.isfile(output_directory+'ms_on_rank_'+str(my_rank)+'.npy') and \
-				os.path.isfile(output_directory+'qs_on_rank_'+str(my_rank)+'.npy'):
-
-				local_Us = np.load(output_directory+'Us_on_rank_'+str(my_rank)+'.npy')
-				local_sigmas = np.load(output_directory+'sigmas_on_rank_'+str(my_rank)+'.npy')
-				local_Vs = np.load(output_directory+'Vs_on_rank_'+str(my_rank)+'.npy')
-				local_ms = np.load(output_directory+'ms_on_rank_'+str(my_rank)+'.npy')
-				local_qs = np.load(output_directory+'qs_on_rank_'+str(my_rank)+'.npy')
+			if os.path.isfile(output_directory+'Us_on_proc_'+str(proc_id)+'.npy') and \
+				os.path.isfile(output_directory+'sigmas_on_proc_'+str(proc_id)+'.npy') and \
+				os.path.isfile(output_directory+'Vs_on_proc_'+str(proc_id)+'.npy') and \
+				os.path.isfile(output_directory+'ms_on_proc_'+str(proc_id)+'.npy') and \
+				os.path.isfile(output_directory+'qs_on_proc_'+str(proc_id)+'.npy'):
+
+				local_Us = np.load(output_directory+'Us_on_proc_'+str(proc_id)+'.npy')
+				local_sigmas = np.load(output_directory+'sigmas_on_proc_'+str(proc_id)+'.npy')
+				local_Vs = np.load(output_directory+'Vs_on_proc_'+str(proc_id)+'.npy')
+				local_ms = np.load(output_directory+'ms_on_proc_'+str(proc_id)+'.npy')
+				local_qs = np.load(output_directory+'qs_on_proc_'+str(proc_id)+'.npy')
 
 				last_datum_generated = min(local_Us.shape[0],local_sigmas.shape[0],local_Vs.shape[0],\
 											local_ms.shape[0],local_qs.shape[0])
@@ -691,10 +889,21 @@ def _construct_low_rank_Jacobians_batched(self,check_for_data = True):
 					local_ms = local_ms[:last_datum_generated,:]
 				if local_qs.shape[0] > last_datum_generated:
 					local_qs = local_qs[:last_datum_generated,:]
-
+			if (self.control_distribution is not None) and os.path.isfile(output_directory+'zs_on_proc_'+str(proc_id)+'.npy'):
+				local_zs = np.load(output_directory+'zs_on_proc_'+str(proc_id)+'.npy')
+				last_z_generated = local_zs.shape[0]
+				if last_z_generated < last_datum_generated:
+					last_datum_generated = last_z_generated
+					# Update slices
+					local_Us = local_Us[:last_datum_generated,:,:]
+					local_sigmas = local_sigmas[:last_datum_generated,:]
+					local_Vs = local_Vs[:last_datum_generated,:,:]
+					local_ms = local_ms[:last_datum_generated,:]
+					local_qs = local_qs[:last_datum_generated,:]
 
 		# Generate the input output pairs that correspond to the 
 		assert len(self.ms) == self.parameters['samples_per_process']
+		assert len(self.zs) == self.parameters['samples_per_process']
 		assert len(self.us) == self.parameters['samples_per_process']
 		# If the us are updated in place then create a method for the observable that just applies B to u
 		t0_mq = time.time()
@@ -705,21 +914,26 @@ def _construct_low_rank_Jacobians_batched(self,check_for_data = True):
 				print('Saving input output data pair '+str(i))
 
 			local_ms = np.concatenate((local_ms,np.expand_dims(self.ms[i].get_local(),0)))
+			
 			qi = self.observables[i].evalu(self.us[i]).get_local()
 			local_qs = np.concatenate((local_qs,np.expand_dims(qi,0)))
-			np.save(output_directory+'ms_on_rank_'+str(my_rank)+'.npy',np.array(local_ms))
-			np.save(output_directory+'qs_on_rank_'+str(my_rank)+'.npy',np.array(local_qs))
+			np.save(output_directory+'ms_on_proc_'+str(proc_id)+'.npy',np.array(local_ms))
+			np.save(output_directory+'qs_on_proc_'+str(proc_id)+'.npy',np.array(local_qs))
+			if self.control_distribution is not None:
+				local_zs = np.concatentate((local_zs,np.expand_dims(self.zs[i].get_local(),0)))
+				np.save(output_directory+'zs_on_proc_'+str(proc_id)+'.npy',np.array(local_zs))
+
 			if self.parameters['verbose']:
 				print('On datum saved every ',(time.time() -t0_mq)/(i - last_datum_generated+1),' s, on average.')
 
 
 		# Initialize randomized Omega
-		input_vector = dl.Vector(self.mesh_constructor_comm)
-		self.Js[0].init_vector(input_vector,1)
-		Omega = MultiVector(input_vector,rank + self.parameters['oversampling'])
+		parameter_vector = dl.Vector(self.mesh_constructor_comm)
+		self.Js[0].init_vector(parameter_vector,1)
+		Omega = hp.MultiVector(parameter_vector,rank + self.parameters['oversampling'])
 		# Omega does not need to be communicated across processes in this case
 		# like with the global reduction collectives
-		parRandom.normal(1.,Omega)
+		hp.parRandom.normal(1.,Omega)
 
 		t0 = time.time()
 		# I think this is all hard coded for a single serial mesh, check if 
@@ -735,18 +949,16 @@ def _construct_low_rank_Jacobians_batched(self,check_for_data = True):
 			# Reusing Omega for each randomized pass, this shouldn't be an issue,
 			# but one could resample at each iteration
 
-			U, sigma, V = accuracyEnhancedSVD(self.Js[i],Omega,rank,s=1)
+			U, sigma, V = hp.accuracyEnhancedSVD(self.Js[i],Omega,rank,s=1)
 
 			local_Us = np.concatenate((local_Us,np.expand_dims(mv_to_dense(U),0)))
 			local_sigmas = np.concatenate((local_sigmas,np.expand_dims(sigma,0)))
 			local_Vs = np.concatenate((local_Vs,np.expand_dims(mv_to_dense(V),0)))
 
-
-
 			if self.mesh_constructor_comm.rank == 0:
-				np.save(output_directory+'Us_on_rank_'+str(my_rank)+'.npy',np.array(local_Us))
-				np.save(output_directory+'sigmas_on_rank_'+str(my_rank)+'.npy',np.array(local_sigmas))
-				np.save(output_directory+'Vs_on_rank_'+str(my_rank)+'.npy',np.array(local_Vs))
+				np.save(output_directory+'Us_on_proc_'+str(proc_id)+'.npy',np.array(local_Us))
+				np.save(output_directory+'sigmas_on_proc_'+str(proc_id)+'.npy',np.array(local_sigmas))
+				np.save(output_directory+'Vs_on_proc_'+str(proc_id)+'.npy',np.array(local_Vs))
 			if self.parameters['verbose']:
 				print('On Jacobian datum generated every ',(time.time() -t0)/(i - last_datum_generated+1),' s, on average.')
 
@@ -792,16 +1004,16 @@ def test_errors(self,test_input = True, test_output = False, ranks = [None],cut_
 			global_std_rel_errors_input = np.zeros_like(ranks,dtype = np.float64)
 
 			# Naive test on output space
-			LocalParameters = MultiVector(self.ms[0],self.parameters['error_test_samples'])
+			LocalParameters = hp.MultiVector(self.ms[0],self.parameters['error_test_samples'])
 			LocalParameters.zero()
 			# Generate samples
 			for i in range(self.parameters['error_test_samples']):
 				t0 = time.time()
-				parRandom.normal(1,self.noise)
+				hp.parRandom.normal(1,self.noise)
 				self.prior.sample(self.noise,LocalParameters[i])
 
-			LocalErrors = MultiVector(self.ms[0],self.parameters['error_test_samples'])
-			projection_vector = self.observable.generate_vector(PARAMETER) 
+			LocalErrors = hp.MultiVector(self.ms[0],self.parameters['error_test_samples'])
+			projection_vector = self.observable.generate_vector(hp.PARAMETER) 
 
 			for rank_index,rank in enumerate(ranks):
 				LocalErrors.zero()
@@ -809,7 +1021,7 @@ def test_errors(self,test_input = True, test_output = False, ranks = [None],cut_
 					V_GN = self.V_GN
 					d_GN = self.d_GN
 				else:
-					V_GN = MultiVector(self.V_GN[0],rank)
+					V_GN = hp.MultiVector(self.V_GN[0],rank)
 					d_GN = self.d_GN[0:rank]
 					for i in range(rank):
 						V_GN[i].axpy(1.,self.V_GN[i])
@@ -817,7 +1029,7 @@ def test_errors(self,test_input = True, test_output = False, ranks = [None],cut_
 				if self.prior_preconditioned:
 					InputProjectorOperator = PriorPreconditionedProjector(V_GN,self.prior.R, input_init_vector_lambda)
 				else:
-					InputProjectorOperator = LowRankOperator(np.ones_like(d_GN),V_GN, input_init_vector_lambda)
+					InputProjectorOperator = hp.LowRankOperator(np.ones_like(d_GN),V_GN, input_init_vector_lambda)
 			
 				rel_errors = np.zeros(LocalErrors.nvec())
 				for i in range(LocalErrors.nvec()):
@@ -837,7 +1049,7 @@ def test_errors(self,test_input = True, test_output = False, ranks = [None],cut_
 					print('Naive global average relative error input = ',global_avg_rel_errors_input[rank_index],' for rank ',rank)
 
 			# Double Loop MC Error test does not work when prior preconditioning is used.
-			# This will be fixed soon.
+			# This will be fixed soon. Or someday. Someday soon! :+)
 			if False:
 				if self.d_GN is None:
 					if self.mesh_constructor_comm.rank == 0:
@@ -864,21 +1076,21 @@ def test_errors(self,test_input = True, test_output = False, ranks = [None],cut_
 				# Instantiate input and output data arrays
 				observable_vector = dl.Vector(self.mesh_constructor_comm)
 				self.observable.init_vector(observable_vector,dim = 0)
-				LocalObservables = MultiVector(observable_vector,self.parameters['error_test_samples'])
+				LocalObservables = hp.MultiVector(observable_vector,self.parameters['error_test_samples'])
 				LocalObservables.zero()
-				LocalParameters = MultiVector(self.ms[0],self.parameters['error_test_samples'])
+				LocalParameters = hp.MultiVector(self.ms[0],self.parameters['error_test_samples'])
 				LocalParameters.zero()
 				# Generate samples
 				for i in range(self.parameters['error_test_samples']):
 					t0 = time.time()
-					parRandom.normal(1,self.noise)
+					hp.parRandom.normal(1,self.noise)
 					self.prior.sample(self.noise,LocalParameters[i])
 					LocalObservables[i].axpy(1.,self.observable.eval(LocalParameters[i],setLinearizationPoint = True))
 					if self.parameters['verbose'] and (self.mesh_constructor_comm.rank == 0):
 						print('Generating local observable ',i,' for input error test took',time.time() -t0, 's')
 
 				# Instantiate array for errors
-				LocalErrors = MultiVector(observable_vector,self.parameters['error_test_samples'])
+				LocalErrors = hp.MultiVector(observable_vector,self.parameters['error_test_samples'])
 
 				m_r = dl.Vector(self.mesh_constructor_comm)
 				self.observable.init_vector(m_r,dim = 1)
@@ -896,7 +1108,7 @@ def test_errors(self,test_input = True, test_output = False, ranks = [None],cut_
 						V_GN = self.V_GN
 						d_GN = self.d_GN
 					else:
-						V_GN = MultiVector(self.V_GN[0],rank)
+						V_GN = hp.MultiVector(self.V_GN[0],rank)
 						V_GN.zero()
 						d_GN = self.d_GN[0:rank]
 						for i in range(rank):
@@ -905,7 +1117,7 @@ def test_errors(self,test_input = True, test_output = False, ranks = [None],cut_
 					if self.prior_preconditioned:
 						InputProjectorOperator = PriorPreconditionedProjector(V_GN,self.prior.R, input_init_vector_lambda)
 					else:
-						InputProjectorOperator = LowRankOperator(np.ones_like(d_GN),V_GN, input_init_vector_lambda)
+						InputProjectorOperator = hp.LowRankOperator(np.ones_like(d_GN),V_GN, input_init_vector_lambda)
 					print('Constructed projection operator for rank ',rank)
 					rel_errors = np.zeros(LocalErrors.nvec())
 
@@ -922,7 +1134,7 @@ def test_errors(self,test_input = True, test_output = False, ranks = [None],cut_
 							t0 = time.time()
 							y.zero()
 							y_r.zero()
-							parRandom.normal(1,self.noise)
+							hp.parRandom.normal(1,self.noise)
 							self.prior.sample(self.noise,y)
 							InputProjectorOperator.mult(y,y_r)
 							y.axpy(-1., y_r)
@@ -980,19 +1192,30 @@ def test_errors(self,test_input = True, test_output = False, ranks = [None],cut_
 			# Naive test on output space
 			observable_vector = dl.Vector(self.mesh_constructor_comm)
 			self.observable.init_vector(observable_vector,dim = 0)
-			LocalObservables = MultiVector(observable_vector,self.parameters['error_test_samples'])
+			LocalObservables = hp.MultiVector(observable_vector,self.parameters['error_test_samples'])
 			LocalObservables.zero()
 			for i in range(LocalObservables.nvec()):
 				t0 = time.time()
-				parRandom.normal(1,self.noise)
+				hp.parRandom.normal(1,self.noise)
 				self.prior.sample(self.noise,self.ms[0])
-				x = [self.us[0],self.ms[0],None]
+				if self.control_distribution is not None:
+					assert self.zs is not None
+					if type(self.zs) is list:
+						assert self.zs[0] is not None
+						self.control_distribution.sample(self.zs[0])
+						x = [self.us[0],self.ms[0],None,self.zs[0]]
+					else:
+						assert self.z is not None
+						self.control_distribution.sample(self.z)
+						x = [self.us[0],self.ms[0],None,self.z]
+				else:
+					x = [self.us[0],self.ms[0],None]
 				self.observable.setLinearizationPoint(x)
 				LocalObservables[i].axpy(1.,self.observable.eval(self.ms[0]))
 				if self.parameters['verbose'] and (self.mesh_constructor_comm.rank == 0):
 					print('Generating local observable ',i, ' for output error test')
 
-			LocalErrors = MultiVector(observable_vector,self.parameters['error_test_samples'])
+			LocalErrors = hp.MultiVector(observable_vector,self.parameters['error_test_samples'])
 			projection_vector = dl.Vector(self.mesh_constructor_comm)
 			self.observable.init_vector(projection_vector,dim = 0)
 
@@ -1002,12 +1225,12 @@ def test_errors(self,test_input = True, test_output = False, ranks = [None],cut_
 					U_NG = self.U_NG
 					d_NG = self.d_NG
 				else:
-					U_NG = MultiVector(self.U_NG[0],rank)
+					U_NG = hp.MultiVector(self.U_NG[0],rank)
 					d_NG = self.d_NG[0:rank]
 					for i in range(rank):
 						U_NG[i].axpy(1.,self.U_NG[i])
 				output_init_vector_lambda = lambda x, dim: self.observable.init_vector(x,dim = 0)
-				OutputProjectorOperator = LowRankOperator(np.ones_like(d_NG),U_NG,output_init_vector_lambda)
+				OutputProjectorOperator = hp.LowRankOperator(np.ones_like(d_NG),U_NG,output_init_vector_lambda)
 			
 				rel_errors = np.zeros(LocalErrors.nvec())
 				for i in range(LocalErrors.nvec()):
diff --git a/hippyflow/modeling/blockVector.py b/hippyflow/modeling/blockVector.py
index 30e928b..bfc750a 100644
--- a/hippyflow/modeling/blockVector.py
+++ b/hippyflow/modeling/blockVector.py
@@ -15,7 +15,7 @@
 
 
 import dolfin as dl
-from hippylib import *
+import hippylib as hp 
 
 class BlockVector(object):
     """
@@ -58,7 +58,7 @@ def randn_perturb(self,std_dev):
         to each of the snapshots.
         """
         for d in self.data:
-            parRandom.normal_perturb(std_dev, d)
+            hp.parRandom.normal_perturb(std_dev, d)
 
     
     def axpy(self, a, other):
@@ -92,5 +92,5 @@ def copy(self):
     
     def export(self,Vh, fid, xname):
         for xi in self.data:
-            xfun = vector2Function(xi,Vh, name=xname)
+            xfun = hp.vector2Function(xi,Vh, name=xname)
             fid << xfun
diff --git a/hippyflow/modeling/cMinimization.py b/hippyflow/modeling/cMinimization.py
index 552e804..2168fec 100644
--- a/hippyflow/modeling/cMinimization.py
+++ b/hippyflow/modeling/cMinimization.py
@@ -19,7 +19,6 @@
 import math
 
 import sys, os
-sys.path.append( os.environ.get('HIPPYLIB_BASE_DIR'))
 import hippylib as hp
 
 
diff --git a/hippyflow/modeling/controlJacobian.py b/hippyflow/modeling/controlJacobian.py
new file mode 100644
index 0000000..9c435dc
--- /dev/null
+++ b/hippyflow/modeling/controlJacobian.py
@@ -0,0 +1,95 @@
+# Copyright (c) 2020-2022, The University of Texas at Austin 
+# & Washington University in St. Louis.
+#
+# All Rights reserved.
+# See file COPYRIGHT for details.
+#
+# This file is part of the hIPPYflow package. For more information see
+# https://github.com/hippylib/hippyflow/
+#
+# hIPPYflow 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) version 2.0 dated June 1991.
+
+
+import dolfin as dl
+import numpy as np
+import hippylib as hp
+from .jacobian import Jacobian
+
+CONTROL = 3
+
+class ObservableControlJacobian(Jacobian):
+    """
+    This class implements matrix free application of the Jacobian operator.
+    The constructor takes the following parameters:
+    - :code:`model`:               the object which contains the description of the problem.
+    
+    Type :code:`help(modelTemplate)` for more information on which methods model should implement.
+    """
+    def __init__(self, observable):
+        """
+        Construct the Observable Jacobian operator
+        """
+        self.observable = observable
+
+        assert hasattr(self.observable,'applyCz')
+        assert hasattr(self.observable,'applyCzt')
+
+        self.ncalls = 0
+        
+        self.rhs_fwd = observable.generate_vector(hp.STATE)
+        self.rhs_adj = observable.generate_vector(hp.ADJOINT)
+        self.rhs_adj2 = observable.generate_vector(hp.ADJOINT)
+        self.uhat    = observable.generate_vector(hp.STATE)
+        self.phat    = observable.generate_vector(hp.ADJOINT)
+        self.yhelp = observable.generate_vector(CONTROL)
+
+        self.Bu = dl.Vector(self.mpi_comm())
+        self.observable.B.init_vector(self.Bu,0)
+        self.Ctphat = observable.generate_vector(CONTROL)
+        self.shape = (self.Bu.get_local().shape[0],self.yhelp.get_local().shape[0])
+
+    def mpi_comm(self):
+        return self.observable.B.mpi_comm()
+    
+    def init_vector(self, x, dim):
+        """
+        Reshape the Vector :code:`x` so that it is compatible with the reduced Hessian
+        operator.
+        Parameters:
+        - :code:`x`: the vector to reshape.
+        - :code:`dim`: if 0 then :code:`x` will be made compatible with the range of the Jacobian, if 1 then :code:`x` will be made compatible with the domain of the Jacobian.
+               
+         """
+        if dim == 0:
+            self.observable.init_vector(x,0)
+        elif dim == 1:
+            # CONTROL = 3
+            self.observable.init_vector(x,3)
+
+
+        
+    def mult(self,x,y):
+        """
+        Apply the Jacobian :code:`x`. Return the result in :code:`y`.
+        Implemented for dl.Vector
+        """
+        self.observable.applyCz(x, self.rhs_fwd)
+        self.observable.solveFwdIncremental(self.uhat, self.rhs_fwd)
+        assert hasattr(self.observable,'applyB'), 'LinearObservable must have attribute applyB'
+        self.observable.applyB(self.uhat,y)
+        y *= -1.0
+        self.ncalls += 1
+
+    def transpmult(self,x,y):
+        """
+        Apply the Jacobian transpose :code:`x`. Return the result in :code:`y`.
+        Implemented for dl.Vector
+        """
+        assert hasattr(self.observable,'applyBt'), 'LinearObservable must have attribute applyBt'
+        self.observable.applyBt(x,self.rhs_adj)
+        self.observable.solveAdjIncremental(self.phat, self.rhs_adj)
+        self.observable.applyCzt(self.phat, y)
+        y *= -1.0
+        self.ncalls += 1
diff --git a/hippyflow/modeling/hippylibModelWrapper.py b/hippyflow/modeling/hippylibModelWrapper.py
index efee41d..16318c2 100644
--- a/hippyflow/modeling/hippylibModelWrapper.py
+++ b/hippyflow/modeling/hippylibModelWrapper.py
@@ -113,7 +113,7 @@ def setLinearizationPoint(self, x):
 			- or partially precompute the block of the hessian (if feasible)
 		"""
 		x[hp.ADJOINT] = self.model.generate_vector(hp.ADJOINT)
-		self.problem.setLinearizationPoint(x, True)
+		self.model.problem.setLinearizationPoint(x, True)
 
 
 	def evalVariationalGradient(self,x,u0 = None,setLinearizationPoint = False,misfit_only = True):
@@ -319,6 +319,8 @@ def evalJacobian(self,x,u0 = None,linearizationPointSet=False):
 
 
 	def samplePrior(self):
+		"""
+		"""
 		if self.noise_help is None:
 			self.noise_help = dl.Vector()
 			self.model.prior.init_vector(self.noise_help,"noise")
@@ -335,7 +337,7 @@ def samplePrior(self):
 
 		return self.sample_help
 
-	def setUpInverseProblem(self):
+	def setUpInverseProblem(self,mtrue_expr = None):
 		"""
 		"""
 		assert self.settings['rel_noise'] is not None
@@ -343,9 +345,16 @@ def setUpInverseProblem(self):
 		if self.mtrue is None:
 			self.mtrue = dl.Vector()
 			self.model.prior.init_vector(self.mtrue,0)
-		self.mtrue.zero()
-		# self.mtrue.axpy(1.0,self.samplePrior())
-		self.mtrue.set_local(self.samplePrior().get_local())
+
+		if mtrue_expr is None:
+			self.mtrue.zero()
+			# self.mtrue.axpy(1.0,self.samplePrior())
+			self.mtrue.set_local(self.samplePrior().get_local())
+		else:
+			mtrue_func = dl.Function(self.model.problem.Vh[hp.PARAMETER])
+			mtrue_func.interpolate(mtrue_expr)
+			self.mtrue.zero()
+			self.mtrue.axpy(1.,mtrue_func.vector())
 
 		self.model.misfit.d.zero()
 
diff --git a/hippyflow/modeling/jacobian.py b/hippyflow/modeling/jacobian.py
index ee298e4..573b6f6 100644
--- a/hippyflow/modeling/jacobian.py
+++ b/hippyflow/modeling/jacobian.py
@@ -12,7 +12,7 @@
 # Software Foundation) version 2.0 dated June 1991.
 
 
-from hippylib import *
+import hippylib as hp 
 import dolfin as dl
 import numpy as np
        
@@ -76,18 +76,18 @@ def __init__(self, observable):
 
         self.ncalls = 0
         
-        self.rhs_fwd = observable.generate_vector(STATE)
-        self.rhs_adj = observable.generate_vector(ADJOINT)
-        self.rhs_adj2 = observable.generate_vector(ADJOINT)
-        self.uhat    = observable.generate_vector(STATE)
-        self.phat    = observable.generate_vector(ADJOINT)
-        self.yhelp = observable.generate_vector(PARAMETER)
+        self.rhs_fwd = observable.generate_vector(hp.STATE)
+        self.rhs_adj = observable.generate_vector(hp.ADJOINT)
+        self.rhs_adj2 = observable.generate_vector(hp.ADJOINT)
+        self.uhat    = observable.generate_vector(hp.STATE)
+        self.phat    = observable.generate_vector(hp.ADJOINT)
+        self.yhelp = observable.generate_vector(hp.PARAMETER)
 
 
 
         self.Bu = dl.Vector(self.mpi_comm())
         self.observable.B.init_vector(self.Bu,0)
-        self.Ctphat = observable.generate_vector(PARAMETER)
+        self.Ctphat = observable.generate_vector(hp.PARAMETER)
         self.shape = (self.Bu.get_local().shape[0],self.yhelp.get_local().shape[0])
 
     def mpi_comm(self):
diff --git a/hippyflow/modeling/maternPrior.py b/hippyflow/modeling/maternPrior.py
index 983929d..1362c15 100644
--- a/hippyflow/modeling/maternPrior.py
+++ b/hippyflow/modeling/maternPrior.py
@@ -14,28 +14,28 @@
 
 import dolfin as dl
 import numpy as np
-from hippylib import *
+import hippylib as hp
 
 def BiLaplacian2D(Vh_parameter,gamma = 0.1,delta = 0.1, theta0 = 2.0, theta1 = 0.5,\
 			 alpha = np.pi/4,mean = None,robin_bc = False):
 	"""
 	Return 2D BiLaplacian prior given function space and coefficients for Matern covariance
 	"""
-	anis_diff = dl.CompiledExpression(ExpressionModule.AnisTensor2D(), degree = 1)
+	anis_diff = dl.CompiledExpression(hp.ExpressionModule.AnisTensor2D(), degree = 1)
 	anis_diff.theta0 = theta0
 	anis_diff.theta1 = theta1
 	anis_diff.alpha = alpha	
 
-	return BiLaplacianPrior(Vh_parameter, gamma, delta, anis_diff,mean = mean,robin_bc = robin_bc)
+	return hp.BiLaplacianPrior(Vh_parameter, gamma, delta, anis_diff,mean = mean,robin_bc = robin_bc)
 
 
 def Laplacian2D(Vh_parameter,gamma = 0.1,delta = 0.1, theta0 = 2.0, theta1 = 0.5, alpha = np.pi/4,mean = None):
 	"""
 	Return 2D Laplacian prior given function space and coefficients for Matern covariance
 	"""
-	anis_diff = dl.CompiledExpression(ExpressionModule.AnisTensor2D(), degree = 1)
+	anis_diff = dl.CompiledExpression(hp.ExpressionModule.AnisTensor2D(), degree = 1)
 	anis_diff.theta0 = theta0
 	anis_diff.theta1 = theta1
 	anis_diff.alpha = alpha	
 
-	return LaplacianPrior(Vh_parameter, gamma, delta, mean = mean)
\ No newline at end of file
+	return hp.LaplacianPrior(Vh_parameter, gamma, delta, mean = mean)
\ No newline at end of file
diff --git a/hippyflow/modeling/multiPDEProblem.py b/hippyflow/modeling/multiPDEProblem.py
index 4e0e64d..791abca 100644
--- a/hippyflow/modeling/multiPDEProblem.py
+++ b/hippyflow/modeling/multiPDEProblem.py
@@ -13,12 +13,12 @@
 
 
 import dolfin as dl
-from hippylib import *
+import hippylib as hp 
 from .blockVector import BlockVector
 
 
 
-class MultiPDEProblem(PDEProblem):
+class MultiPDEProblem(hp.PDEProblem):
     def __init__(self, problems):
         self.Vh = problems[0].Vh
         self.problems = problems
@@ -28,7 +28,7 @@ def generate_state(self):
         """ 
         Return a vector in the shape of the state 
         """
-        return BlockVector(self.Vh[STATE], self.n_problems)
+        return BlockVector(self.Vh[hp.STATE], self.n_problems)
         
     def generate_parameter(self):
         return self.problems[0].generate_parameter()
@@ -121,9 +121,9 @@ def apply_ij(self,i,j, dir, out):
 
         """
         out.zero()
-        if i == PARAMETER:
+        if i == hp.PARAMETER:
             tmp = self.generate_parameter()
-            if j == PARAMETER:
+            if j == hp.PARAMETER:
                 for k in range(self.n_problems):
                     self.problems[k].apply_ij(i,j, dir,tmp)
                     out.axpy(1., tmp)
@@ -133,7 +133,7 @@ def apply_ij(self,i,j, dir, out):
                     out.axpy(1., tmp)
         else:
             assert type(out) is BlockVector
-            if j == PARAMETER:
+            if j == hp.PARAMETER:
                 for k in range(self.n_problems):
                     self.problems[k].apply_ij(i,j, dir,out.data[k] )
             else:
diff --git a/hippyflow/modeling/multiStateLinearObservable.py b/hippyflow/modeling/multiStateLinearObservable.py
index 287aed1..9d60978 100644
--- a/hippyflow/modeling/multiStateLinearObservable.py
+++ b/hippyflow/modeling/multiStateLinearObservable.py
@@ -11,7 +11,7 @@
 # terms of the GNU General Public License (as published by the Free
 # Software Foundation) version 2.0 dated June 1991.
 
-from hippylib import *
+import hippylib as hp 
 import dolfin as dl
 import numpy as np
       
@@ -62,11 +62,11 @@ def generate_vector(self, component = "ALL"):
 			x = [self.problem.generate_state(),
 				 self.problem.generate_parameter(),
 				 self.problem.generate_state()]
-		elif component == STATE:
+		elif component == hp.STATE:
 			x = self.problem.generate_state()
-		elif component == PARAMETER:
+		elif component == hp.PARAMETER:
 			x = self.problem.generate_parameter()
-		elif component == ADJOINT:
+		elif component == hp.ADJOINT:
 			x = self.problem.generate_state()
 
 		return x
@@ -153,7 +153,7 @@ def setLinearizationPoint(self, x):
 			- simply store a copy of x and evaluate action of blocks of the Hessian on the fly
 			- or partially precompute the block of the hessian (if feasible)
 		"""
-		x[ADJOINT] = self.problem.problems[0].generate_state()
+		x[hp.ADJOINT] = self.problem.problems[0].generate_state()
 		self.problem.setLinearizationPoint(x, True)
 
 		
@@ -197,7 +197,7 @@ def applyC(self, dm, out):
 			
 		.. note:: This routine assumes that :code:`out` has the correct shape.
 		"""
-		self.problem.apply_ij(ADJOINT,PARAMETER, dm, out)
+		self.problem.apply_ij(hp.ADJOINT,hp.PARAMETER, dm, out)
 	
 	def applyCt(self, dp, out):
 		"""
@@ -209,4 +209,4 @@ def applyCt(self, dp, out):
 			
 		..note:: This routine assumes that :code:`out` has the correct shape.
 		"""
-		self.problem.apply_ij(PARAMETER,ADJOINT, dp, out)
+		self.problem.apply_ij(hp.PARAMETER,hp.ADJOINT, dp, out)
diff --git a/hippyflow/modeling/observable.py b/hippyflow/modeling/observable.py
index 63ad506..97dc65a 100644
--- a/hippyflow/modeling/observable.py
+++ b/hippyflow/modeling/observable.py
@@ -11,7 +11,7 @@
 # terms of the GNU General Public License (as published by the Free
 # Software Foundation) version 2.0 dated June 1991.
 
-from hippylib import *
+import hippylib as hp 
 import dolfin as dl
 import numpy as np
 
@@ -76,6 +76,7 @@ def __init__(self, problem, B):
 			- prior: the prior 
 		"""
 		self.problem = problem
+		self.is_control_problem = hasattr(self.problem,'Cz')
 		self.B = B
 		# self.Bu = dl.Vector(self.B.mpi_comm())
 		# self.B.init_vector(self.Bu,0)
@@ -114,14 +115,14 @@ def generate_vector(self, component = "ALL"):
 				x = [self.problem.generate_state(),
 					 self.problem.generate_parameter(),
 					 self.problem.generate_state()]
-		elif component == STATE:
+		elif component == hp.STATE:
 			x = self.problem.generate_state()
-		elif component == PARAMETER:
+		elif component == hp.PARAMETER:
 			x = self.problem.generate_parameter()
-		elif component == ADJOINT:
+		elif component == hp.ADJOINT:
 			x = self.problem.generate_state()
-		elif component == 3:
-			assert hasattr(self.problem,'Cz'), 'Assuming it is control problem'
+		elif component == CONTROL:
+			assert self.is_control_problem, 'Assuming it is a control problem'
 			# 3 denotes a control variable needs to be generated
 			x = self.problem.generate_control()
 
@@ -144,6 +145,7 @@ def init_vector(self, x, dim):
 			# self.model.init_parameter(x)
 			self.problem.C.init_vector(x,1)
 		elif dim == 3:
+			assert self.is_control_problem, 'Assuming it is a control problem'
 			self.problem.Cz.init_vector(x,1)
 		else: 
 			raise
@@ -158,7 +160,7 @@ def init_parameter(self, m):
 		self.problem.init_parameter(m)
 
 			
-	def eval(self, m, u0 = None,setLinearizationPoint = False):
+	def eval(self, m, u0 = None, z = None, setLinearizationPoint = False):
 		"""
 		Given the input parameter :code:`m` solve for the state field $u(m)$, and evaluate 
 		the linear state observable $Bu(m)$
@@ -168,7 +170,11 @@ def eval(self, m, u0 = None,setLinearizationPoint = False):
 		"""
 		if u0 is None:
 			u0 = self.problem.generate_state()
-		x = [u0, m, None]
+		if self.is_control_problem:
+			assert z is not None
+			x = [u0,m,None,z]
+		else:
+			x = [u0, m, None]
 		self.problem.solveFwd(u0,x)
 		if setLinearizationPoint:
 			self.setLinearizationPoint(x)
@@ -218,7 +224,7 @@ def setLinearizationPoint(self, x):
 			- simply store a copy of x and evaluate action of blocks of the Hessian on the fly
 			- or partially precompute the block of the hessian (if feasible)
 		"""
-		x[ADJOINT] = self.problem.generate_state()
+		x[hp.ADJOINT] = self.problem.generate_state()
 		self.problem.setLinearizationPoint(x, True)
 
 		
@@ -263,7 +269,7 @@ def applyC(self, dm, out):
 			
 		.. note:: This routine assumes that :code:`out` has the correct shape.
 		"""
-		self.problem.apply_ij(ADJOINT,PARAMETER, dm, out)
+		self.problem.apply_ij(hp.ADJOINT,hp.PARAMETER, dm, out)
 	
 	def applyCt(self, dp, out):
 		"""
@@ -275,7 +281,7 @@ def applyCt(self, dp, out):
 			
 		..note:: This routine assumes that :code:`out` has the correct shape.
 		"""
-		self.problem.apply_ij(PARAMETER,ADJOINT, dp, out)
+		self.problem.apply_ij(hp.PARAMETER,hp.ADJOINT, dp, out)
 
 	def applyCz(self, dz, out):
 		"""
@@ -289,7 +295,7 @@ def applyCz(self, dz, out):
 			
 		.. note:: This routine assumes that :code:`out` has the correct shape.
 		"""
-		self.problem.apply_ij(ADJOINT,CONTROL, dz, out)
+		self.problem.apply_ij(hp.ADJOINT,CONTROL, dz, out)
 	
 	def applyCzt(self, dp, out):
 		"""
@@ -301,7 +307,7 @@ def applyCzt(self, dp, out):
 			
 		..note:: This routine assumes that :code:`out` has the correct shape.
 		"""
-		self.problem.apply_ij(CONTROL,ADJOINT, dp, out)
+		self.problem.apply_ij(CONTROL,hp.ADJOINT, dp, out)
 
 def hippylibModelLinearStateObservable(model):
 	"""
diff --git a/hippyflow/modeling/operatorWrappers.py b/hippyflow/modeling/operatorWrappers.py
index 45641e0..6140de1 100644
--- a/hippyflow/modeling/operatorWrappers.py
+++ b/hippyflow/modeling/operatorWrappers.py
@@ -51,4 +51,61 @@ def transpmult(self,x,y):
 		y.zero()
 		y.set_local(self.matrix.T@x.get_local())
 
+
+class MeanJTJfromDataOperator:
+	"""
+	"""
+	def __init__(self, J, prior):
+		"""
+		"""
+
+		# Assumes J.shape = (ndata,rank,dM)
+		self._J = J
+		self.ndata, self.r, self.dM = self.J.shape
+		self._prior = prior
+
+		if hasattr(self.prior, "R"):
+			self.init_vector_lambda = lambda x,dim: prior.R.init_vector(x,dim)
+		else:
+			self.init_vector_lambda = lambda x,dim: prior.Hlr.init_vector(x,dim)
+
+	@property
+	def J(self):
+		return self._J
+
+	@property
+	def prior(self):
+		return self._prior
+	
+	def init_vector(self,x,dim):
+		"""
+		Initialize :code:`x` to be compatible with the range (:code:`dim=0`) or domain (:code:`dim=1`) of :code:`A`.
+		"""
+		assert init_vector_lambda is not None
+		self.init_vector_lambda(x,dim)
+
+	def mult(self,x,y):
+		"""
+		Compute :math:`y = mean(JTJ)x `
+		"""
+		x_np = x.get_local()
+		# print('x_np.shape = ',x_np.shape)
+		X_np = np.tile(x_np,(self.ndata,1))
+		# print('X_np.shape = ',X_np.shape)
+		assert X_np.shape == (self.ndata,self.dM)
+		JX_np = np.einsum('ijk,ik->ij',self.J,X_np)
+		# print('PhiTJX_np.shape = ',JX_np.shape)
+		JTJX_np = np.einsum('ijk,ij->ik',self.J,JX_np)
+		# print('JTPhiPhiTJX_np.shape = ',JTJX_np.shape)
+		y.set_local(np.mean(JTJX_np,axis = 0))
+
+	def transpmult(self,x,y):
+		"""
+		Compute :math:`y = mean(JTJ)x `
+		JTJ is naturally self-adjoint so tranpsmult is mult.
+		"""
+		return self.mult(x,y)
+
+
+
 		
\ No newline at end of file
diff --git a/hippyflow/modeling/priorPreconditionedProjector.py b/hippyflow/modeling/priorPreconditionedProjector.py
index 7c752c9..9347cb6 100644
--- a/hippyflow/modeling/priorPreconditionedProjector.py
+++ b/hippyflow/modeling/priorPreconditionedProjector.py
@@ -11,7 +11,7 @@
 # terms of the GNU General Public License (as published by the Free
 # Software Foundation) version 2.0 dated June 1991.
 
-from hippylib import *
+import hippylib as hp 
 import dolfin as dl
 import numpy as np
 
diff --git a/hippyflow/utilities/plotting.py b/hippyflow/utilities/plotting.py
index bda736d..49e63c7 100644
--- a/hippyflow/utilities/plotting.py
+++ b/hippyflow/utilities/plotting.py
@@ -14,10 +14,13 @@
 from matplotlib.backends.backend_pdf import PdfPages
 from matplotlib import rc
 from pylab import *
-plt.rc('text', usetex=True)
-plt.rc('font', family='serif')
-plt.rc('text.latex', preamble=r'\usepackage{amsfonts}')
-
+try:
+    plt.rc('text', usetex=True)
+    plt.rc('font', family='serif')
+    plt.rc('text.latex', preamble=r'\usepackage{amsfonts}')
+except:
+    print('Error loading latex, will not be used in plots')
+    pass
 
 def spectrum_plot(lambdas, axis_label = ['i','$\lambda$','Spectrum'],\
                  ylims = None, out_name = None):