Rbm/source/Rbm.cpp

/*
 * To change this license header, choose License Headers in Project Properties.
 * To change this template file, choose Tools | Templates
 * and open the template in the editor.
 */

/* 
 * File:   Rbm.cpp
 * Author: jens
 * 
 * Created on 21. Oktober 2019, 21:28
 */

#include <cassert>
#include "Rbm.hpp"
#include "matutils.hpp"

using namespace Matutils;

#define USE_CD_HINTON 0

Rbm::Rbm(size_t numVisible, size_t numHidden)
: m_params()
, m_whv(numVisible, numHidden)
, m_bh(1, numHidden)
, m_bv(1, numVisible)
{
	assert(numVisible > 0);
	assert(numHidden > 0);
}

Rbm::Rbm(const Rbm& orig)
: m_params(orig.m_params)
, m_whv(orig.m_whv)
, m_bh(orig.m_bh)
, m_bv(orig.m_bv)
{
}

Rbm::~Rbm()
{
}

size_t Rbm::numHidden() const
{
	return m_bh.size();
}

size_t Rbm::numVisible() const
{
	return m_bv.size();
}

Rbm::Params& Rbm::params()
{
	return m_params;
}

arma::mat Rbm::rms_error(arma::mat diffErr)
{
	arma::mat diffErr_squared = diffErr % diffErr;
	return arma::sum(diffErr_squared, 1) * 1.0 / diffErr_squared.n_cols;
}

double Rbm::rms_error_accu(arma::mat diffErr)
{
	arma::mat diffErr_squared = diffErr % diffErr;
	return arma::accu(diffErr_squared) / diffErr_squared.n_elem;
}

arma::mat Rbm::toHiddenProbs(const arma::mat& visible) const
{
	arma::mat state = v_to_h(visible);
	if (m_params.doGaussianHidden)
	{
		return state;
	}
	return prob(state);
}

arma::mat Rbm::toVisibleProbs(const arma::mat& hidden) const
{
	arma::mat state = h_to_v(hidden);
	if (m_params.doGaussianVisible)
	{
		return state;
	}
	return prob(state);
}

arma::mat Rbm::sampleVisible(const arma::mat& visible) const
{
	return m_params.doGaussianVisible ? sample_gaussian(visible) : sample(visible);
}

arma::mat Rbm::sampleHidden(const arma::mat& hidden) const
{
	return m_params.doGaussianHidden ? sample_gaussian(hidden) : sample(hidden);
}

void Rbm::weightsAssign(const arma::mat& w, const arma::mat& bh, const arma::mat& bv)
{
	m_whv.submat(0, 0, w.n_rows - 1, w.n_cols - 1) = w;
	m_bh.submat(0, 0, bh.n_rows - 1, bh.n_cols - 1) = bh;
	m_bv.submat(0, 0, bv.n_rows - 1, bv.n_cols - 1) = bv;
}

void Rbm::weightsInit(double stddev, double mu)
{
	m_whv = uniform(m_whv, stddev, mu);
	m_bh = uniform(m_bh, 0, 0);
	m_bv = uniform(m_bv, 0, 0);
}

void Rbm::fromJson(Json::Value rbm)
{
	std::cout << "Importing Rbm" << std::endl;
	m_params.fromJson(rbm["params"]);
}

Json::Value Rbm::toJson() const
{
	std::cout << "Exporting Rbm" << std::endl;
	Json::Value rbm;
	rbm["params"] = m_params.toJson();
	return rbm;
}

void Rbm::gibbs_vh(arma::mat &v_probs, arma::mat &h_probs) const
{
	for (int i=0; i < m_params.numGibbs; i++)
	{
		// Create hidden representation given v
		h_probs = toHiddenProbs(v_probs);

		// Create visible reconstruction (a fantasy...) given hid
		v_probs = toVisibleProbs(h_probs);
	}
}

void Rbm::gibbs_hv(arma::mat &h_probs, arma::mat &v_probs) const
{
	for (int i=0; i < m_params.numGibbs; i++)
	{
		// Create visible reconstruction (a fantasy...) given hid
		v_probs = toVisibleProbs(h_probs);

		// Create hidden representation given v
		h_probs = toHiddenProbs(v_probs);
	}
}

void Rbm::cd(arma::mat const &v_data, arma::mat &dw, arma::mat &dbh, arma::mat &dbv)
{
#if USE_CD_HINTON
	{
		cd_hinton(v_data, dw, dbh, dbv);
	}
#else
	{
		cd_jens(v_data, dw, dbh, dbv);
	}
#endif
}

