void apply(cube<real>& v, real bias) noexcept override
 {
     real* d = v.data();
     size_t  n = v.num_elements();
     for ( size_t i = 0; i < n; ++i )
         d[i] = f_(d[i] + bias);
 }

//HELIX ANIMATION 
uint8_t animation::helix(cube &c){
    float X = 0;
    float Y = 0;
    float Z = 0;
    
    slow++;
    if(slow < 100){
        return 0;
    }
    slow = 0;
    
    c.clearVoxels();
    
    //use my fancy pants sine function
    for(uint8_t i = 0; i < 3; i++){
        for(uint8_t z = 0; z < 4; z++){
            Z = z*52;
            X = get_sinA(Z + phase + 18*i);
            Y = get_sinA(Z + phase + 90 + 18*i);
            X = (X+1)*4;
            Y = (Y+1)*4;
            c.setVoxel((uint8_t)X, (uint8_t)Y, z);
            c.setVoxel((uint8_t)(8-X), (uint8_t)(8-Y), z+4);
        }
    }
    
    //increment the phase 
    phase+=18;
    
    if(phase > 360){
        phase -= 360;
    }
    
    return 1;
}

cube barrel_distort_rgb(const cube &F, double offset_x) {
  cube G(F.n_rows, F.n_cols, F.n_slices);
  G.slice(0) = barrel_distort(F.slice(0), offset_x);
  G.slice(1) = barrel_distort(F.slice(1), offset_x);
  G.slice(2) = barrel_distort(F.slice(2), offset_x);
  return G;
}

SEXP move_B(const mat& y, cube& B, const mat& kappa_star, const field<mat>& C,
            //const field<sp_mat>& C,
            mat& gamma, const mat& D, const ucolvec& s, double tau_e)
{
     BEGIN_RCPP
     // N x T matrix of standardized counts, y
     // B = (B_1,...,B_K), where B_k is N x T matrix for iGMRF term k
     // N x T matrix, gamma = sum_{k=1^K}(B_k)
     // K x T, D, where row k contains T x 1 diagonal elements of Q_k
     // sample cluster assignments, s(1), ..., s(N)  
     // Q = (Q_1,...,Q_K), where Q_k is a T x T de-scaled iGMRF precision matrix
     // C = (C_1,...,C_K), where C_k = D_k^-1 * Omega_k, 
     // where Omega_k is the T x T adjacency matrix for iGMRF term, k
     // D is a K x T matrix where row k contains T diagonal elements of Q_k
     // K x M matrix, kappa_star records locations for each iGMRF term 
     int K     = B.n_slices;
     int N     = y.n_rows;
     int T     = D.n_cols;
     colvec bbar_ki(T); bbar_ki.zeros();
     rowvec gammatilde_ki(T); gammatilde_ki.zeros();
     rowvec ytilde_ki(T); ytilde_ki.zeros();
     rowvec d_k(T);
     vec zro(T); zro.zeros();
     double e_ij, phi_ij, bkij, h_ij;
     int k = 0, i = 0, j = 0; /* loop variables */
     for( k = 0; k < K; k++ ) /* over iGMRF terms */
     {
          d_k                           = D.row(k);
          for( i = 0; i < N; i++ ) /* over units */
          {
               /* take out T x 1, b_ki, to be sampled from gamma */
               //gammatilde_ki        = gamma.row(i) - B.slice(k).row(i); /* 1 x T */
               //ytilde_ki            = y.row(i) - gammatilde_ki;
               gammatilde_ki            = gamma.row(i); /* 1 x T */
               // mean of univariate iGMRF, b_kij = 1/d_kj * (omega_kj(-j) * b_ki(-j))
               bbar_ki                  = C(k,0) * B.slice(k).row(i).t(); /* T x 1 */
               for( j = 0; j < T; j++ ) /* over time points */
               {
                    gammatilde_ki(j)    -= B.slice(k)(i,j);
                    ytilde_ki(j)        = y(i,j) - gammatilde_ki(j);
                    // mean of univariate iGMRF, b_kij = 1/d_kj * (omega_kj(-j) * b_ki(-j))
                    B.slice(k)(i,j)     = 0;
                    //bbar_ki(j)          = dot( C(k,0).row(j), B.slice(k).row(i) );
                    e_ij                = tau_e*ytilde_ki(j) 
                                             + d_k(j)*kappa_star(k,s(i)) * bbar_ki(j);
                    phi_ij              = tau_e + d_k(j)*kappa_star(k,s(i));
                    h_ij                = (e_ij / phi_ij);
                    bkij                = rnorm( 1, (h_ij), sqrt(1/phi_ij) )[0];
                    B.slice(k)(i,j)     = bkij;
                    /* put back new sampled values for b_kij */
                    gammatilde_ki(j)    += B.slice(k)(i,j);
               } /* end loop j over time points */
               /* put back new sampled values for b_ki */
               //gamma.row(i)         = gammatilde_ki + B.slice(k).row(i); 
               gamma.row(i)        = gammatilde_ki;
          } /* end loop i over units */
     } /* end loop k over iGMRF terms */
     END_RCPP
} /* end function ustep to sample B_1,...,B_K */

