Repository: IBM/aihwkit (language: python)

Issue:
Decay weight value to arbitrary point 
## Description and motivation

Currently, the `lifetime` parameter decays the weight value exponentially to 0. However, for some materials the decay might not go to strictly zero.  

## Proposed solution

Add a "decay_to" parameter that sets the decay to value for each device. 

## Alternatives and other information

<!-- A clear and concise description of any alternative solutions or features you've considered. -->



Discussion / hints:

Fix patch:
diff --git a/src/aihwkit/simulator/configs/devices.py b/src/aihwkit/simulator/configs/devices.py
index 003fefd0..92ac933c 100644
--- a/src/aihwkit/simulator/configs/devices.py
+++ b/src/aihwkit/simulator/configs/devices.py
@@ -86,7 +86,7 @@ class PulsedDevice(_PrintableMixin):
 
     **Decay**:
 
-    .. math:: w_{ij} \leftarrow w_{ij}\,(1-\alpha_\text{decay}\delta_{ij})
+    .. math:: w_{ij} \leftarrow (w_{ij} - b_{ij})\,(1-\alpha_\text{decay}\delta_{ij}) + b_{ij}
 
     Weight decay is only activated by inserting a specific call to
     :meth:`~aihwkit.simulator.tiles.base.Base.decay_weights`, which is
@@ -97,6 +97,16 @@ class PulsedDevice(_PrintableMixin):
     initialization. :math:`\alpha_\text{decay}` is a scaling factor
     that can be given during run-time.
 
+    The bias :math:`b_{ij}` is given by the reset bias and which is
+    determined by the parameter ``reset`` (mean value) and
+    ``reset_dtod``  (device-to-device variability). Thus
+
+    .. math:: b_{ij} = \mu_\text{reset} \left(1 + \sigma_\text{reset-dtod}\xi\right)
+
+    Note that the reset bias is also applied in case the device is
+    reset (see above).
+
+
     **Diffusion**:
 
     Similar to the decay, diffusion is only activated by inserting a
@@ -125,6 +135,7 @@ class PulsedDevice(_PrintableMixin):
         mini-batch but requires an explicit call to
         :meth:`~aihwkit.simulator.tiles.base.Base.drift_weights` each
         time the drift should be applied.