void Rbm::cd_hinton_hid_linear(arma::mat const &v_data, arma::mat &dw, arma::mat &dbh, arma::mat &dbv)
{
	// Start positive phase	
	arma::mat poshidprobs = v_to_h(v_data);
	arma::mat posprods = v_data.t() * poshidprobs;
	arma::mat poshidact = arma::sum(poshidprobs);
	arma::mat posvisact = arma::sum(v_data);
	
	// End of positive phase	
	arma::mat poshidstates = poshidprobs + arma::randn(arma::size(poshidprobs));

	// Start negative phase	
	arma::mat negdata = prob(h_to_v(poshidstates));
	arma::mat neghidprobs = v_to_h(negdata); 
	arma::mat negprods = negdata.t() * neghidprobs; 
	arma::mat neghidact = arma::sum(neghidprobs);
	arma::mat negvisact = arma::sum(negdata);

	// Update weight deltas
	dw = posprods - negprods;
	dbv = posvisact - negvisact;
	dbh = poshidact - neghidact;
}

void Rbm::cd_hinton(arma::mat const &v_data, arma::mat &dw, arma::mat &dbh, arma::mat &dbv)
{
	// Start positive phase	
	arma::mat poshidprobs = toHiddenProbs(v_data);
	arma::mat posprods;
	arma::mat poshidact = arma::sum(poshidprobs);
	arma::mat posvisact = arma::sum(v_data);
	
	// End of positive phase
	arma::mat poshidstates;

	if (m_params.doGaussianHidden)
	{
		if (m_params.gibbsDoSampleHidden)
		{
			poshidstates = poshidprobs + arma::randn(arma::size(poshidprobs));
		}
		else
		{
			poshidstates = poshidprobs;
		}
	}
	else
	{
		poshidstates = sample(poshidprobs);
	}

	if (m_params.doRaoBlackwell)
	{
		posprods = v_data.t() * poshidprobs;
	}
	else
	{
		posprods = v_data.t() * sampleHidden(poshidprobs);
	}

	arma::mat negprods;
	arma::mat neghidact;
	arma::mat negvisact;
	for (int i=0; i < m_params.numGibbs; i++)
	{
		// Start negative phase
		arma::mat negdata;
		if (m_params.gibbsDoSampleHidden)
		{
			negdata = toVisibleProbs(sampleHidden(poshidstates));
		}
		else
		{
			negdata = toVisibleProbs(poshidstates);
		}
		arma::mat neghidprobs;
		if (m_params.gibbsDoSampleVisible)
		{
			neghidprobs = toHiddenProbs(sampleVisible(negdata));
		}
		else
		{
			neghidprobs = toHiddenProbs(negdata);
		}
		negprods = negdata.t() * neghidprobs; 
		neghidact = arma::sum(neghidprobs);
		negvisact = arma::sum(negdata);

		poshidprobs = neghidprobs;
	}
	// Update weight deltas
	dw = posprods - negprods;
	dbv = posvisact - negvisact;
	dbh = poshidact - neghidact;
}

void Rbm::cd_jens(arma::mat const &v_states, arma::mat &dw, arma::mat &dbh, arma::mat &dbv)
{
	arma::mat v_probs(v_states);
	arma::mat h_probs = toHiddenProbs(v_states);
	arma::mat h_states;

	// Sample hidden. Rao-Blackwellization uses the mean for the positive-phase
	// gradient (lower variance) regardless of unit type; previously a Gaussian
	// hidden unit always got the noisy sample here even with doRaoBlackwell
	// set, ignoring the flag and adding avoidable gradient variance.
	if (m_params.doRaoBlackwell)
	{
		h_states = h_probs;
	}
	else
	{
		h_states = sampleHidden(h_probs);
	}
	
	// Update weights (positive phase)
	dw = v_states.t() * h_states;
	dbv = sum(v_states, 0);
	dbh = sum(h_states, 0);

	// Gibbs sampling with training params
	for (int i=0; i < m_params.numGibbs; i++)
	{
		// Create visible reconstruction (a fantasy...) given hid
		if (m_params.gibbsDoSampleHidden)
		{
			v_probs = toVisibleProbs(sampleHidden(h_probs));
		}
		else
		{
			v_probs = toVisibleProbs(h_probs);
		}

		// Create hidden representation given v
		if (m_params.gibbsDoSampleVisible)
		{
			h_probs = toHiddenProbs(sampleVisible(v_probs));
		}
		else
		{
			h_probs = toHiddenProbs(v_probs);
		}
	}

	// Update weights (negative phase)
	dw -= v_probs.t() * h_probs;
	dbv -= sum(v_probs, 0);
	dbh -= sum(h_probs, 0);
}

