ebm.cpp

Gaussian Mixture Algorithm

{
	idx_dtanh(in->x, in->ddx);
	idx_mul(in->ddx, in->ddx, in->ddx);
	idx_mul(in->ddx, out->ddx, in->ddx);
}

void tanh_module::forget(forget_param_linear& fp)
{
}

void tanh_module::normalize()
{
}

////////////////////////////////////////////////////////////////

addc_module::addc_module(parameter *p, intg size)
{
	bias = new state_idx(p,size);
}

addc_module::~addc_module()
{
	delete bias;
}

void addc_module::fprop(state_idx* in, state_idx* out)
{
	out->resize(bias->x.dim(0));
	idx_add(in->x,bias->x,out->x);
}

void addc_module::bprop(state_idx* in, state_idx* out)
{
	idx_copy(out->dx,in->dx);
	idx_copy(out->dx,bias->dx);
}

void addc_module::bbprop(state_idx* in, state_idx* out)
{
	idx_copy(out->ddx,in->ddx);
	idx_copy(out->ddx,bias->ddx);
}

void addc_module::forget(forget_param_linear& fp)
{
	idx_clear(bias->x);
}

void addc_module::normalize()
{
}
////////////////////////////////////////////////////////////////

nn_layer_full::nn_layer_full(parameter *p, intg ninputs, intg noutputs)
{
	linear = new linear_module(p,ninputs,noutputs);
	bias = new state_idx(p,noutputs);
	sum = new state_idx(noutputs);
	sigmoid = new tanh_module();
}

nn_layer_full::~nn_layer_full()
{
	delete sigmoid;
	delete sum;
	delete bias;
	delete linear;
}

void nn_layer_full::fprop(state_idx *in, state_idx *out)
{
	out->resize(bias->x.dim(0));
	linear->fprop(in, sum);
	idx_add(sum->x, bias->x, sum->x);
	sigmoid->fprop(sum, out);
}

void nn_layer_full::bprop(state_idx *in, state_idx *out)
{
	sigmoid->bprop(sum, out);
	idx_copy(sum->dx, bias->dx);
	linear->bprop(in, sum);
}

void nn_layer_full::bbprop(state_idx *in, state_idx *out)
{
	sigmoid->bbprop(sum, out);
	idx_copy(sum->ddx, bias->ddx);
	linear->bbprop(in, sum);
}

void nn_layer_full::forget(forget_param_linear &fp)
{
	linear->forget(fp);
	idx_clear(bias->x);
}

////////////////////////////////////////////////////////////////

f_layer::f_layer(parameter *p, intg tin, intg tout, intg si, intg sj,
		module_1_1<state_idx,state_idx> *sq)
{
	weight = new state_idx(p, tout, tin);
	bias = new state_idx(p, tout);
	sum = new state_idx(tout, si, sj);
	squash = sq;
}

f_layer::~f_layer()
{
	delete weight;
	delete bias;
	delete sum;
}

void f_layer::forget(forget_param_linear &fp)
{
	idx_clear(bias->x);
	double z = fp.value / pow(weight->x.dim(1), fp.exponent);
	if(!drand_ini) printf("You have not initialized random sequence. Please call init_drand() before using this function !\n");
	idx_aloop1(w,weight->x,double)
	{	*w = drand(z);}
}

void f_layer::fprop(state_idx *in, state_idx *out)
{
	intg inx_d1 = in->x.dim(1);
	intg inx_d2 = in->x.dim(2);
	intg ws = weight->x.dim(0);
	// resize sum and output
	sum->resize(ws, inx_d1, inx_d2);
	out->resize(ws, inx_d1, inx_d2);
	// main matrix multiplication
	{
		int tr[] = { 2, 1, 0 };
		Idx<double> inx(in->x.transpose(tr));
		Idx<double> outx(sum->x.transpose(tr));
		// loop over spatial dimensions
		idx_bloop2(linx,inx,double, loutx,outx,double)
		{
			idx_bloop2(llinx,linx,double, lloutx,loutx,double)
			{
				// fprintf(stdout,"f_layer::fprop\n");
				// weight->x.pretty(stdout);
				// weight->x.fdump(stdout);
				// llinx.pretty(stdout);
				// llinx.fdump(stdout);
				// lloutx.pretty(stdout);
				// multiply weight matrix by input
				idx_m2dotm1(weight->x, llinx, lloutx);
				// lloutx.pretty(stdout);
				// lloutx.fdump(stdout);
			}
		}
	}
	{ // add bias
		// fprintf(stdout,"f_layer::fprop adding bias\n");
		// bias->x.pretty(stdout);
		// bias->x.fdump(stdout);
		idx_bloop2(sumx,sum->x,double, biasx,bias->x,double)
		{
			idx_addc(sumx, biasx.get(), sumx);
		}
		// sum->x.pretty(stdout);
		// sum->x.fdump(stdout);

	}
	// call squashing function
	squash->fprop(sum, out);
}

