Files
Rbm/source/Rbm.cpp
T
2022-01-08 13:30:26 +00:00

433 lines
10 KiB
C++

/*
* 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"
#define RBM_TRAIN_FLAT 0
Rbm::Rbm(size_t numVisible, size_t numHidden, size_t numContext)
: m_params()
, m_whv(numVisible, numHidden)
, m_whc(numContext, numHidden)
, m_bhv(1, numHidden)
, m_bhc(1, numHidden)
, m_bv(1, numVisible)
, m_bc(1, numContext)
{
assert(numVisible > 0);
assert(numHidden > 0);
if (numContext)
{
assert(numContext == numHidden);
}
Noise_Init(&m_noise, 0x32727155);
}
Rbm::Rbm(const Rbm& orig)
: m_params(orig.m_params)
, m_whv(orig.m_whv)
, m_whc(orig.m_whc)
, m_bhv(orig.m_bhv)
, m_bhc(orig.m_bhc)
, m_bv(orig.m_bv)
, m_bc(orig.m_bc)
{
}
Rbm::~Rbm()
{
Noise_Free(&m_noise);
}
void Rbm::weightsInit(double stddev, double mu)
{
uniform(m_whv, stddev, mu);
uniform(m_whc, 0.01*stddev, mu);
uniform(m_bhv, stddev, mu);
uniform(m_bhc, 0.01*stddev, mu);
uniform(m_bv, stddev, mu);
uniform(m_bc, stddev, mu);
}
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_hv(arma::mat &hv_probs, arma::mat &v_probs)
{
for (int gibbs=0; gibbs < m_params.numGibbs; gibbs++)
{
// Create visible reconstruction (a fantasy...) given hid
if (m_params.gibbsDoSampleHidden)
{
v_probs = prob(h_to_v(sample(hv_probs)));
}
else
{
v_probs = prob(h_to_v(hv_probs));
}
// Create hidden representation given v
if (m_params.gibbsDoSampleVisible)
{
hv_probs = prob(v_to_h(sample(v_probs)));
}
else
{
hv_probs = prob(v_to_h(v_probs));
}
}
}
void Rbm::gibbs_hc(arma::mat &hc_probs, arma::mat &c_probs)
{
for (int gibbs=0; gibbs < m_params.numGibbs; gibbs++)
{
// Create visible reconstruction (a fantasy...) given hid
if (m_params.gibbsDoSampleHidden)
{
c_probs = prob(h_to_c(sample(hc_probs)));
}
else
{
c_probs = prob(h_to_c(hc_probs));
}
// Create hidden representation given v
if (m_params.gibbsDoSampleVisible)
{
hc_probs = prob(c_to_h(sample(c_probs)));
}
else
{
hc_probs = prob(c_to_h(c_probs));
}
}
}
void Rbm::weightUpdate_hv(arma::mat const &v_states, arma::mat &dw, arma::mat &dbhv, arma::mat &dbv)
{
arma::mat v_probs(v_states);
arma::mat h_states = v_to_h(v_states);
arma::mat h_probs = prob(h_states);
// Sample hidden
if (m_params.doRaoBlackwell)
{
h_states = h_probs;
}
else
{
h_states = sample(h_probs);
}
// Update weights (positive phase)
dw = v_states.t() * h_states;
dbv = sum(v_states, 0);
dbhv = sum(h_states, 0);
gibbs_hv(h_probs, v_probs);
// Update weights (negative phase)
dw -= v_probs.t() * h_probs;
dbv -= sum(v_probs, 0);
dbhv -= sum(h_probs, 0);
}
void Rbm::weightUpdate_hc(arma::mat const &c_states, arma::mat &dw, arma::mat &dbhc, arma::mat &dbc)
{
arma::mat c_probs(c_states);
arma::mat h_states = c_to_h(c_states);
arma::mat h_probs = prob(h_states);
// Sample hidden
if (m_params.doRaoBlackwell)
{
h_states = h_probs;
}
else
{
h_states = sample(h_probs);
}
// Update weights (positive phase)
dw = c_states.t() * h_states;
dbc = sum(c_states, 0);
dbhc = sum(h_states, 0);
gibbs_hc(h_probs, c_probs);
// Update weights (negative phase)
dw -= c_probs.t() * h_probs;
dbc -= sum(c_probs, 0);
dbhc -= sum(h_probs, 0);
}
void Rbm::train(const arma::mat& batch, IListener* pListener)
{
Status status;
double dProgress = 100.0/(batch.n_rows*m_params.numEpochs);
double progress = 0;
int lastProgress = -100;
int batchRowIndex = 0;
arma::mat grad_bias_v(arma::zeros(1, m_bv.n_cols));
arma::mat grad_bias_c(arma::zeros(1, m_bc.n_cols));
arma::mat grad_bias_hv(arma::zeros(1, m_bhv.n_cols));
arma::mat grad_bias_hc(arma::zeros(1, m_bhc.n_cols));
arma::mat grad_weight_hv(arma::zeros(m_whv.n_rows, m_whv.n_cols));
arma::mat grad_weight_hc(arma::zeros(m_whc.n_rows, m_whc.n_cols));
arma::mat momentum_whv = arma::zeros(m_whv.n_rows, m_whv.n_cols);
arma::mat momentum_whc = arma::zeros(m_whc.n_rows, m_whc.n_cols);
arma::mat momentum_bias_v(arma::zeros(1, m_bv.n_cols));
arma::mat momentum_bias_c(arma::zeros(1, m_bc.n_cols));
arma::mat momentum_bias_hv(arma::zeros(1, m_bhv.n_cols));
arma::mat momentum_bias_hc(arma::zeros(1, m_bhc.n_cols));
arma::mat penalty_weights = arma::zeros(m_whv.n_rows, m_whv.n_cols);
int trainingSizeRemain = batch.n_rows;
bool shouldAbort = false;
while (trainingSizeRemain && !shouldAbort)
{
int miniBatchSizeActual = std::min(m_params.miniBatchSize, trainingSizeRemain);
arma::mat miniBatch = batch.rows(batchRowIndex, batchRowIndex+miniBatchSizeActual-1);
trainingSizeRemain -= miniBatchSizeActual;
batchRowIndex += miniBatchSizeActual;
int scaler = std::min(m_params.miniBatchSize, (int)batch.n_rows);
double learning_rate = m_params.learningRate/scaler;
double weight_decay = m_params.weightDecay/scaler;
#if RBM_TRAIN_FLAT
arma::mat h_probs(miniBatchSizeActual, m_bhv.n_cols);
arma::mat v_probs(miniBatchSizeActual, m_bv.n_cols);
arma::mat hid_states(miniBatchSizeActual, m_bhv.n_cols);
#endif
arma::mat v_states(miniBatch);
arma::mat c_states(arma::zeros(miniBatchSizeActual, m_bc.n_cols));
// Create hidden layer base on training data
if (m_params.doSampleBatch)
{
// When the hidden units are being driven by data, always use stochastic binary states
v_states = sample(miniBatch);
}
for (int epoch=0; epoch < m_params.numEpochs; epoch++)
{
#if RBM_TRAIN_FLAT
// Sample hidden
h_probs = prob(v_to_h(v_states));
if (m_params.doRaoBlackwell)
{
hid_states = h_probs;
}
else
{
hid_states = sample(h_probs);
}
// Update weights (positive phase)
grad_weight = v_states.t() * hid_states;
grad_bias_v = sum(v_states, 0);
grad_bias_hv = sum(hid_states, 0);
for (int gibbs_hv=0; gibbs_hv < m_params.numGibbs; gibbs_hv++)
{
// Create visible reconstruction (a fantasy...) given hid
if (m_params.gibbsDoSampleHidden)
{
v_probs = prob(h_to_v(sample(h_probs)));
}
else
{
v_probs = prob(h_to_v(h_probs));
}
// Create hidden representation given v
if (m_params.gibbsDoSampleVisible)
{
hid_states = v_to_h(sample(v_probs));
}
else
{
hid_states = v_to_h(v_probs);
}
h_probs = prob(hid_states);
}
// Update weights (negative phase)
grad_weight -= v_probs.t() * h_probs;
grad_bias_v -= sum(v_probs, 0);
grad_bias_hv -= sum(h_probs, 0);
#else
if (m_bc.n_cols == 0)
{
weightUpdate_hv(v_states, grad_weight_hv, grad_bias_hv, grad_bias_v);
}
else
{
weightUpdate_hv(v_states, grad_weight_hv, grad_bias_hv, grad_bias_v);
weightUpdate_hc(c_states, grad_weight_hc, grad_bias_hc, grad_bias_c);
}
#endif
penalty_weights = weight_decay*arma::sign(m_whv);
status.L1 = accu(abs(m_whv));
status.L2 = accu(m_whv % m_whv);
momentum_bias_v = m_params.momentum*momentum_bias_v + grad_bias_v;
momentum_bias_hv = m_params.momentum*momentum_bias_hv + grad_bias_hv;
momentum_bias_c = m_params.momentum*momentum_bias_c + grad_bias_c;
momentum_bias_hc = m_params.momentum*momentum_bias_hc + grad_bias_hc;
momentum_whv = m_params.momentum*momentum_whv + grad_weight_hv - status.L2*penalty_weights;
momentum_whc = m_params.momentum*momentum_whc + grad_weight_hc;
m_bv += learning_rate*momentum_bias_v;
m_bc += learning_rate*momentum_bias_c;
m_bhv += learning_rate*momentum_bias_hv;
m_bhc += learning_rate*momentum_bias_hc;
m_whv += learning_rate*momentum_whv;
m_whc += learning_rate*momentum_whc;
progress += dProgress*miniBatchSizeActual;
status.progress = (int)(progress + 0.5);
if (status.progress != lastProgress)
{
lastProgress = status.progress;
// Calculate error
arma::mat diffErr = miniBatch - prob(h_to_v(prob(v_to_h(v_states))));
arma::mat diffErr_squared = diffErr % diffErr;
status.err = accu(diffErr_squared)/diffErr_squared.n_elem;
if (pListener)
{
if(!pListener->onProgress(this, status))
{
shouldAbort = true;
break;
}
}
}
} // Number of epochs
} // number of mini batches
arma::mat diffErr = batch - prob(h_to_v(prob(v_to_h(batch))));
arma::mat diffErr_squared = diffErr % diffErr;
status.err_total = accu(diffErr_squared)/diffErr_squared.n_elem;
if (pListener)
{
pListener->onProgress(this, status);
}
}
arma::mat Rbm::prob(const arma::mat &src)
{
return 1 / (1 + (arma::exp(-src)));
}
arma::mat Rbm::sample(const arma::mat &src)
{
arma::mat dst = src;
uniform(dst);
for (size_t i=0; i < src.n_rows; i++)
{
for (size_t j=0; j < src.n_cols; j++)
{
dst(i, j) = src(i, j) >= dst(i, j);
}
}
return dst;
}
arma::mat Rbm::v_to_h(const arma::mat &visible) const
{
return visible * m_whv + arma::repmat(m_bhv, 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);
}
arma::mat Rbm::c_to_h(const arma::mat &context) const
{
return context * m_whc + arma::repmat(m_bhc, context.n_rows, 1);
}
arma::mat Rbm::h_to_c(const arma::mat &hidden) const
{
return hidden * m_whc.t() + arma::repmat(m_bc, hidden.n_rows, 1);
}
arma::mat Rbm::normalize(const arma::mat& src)
{
double mean = arma::accu(src)/src.n_elem;
arma::mat x = src - mean;
arma::mat x2 = x % x;
double stddev = sqrt(arma::accu(x2)/x2.n_elem);
std::cout << "mean" << " : " << std::endl << mean << std::endl;
std::cout << "stddev" << ": " << std::endl << stddev << std::endl;
return x/stddev;
}
void Rbm::uniform(arma::mat& srcDst, double stdDev, double mu)
{
#if 1
for (size_t i=0; i < srcDst.n_rows; i++)
{
for (size_t j=0; j < srcDst.n_cols; j++)
{
srcDst(i, j) = stdDev*(Noise_Uniform(&m_noise) + mu - 0.5);
}
}
#else
srcDst = stdDev*(arma::randu(srcDst.n_rows, srcDst.n_cols) + mu - 0.5);
#endif
}
const arma::mat& Rbm::w() const
{
return m_whv;
}
const arma::mat& Rbm::bv() const
{
return m_bv;
}
const arma::mat& Rbm::bh() const
{
return m_bhv;
}