void Rbm::train(arma::mat const &batch, IListener* pListener)
{
	Status status;
	double dProgress = 100.0/(batch.n_rows*m_params.numEpochs);
	double progress = 0;
	int lastProgress = -100;

	arma::mat dbv(arma::zeros(1, m_bv.n_cols));
	arma::mat dbh(arma::zeros(1, m_bh.n_cols));
	arma::mat dwhv(arma::zeros(m_whv.n_rows, m_whv.n_cols));
	arma::mat inc_whv = arma::zeros(m_whv.n_rows, m_whv.n_cols);
	arma::mat inc_bv(arma::zeros(1, m_bv.n_cols));
	arma::mat inc_bh(arma::zeros(1, m_bh.n_cols));

	bool shouldAbort = false;
	// One epoch = one full pass over the (freshly reshuffled) data, mini-batch
	// by mini-batch. Epochs must be the outer loop: nesting mini-batches on
	// the outside would run all numEpochs steps on one fixed chunk before ever
	// moving to the next, biasing training toward whatever data comes last.
	for (int epoch=0; epoch < m_params.numEpochs && !shouldAbort; epoch++)
	{
		arma::mat shuffledBatch = arma::shuffle(batch);
		int trainingSizeRemain = shuffledBatch.n_rows;
		int batchRowIndex = 0;

		while (trainingSizeRemain && !shouldAbort)
		{
			int miniBatchSizeActual = std::min(m_params.miniBatchSize, trainingSizeRemain);
			arma::mat miniBatch = shuffledBatch.rows(batchRowIndex, batchRowIndex+miniBatchSizeActual-1);
			trainingSizeRemain -= miniBatchSizeActual;
			batchRowIndex += miniBatchSizeActual;
			int numcases = std::min(m_params.miniBatchSize, (int)batch.n_rows);

			arma::mat v_states(miniBatch);

			// Create hidden layer base on training data
			if (m_params.doSampleBatch)
			{
				// Sample the visible batch (respects doGaussianVisible)
				v_states = sampleVisible(miniBatch);
			}

			// Contrastive divergence learning: calculate gradients
			cd(v_states, dwhv, dbh, dbv);

			// Adjust weight and biases
			inc_bv = m_params.momentum*inc_bv + m_params.learningRate/numcases*dbv;
			inc_bh = m_params.momentum*inc_bh + m_params.learningRate/numcases*dbh;
			inc_whv = m_params.momentum*inc_whv + m_params.learningRate*dwhv/numcases;

			m_bv += inc_bv;
			m_bh += inc_bh;
			m_whv += inc_whv;

			// L2 and L1 applied directly to the weights, outside the momentum
			// recursion above -- folding either into inc_whv would let momentum
			// geometrically amplify its effective strength to roughly
			// lambda/(1-momentum), well past the configured value. Neither
			// touches the biases, matching pyRBM's state_adjust().
			if (m_params.l2Lambda != 0.0)
			{
				m_whv -= m_params.learningRate*m_params.l2Lambda*m_whv;
			}
			if (m_params.l1Lambda != 0.0)
			{
				// Subgradient of lambda*||W||_1 is lambda*sign(W) (sign(0)=0).
				m_whv -= m_params.learningRate*m_params.l1Lambda*arma::sign(m_whv);
			}

			progress += dProgress*miniBatchSizeActual;
			status.progress = (int)(progress + 0.5);

			// Update status
			if (status.progress != lastProgress)
			{
				lastProgress = status.progress;

				// Calculate error
				status.err = rms_error_accu(miniBatch - toVisibleProbs(toHiddenProbs(v_states)));
				status.L1 = arma::accu(arma::abs(m_whv));
				if (pListener)
				{
					if(!pListener->onProgress(this, status))
					{
						shouldAbort = true;
						break;
					}
				}
			}

		} // mini batches

	} // epochs

	// Update final status
	status.err_total = rms_error_accu(batch - toVisibleProbs(toHiddenProbs(batch)));
	status.L1 = arma::accu(arma::abs(m_whv));

	if (pListener)
	{
		pListener->onProgress(this, status);
	}
}

arma::mat Rbm::v_to_h(const arma::mat &visible) const
{
	return visible * m_whv + arma::repmat(m_bh, visible.n_rows, 1);
}

arma::mat Rbm::h_to_v(const arma::mat &hidden) const
{
	return hidden * m_whv.t() + arma::repmat(m_bv, hidden.n_rows, 1);
}

const arma::mat& Rbm::whv() const
{
	return m_whv;
}

const arma::mat& Rbm::bv() const
{
	return m_bv;
}

const arma::mat& Rbm::bh() const
{
	return m_bh;
}