void f_layer::bprop(state_idx *in, state_idx *out)
{
	// backprop through squasher
	squash->bprop(sum, out);
	// compute gradient of bias
	{
		idx_bloop2(lha,sum->dx,double, lb,bias->dx,double)
		{	*(lb.ptr()) += idx_sum(lha);}
	}
	// backprop through weight matrix
	int tr[] = { 2, 1, 0 };
	Idx<double> inx(in->x.transpose(tr));
	Idx<double> indx(in->dx.transpose(tr));
	Idx<double> outdx(sum->dx.transpose(tr));
	Idx<double> tkerx(weight->x.transpose(0, 1));
	{ idx_bloop3(linx,inx,double, lindx,indx,double, loutdx,outdx,double) {
		{ idx_bloop3(llinx,linx,double, llindx,lindx,double, lloutdx,loutdx,double) {
			idx_m1extm1acc(lloutdx,llinx,weight->dx);
			idx_m2dotm1(tkerx,lloutdx,llindx);
			// direct call because idxm2dotm1 didn't use to work on tranposed matrices
			// cblas_dgemv(CblasRowMajor, CblasTrans, weight->x.dim(0), weight->x.dim(1),
			//             1.0, weight->x.idx_ptr(), weight->x.mod(0),
			//             lloutdx.idx_ptr(), lloutdx.mod(0),
			//             0.0, llindx.idx_ptr(), llindx.mod(0));
		}}
	}}
}

void f_layer::bbprop(state_idx *in, state_idx *out)
{
	// backprop through squasher
	squash->bbprop(sum, out);
	// compute gradient of bias
	{	idx_bloop2(lha, sum->ddx, double, lb, bias->ddx, double) {
		idx_sumacc(lha, lb); }}
	// backprop through weight matrix
	int tr[] = { 2, 1, 0 };
	Idx<double> inx(in->x.transpose(tr));
	Idx<double> indx(in->ddx.transpose(tr));
	Idx<double> outdx(sum->ddx.transpose(tr));
	Idx<double> tkerx(weight->x.transpose(1, 0));
	idx_bloop3(linx,inx,double, lindx,indx,double, loutdx,outdx,double)
	{
		idx_bloop3(llinx,linx,double, llindx,lindx,double, lloutdx,loutdx,double)
		{
			idx_m1squextm1acc(lloutdx,llinx,weight->ddx);
			idx_m2squdotm1(tkerx,lloutdx,llindx);
		}
	}
}

////////////////////////////////////////////////////////////////

c_layer::c_layer(parameter *p, intg ki, intg kj, intg ri, intg rj, Idx<intg> *tbl,
		intg thick, intg si, intg sj, module_1_1<state_idx,state_idx> *sqsh)
{
	thickness = thick;
	stridei = ri;
	stridej = rj;
	kernel = new state_idx(p, tbl->dim(0), ki, kj);
	table = tbl;
	bias = new state_idx(p, thick);
	sum = new state_idx(thick, si, sj);
	squash = sqsh;
}

c_layer::~c_layer()
{
	delete kernel;
	delete bias;
	delete sum;
}

void c_layer::set_stride(intg ri, intg rj)
{
	stridei = ri;
	stridej = rj;
}

void c_layer::forget(forget_param_linear &fp)
{
	idx_clear(bias->x);
	Idx<double> kx(kernel->x);
	intg vsize = kx.dim(1);
	intg hsize = kx.dim(2);
	Idx<intg> ts(table->select(1, 1));
	Idx<int> fanin(1 + idx_max(ts));
	if(!drand_ini) printf("You have not initialized random sequence. Please call init_drand() before using this function !\n");
  idx_clear(fanin);
	{ idx_bloop1(tab, *table, intg)	{
			fanin.set(1 + fanin.get(tab.get(1)), tab.get(1)); }}
	{ idx_bloop2(tab, *table, intg, x, kx, double) {
		double s = fp.value / pow((vsize * hsize * fanin.get(tab.get(1))), fp.exponent);
		{	idx_bloop1(lx, x, double)	{
			{	idx_bloop1(llx, lx, double) {
				double n = drand(-s, s);
				llx.set(n);
			}}
		}}
	}}
}