+
     """
 
     bindings_class: ClassVar[Type] = devices.PulsedResistiveDeviceParameter
@@ -186,7 +197,7 @@ class PulsedDevice(_PrintableMixin):
     """No up-down differences and device-to-device variability in the bounds
     for the devices in the bias row."""
 
-    reset: float = 0.01
+    reset: float = 0.0
     """The reset values and spread per cross-point ``ij`` when using reset
     functionality of the device."""
 
diff --git a/src/rpucuda/cuda/cuda_math_util.cu b/src/rpucuda/cuda/cuda_math_util.cu
index 59ed86ff..c18afa54 100644
--- a/src/rpucuda/cuda/cuda_math_util.cu
+++ b/src/rpucuda/cuda/cuda_math_util.cu
@@ -243,11 +243,7 @@ void ger<double>(
 // W += A
 template <typename T> __global__ void kernelElemAdd(T *dev_W, const int size, const T *dev_A) {
 
-  volatile unsigned int tid = blockDim.x * blockIdx.x + threadIdx.x;
-
-  if (tid < size) {
-    dev_W[tid] += dev_A[tid];
-  }
+  RPU_CUDA_1D_KERNEL_LOOP(idx, size) { dev_W[idx] += dev_A[idx]; }
 }
 template <typename T>
 void elemadd(const CudaContext *context, T *dev_W, const int size, const T *dev_A) {
@@ -264,11 +260,9 @@ template void elemadd<double>(const CudaContext *, double *, const int, const do
 // W = W.*W
 template <typename T> __global__ void kernelElemPow2(T *dev_W, const int size) {
 
-  volatile unsigned int tid = blockDim.x * blockIdx.x + threadIdx.x;
-
-  if (tid < size) {
-    T x = dev_W[tid];
-    dev_W[tid] = x * x;
+  RPU_CUDA_1D_KERNEL_LOOP(idx, size) {
+    T x = dev_W[idx];
+    dev_W[idx] = x * x;
   }
 }
 template <typename T> void elempow2(const CudaContext *context, T *dev_W, const int size) {
@@ -285,11 +279,7 @@ template void elempow2<double>(const CudaContext *, double *, const int);
 // V = abs(W )
 template <typename T> __global__ void kernelElemAbs(T *dev_V, const T *dev_W, const int size) {
 
-  volatile unsigned int tid = blockDim.x * blockIdx.x + threadIdx.x;
-
-  if (tid < size) {
-    dev_V[tid] = fabs(dev_W[tid]);
-  }
+  RPU_CUDA_1D_KERNEL_LOOP(idx, size) { dev_V[idx] = fabs(dev_W[idx]); }
 }
 template <typename T>
 void elemabs(const CudaContext *context, T *dev_V, const T *dev_W, const int size) {
@@ -307,16 +297,13 @@ template void elemabs<double>(const CudaContext *, double *, const double *, con
 template <typename T, typename T_A>
 __global__ void kernelElemAddScale(T *dev_W, const int size, const T_A *dev_A, const T beta) {
 
-  volatile unsigned int tid = blockDim.x * blockIdx.x + threadIdx.x;
-
   T b = beta;
-
-  if (tid < size) {
-    T w = dev_W[tid];
-    T a = dev_A[tid];
+  RPU_CUDA_1D_KERNEL_LOOP(idx, size) {
+    T w = dev_W[idx];
+    T a = dev_A[idx];
 
     w += a * b;
-    dev_W[tid] = w;
+    dev_W[idx] = w;
   }
 }
 template <typename T, typename T_A>
@@ -340,15 +327,13 @@ elemaddscale<double, float>(const CudaContext *, double *, const int, const floa
 template <typename T>
 __global__ void kernelElemAddScale(T *dev_W, const int size, const T *dev_A, const T *dev_B) {
 
-  volatile unsigned int tid = blockDim.x * blockIdx.x + threadIdx.x;
-
-  if (tid < size) {
-    T w = dev_W[tid];
-    T a = dev_A[tid];
-    T b = dev_B[tid];
+  RPU_CUDA_1D_KERNEL_LOOP(idx, size) {
+    T w = dev_W[idx];
+    T a = dev_A[idx];
+    T b = dev_B[idx];
 
     w += a * b;
-    dev_W[tid] = w;
+    dev_W[idx] = w;
   }
 }
 template <typename T>
@@ -371,20 +356,18 @@ template <typename T, typename T_A>
 __global__ void
 kernelElemASB02(T *dev_W, const int size, const T_A *dev_A, const T *dev_B, float *dev_4params) {
 
-  volatile unsigned int tid = blockDim.x * blockIdx.x + threadIdx.x;
-
-  if (tid < size) {
-    T w = dev_W[tid];
-    T a = dev_A[tid];
-    T b = dev_B[tid];
-    float4 parij = reinterpret_cast<float4 *>(dev_4params)[tid];
+  RPU_CUDA_1D_KERNEL_LOOP(idx, size) {
+    T w = dev_W[idx];
+    T a = dev_A[idx];
+    T b = dev_B[idx];
+    float4 parij = reinterpret_cast<float4 *>(dev_4params)[idx];
 
     w += a * b;
     // check bounds
     w = (w > parij.z) ? parij.z : w;
     w = (w < parij.x) ? parij.x : w;
 
-    dev_W[tid] = w;
+    dev_W[idx] = w;
   }
 }
 template <typename T, typename T_A>
@@ -411,49 +394,55 @@ template void elemasb02<double, float>(
 
 // sat(W *= A)
 template <typename T>
-__global__ void kernelElemScale(T *dev_W, const int size, const T *dev_A, float *dev_4params) {
-
-  volatile unsigned int tid = blockDim.x * blockIdx.x + threadIdx.x;
+__global__ void
+kernelElemScale(T *dev_W, const int size, const T *dev_A, float *dev_4params, const T *dev_shift) {
 
-  if (tid < size) {
-    T w = dev_W[tid];
-    T a = dev_A[tid];
-    float4 parij = reinterpret_cast<float4 *>(dev_4params)[tid];
+  bool with_shift = dev_shift != nullptr;
+  RPU_CUDA_1D_KERNEL_LOOP(idx, size) {
+    T w = dev_W[idx];
+    T a = dev_A[idx];
+    T s = with_shift ? dev_shift[idx] : 0.0;
+    float4 parij = reinterpret_cast<float4 *>(dev_4params)[idx];
 
-    w *= a;
+    w = (w - s) * a + s;
     // check bounds
     w = (w > parij.z) ? parij.z : w;
     w = (w < parij.x) ? parij.x : w;
 
-    dev_W[tid] = w;
+    dev_W[idx] = w;
   }
 }
 template <typename T>
 void elemscale(
-    const CudaContext *context, T *dev_W, const int size, const T *dev_A, float *dev_4params) {
+    const CudaContext *context,
+    T *dev_W,
+    const int size,
+    const T *dev_A,
+    float *dev_4params,
+    const T *dev_shift) {
 
   int nthreads = context->getNThreads();
   int nblocks = context->getNBlocks(size, nthreads);
   kernelElemScale<T>
-      <<<nblocks, nthreads, 0, context->getStream()>>>(dev_W, size, dev_A, dev_4params);
+      <<<nblocks, nthreads, 0, context->getStream()>>>(dev_W, size, dev_A, dev_4params, dev_shift);
 }
-template void elemscale<float>(const CudaContext *, float *, const int, const float *, float *);
+template void
+elemscale<float>(const CudaContext *, float *, const int, const float *, float *, const float *);
 #ifdef RPU_USE_DOUBLE
-template void elemscale<double>(const CudaContext *, double *, const int, const double *, float *);
+template void elemscale<double>(
+    const CudaContext *, double *, const int, const double *, float *, const double *);
 #endif
 
 // sat(W)
 template <typename T> __global__ void kernelElemSat(T *dev_W, const int size, float *dev_4params) {
 
-  volatile unsigned int tid = blockDim.x * blockIdx.x + threadIdx.x;
-
-  if (tid < size) {
-    T w = dev_W[tid];
-    float4 parij = reinterpret_cast<float4 *>(dev_4params)[tid];
+  RPU_CUDA_1D_KERNEL_LOOP(idx, size) {
+    T w = dev_W[idx];
+    float4 parij = reinterpret_cast<float4 *>(dev_4params)[idx];
     // check bounds
     w = (w > parij.z) ? parij.z : w;
     w = (w < parij.x) ? parij.x : w;
-    dev_W[tid] = w;
+    dev_W[idx] = w;
   }
 }
 template <typename T>
@@ -471,22 +460,29 @@ template void elemsat<double>(const CudaContext *, double *, const int, float *)
 // sat(W *= 1+(A-1)*alpha)
 template <typename T>
 __global__ void kernelElemScaleAlpha(
-    T *dev_W, const int size, const T *dev_A, float *dev_4params, const T alpha_in) {
+    T *dev_W,
+    const int size,
+    const T *dev_A,
+    float *dev_4params,
+    const T alpha_in,
+    const T *dev_shift) {
 
-  volatile unsigned int tid = blockDim.x * blockIdx.x + threadIdx.x;
   volatile T alpha = alpha_in;
+  bool with_shift = dev_shift != nullptr;
+  RPU_CUDA_1D_KERNEL_LOOP(idx, size) {
+    T w = dev_W[idx];
+    T a = dev_A[idx];
+    T s = with_shift ? dev_shift[idx] : 0.0;
+    float4 parij = reinterpret_cast<float4 *>(dev_4params)[idx];
 
-  if (tid < size) {
-    T w = dev_W[tid];
-    T a = dev_A[tid];
-    float4 parij = reinterpret_cast<float4 *>(dev_4params)[tid];
+    T scale = 1.0 + alpha * (a - 1.0);
+    w = (w - s) * scale + s;
 
-    w *= 1.0 + alpha * (a - 1.0);
     // check bounds
     w = (w > parij.z) ? parij.z : w;
     w = (w < parij.x) ? parij.x : w;
 
-    dev_W[tid] = w;
+    dev_W[idx] = w;
   }
 }
 template <typename T>
@@ -496,30 +492,35 @@ void elemscalealpha(
     const int size,
     const T *dev_A,
     float *dev_4params,
-    const T alpha) {
+    const T alpha,
+    const T *dev_shift) {
 
   int nthreads = context->getNThreads();
   int nblocks = context->getNBlocks(size, nthreads);
-  kernelElemScaleAlpha<T>
-      <<<nblocks, nthreads, 0, context->getStream()>>>(dev_W, size, dev_A, dev_4params, alpha);
+  kernelElemScaleAlpha<T><<<nblocks, nthreads, 0, context->getStream()>>>(
+      dev_W, size, dev_A, dev_4params, alpha, dev_shift);
 }
-template void
-elemscalealpha<float>(const CudaContext *, float *, const int, const float *, float *, const float);
+template void elemscalealpha<float>(
+    const CudaContext *, float *, const int, const float *, float *, const float, const float *);
 #ifdef RPU_USE_DOUBLE
 template void elemscalealpha<double>(
-    const CudaContext *, double *, const int, const double *, float *, const double);
+    const CudaContext *,
+    double *,
+    const int,
+    const double *,
+    float *,
+    const double,
+    const double *);
 #endif
 
 // W += A, A = W
 template <typename T> __global__ void kernelElemAddCopy(T *dev_W, T *dev_A, const int size) {
 
-  volatile unsigned int tid = blockDim.x * blockIdx.x + threadIdx.x;
-
-  if (tid < size) {
-    T w = dev_W[tid];
-    w += dev_A[tid];
-    dev_W[tid] = w;
-    dev_A[tid] = w;
+  RPU_CUDA_1D_KERNEL_LOOP(idx, size) {
+    T w = dev_W[idx];
+    w += dev_A[idx];
+    dev_W[idx] = w;
+    dev_A[idx] = w;
   }
 }
 template <typename T>
@@ -538,19 +539,17 @@ template void elemaddcopy<double>(const CudaContext *, double *, double *, const
 template <typename T>
 __global__ void kernelElemAddCopySat(T *dev_W, T *dev_A, const int size, const float *dev_4params) {
 
-  volatile unsigned int tid = blockDim.x * blockIdx.x + threadIdx.x;
-
-  if (tid < size) {
-    T w = dev_W[tid];
-    T a = dev_A[tid];
-    const float4 parij = reinterpret_cast<const float4 *>(dev_4params)[tid];
+  RPU_CUDA_1D_KERNEL_LOOP(idx, size) {
+    T w = dev_W[idx];
+    T a = dev_A[idx];
+    const float4 parij = reinterpret_cast<const float4 *>(dev_4params)[idx];
     w += a;
     // check bounds
     w = (w > parij.z) ? parij.z : w;
     w = (w < parij.x) ? parij.x : w;
     a = w;
-    dev_W[tid] = w;
-    dev_A[tid] = a;
+    dev_W[idx] = w;
+    dev_A[idx] = a;
   }
 }
 template <typename T>
@@ -581,22 +580,20 @@ __global__ void kernelElemResetSat(
     const T thres,
     const float *dev_4params) {
 
-  volatile unsigned int tid = blockDim.x * blockIdx.x + threadIdx.x;
-
-  if (tid < size) {
+  RPU_CUDA_1D_KERNEL_LOOP(idx, size) {
     T th = thres;
-    T p = dev_P[tid];
-    T a = dev_A[tid]; // load a bit more ?
-    T b = dev_B[tid];
-    T w = dev_W[tid];
-    const float4 parij = reinterpret_cast<const float4 *>(dev_4params)[tid];
+    T p = dev_P[idx];
+    T a = dev_A[idx]; // load a bit more ?
+    T b = dev_B[idx];
+    T w = dev_W[idx];
+    const float4 parij = reinterpret_cast<const float4 *>(dev_4params)[idx];
 
     if (p < th) {
       w = a + b;
       // check bounds [only those that changed]
       w = (w > parij.z) ? parij.z : w;
       w = (w < parij.x) ? parij.x : w;
-      dev_W[tid] = w;
+      dev_W[idx] = w;
     }
   }
 }
@@ -713,13 +710,11 @@ template void elemresetsatmsk<double>(
 template <typename T>
 __global__ void kernelElemSubCopy(T *dev_W, T *dev_A, const int size, const T scale) {
 
-  volatile unsigned int tid = blockDim.x * blockIdx.x + threadIdx.x;
-
-  if (tid < size) {
-    T new_w = dev_W[tid];
-    T old_w = dev_A[tid];
-    dev_W[tid] = old_w;
-    dev_A[tid] = scale * (new_w - old_w); // dw
+  RPU_CUDA_1D_KERNEL_LOOP(idx, size) {
+    T new_w = dev_W[idx];
+    T old_w = dev_A[idx];
+    dev_W[idx] = old_w;
+    dev_A[idx] = scale * (new_w - old_w); // dw
   }
 }
 template <typename T>
@@ -736,10 +731,8 @@ template void elemsubcopy<double>(const CudaContext *, double *, double *, const
 
 // set all elements to a
 template <typename T> __global__ void kernelSetConstAlpha(T *values, int size, T alpha) {
-  int tid = blockDim.x * blockIdx.x + threadIdx.x;
-  if (tid < size) {
-    values[tid] = alpha;
-  }
+
+  RPU_CUDA_1D_KERNEL_LOOP(idx, size) { values[idx] = alpha; }
 }
 
 template <typename T>
@@ -759,10 +752,8 @@ template void elemconst<char>(const CudaContext *, char *, const int, const char
 
 // w = max(min(w,|a|),-|a|)
 template <typename T> __global__ void kernelAClip(T *values, int size, T abs_a) {
-  int tid = blockDim.x * blockIdx.x + threadIdx.x;
-  if (tid < size) {
-    values[tid] = MIN(MAX(values[tid], -abs_a), abs_a);
-  }
+
+  RPU_CUDA_1D_KERNEL_LOOP(idx, size) { values[idx] = MIN(MAX(values[idx], -abs_a), abs_a); }
 }
 
 template <typename T> void aclip(const CudaContext *context, T *W, const int size, const T a) {
@@ -778,10 +769,8 @@ template void aclip<double>(const CudaContext *, double *, const int, const doub
 
 // w = max(w,a)
 template <typename T> __global__ void kernelElemmax(T *values, int size, T a) {
-  int tid = blockDim.x * blockIdx.x + threadIdx.x;
-  if (tid < size) {
-    values[tid] = MAX(values[tid], a);
-  }
+
+  RPU_CUDA_1D_KERNEL_LOOP(idx, size) { values[idx] = MAX(values[idx], a); }
 }
 
 template <typename T> void elemmax(const CudaContext *context, T *W, const int size, const T a) {
@@ -797,10 +786,8 @@ template void elemmax<double>(const CudaContext *, double *, const int, const do
 
 // w = min(w,a)
 template <typename T> __global__ void kernelElemmin(T *values, int size, T a) {
-  int tid = blockDim.x * blockIdx.x + threadIdx.x;
-  if (tid < size) {
-    values[tid] = MIN(values[tid], a);
-  }
+
+  RPU_CUDA_1D_KERNEL_LOOP(idx, size) { values[idx] = MIN(values[idx], a); }
 }
 
 template <typename T> void elemmin(const CudaContext *context, T *W, const int size, const T a) {
@@ -816,10 +803,10 @@ template void elemmin<double>(const CudaContext *, double *, const int, const do
 
 // w = w<a?0:w elementwise
 template <typename T> __global__ void kernelElemSetBelowZero(T *values, int size, T a) {
-  int tid = blockDim.x * blockIdx.x + threadIdx.x;
-  if (tid < size) {
-    T v = values[tid];
-    values[tid] = v < a ? (T)0.0 : v;
+
+  RPU_CUDA_1D_KERNEL_LOOP(idx, size) {
+    T v = values[idx];
+    values[idx] = v < a ? (T)0.0 : v;
   }
 }
 
@@ -840,30 +827,21 @@ template void elemsetbelowzero<double>(const CudaContext *, double *, const int,
 template <typename T>
 __global__ void kernelElemAdd(T *dev_W, const int size, const T *dev_A, const T *dev_B) {
 
-  volatile unsigned int tid = blockDim.x * blockIdx.x + threadIdx.x;
-  if (tid < size) {
-    dev_W[tid] = dev_A[tid] + dev_B[tid];
-  }
+  RPU_CUDA_1D_KERNEL_LOOP(idx, size) { dev_W[idx] = dev_A[idx] + dev_B[idx]; }
 }
 template <typename T>
 __global__ void kernelElemSub(T *dev_W, const int size, const T *dev_A, const T *dev_B) {
 
-  volatile unsigned int tid = blockDim.x * blockIdx.x + threadIdx.x;
-  if (tid < size) {
-    dev_W[tid] = dev_A[tid] - dev_B[tid];
-  }
+  RPU_CUDA_1D_KERNEL_LOOP(idx, size) { dev_W[idx] = dev_A[idx] - dev_B[idx]; }
 }
 
 template <typename T>
 __global__ void kernelElemWeightedSum(
     T *dev_W, const int size, const T *dev_A, const T a_in, const T *dev_B, const T b_in) {
 
-  volatile unsigned int tid = blockDim.x * blockIdx.x + threadIdx.x;
   T a = a_in;
   T b = b_in;
-  if (tid < size) {
-    dev_W[tid] = a * dev_A[tid] + b * dev_B[tid];
-  }
+  RPU_CUDA_1D_KERNEL_LOOP(idx, size) { dev_W[idx] = a * dev_A[idx] + b * dev_B[idx]; }
 }
 
 template <typename T>
@@ -912,16 +890,14 @@ template void elemweightedsum<double>(
 template <typename T>
 __global__ void kernelElemAverage(T *dev_W, const int size, T **dev_Ms, const int m) {
 
-  volatile unsigned int tid = blockDim.x * blockIdx.x + threadIdx.x;
-
-  if (tid < size) {
+  RPU_CUDA_1D_KERNEL_LOOP(idx, size) {
     T w = 0;
 
     for (int i = 0; i < m; i++) {
       const T *w1 = dev_Ms[i];
-      w += w1[tid];
+      w += w1[idx];
     }
-    dev_W[tid] = w / (T)m;
+    dev_W[idx] = w / (T)m;
   }
 }
 
@@ -1098,11 +1074,8 @@ template void copyWithoutBias<double>(
 template <typename OutputIteratorT, typename InputIteratorT>
 __global__ void kernelCopyWithIterator(
     OutputIteratorT out_tensor, const InputIteratorT in_tensor, const int sz_all) {
-  int tid = blockDim.x * blockIdx.x + threadIdx.x;
 
-  if (tid < sz_all) {
-    out_tensor[tid] = in_tensor[tid];
-  }
+  RPU_CUDA_1D_KERNEL_LOOP(idx, sz_all) { out_tensor[idx] = in_tensor[idx]; }
 }
 
 template <typename OutputIteratorT, typename InputIteratorT>
diff --git a/src/rpucuda/cuda/cuda_math_util.h b/src/rpucuda/cuda/cuda_math_util.h
index 271b4f53..44732c3c 100644
--- a/src/rpucuda/cuda/cuda_math_util.h
+++ b/src/rpucuda/cuda/cuda_math_util.h
@@ -121,15 +121,21 @@ void elemasb02(
     const T *B,
     float *dev_4params); // bounds in [0,2] // 4params and 2params always float !
 
-// sat(W *= A)
-template <typename T>
-void elemscale(const CudaContext *context, T *W, const int size, const T *A, float *dev_4params);
-
 // sat(W)
 template <typename T>
 void elemsat(const CudaContext *context, T *W, const int size, float *dev_4params);
 
-// sat(W *= 1+alpha*(A-1))
+// sat(W = (W - S) * A + S)
+template <typename T>
+void elemscale(
+    const CudaContext *context,
+    T *W,
+    const int size,
+    const T *A,
+    float *dev_4params,
+    const T *dev_shift = nullptr);
+
+// sat(W = (W - S) (1+(A-1)*alpha) + S)
 template <typename T>
 void elemscalealpha(
     const CudaContext *context,
@@ -137,7 +143,8 @@ void elemscalealpha(
     const int size,
     const T *A,
     float *dev_4params,
-    const T alpha);
+    const T alpha,
+    const T *dev_shift = nullptr);
 
 // W += A, A = W
 template <typename T> void elemaddcopy(const CudaContext *context, T *W, T *A, const int size);
diff --git a/src/rpucuda/cuda/rpucuda_pulsed.h b/src/rpucuda/cuda/rpucuda_pulsed.h
index 8d960f00..567ea523 100644
--- a/src/rpucuda/cuda/rpucuda_pulsed.h
+++ b/src/rpucuda/cuda/rpucuda_pulsed.h
@@ -201,9 +201,6 @@ template <typename T> class RPUCudaPulsed : public RPUCudaSimple<T> {
 private:
   PulsedMetaParameter<T> par_;
 
-  std::unique_ptr<CudaArray<T>> dev_decay_scale_ = nullptr;
-  std::unique_ptr<CudaArray<T>> dev_diffusion_rate_ = nullptr;
-
   // forward
   std::unique_ptr<CudaArray<T>> dev_f_x_vector_inc1_ = nullptr;
   std::unique_ptr<CudaArray<T>> dev_f_d_vector_inc1_ = nullptr;
diff --git a/src/rpucuda/cuda/rpucuda_pulsed_device.cu b/src/rpucuda/cuda/rpucuda_pulsed_device.cu
index 4ccf4909..8c0268d2 100644
--- a/src/rpucuda/cuda/rpucuda_pulsed_device.cu
+++ b/src/rpucuda/cuda/rpucuda_pulsed_device.cu
@@ -106,36 +106,49 @@ void PulsedRPUDeviceCuda<T>::populateFrom(const AbstractRPUDevice<T> &rpu_device
 
   PulsedRPUDeviceCudaBase<T>::populateFrom(rpu_device_in);
 
-  PulsedDPStruc<T> **sup = rpu_device.getDPStruc();
-
   // copy RPU to device variables
   float *tmp = new float[4 * size];
+
   T *tmp_ds = new T[size];
   T *tmp_df = new T[size];
   T *tmp_rb = new T[size];
   T *tmp_pw = new T[size];
 
+  T *mn = rpu_device.getMinBound()[0];
+  T *mx = rpu_device.getMaxBound()[0];
+  T *su = rpu_device.getScaleUp()[0];
+  T *sd = rpu_device.getScaleDown()[0];
+
   T *ds = rpu_device.getDecayScale()[0];
   T *df = rpu_device.getDiffusionRate()[0];
   T *rb = rpu_device.getResetBias()[0];
   T *pw = rpu_device.getPersistentWeights()[0];
+  bool with_diffusion = false;
+  bool with_reset_bias = false;
 
   for (int i = 0; i < d_size; ++i) {
     for (int j = 0; j < x_size; ++j) {
 
-      int k = j * (d_size * 4) +
-              4 * i; // transposed: col major required by cuBLAS .. linear arangmenet for now
-      tmp[k] = sup[i][j].min_bound;
-      tmp[k + 1] = sup[i][j].scale_down;
-      tmp[k + 2] = sup[i][j].max_bound;
-      tmp[k + 3] = sup[i][j].scale_up;
-
       int l_t = j * (d_size) + i;
       int l = i * (x_size) + j;
+      // transposed: col major required by cuBLAS .. linear arangmenet for now
+      int k = j * (d_size * 4) + 4 * i;
+      tmp[k] = mn[l];
+      tmp[k + 1] = sd[l];
+      tmp[k + 2] = mx[l];
+      tmp[k + 3] = su[l];
+
       tmp_ds[l_t] = ds[l];
       tmp_df[l_t] = df[l];
       tmp_rb[l_t] = rb[l];
       tmp_pw[l_t] = pw[l];
+
+      if (df[l] != 0.0) {
+        with_diffusion = true;
+      }
+      if (rb[l] != 0.0) {
+        with_reset_bias = true;
+      }
     }
   }
 
@@ -143,18 +156,17 @@ void PulsedRPUDeviceCuda<T>::populateFrom(const AbstractRPUDevice<T> &rpu_device
   dev_decay_scale_->assign(tmp_ds);
 
   // other parameters (on the fly)
-  const auto &par = getPar();
-  if (par.diffusion > 0.0) {
+  if (with_diffusion) {
     dev_diffusion_rate_ = RPU::make_unique<CudaArray<T>>(this->context_, size);
     dev_diffusion_rate_->assign(tmp_df);
   }
 
-  if (par.needsResetBias()) {
+  if (with_reset_bias) {
     dev_reset_bias_ = RPU::make_unique<CudaArray<T>>(this->context_, size);
     dev_reset_bias_->assign(tmp_rb);
   }
 
-  if (par.usesPersistentWeight()) {
+  if (getPar().usesPersistentWeight()) {
     dev_persistent_weights_ = RPU::make_unique<CudaArray<T>>(this->context_, size);
     dev_persistent_weights_->assign(tmp_pw);
   }
@@ -186,18 +198,22 @@ void PulsedRPUDeviceCuda<T>::decayWeights(T *weights, T alpha, bool bias_no_deca
 
   RPU::math::elemscalealpha<T>(
       this->context_, w, bias_no_decay ? MAX(this->size_ - this->d_size_, 0) : this->size_,
-      dev_decay_scale_->getData(), dev_4params_->getData(), alpha);
+      dev_decay_scale_->getData(), dev_4params_->getData(), alpha,
+      dev_reset_bias_ != nullptr ? dev_reset_bias_->getData() : nullptr);
 
   applyUpdateWriteNoise(weights);
 }
 
 template <typename T> void PulsedRPUDeviceCuda<T>::decayWeights(T *weights, bool bias_no_decay) {
 
-  T *w = getPar().usesPersistentWeight() ? dev_persistent_weights_->getData() : weights;
+  const auto &par = getPar();
+
+  T *w = par.usesPersistentWeight() ? dev_persistent_weights_->getData() : weights;
 
   RPU::math::elemscale<T>(
       this->context_, w, bias_no_decay ? MAX(this->size_ - this->d_size_, 0) : this->size_,
-      dev_decay_scale_->getData(), dev_4params_->getData());
+      dev_decay_scale_->getData(), dev_4params_->getData(),
+      dev_reset_bias_ != nullptr ? dev_reset_bias_->getData() : nullptr);
 
   applyUpdateWriteNoise(weights);
 }
diff --git a/src/rpucuda/rpu_linearstep_device.cpp b/src/rpucuda/rpu_linearstep_device.cpp
index 2e4030f8..12e98fda 100644
--- a/src/rpucuda/rpu_linearstep_device.cpp
+++ b/src/rpucuda/rpu_linearstep_device.cpp
@@ -32,8 +32,6 @@ void LinearStepRPUDevice<T>::populate(
 
   for (int i = 0; i < this->d_size_; ++i) {
 
-    PulsedDPStruc<T> *s = this->sup_[i];
-
     for (int j = 0; j < this->x_size_; ++j) {
 
       T diff_slope_at_bound_up = par.ls_decrease_up + par.ls_decrease_up_dtod * rng->sampleGauss();
@@ -52,12 +50,14 @@ void LinearStepRPUDevice<T>::populate(
            bound, which does not make sense both slopes are negative
            (sign of scale_up/scale_down here both positive and later
            corrected in update rule) */
-        w_slope_up_[i][j] = -diff_slope_at_bound_up * s[j].scale_up / par.w_max;
-        w_slope_down_[i][j] = -diff_slope_at_bound_down * s[j].scale_down / par.w_min;
+        w_slope_up_[i][j] = -diff_slope_at_bound_up * this->w_scale_up_[i][j] / par.w_max;
+        w_slope_down_[i][j] = -diff_slope_at_bound_down * this->w_scale_down_[i][j] / par.w_min;
       } else {
         /* In this case slope depends on the bound*/
-        w_slope_up_[i][j] = -diff_slope_at_bound_up * s[j].scale_up / s[j].max_bound;
-        w_slope_down_[i][j] = -diff_slope_at_bound_down * s[j].scale_down / s[j].min_bound;
+        w_slope_up_[i][j] =
+            -diff_slope_at_bound_up * this->w_scale_up_[i][j] / this->w_max_bound_[i][j];
+        w_slope_down_[i][j] =
+            -diff_slope_at_bound_down * this->w_scale_down_[i][j] / this->w_min_bound_[i][j];
       }
     }
   }
@@ -78,14 +78,14 @@ template <typename T> void LinearStepRPUDevice<T>::printDP(int x_count, int d_co
     for (int j = 0; j < x_count; ++j) {
       std::cout.precision(5);
       std::cout << i << "," << j << ": ";
-      std::cout << "[<" << this->sup_[i][j].max_bound << ",";
-      std::cout << this->sup_[i][j].min_bound << ">,<";
-      std::cout << this->sup_[i][j].scale_up << ",";
-      std::cout << this->sup_[i][j].scale_down << ">,<";
+      std::cout << "[<" << this->w_max_bound_[i][j] << ",";
+      std::cout << this->w_min_bound_[i][j] << ">,<";
+      std::cout << this->w_scale_up_[i][j] << ",";
+      std::cout << this->w_scale_down_[i][j] << ">,<";
       std::cout << w_slope_up_[i][j] << ",";
       std::cout << w_slope_down_[i][j] << ">]";
       std::cout.precision(10);
-      std::cout << this->sup_[i][j].decay_scale << ", ";
+      std::cout << this->w_decay_scale_[i][j] << ", ";
       std::cout.precision(6);
       std::cout << this->w_diffusion_rate_[i][j] << ", ";
       std::cout << this->w_reset_bias_[i][j];
diff --git a/src/rpucuda/rpu_powstep_device.cpp b/src/rpucuda/rpu_powstep_device.cpp
index ce880357..3a1bb6fd 100644
--- a/src/rpucuda/rpu_powstep_device.cpp
+++ b/src/rpucuda/rpu_powstep_device.cpp
@@ -70,14 +70,14 @@ template <typename T> void PowStepRPUDevice<T>::printDP(int x_count, int d_count
     for (int j = 0; j < x_count; ++j) {
       std::cout.precision(5);
       std::cout << i << "," << j << ": ";
-      std::cout << "[<" << this->sup_[i][j].max_bound << ",";
-      std::cout << this->sup_[i][j].min_bound << ">,<";
-      std::cout << this->sup_[i][j].scale_up << ",";
-      std::cout << this->sup_[i][j].scale_down << ">,<";
+      std::cout << "[<" << this->w_max_bound_[i][j] << ",";
+      std::cout << this->w_min_bound_[i][j] << ">,<";
+      std::cout << this->w_scale_up_[i][j] << ",";
+      std::cout << this->w_scale_down_[i][j] << ">,<";
       std::cout << w_gamma_up_[i][j] << ",";
       std::cout << w_gamma_down_[i][j] << ">]";
       std::cout.precision(10);
-      std::cout << this->sup_[i][j].decay_scale << ", ";
+      std::cout << this->w_decay_scale_[i][j] << ", ";
       std::cout.precision(6);
       std::cout << this->w_diffusion_rate_[i][j] << ", ";
       std::cout << this->w_reset_bias_[i][j];
diff --git a/src/rpucuda/rpu_pulsed_device.cpp b/src/rpucuda/rpu_pulsed_device.cpp
index 4d6d38c8..7cd3eafa 100644
--- a/src/rpucuda/rpu_pulsed_device.cpp
+++ b/src/rpucuda/rpu_pulsed_device.cpp
@@ -91,13 +91,14 @@ template <typename T> void PulsedRPUDevice<T>::allocateContainers() {
   int d_sz = this->d_size_;
   int x_sz = this->x_size_;
 
-  sup_ = Array_2D_Get<PulsedDPStruc<T>>(d_sz, x_sz);
-
   w_max_bound_ = Array_2D_Get<T>(d_sz, x_sz);
   w_min_bound_ = Array_2D_Get<T>(d_sz, x_sz);
+
   w_scale_up_ = Array_2D_Get<T>(d_sz, x_sz);
   w_scale_down_ = Array_2D_Get<T>(d_sz, x_sz);
+
   w_decay_scale_ = Array_2D_Get<T>(d_sz, x_sz);
+
   w_diffusion_rate_ = Array_2D_Get<T>(d_sz, x_sz);
   w_reset_bias_ = Array_2D_Get<T>(d_sz, x_sz);
   w_persistent_ = Array_2D_Get<T>(d_sz, x_sz);
@@ -113,14 +114,6 @@ template <typename T> void PulsedRPUDevice<T>::allocateContainers() {
       w_diffusion_rate_[i][j] = (T)0.0;
       w_reset_bias_[i][j] = (T)0.0;
       w_persistent_[i][j] = (T)0.0;
-
-      sup_[i][j].max_bound = std::numeric_limits<T>::max();
-      sup_[i][j].min_bound = std::numeric_limits<T>::min();
-      sup_[i][j].scale_up = (T)0.0;
-      sup_[i][j].scale_down = (T)0.0;
-      sup_[i][j].diffusion_rate = (T)0.0;
-      sup_[i][j].decay_scale = (T)1.0;
-      sup_[i][j].reset_bias = (T)0.0;
     }
   }
   containers_allocated_ = true;
@@ -130,8 +123,6 @@ template <typename T> void PulsedRPUDevice<T>::freeContainers() {
 
   if (containers_allocated_) {
 
-    Array_2D_Free<PulsedDPStruc<T>>(sup_);
-
     Array_2D_Free<T>(w_max_bound_);
     Array_2D_Free<T>(w_min_bound_);
     Array_2D_Free<T>(w_scale_up_);
@@ -172,7 +163,6 @@ PulsedRPUDevice<T>::PulsedRPUDevice(const PulsedRPUDevice<T> &other)
 
   for (int j = 0; j < this->x_size_; ++j) {
     for (int i = 0; i < this->d_size_; ++i) {
-      sup_[i][j] = other.sup_[i][j];
 
       w_scale_up_[i][j] = other.w_scale_up_[i][j];
       w_scale_down_[i][j] = other.w_scale_down_[i][j];
@@ -184,6 +174,8 @@ PulsedRPUDevice<T>::PulsedRPUDevice(const PulsedRPUDevice<T> &other)
       w_diffusion_rate_[i][j] = other.w_diffusion_rate_[i][j];
 
       w_reset_bias_[i][j] = other.w_reset_bias_[i][j];
+
+      w_persistent_[i][j] = other.w_persistent_[i][j];
     }
   }
 }
@@ -212,7 +204,6 @@ PulsedRPUDevice<T> &PulsedRPUDevice<T>::operator=(PulsedRPUDevice<T> &&other) {
   containers_allocated_ = other.containers_allocated_;
 
   // pointers
-  sup_ = other.sup_;
   w_scale_up_ = other.w_scale_up_;
   w_scale_down_ = other.w_scale_down_;
 
@@ -225,7 +216,6 @@ PulsedRPUDevice<T> &PulsedRPUDevice<T>::operator=(PulsedRPUDevice<T> &&other) {
   w_persistent_ = other.w_persistent_;
 
   // set pointers to null
-  other.sup_ = nullptr;
   other.w_scale_up_ = nullptr;
   other.w_scale_down_ = nullptr;
 
@@ -324,16 +314,6 @@ void PulsedRPUDevice<T>::setDeviceParameter(T **out_weights, const std::vector<T
     }
 
     dw_min += (fabs(w_scale_up_[0][i]) + fabs(w_scale_down_[0][i])) / (T)2.0;
-
-    // copy to sup
-    sup_[0][i].max_bound = w_max_bound_[0][i];
-    sup_[0][i].min_bound = w_min_bound_[0][i];
-    sup_[0][i].scale_up = w_scale_up_[0][i];
-    sup_[0][i].scale_down = w_scale_down_[0][i];
-
-    sup_[0][i].decay_scale = w_decay_scale_[0][i];
-    sup_[0][i].diffusion_rate = w_diffusion_rate_[0][i];
-    sup_[0][i].reset_bias = w_reset_bias_[0][i];
   }
 
   if (!getPar().legacy_params && this->hasWDrifter()) {
@@ -363,11 +343,12 @@ template <typename T> void PulsedRPUDevice<T>::decayWeights(T **weights, bool bi
   T *wd = w_decay_scale_[0];
   T *max_bound = w_max_bound_[0];
   T *min_bound = w_min_bound_[0];
+  T *b = w_reset_bias_[0];
 
   if (!bias_no_decay) {
     PRAGMA_SIMD
     for (int i = 0; i < this->size_; ++i) {
-      w[i] *= wd[i];
+      w[i] = (w[i] - b[i]) * wd[i] + b[i];
       w[i] = MIN(w[i], max_bound[i]);
       w[i] = MAX(w[i], min_bound[i]);
     }
@@ -375,7 +356,8 @@ template <typename T> void PulsedRPUDevice<T>::decayWeights(T **weights, bool bi
     const int last_col = this->x_size_ - 1; // x-major (ie row major)
     PRAGMA_SIMD
     for (int i = 0; i < this->size_; ++i) {
-      w[i] *= (i % this->x_size_ == last_col) ? (T)1.0 : wd[i];
+      T s = (i % this->x_size_ == last_col) ? (T)1.0 : wd[i];
+      w[i] = (w[i] - b[i]) * s + b[i];
       w[i] = MIN(w[i], max_bound[i]);
       w[i] = MAX(w[i], min_bound[i]);
     }
@@ -391,11 +373,13 @@ void PulsedRPUDevice<T>::decayWeights(T **weights, T alpha, bool bias_no_decay)
   T *wd = w_decay_scale_[0];
   T *max_bound = w_max_bound_[0];
   T *min_bound = w_min_bound_[0];
+  T *b = w_reset_bias_[0];
 
   if (!bias_no_decay) {
     PRAGMA_SIMD
     for (int i = 0; i < this->size_; ++i) {
-      w[i] *= 1 + alpha * (wd[i] - 1);
+      T s = 1 + alpha * (wd[i] - 1);
+      w[i] = (w[i] - b[i]) * s + b[i];
       w[i] = MIN(w[i], max_bound[i]);
       w[i] = MAX(w[i], min_bound[i]);
     }
@@ -403,7 +387,8 @@ void PulsedRPUDevice<T>::decayWeights(T **weights, T alpha, bool bias_no_decay)
     const int last_col = this->x_size_ - 1; // x-major (ie row major)
     PRAGMA_SIMD
     for (int i = 0; i < this->size_; ++i) {
-      w[i] *= (i % this->x_size_ == last_col) ? (T)1.0 : (1 + alpha * (wd[i] - 1));
+      T s = (i % this->x_size_ == last_col) ? (T)1.0 : (1 + alpha * (wd[i] - 1));
+      w[i] = (w[i] - b[i]) * s + b[i];
       w[i] = MIN(w[i], max_bound[i]);
       w[i] = MAX(w[i], min_bound[i]);
     }
@@ -518,24 +503,16 @@ void PulsedRPUDevice<T>::copyInvertDeviceParameter(const PulsedRPUDeviceBase<T>
     RPU_FATAL("Expect RPU Pulsed device");
   }
 
-  PulsedDPStruc<T> **sup = rpu->getDPStruc();
-
   for (int j = 0; j < this->x_size_; ++j) {
     for (int i = 0; i < this->d_size_; ++i) {
 
       // scaleup/down both have same sign
-      sup_[i][j].scale_up = sup[i][j].scale_down;
-      w_scale_up_[i][j] = sup[i][j].scale_down;
-
-      sup_[i][j].scale_down = sup[i][j].scale_up;
-      w_scale_down_[i][j] = sup[i][j].scale_up;
+      std::swap(w_scale_up_[i][j], w_scale_down_[i][j]);
 
       // min max have sign. mirror
-      sup_[i][j].min_bound = -sup[i][j].max_bound;
-      sup_[i][j].max_bound = -sup[i][j].min_bound;
-
-      w_max_bound_[i][j] = -sup[i][j].min_bound;
-      w_min_bound_[i][j] = -sup[i][j].max_bound;
+      T b = w_max_bound_[i][j];
+      w_max_bound_[i][j] = -w_min_bound_[i][j];
+      w_min_bound_[i][j] = -b;
     }
   }
 }
@@ -605,89 +582,75 @@ void PulsedRPUDevice<T>::populate(const PulsedRPUDeviceMetaParameter<T> &p, Real
   for (int j = 0; j < this->x_size_; ++j) {
     for (int i = 0; i < this->d_size_; ++i) {
 
-      sup_[i][j].max_bound = par.w_max * (1 + par.w_max_dtod * rng->sampleGauss());
-      sup_[i][j].min_bound = par.w_min * (1 + par.w_min_dtod * rng->sampleGauss());
+      w_max_bound_[i][j] = par.w_max * (1 + par.w_max_dtod * rng->sampleGauss());
+      w_min_bound_[i][j] = par.w_min * (1 + par.w_min_dtod * rng->sampleGauss());
 
       T gain = ((T)1.0 + gain_std * rng->sampleGauss());
       T r = up_down_std * rng->sampleGauss();
 
-      sup_[i][j].scale_up = (up_bias + gain + r) * par.dw_min; // to reduce mults in updates
-      sup_[i][j].scale_down = (down_bias + gain - r) * par.dw_min;
+      w_scale_up_[i][j] = (up_bias + gain + r) * par.dw_min; // to reduce mults in updates
+      w_scale_down_[i][j] = (down_bias + gain - r) * par.dw_min;
 
       // enforce consistency
       if (par.enforce_consistency) {
-        sup_[i][j].scale_up = fabs(sup_[i][j].scale_up);
-        sup_[i][j].scale_down = fabs(sup_[i][j].scale_down);
+        w_scale_up_[i][j] = fabs(w_scale_up_[i][j]);
+        w_scale_down_[i][j] = fabs(w_scale_down_[i][j]);
 
-        if (sup_[i][j].min_bound > sup_[i][j].max_bound) {
-          std::swap(sup_[i][j].min_bound, sup_[i][j].max_bound);
+        if (w_min_bound_[i][j] > w_max_bound_[i][j]) {
+          std::swap(w_min_bound_[i][j], w_max_bound_[i][j]);
         }
-        sup_[i][j].max_bound = fabs(sup_[i][j].max_bound);
-        sup_[i][j].min_bound = -fabs(sup_[i][j].min_bound);
+        w_max_bound_[i][j] = fabs(w_max_bound_[i][j]);
+        w_min_bound_[i][j] = -fabs(w_min_bound_[i][j]);
       } else {
         // "turn off" weight if max<min
-        if (sup_[i][j].min_bound > sup_[i][j].max_bound) {
-          T m = sup_[i][j].max_bound + (sup_[i][j].min_bound - sup_[i][j].max_bound) / ((T)2.0);
-          sup_[i][j].max_bound = m;
-          sup_[i][j].min_bound = m;
+        if (w_min_bound_[i][j] > w_max_bound_[i][j]) {
+          T m = w_max_bound_[i][j] + (w_min_bound_[i][j] - w_max_bound_[i][j]) / ((T)2.0);
+          w_max_bound_[i][j] = m;
+          w_min_bound_[i][j] = m;
         }
       }
 
       // corrupt devices
       if (par.corrupt_devices_prob > rng->sampleUniform()) {
         // stuck somewhere in min_max
-        T mn =
-            MAX(MIN(sup_[i][j].max_bound, sup_[i][j].min_bound), -fabs(par.corrupt_devices_range));
-        T mx =
-            MIN(MAX(sup_[i][j].max_bound, sup_[i][j].min_bound), fabs(par.corrupt_devices_range));
+        T mn = MAX(MIN(w_max_bound_[i][j], w_min_bound_[i][j]), -fabs(par.corrupt_devices_range));
+        T mx = MIN(MAX(w_max_bound_[i][j], w_min_bound_[i][j]), fabs(par.corrupt_devices_range));
 
         T value = mn + (mx - mn) * rng->sampleUniform();
-        sup_[i][j].max_bound = value;
-        sup_[i][j].min_bound = value;
-        sup_[i][j].scale_up = (T)0.0;
-        sup_[i][j].scale_down = (T)0.0;
+        w_max_bound_[i][j] = value;
+        w_min_bound_[i][j] = value;
+        w_scale_up_[i][j] = (T)0.0;
+        w_scale_down_[i][j] = (T)0.0;
       }
 
       // perfect bias
       if ((par.perfect_bias) && (j == this->x_size_ - 1)) {
-        sup_[i][j].scale_up = par.dw_min;
-        sup_[i][j].scale_down = par.dw_min;
-        sup_[i][j].min_bound = 100 * par.w_min; // essentially no bound
-        sup_[i][j].max_bound = 100 * par.w_max; // essentially no bound
+        w_scale_up_[i][j] = par.dw_min;
+        w_scale_down_[i][j] = par.dw_min;
+        w_min_bound_[i][j] = 100 * par.w_min; // essentially no bound
+        w_max_bound_[i][j] = 100 * par.w_max; // essentially no bound
       }
 
       //--------------------
       // diffusion
       {
         T t = fabs(par.diffusion * (1 + par.diffusion_dtod * rng->sampleGauss()));
-        sup_[i][j].diffusion_rate = t;
+        w_diffusion_rate_[i][j] = t;
       }
 
       //--------------------
       // reset
       { // additive dtod
         T t = par.reset + par.reset_dtod * rng->sampleGauss();
-        sup_[i][j].reset_bias = t;
+        w_reset_bias_[i][j] = t;
       }
 
       //--------------------
       // decay
       {
         T t = par.lifetime * ((T)1.0 + par.lifetime_dtod * rng->sampleGauss());
-        sup_[i][j].decay_scale = t > 1.0 ? (T)((T)1. - ((T)1. / t)) : (T)0.0;
+        w_decay_scale_[i][j] = t > 1.0 ? (T)((T)1. - ((T)1. / t)) : (T)0.0;
       }
-
-      //--------------------
-      // copy for for linear access
-      w_max_bound_[i][j] = sup_[i][j].max_bound;
-      w_min_bound_[i][j] = sup_[i][j].min_bound;
-      w_scale_up_[i][j] = sup_[i][j].scale_up;
-      w_scale_down_[i][j] = sup_[i][j].scale_down;
-
-      w_decay_scale_[i][j] = sup_[i][j].decay_scale;
-      w_diffusion_rate_[i][j] = sup_[i][j].diffusion_rate;
-
-      w_reset_bias_[i][j] = sup_[i][j].reset_bias;
     }
   }
 }
@@ -707,12 +670,12 @@ template <typename T> void PulsedRPUDevice<T>::printDP(int x_count, int d_count)
 
   for (int i = 0; i < d_count1; ++i) {
     for (int j = 0; j < x_count1; ++j) {
-      std::cout << "[<" << sup_[i][j].max_bound << ", ";
-      std::cout << sup_[i][j].min_bound << ">, <";
-      std::cout << sup_[i][j].scale_up << ", ";
-      std::cout << sup_[i][j].scale_down << "> ";
+      std::cout << "[<" << w_max_bound_[i][j] << ", ";
+      std::cout << w_min_bound_[i][j] << ">, <";
+      std::cout << w_scale_up_[i][j] << ", ";
+      std::cout << w_scale_down_[i][j] << "> ";
       std::cout.precision(10);
-      std::cout << sup_[i][j].decay_scale << ", ";
+      std::cout << w_decay_scale_[i][j] << ", ";
       std::cout.precision(6);
       std::cout << w_diffusion_rate_[i][j] << ", ";
       std::cout << w_reset_bias_[i][j];
diff --git a/src/rpucuda/rpu_pulsed_device.h b/src/rpucuda/rpu_pulsed_device.h
index 09d33636..0884e5f2 100644
--- a/src/rpucuda/rpu_pulsed_device.h
+++ b/src/rpucuda/rpu_pulsed_device.h
@@ -101,7 +101,6 @@ template <typename T> struct PulsedRPUDeviceMetaParameter : PulsedRPUDeviceMetaP
   virtual bool implementsWriteNoise() const { return false; }; // needs to be activated in derived
   inline bool usesPersistentWeight() const { return write_noise_std > 0; };
   inline T getScaledWriteNoise() const { return write_noise_std * this->dw_min; };
-  inline bool needsResetBias() const { return reset > 0 || reset_dtod > 0; };
 
   void initialize() override {
     PulsedRPUDeviceMetaParameterBase<T>::initialize();
@@ -110,7 +109,6 @@ template <typename T> struct PulsedRPUDeviceMetaParameter : PulsedRPUDeviceMetaP
     }
     reset_dtod = MAX(reset_dtod, (T)0.0);
     reset_std = MAX(reset_std, (T)0.0);
-    reset = MAX(reset, (T)0.0);
   };
 };
 
@@ -176,16 +174,6 @@ template <typename T> class PulsedRPUDeviceBase : public SimpleRPUDevice<T> {
   T weight_granularity_ = 0.0;
 };
 
-template <typename T> struct PulsedDPStruc {
-  T max_bound = (T)0.0;
-  T min_bound = (T)0.0;
-  T scale_up = (T)0.0;
-  T scale_down = (T)0.0;
-  T decay_scale = (T)0.0;
-  T diffusion_rate = (T)0.0;
-  T reset_bias = (T)0.0;
-};
-
 template <typename T> class PulsedRPUDevice : public PulsedRPUDeviceBase<T> {
 
 public:
@@ -215,7 +203,6 @@ template <typename T> class PulsedRPUDevice : public PulsedRPUDeviceBase<T> {
     swap(a.w_diffusion_rate_, b.w_diffusion_rate_);
     swap(a.w_persistent_, b.w_persistent_);
     swap(a.w_reset_bias_, b.w_reset_bias_);
-    swap(a.sup_, b.sup_);
 
     swap(a.containers_allocated_, b.containers_allocated_);
   }
@@ -235,7 +222,6 @@ template <typename T> class PulsedRPUDevice : public PulsedRPUDeviceBase<T> {
   inline T **getResetBias() const { return w_reset_bias_; };
   inline T **getScaleUp() const { return w_scale_up_; };
   inline T **getScaleDown() const { return w_scale_down_; };
-  inline PulsedDPStruc<T> **getDPStruc() const { return sup_; };
   PulsedRPUDeviceMetaParameter<T> &getPar() const override {
     return static_cast<PulsedRPUDeviceMetaParameter<T> &>(SimpleRPUDevice<T>::getPar());
   };
@@ -259,8 +245,6 @@ template <typename T> class PulsedRPUDevice : public PulsedRPUDeviceBase<T> {
 protected:
   void populate(const PulsedRPUDeviceMetaParameter<T> &par, RealWorldRNG<T> *rng);
 
-  PulsedDPStruc<T> **sup_ = nullptr;
-
   T **w_max_bound_ = nullptr;
   T **w_min_bound_ = nullptr;
   T **w_scale_up_ = nullptr;


Test patch:
diff --git a/src/rpucuda/cuda/rpucuda_pulsed_device_test.cpp b/src/rpucuda/cuda/rpucuda_pulsed_device_test.cpp
index 55fe3285..e90943e9 100644
--- a/src/rpucuda/cuda/rpucuda_pulsed_device_test.cpp
+++ b/src/rpucuda/cuda/rpucuda_pulsed_device_test.cpp
@@ -17,8 +17,8 @@
 #include "rpu_pulsed.h"
 #include "rpucuda_constantstep_device.h"
 #include "rpucuda_expstep_device.h"
-#include "rpucuda_hidden_device.h"
 #include "rpucuda_linearstep_device.h"
+#include "rpucuda_powstep_device.h"
 #include "rpucuda_pulsed.h"
 #include "rpucuda_pulsed_device.h"
 #include "utility_functions.h"
@@ -245,10 +245,11 @@ template <typename DeviceParT> class RPUDeviceTestFixture : public ::testing::Te
 
 // types
 typedef ::testing::Types<
-    HiddenStepRPUDeviceMetaParameter<num_t>,
     LinearStepRPUDeviceMetaParameter<num_t>,
     ExpStepRPUDeviceMetaParameter<num_t>,
+    PowStepRPUDeviceMetaParameter<num_t>,
     ConstantStepRPUDeviceMetaParameter<num_t>>
+
     MetaPar;
 
 TYPED_TEST_CASE(RPUDeviceTestFixture, MetaPar);
diff --git a/tests/test_specific_tiles.py b/tests/test_specific_tiles.py
index 4bd6f29a..c81866ec 100644
--- a/tests/test_specific_tiles.py
+++ b/tests/test_specific_tiles.py
@@ -82,7 +82,6 @@ def test_hidden_parameter_setting(self):
 
             model.analog_tile.set_hidden_parameters(params)
 
-            print(model.analog_tile.tile)
             weight, bias = model.get_weights()
 
             # should be
@@ -99,9 +98,11 @@ def test_decay(self):
 
         lifetime = 100.  # initial setting (needs to be larger 1)
         gamma = 0.1
+        reset_bias = 0.1  # decay shift
         rpu_config = self.get_transfer_compound(gamma=gamma,
                                                 lifetime=lifetime,
-                                                lifetime_dtod=0.0)
+                                                lifetime_dtod=0.0,
+                                                reset=reset_bias)
 
         model = self.get_layer(in_features=2, out_features=1, rpu_config=rpu_config)
 
@@ -145,10 +146,10 @@ def test_decay(self):
         weight, bias = model.get_weights()
 
         # reference values
-        a = a * pow(a_dcy, epochs)
-        b = b * pow(b_dcy, epochs)
-        c = c * pow(c_dcy, epochs)
-        d = d * pow(d_dcy, epochs)
+        a = (a - reset_bias) * pow(a_dcy, epochs) + reset_bias
+        b = (b - reset_bias) * pow(b_dcy, epochs) + reset_bias
+        c = (c - reset_bias) * pow(c_dcy, epochs) + reset_bias
+        d = (d - reset_bias) * pow(d_dcy, epochs) + reset_bias
 
         if self.digital_bias:
             self.assertEqual(bias[0], 0.0)