<commit_before>"""
An implementation of overcomplete ISA using Gaussian scale mixtures.
"""

__license__ = 'MIT License <http://www.opensource.org/licenses/mit-license.php>'
__author__ = '<NAME> <<EMAIL>>'
__docformat__ = 'epytext'

from distribution import Distribution
from numpy import *
from numpy import min, max, round
from numpy.random import randint, randn, rand, logseries, permutation
from numpy.linalg import svd, pinv, inv, det, slogdet, cholesky
from scipy.linalg import solve
from scipy.optimize import fmin_l_bfgs_b, fmin_cg, check_grad
from scipy.stats import laplace, t
from tools import gaborf, mapp, logmeanexp, asshmarray, sqrtmi, sqrtm
from gsm import GSM
from copy import deepcopy

class ISA(Distribution):
	"""
	An implementation of overcomplete ISA using Gaussian scale mixtures.
	"""

	def __init__(self, num_visibles, num_hiddens=None, ssize=1, noise=False):
		"""
		@type  num_visibles: integer
		@param num_visibles: data dimensionality

		@type  num_hiddens: integer
		@param num_hiddens: number of hidden units

		@type  ssize: integer
		@param ssize: subspace dimensionality

		@type  noise: bool
		@param noise: add additional hidden units for noise
		"""

		if num_hiddens is None:
			num_hiddens = num_visibles

		self.dim = num_visibles
		self.num_visibles = num_visibles
		self.num_hiddens = num_hiddens

		# random linear feature 
		self.A = randn(self.num_visibles, self.num_hiddens) / 10.

		# subspace densities
		self.subspaces = [
			GSM(ssize) for _ in range(int(num_hiddens) / int(ssize))]

		if mod(num_hiddens, ssize) > 0:
			self.subspaces.append(GSM(mod(num_hiddens, ssize)))

		self._noise = False
		self.noise = noise



	def initialize(self, X=None, method='data'):
		"""
		Initializes parameter values with more sensible values.

		@type  X: array_like
		@param X: data points stored in columns

		@type  method: string
		@param method: type of initialization ('data', 'gabor' or 'random')
		"""

		if method.lower() == 'data':
			# initialize features with data points
			if X is not None:
				# center data
				self.A = X[:, randint(X.shape[1], size=self.num_hiddens)] / sqrt(self.num_hiddens)
				self.A /= sqrt(sum(square(self.A), 0))

		elif method.lower() == 'gabor':
			# initialize features with Gabor filters
			if self.subspaces[0].dim > 1 and not mod(self.num_hiddens, 2):
				for i in range(self.num_hiddens / 2):
					G = gaborf(self.num_visibles)
					self.A[:, 2 * i] = real(G)
					self.A[:, 2 * i + 1] = imag(G)
			else:
				for i in range(len(self.subspaces)):
					self.A[:, i] = gaborf(self.num_visibles, complex=False)

		elif method.lower() == 'random':
			# initialize with Gaussian white noise
			self.A = randn(num_visibles, num_hiddens)

		elif method.lower() == 'laplace' or method.lower() == 'student':
			if method.lower() == 'laplace':
				# approximate multivariate Laplace with GSM
				samples = randn(self.subspaces[0].dim, 10000)
				samples = samples / sqrt(sum(square(samples), 0))
				samples = laplace.rvs(size=[1, 10000]) * samples
			else:
				samples = randn(self.subspaces[0].dim, 10000)
				samples = samples / sqrt(sum(square(samples), 0))
				samples = t.rvs(2., size=[1, 10000]) * samples

			if self.noise:
				# ignore first subspace
				gsm = GSM(self.subspaces[1].dim, self.subspaces[1].num_scales)
				gsm.train(samples, max_iter=200, tol=1e-8)

				for m in self.subspaces[1:]:
					m.scales = gsm.scales.copy()
			else:
				gsm = GSM(self.subspaces[0].dim, self.subspaces[0].num_scales)
				gsm.train(samples, max_iter=200, tol=1e-8)

				for m in self.subspaces:
					m.scales = gsm.scales.copy()

		else:
			raise ValueError('Unknown initialization method \'{0}\'.'.format(method))



	def train(self, X, method=('sgd', {}), sampling_method=('gibbs', {}), **kwargs):
		"""
		Optimizes model parameters with respect to the log-likelihoood.

		@type  max_iter: integer
		@param max_iter: maximum number of iterations through the dataset

		@type  method: tuple
		@param method: optimization method used to optimize filters

		@type  adative: bool
		@param adaptive: automatically adjust step width when using SGD (default: True)

		@type  train_prior: bool
		@param train_prior: whether or not to optimize the marginal distributions (default: True)

		@type  train_subspaces: bool
		@param train_subspaces: automatically learn subspace sizes (default: False)

		@type  orthogonalize: bool
		@param orthogonalize: after each step, orthogonalize rows of feature matrix (default: False)

		@type  sampling_method: tuple
		@param sampling_method: method and parameters to generate hidden representations

		@type  init_sampling_steps: integer
		@param init_sampling_steps: number of steps used to initialize persistent samples (default: 0)
		"""

		max_iter = kwargs.get('max_iter', 100)
		adaptive = kwargs.get('adaptive', True) 
		train_prior = kwargs.get('train_prior', True)
		train_subspaces = kwargs.get('train_subspaces', False)
		orthogonalize = kwargs.get('orthogonalize', False)
		persistent = kwargs.get('persistent', True)
		init_sampling_steps = kwargs.get('init_sampling_steps', 0)

		if Distribution.VERBOSITY > 0:
			if self.num_hiddens > self.num_visibles:
				print 0
			else:
				print 0, self.evaluate(X)

		if isinstance(method, str):
			method = (method, {})

		if isinstance(sampling_method, str):
			sampling_method = (sampling_method, {})

		if persistent and init_sampling_steps:
			# initialize samples
			sampling_method[1]['Z'] = self.sample_nullspace(X, 
				method=(sampling_method[0], 
					dict(sampling_method[1], num_steps=init_sampling_steps)))

		if adaptive and 'step_width' not in method[1]:
			method[1]['step_width'] = 0.001

		for i in range(max_iter):
			# complete data (E)
			Y = self.sample_posterior(X, method=sampling_method)

			if train_prior:
				# optimize parameters of the prior (M)
				self.train_prior(Y)

			if train_subspaces:
				# learn subspaces (M)
				Y = self.train_subspaces(Y)

			if not orthogonalize and train_prior or train_subspaces:
				# normalize variances of marginals
				self.normalize_prior()

			if persistent:
				# initializes samples in next iteration
				sampling_method[1]['Z'] = dot(self.nullspace_basis(), Y)

			# optimize linear features (M)
			if method[0].lower() == 'sgd':
				improved = self.train_sgd(Y, **method[1])

				if adaptive:
					# adjust learning rate
					method[1]['step_width'] *= 1.1 if improved else 0.5

			elif method[0].lower() == 'lbfgs':
				self.train_lbfgs(Y, **method[1])

			else:
				raise ValueError('Unknown training method \'{0}\'.'.format(method[0]))

			if orthogonalize:
				# normalize feature matrix
				self.orthogonalize()

			if Distribution.VERBOSITY > 0:
				if self.num_hiddens > self.num_visibles:
					print i + 1
				else:
					print i + 1, self.evaluate(X)



	def train_prior(self, Y, **kwargs):
		"""
		Optimize parameters of the marginal distribution over the hidden variables.
		The parameters are fit to maximize the average log-likelihood of the
		columns in `Y`.

		@type  Y: array_like
		@param Y: hidden states
		"""

		offset = [0]
		for model in self.subspaces:
			offset.append(offset[-1] + model.dim)