void c_layer::fprop(state_idx *in, state_idx *out)
{
	intg ki = kernel->x.dim(1);
	intg kj = kernel->x.dim(2);
	intg sini = in->x.dim(1);
	intg sinj = in->x.dim(2);

	if (((sini - (ki - stridei)) % stridei != 0) || ((sinj - (kj - stridej))
			% stridej != 0))ylerror("inconsistent input size, kernel size, and subsampling ratio.");
	if ((stridei != 1) || (stridej != 1))ylerror("stride > 1 not implemented yet.");
	Idx<double> uuin(in->x.unfold(1, ki, stridei));
	uuin = uuin.unfold(2, kj, stridej);
	Idx<double> lki(kernel->x.dim(1), kernel->x.dim(2));
	// resize output if necessary
	sum->resize(thickness, uuin.dim(1), uuin.dim(2));
	out->resize(thickness, uuin.dim(1), uuin.dim(2));
	idx_clear(sum->x);
	// generic convolution
	{	idx_bloop2(lk, kernel->x, double, lt, *table, intg)	{
		Idx<double> suin(uuin.select(0, lt->get(0)));
		Idx<double> sout((sum->x).select(0, lt->get(1)));
		idx_m4dotm2acc(suin, lk, sout);
	}}
	// add bias
	{	idx_bloop3(sumx, sum->x, double, biasx, bias->x, double,
			outx, out->x, double)	{
		idx_addc(sumx, biasx.get(), sumx);
	}}
	// call squashing function
	squash->fprop(sum, out);
}

void c_layer::bprop(state_idx *in, state_idx *out)
{
	// backprop gradient through squasher
	squash->bprop(sum, out);
	// compute gradient of bias
	{	idx_bloop2(lha, sum->dx, double, lb, bias->dx, double) {
			idx_sumacc(lha, lb); }}
	// backprop through convolution
	idx_clear(in->dx);
	Idx<double> uuin(in->dx.unfold(1, (kernel->dx).dim(1), stridei));
	uuin = uuin.unfold(2, (kernel->dx).dim(2), stridej);
	Idx<double> uuinf(in->x.unfold(1, (kernel->dx).dim(1), stridei));
	uuinf = uuinf.unfold(2, (kernel->dx).dim(2), stridej);
	int transp[5] = { 0, 3, 4, 1, 2 };
	Idx<double> borp(uuinf.transpose(transp));
	{ idx_bloop3 (lk, kernel->dx, double, lkf, kernel->x, double, lt, *table, intg) {
			intg islice = lt.get(0);
			Idx<double> suin(uuin.select(0, islice));
			Idx<double> sborp(borp.select(0, islice));
			Idx<double> sout((sum->dx).select(0, lt.get(1)));
			// backward convolution
			idx_m2extm2acc(sout, lkf, suin);
			// compute gradient for kernel
			idx_m4dotm2acc(sborp, sout, lk);
	}}
}

void c_layer::bbprop(state_idx *in, state_idx *out)
{
	// backprop gradient through squasher
	squash->bbprop(sum, out);
	// compute gradient of bias
	{ idx_bloop2(lha, sum->ddx, double, lb, bias->ddx, double) {
			idx_sumacc(lha, lb); }}
	// backprop through convolution
	idx_clear(in->ddx);
	Idx<double> uuin(in->ddx.unfold(1, (kernel->ddx).dim(1), stridei));
	uuin = uuin.unfold(2, (kernel->ddx).dim(2), stridej);
	Idx<double> uuinf(in->x.unfold(1, (kernel->ddx).dim(1), stridei));
	uuinf = uuinf.unfold(2, (kernel->ddx).dim(2), stridej);
	int transp[5] = { 0, 3, 4, 1, 2 };
	Idx<double> borp(uuinf.transpose(transp));
	{	idx_bloop3 (lk, kernel->ddx, double, lkf, kernel->x, double, lt, *table, intg) {
			intg islice = lt.get(0);
			Idx<double> suin(uuin.select(0, islice));
			Idx<double> sborp(borp.select(0, islice));
			Idx<double> sout((sum->ddx).select(0, lt.get(1)));
			// backward convolution
			idx_m2squextm2acc(sout, lkf, suin);
			// compute gradient for kernel
			idx_m4squdotm2acc(sborp, sout, lk);
	}}
}