inline bool equal(const cube<T>& c1, const cube<T>& c2)
{
    if ( c1.n_elem != c2.n_elem )
    {
        return false;
    }

    return std::equal(c1.memptr(), c1.memptr() + c1.n_elem, c2.memptr());
}

 typename std::enable_if<has_public_member_grad<T>::value>::type
 apply_grad(cube<real>& g, const cube<real>& f, const T& fn) noexcept
 {
     real* gp = g.data();
     const real* fp = f.data();
     size_t  n = g.num_elements();
     for ( size_t i = 0; i < n; ++i )
         gp[i] *= fn.grad(fp[i]);
 }

inline void pairwise_div( cube<T>& a, const cube<T>& b )
{
    ZI_ASSERT(a.n_elem==b.n_elem);
    T* ap = a.memptr();
    const T* bp = b.memptr();

    for ( size_t i = 0; i < a.n_elem; ++i )
        ap[i] /= bp[i];
}

inline void convolve_constant( cube<T> const & a,
                               identity_t<T> b,
                               cube<T> & r) noexcept
{
    ZI_ASSERT(size(a)==size(r));

    T const * ap = a.data();
    T * rp = r.data();

    for ( size_t i = 0; i < r.num_elements(); ++i )
        rp[i] = ap[i] * b;
}

inline T convolve_constant_flipped( cube<T> const & a,
                                    cube<T> const & b ) noexcept
{
    ZI_ASSERT(size(a)==size(b));

    T r = 0;
    T const * ap = a.data();
    T const * bp = b.data();

    for ( size_t i = 0; i < a.num_elements(); ++i )
        r += ap[i] * bp[i];

    return r;
}

/**********************************************************************
*   MSE Class
***********************************************************************/
double MSE::cost(cube pred, cube y){
    dbg_assert(y.n_cols == 1 && y.n_slices == 1);
    pred.reshape(pred.n_elem,1,1);
    y.reshape(y.n_elem,1,1);
    
    double retn = 0.0f;
    for(int i=0; i < y.n_elem; i++) {
        retn += pow(pred(i,0,0) - y(i,0,0),2);
    }
    retn /= y.n_elem;
    retn *= 0.5f;
    //return(retn);
    return(0.5 * mean(mean(square(pred.slice(0) - y.slice(0)))));        
}

//RANDOM LIGHTS
uint8_t animation::randomExpand(cube &c){
    uint16_t count = 0;
	uint8_t rX, rY, rZ;
    
    slow++;
    if(slow < 20){
        return 0;
    }
    slow = 0;
    
    randomSeed(analogRead(0));
    rX = rand()%8+0;
    rY = rand()%8+0;
    rZ = rand()%8+0;
	
    //find an empty voxel
    while(c.getVoxel(rX, rY, rZ) == 1 && dir == 1){
        rX = random(0, 8);
        rY = random(0, 8);
        rZ = random(0, 8);
        count++;
        
        if(count > 200){
            dir *= -1;
        }
    }
    count = 0;

    //find a full voxel
    while(c.getVoxel(rX, rY, rZ) == 0 && dir == -1){
        rX = random(0, 8);
        rY = random(0, 8);
        rZ = random(0, 8);
        count++;

        if(count > 200){
            dir *= -1;
        }
    }

    //fill or clear the voxel found
    if(dir == 1){
        c.setVoxel(rX, rY, rZ);
    }else{
        c.clearVoxel(rX, rY, rZ);
    }
    
    return 1;
}