<commit_msg>Changed heuristic for initializing with data points. Added callback function to training method.<commit_after>"""
An implementation of overcomplete ISA using Gaussian scale mixtures.
"""

__license__ = 'MIT License <http://www.opensource.org/licenses/mit-license.php>'
__author__ = '<NAME> <<EMAIL>>'
__docformat__ = 'epytext'

from distribution import Distribution
from numpy import *
from numpy import min, max, round
from numpy.random import randint, randn, rand, logseries, permutation
from numpy.linalg import svd, pinv, inv, det, slogdet, cholesky, eig
from scipy.linalg import solve
from scipy.optimize import fmin_l_bfgs_b, fmin_cg, check_grad
from scipy.stats import laplace, t
from tools import gaborf, mapp, logmeanexp, asshmarray, sqrtmi, sqrtm
from gsm import GSM
from copy import deepcopy

class ISA(Distribution):
	"""
	An implementation of overcomplete ISA using Gaussian scale mixtures.
	"""

	def __init__(self, num_visibles, num_hiddens=None, ssize=1, noise=False):
		"""
		@type  num_visibles: integer
		@param num_visibles: data dimensionality

		@type  num_hiddens: integer
		@param num_hiddens: number of hidden units

		@type  ssize: integer
		@param ssize: subspace dimensionality

		@type  noise: bool
		@param noise: add additional hidden units for noise
		"""

		if num_hiddens is None:
			num_hiddens = num_visibles

		self.dim = num_visibles
		self.num_visibles = num_visibles
		self.num_hiddens = num_hiddens

		# random linear feature 
		self.A = randn(self.num_visibles, self.num_hiddens) / 10.

		# subspace densities
		self.subspaces = [
			GSM(ssize) for _ in range(int(num_hiddens) / int(ssize))]

		if mod(num_hiddens, ssize) > 0:
			self.subspaces.append(GSM(mod(num_hiddens, ssize)))

		self._noise = False
		self.noise = noise



	def initialize(self, X=None, method='data'):
		"""
		Initializes parameter values with more sensible values.

		@type  X: array_like
		@param X: data points stored in columns

		@type  method: string
		@param method: type of initialization ('data', 'gabor' or 'random')
		"""

		if self.noise:
			L = self.A[:, :self.num_visibles]

		if method.lower() == 'data':
			# initialize features with data points
			if X is not None:
				if X.shape[1] < self.num_hiddens:
					raise ValueError('Number of data points to small.')

				elif X.shape[1] / 4 < self.num_hiddens:
					self.A = X[:, permutation(X.shape[1])[:self.num_hiddens]]

				else:
					# whitening matrix
					val, vec = eig(cov(X))
					M = dot(diag(1. / val), vec.T)

					# sort by norm in whitened space
					indices = argsort(sqrt(sum(square(dot(M, X)), 0)))

					# pick random subset
					indices = indices[permutation(min([X.shape[1] / 4, 10000]))[:self.num_hiddens]]

					self.A = X[:, indices]

				self.orthogonalize()

		elif method.lower() == 'gabor':
			# initialize features with Gabor filters
			if self.subspaces[0].dim > 1 and not mod(self.num_hiddens, 2):
				for i in range(self.num_hiddens / 2):
					G = gaborf(self.num_visibles)
					self.A[:, 2 * i] = real(G)
					self.A[:, 2 * i + 1] = imag(G)
			else:
				for i in range(len(self.subspaces)):
					self.A[:, i] = gaborf(self.num_visibles, complex=False)

		elif method.lower() == 'random':
			# initialize with Gaussian white noise
			self.A = randn(num_visibles, num_hiddens)

		elif method.lower() == 'laplace' or method.lower() == 'student':
			if method.lower() == 'laplace':
				# approximate multivariate Laplace with GSM
				samples = randn(self.subspaces[0].dim, 10000)
				samples = samples / sqrt(sum(square(samples), 0))
				samples = laplace.rvs(size=[1, 10000]) * samples
			else:
				samples = randn(self.subspaces[0].dim, 10000)
				samples = samples / sqrt(sum(square(samples), 0))
				samples = t.rvs(2., size=[1, 10000]) * samples

			if self.noise:
				# ignore first subspace
				gsm = GSM(self.subspaces[1].dim, self.subspaces[1].num_scales)
				gsm.train(samples, max_iter=200, tol=1e-8)

				for m in self.subspaces[1:]:
					m.scales = gsm.scales.copy()
			else:
				gsm = GSM(self.subspaces[0].dim, self.subspaces[0].num_scales)
				gsm.train(samples, max_iter=200, tol=1e-8)

				for m in self.subspaces:
					m.scales = gsm.scales.copy()

		else:
			raise ValueError('Unknown initialization method \'{0}\'.'.format(method))

		if self.noise:
			# don't initialize noise covariance
			self.A[:, :self.num_visibles] = L



	def train(self, X, method=('sgd', {}), sampling_method=('gibbs', {}), **kwargs):
		"""
		Optimizes model parameters with respect to the log-likelihoood. If the model is
		overcomplete, expectation maximization (EM) is used.

		The callback function takes two arguments: the model and the current iteration of the
		EM algorithm. It is called before the training starts and then after every iteration.

		@type  max_iter: integer
		@param max_iter: maximum number of iterations through the dataset

		@type  method: tuple
		@param method: optimization method used to optimize filters

		@type  adative: bool
		@param adaptive: automatically adjust step width when using SGD (default: True)

		@type  train_prior: bool
		@param train_prior: whether or not to optimize the marginal distributions (default: True)

		@type  train_subspaces: bool
		@param train_subspaces: automatically learn subspace sizes (default: False)

		@type  orthogonalize: bool
		@param orthogonalize: after each step, orthogonalize rows of feature matrix (default: False)

		@type  sampling_method: tuple
		@param sampling_method: method and parameters to generate hidden representations

		@type  init_sampling_steps: integer
		@param init_sampling_steps: number of steps used to initialize persistent samples (default: 0)

		@type  callback: function
		@param callback: called after every iteration
		"""

		max_iter = kwargs.get('max_iter', 100)
		adaptive = kwargs.get('adaptive', True) 
		train_prior = kwargs.get('train_prior', True)
		train_subspaces = kwargs.get('train_subspaces', False)
		orthogonalize = kwargs.get('orthogonalize', False)
		persistent = kwargs.get('persistent', True)
		init_sampling_steps = kwargs.get('init_sampling_steps', 0)
		callback = kwargs.get('callback', None)

		if Distribution.VERBOSITY > 0:
			if self.num_hiddens > self.num_visibles:
				print 0
			else:
				print 0, self.evaluate(X)

		if isinstance(method, str):
			method = (method, {})

		if isinstance(sampling_method, str):
			sampling_method = (sampling_method, {})

		if persistent and init_sampling_steps:
			# initialize samples
			sampling_method[1]['Z'] = self.sample_nullspace(X, 
				method=(sampling_method[0], 
					dict(sampling_method[1], num_steps=init_sampling_steps)))

		if adaptive and 'step_width' not in method[1]:
			method[1]['step_width'] = 0.001

		if callback:
			callback(self, 0)

		for i in range(max_iter):
			# complete data (E)
			Y = self.sample_posterior(X, method=sampling_method)

			if train_prior:
				# optimize parameters of the prior (M)
				self.train_prior(Y)

			if train_subspaces:
				# learn subspaces (M)
				Y = self.train_subspaces(Y)

			if not orthogonalize and train_prior or train_subspaces:
				# normalize variances of marginals
				self.normalize_prior()

			if persistent:
				# initializes samples in next iteration
				sampling_method[1]['Z'] = dot(self.nullspace_basis(), Y)

			# optimize linear features (M)
			if method[0].lower() == 'sgd':
				improved = self.train_sgd(Y, **method[1])

				if adaptive:
					# adjust learning rate
					method[1]['step_width'] *= 1.1 if improved else 0.5

			elif method[0].lower() == 'lbfgs':
				self.train_lbfgs(Y, **method[1])

			else:
				raise ValueError('Unknown training method \'{0}\'.'.format(method[0]))

			if orthogonalize:
				# normalize feature matrix
				self.orthogonalize()

			if callback:
				callback(self, i + 1)

			if Distribution.VERBOSITY > 0:
				if self.num_hiddens > self.num_visibles:
					print i + 1
				else:
					print i + 1, self.evaluate(X)



	def train_prior(self, Y, **kwargs):
		"""
		Optimize parameters of the marginal distribution over the hidden variables.
		The parameters are fit to maximize the average log-likelihood of the
		columns in `Y`.

		@type  Y: array_like
		@param Y: hidden states
		"""

		offset = [0]
		for model in self.subspaces:
			offset.append(offset[-1] + model.dim)

{'fraction_non_alphanumeric': 0.1152599507325295, 'fraction_numerical': 0.012511344483339817, 'mean_word_length': 3.843642072213501, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 23, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5242146', 'n_tokens_mistral': 6115, 'n_tokens_neox': 5012, 'n_words': 1474}