////////////////////////////////////////////////////////////////

#ifdef USE_IPP

// TODO: Copy IPP in project
// TODO: ipp 64 for doubles?

void c_layer_ipp::fprop (state_idx *in, state_idx *out) {
  intg ki = kernel->x.dim(1);
  intg kj = kernel->x.dim(2);
  intg sini = in->x.dim(1);
  intg sinj = in->x.dim(2);

  if (((sini - (ki - stridei) % stridei) != 0) ||
      ((sinj - (kj - stridej) % stridej) != 0))
    ylerror("inconsistent input size, kernel size, and subsampling ratio.");
  if ((stridei != 1) || (stridej != 1))
    ylerror("stride > 1 not implemented yet.");
  Idx<double> uuin = in->x.unfold(1, ki, stridei);
  uuin = uuin.spec.unfold_inplace(2, kj, stridej);
  Idx<double> lki = Idx<double>(kernel->x.dim(1), kernel->x.dim(2));
  // resize output if necessary
  sum->resize(thickness, uuin.dim(1), uuin.dim(2));
  out->resize(thickness, uuin.dim(1), uuin.dim(2));
  idx_clear(sum->x);
  // generic convolution
  Idx<double> tout = Idx<double>(sum->x.dim(1), sum->x.dim(2));
  { idx_bloop2(lk, kernel->x, double, lt, *table, intg) {
      rev_idx2_tr(*lk, lki);
//      ipp_convolution_float(in->x.select(0, lt.get(0)), lki, tout);
//      ipp_add_float(tout, sum->x.select(0, lt.get(1)));
    }
  }
  // add bias
  { idx_bloop3(sumx, sum->x, double, biasx, bias->x, double,
	       outx, out->x, double) {
//      ipp_addc_nip_float(sumx, biasx.get(), outx);
    }
  }
  // call squashing function
  squash->fprop(sum, out);
}

void c_layer_ipp::bprop (state_idx *in, state_idx *out) {
  // backprop gradient through squasher
  squash->bprop(sum, out);
  // compute gradient of bias
  { idx_bloop2(lha, sum->dx, double, lb, bias->dx, double) {
      idx_sumacc(lha, lb);
    }
  }
  // backprop through convolution
  idx_clear(in->dx);
  /*
  (let* ((ki (idx-dim :kernel:dx 1))	 (kj (idx-dim :kernel:dx 2))
	 (ini (idx-dim :in:dx 1))	 (inj (idx-dim :in:dx 2))
	 (outi (idx-dim :out:dx 1))	 (outj (idx-dim :out:dx 2))
	 (sumi (idx-dim :sum:dx 1))	 (sumj (idx-dim :sum:dx 2))
	 (souti (gbtype-matrix sumi sumj))
	 (tout (gbtype-matrix ki kj)))
	 (idx-bloop ((lk :kernel:dx) (lkf :kernel:x) (lt table))
	 (let* ((islice (lt 0))
	 (sout  (select :sum:dx 0 (lt 1)))
	 )
	 ;; backward convolution
	 (ipp-convolution-full-float sout lkf (select :in:dx 0 islice))
	 ;; compute gradient for kernel
	 (rev-idx2-tr-float sout souti)
	 (ipp-convolution-float (select :in:x 0 islice) souti tout)
	 (ipp-add-float tout lk)
	 )))
	 */
}

#endif

////////////////////////////////////////////////////////////////

Idx<intg> full_table(intg a, intg b) {
  Idx<intg> m(a * b, 2);
  intg p = 0;
  for (intg j = 0; j < b; ++j) {
  	for (intg i = 0; i < a; ++i) {
  		m.set(i, p, 0);
  		m.set(j, p, 1);
  		p++;
  	}
  }
  return m;
}

////////////////////////////////////////////////////////////////////////

s_layer::s_layer(parameter *p, intg ki, intg kj, intg thick, intg si, intg sj,
		module_1_1<state_idx, state_idx> *sqsh) :
	stridei(ki), stridej(kj), squash(sqsh)
{
	coeff = new state_idx(p, thick);
	bias = new state_idx(p, thick);
	sub = new state_idx(thick, si, sj);
	sum = new state_idx(thick, si, sj);
}