cube Utils::conv3(cube body, cube kernel)
{
	if (body.max() <= 0)
	{
		return zeros(body.n_rows, body.n_cols, body.n_slices);
	}
	cube Q(body);
	for (int i = 1; i < body.n_slices - 1; i++)
	{
		Q.slice(i) = conv2(body.slice(i), kernel.slice(1), "same") +
			conv2(body.slice(i - 1), kernel.slice(0), "same") +
			conv2(body.slice(i + 1), kernel.slice(2), "same");
	}
	return Q;
}

    // performs inplace dropout backward
    inline void dropout_backward(cube<real> & g)
    {
        ZI_ASSERT(mask_);

        size_t s = g.num_elements();
        for ( size_t i = 0; i < s; ++i )
        {
            if ( mask_->data()[i] )
                g.data()[i] *= scale();
            else
                g.data()[i]  = static_cast<real>(0);
        }

        // Should we reset mask_ here?
    }

//SINGLE PLANE BOUNCING
uint8_t animation::bouncePlane(cube &c, uint8_t axis){
    slow++;
    if(slow < 200){
        return 0;
    }
    slow = 0;
    
    //the X animation looks the same but actually clears everything in its path
    //rather than clear everything at the start, it makes a simple but cool 
    //transition between some animations
    if(axis != 'X'){
        c.clearVoxels();
    }
    
    switch(axis){
        case 'X':
            for(uint8_t z = 0; z < Zd; z++){
                for(uint8_t y = 0; y < Yd; y++){
                    for(uint8_t x = 0; x < Xd; x++){
                        if(x == pos - dir){
                            c.clearVoxel(x, y, z);
                        }
                        c.setVoxel(pos, y, z);
                    }
                }
            }
            break;
        case 'Y':
            for(uint8_t x = 0; x < Xd; x++){
                for(uint8_t z = 0; z < Zd; z++){
                    c.setVoxel(x, pos, z);
                }
            }
            break;
        case 'Z':
            for(uint8_t x = 0; x < Xd; x++){
                for(uint8_t y = 0; y < Yd; y++){
                    c.setVoxel(x, y, pos);
                }
            }
            break;
    }
    
    //bounce the pos variable between 0 and 7
    bouncePos();
    
    return 1;
}   

cube imresize2(const cube &C, uword m, uword n) {
  cube F(m, n, C.n_slices);
  for (uword k = 0; k < C.n_slices; k++) {
    F.slice(k) = imresize2(C.slice(k), m, n);
  }
  return F;
}

cube noncont_slices(const cube & X, const uvec & index, int r_start, int c_start, int r_stop, int c_stop){
	cube Xsub(r_stop-r_start+1, c_stop-c_start+1, index.n_elem);
	for(unsigned int i=0; i<index.n_elem; i++){
		Xsub.slice(i) = X.slice(index[i]).submat(r_start, c_start, r_stop, c_stop);
	}
	return Xsub; 
}

cube Utils::erode(cube body, cube kernel)
{
	int minx, miny, minz, maxx, maxy, maxz;
	Utils::bounds(body, minx, miny, minz, maxx, maxy, maxz);
	cube result = cube(body);
	int size_x = kernel.n_cols / 2;
	int size_y = kernel.n_rows / 2;
	int size_z = kernel.n_slices / 2;
	maxx -= kernel.n_cols;
	maxy -= kernel.n_rows;
	maxz -= kernel.n_slices;
	for (int i = minx; i <= maxx; i++)
	{
		for (int j = miny; j <= maxy; j++)
		{
			for (int k = minz; k <= maxz; k++)
			{
				cube subcube = body.subcube(j, i, k, size(kernel));
				subcube -= kernel;
				result(j + size_y, i + size_x, k + size_z) = subcube.min() < 0 ? 0 : 1;
			}
		}
	}
	return result;
}

cube conv2(const cube &F, const mat &H) {
  cube G(F.n_rows, F.n_cols, F.n_slices);
  for (uword i = 0; i < F.n_slices; i++) {
    G.slice(i) = conv2(F.slice(i), H);
  }
  return G;
}

inline unique_cube<T> crop_right( const cube<T>& in, const vec3s& s )
{
    auto r = pool<T>::get_unique(s);
    vec3s b = size(in) - s;
    *r = in.subcube(b[0], b[1], b[2], b[0]+s[0]-1, b[1]+s[1]-1, b[2]+s[2]-1);
    return r;
}

//! compute (trial)point-to-(model)point Mahalanobis distance
//! @param[in] index		index of the points (in trial and model) to be compared
//! @param[in] &trial		reference to the trial
//! @param[in] &model		reference to the model
//! @param[in] &variance	reference to the model variance
//! @return 			    Mahalanobis distance between trial-point and model-point
float Classifier::mahalanobisDist(int index,mat &trial,mat &model,cube &variance)
{
	mat difference = trial.col(index) - model.col(index);
	mat distance = (difference.t() * (variance.slice(index)).i()) * difference;

	return distance(0,0);
}