s_layer::~s_layer()
{
	delete coeff;
	delete bias;
	delete sub;
	delete sum;
}

void s_layer::fprop(state_idx *in, state_idx *out)
{
	intg sin_t = in->x.dim(0);
	intg sin_i = in->x.dim(1);
	intg sin_j = in->x.dim(2);
	intg si = sin_i / stridei;
	intg sj = sin_j / stridej;

	if( (sin_i % stridei) != 0 ||
			(sin_j % stridej) != 0)
		ylerror("inconsistent input size and subsampleing ratio");
	sub->resize(sin_t, si, sj);
	sum->resize(sin_t, si, sj);
	out->resize(sin_t, si, sj);
	// 1. subsampling ( coeff * average )
	idx_clear(sub->x);
	{ idx_bloop4(lix, in->x, double, lsx, sub->x, double,
			lcx, coeff->x, double, ltx, sum->x, double) {
		Idx<double> uuin(lix->unfold(1, stridej, stridej));
		uuin = uuin.unfold(0, stridei, stridei);
		{ idx_eloop1(z1, uuin, double) {
			{ idx_eloop1(z2, z1, double) {
				idx_add(z2, lsx, lsx);
			 }
			}
		 }
		}
		idx_dotc(lsx, lcx.get(), ltx);
	 }
	}
	// 2. add bias
	{ idx_bloop3(sumx, sum->x, double, biasx, bias->x, double,
			outx, out->x, double) {
		idx_addc(sumx, biasx.get(), sumx);
	 }
	}
	// 3. call squashing function
	squash->fprop(sum, out);
}

void s_layer::bprop(state_idx *in, state_idx *out)
{
	// 1.
	squash->bprop(sum, out);
	// 2.
	{ idx_bloop2(lha, sum->dx, double, lb, bias->dx, double) {
		idx_sumacc(lha, lb);
	}}
	// 3.
	{ idx_bloop3(lcdx, coeff->dx, double, ltdx, sum->dx, double,
			lsx, sub->x, double) {
		idx_dotacc(lsx, ltdx, lcdx);
	}}
	// 4.
	{ idx_bloop4(lidx, in->dx, double, lsdx, sub->dx, double,
			lcx, coeff->x, double, ltdx2, sum->dx, double) {
		idx_dotc(ltdx2, lcx.get(), lsdx);
		idx_m2oversample(lsdx, stridei, stridej, lidx);
	}}
}

void s_layer::bbprop(state_idx *in, state_idx *out)
{
	// 1.
	squash->bbprop(sum, out);
	// 2.
	{ idx_bloop2(lha, sum->ddx, double, lb, bias->ddx, double) {
		idx_sumacc(lha, lb);
	}}
	// 3.
	{ idx_bloop3(lcdx, coeff->ddx, double, ltdx, sum->ddx, double,
			lsx, sub->x, double) {
		idx_m2squdotm2acc(lsx, ltdx, lcdx);
	}}
	// 4.
	{ idx_bloop4(lidx, in->ddx, double, lsdx, sub->ddx, double,
			lcx, coeff->x, double, ltdx2, sum->ddx, double) {
		double cf = lcx.get();
		idx_dotc(ltdx2, cf * cf, lsdx);
		idx_m2oversample(lsdx, stridei, stridej, lidx);
	}}
}

void s_layer::forget(forget_param_linear &fp) {
	double c = fp.value / pow(stridei * stridej, fp.exponent);
	idx_clear(bias->x);
	idx_fill(coeff->x, c);
}
////////////////////////////////////////////////////////////////////////

logadd_layer::logadd_layer(intg thick, intg si, intg sj) {
  expdist = Idx<double>(thick, si, sj);
  sumexp = Idx<double>(thick);		// scaled partition function
}

void logadd_layer::fprop(state_idx *in, state_idx *out) {
  intg thick = in->x.dim(0);
	intg si = in->x.dim(1);
	intg sj = in->x.dim(2);
	expdist.resize(thick, si, sj);
  out->x.resize(thick);
	if (1 == (si * sj)) {
		// save time and precision if no replication
		Idx<double> inx(in->x.select(2, 0));
		Idx<double> m(inx.select(1, 0));
		Idx<double> ed(expdist.select(2, 0));
		Idx<double> ed1(ed.select(1, 0));
		idx_fill(ed1, 1.0);
		idx_fill(sumexp, 1.0);