Repository: IBM/aihwkit (language: python)

Issue:
Segmentation fault (core dumped)
## Description
Segmentation fault (core dumped) during rebuilding a neural network model.
This is sample code:
```python
def create_model(k):
    model = AnalogSequential(AnalogLinear(784, k), nn.Sigmoid, AnalogLinear(k,10)).cuda()
return model

x = [100,200,300]
for i in x:
    ...
    model = create_model(i)
    ...
```
When I run the above code, after the first iteration, I got the message, "segmentation fault (core dumped)".
If I replace `AnalogSequential` with `nn.Sequential` and `AnalogLinear` with `nn.Linear`, it works. So, I think PyTorch and cuda are not the cause of this issue. In addition, if I run without `cuda()`, it also works.
 
<!-- A clear and concise description of what the bug is. -->

## How to reproduce

<!-- An step by step description of how to reproduce the bug. -->

## Expected behavior

<!-- A clear and concise description of what you expected to happen. -->

## Other information

<!-- Feel free to add more information that can help solving the bug. -->

* Pytorch version: 1.7.1
* Package version: git clone
* OS: Ubuntu 18.04.2
* Python version: 3.8
* Conda version (or N/A): 
* NVIDIA Graphic Driver: 450.51.05
* CUDA Version: 11.0
* GPU: RTX 6000 



Discussion / hints:
Thanks for the detailed report @chaeunl - I'm trying to replicate the issue, as we recently had a segmentation fault appearing in #112 that might have some relation.

Can you confirm the issue can be reproduced constantly (ie. is it always happening consistently at the first  iteration of the loop)? 

Edit: so far I have _not_ been able to reproduce it with:
```python
from torch import nn
from aihwkit.nn import AnalogSequential, AnalogLinear

def create_model(k):
    model = AnalogSequential(
        AnalogLinear(784, k),
        nn.Sigmoid(),
        AnalogLinear(k, 10)).cuda()
    return model

x = [100, 200, 300]
for i in x:
    model = create_model(i)
```
under PyTorch 1.7.1, Python 3.6 / 3.8 / 3.9.
@diego-plan9 , as I simplify my code, I didn't notice that the above code has no problem.
Instead, I found where the error comes from. I think it's about declaring the `TransferCompound` device.

This is a simplified code and I confirmed it raises error:
```python
from torch import nn
from aihwkit.nn import AnalogSequential, AnalogLinear
from aihwkit.simulator.configs import UnitCellRPUConfig
from aihwkit.simulator.configs.devices import LinearStepDevice, TransferCompound

def create_rpu(which_device="LinearStepDevice"):
    if which_device == "LinearStepDevice":
        rpu = UnitCellRPUConfig(device=LinearStepDevice())
    elif which_device == "TransferCompound":
        rpu = UnitCellRPUConfig(device=TransferCompound(unit_cell_devices=[LinearStepDevice(), LinearStepDevice()]))
    else:
        raise ValueError("Undefined Device")
    return rpu

def create_model(k, which_device="LinearStepDevice"):
    model = AnalogSequential(AnalogLinear(784, k, rpu_config=create_rpu(which_device=which_device)),nn.Sigmoid(),AnalogLinear(k, 10, rpu_config=create_rpu(which_device=which_device))).cuda()
    return model

x = [100, 200, 300]
device = "TransferCompound"
for i in x:
    model = create_model(i, which_device=device)
    print("model built")
```
If I change `device` to `"LinearStepDevice"`, it does not raise seg fault. Also, even though using `TransferCompound`, if I build a network without `cuda()` , it works.

*ps: the above code does not give any problem at the time I installed aihwkit through git clone on Nov. last year (it's the first time to install it while we contact through e-mail.). At that time, I use PyTorch 1.8.0 and cuda 11.2 (because of GPU model; RTX 3090). As of now, I am trying with another sever. This one has RTX 6000 and PyTorch 1.7.1 under cuda 11.0 is installed. So, the difference is that we installed aihwkit at different time and use different PyTorch version. And so, while installing, we use different command to specify the cuda architecture: `-DRPU_CUDA_ARCHITECTURES`. I hope the information is helpful to you :)  
Thanks @chaeunl as usual for all the detailed info: I still could not reproduce it but still have some options and configurations to try (and having the exact code and all the pointers help) - will update the issue after exhausting them.
Thank you, @diego-plan9 .
I am also trying the above code in another environment (CUDA/PyTorch are differed from this case).
Hi @chaeunl,
Indeed, I now found an issue that produced a segfault for the CUDA vector device (the base class of the transfer compound). I am not completely sure whether it is the same issue that you observed, but it occurred when only one of mutliple sub devices of the vector device had diffusion or reset defined.  It is fixed in the upcoming releases. Hopefully, that resolves your issue as well.   
 
Thank you so much, @maljoras
Is it merged into github source as of now?  As I clone the code from github directly, I can confirm it now :)

Fix patch:
diff --git a/src/aihwkit/simulator/rpu_base_src/rpu_base_tiles.cpp b/src/aihwkit/simulator/rpu_base_src/rpu_base_tiles.cpp
index b648a519..f343b0c3 100644
--- a/src/aihwkit/simulator/rpu_base_src/rpu_base_tiles.cpp
+++ b/src/aihwkit/simulator/rpu_base_src/rpu_base_tiles.cpp
@@ -679,7 +679,7 @@ void declare_rpu_tiles(py::module &m) {
             CHECK_TORCH_INPUT(d_tensor);
 
             int N = x_input.size(0); // batch
-            int d_image_size = ((d_tensor.numel()/d_tensor.size(0))/d_tensor.size(1));
+            int d_image_size = ((d_tensor.numel() / d_tensor.size(0)) / d_tensor.size(1));
 
             // Call RPU function.
             std::lock_guard<std::mutex> lock(self.mutex_);
@@ -714,7 +714,7 @@ void declare_rpu_tiles(py::module &m) {
             CHECK_TORCH_INPUT(x_tensor);
 
             int N = d_input.size(0); // batch
-            int d_image_size = ((d_input.numel()/d_input.size(0))/d_input.size(1));
+            int d_image_size = ((d_input.numel() / d_input.size(0)) / d_input.size(1));
 
             // Call RPU function.
             std::lock_guard<std::mutex> lock(self.mutex_);
@@ -748,7 +748,7 @@ void declare_rpu_tiles(py::module &m) {
             CHECK_TORCH_INPUT(d_input);
 
             int N = d_input.size(0); // batch
-            int d_image_size = d_input.numel()/(d_input.size(0)*d_input.size(1));
+            int d_image_size = d_input.numel() / (d_input.size(0) * d_input.size(1));
 
             // Call RPU function.
             std::lock_guard<std::mutex> lock(self.mutex_);
diff --git a/src/aihwkit/simulator/rpu_base_src/rpu_base_tiles_cuda.cpp b/src/aihwkit/simulator/rpu_base_src/rpu_base_tiles_cuda.cpp
index cd9fa8ac..2fa9fc83 100644
--- a/src/aihwkit/simulator/rpu_base_src/rpu_base_tiles_cuda.cpp
+++ b/src/aihwkit/simulator/rpu_base_src/rpu_base_tiles_cuda.cpp
@@ -256,7 +256,7 @@ void declare_rpu_tiles_cuda(py::module &m) {
             CHECK_TORCH_CUDA_INPUT(d_tensor);
 
             int N = x_input.size(0); // batch
-            int d_image_size = ((d_tensor.numel()/d_tensor.size(0))/d_tensor.size(1));
+            int d_image_size = ((d_tensor.numel() / d_tensor.size(0)) / d_tensor.size(1));
 
             // Call RPU function.
             self.finishUpdateCalculations();
@@ -267,7 +267,7 @@ void declare_rpu_tiles_cuda(py::module &m) {
                 d_image_size, N, true, is_test);
             return d_tensor;
           },
-           py::arg("x_input"), py::arg("d_tensor"), py::arg("is_test") = false,
+          py::arg("x_input"), py::arg("d_tensor"), py::arg("is_test") = false,
           R"pbdoc(
            Compute the dot product using an index matrix (forward pass).
 
@@ -292,7 +292,7 @@ void declare_rpu_tiles_cuda(py::module &m) {
             CHECK_TORCH_CUDA_INPUT(x_tensor);
 
             int N = d_input.size(0); // batch
-            int d_image_size = ((d_input.numel()/d_input.size(0))/d_input.size(1));
+            int d_image_size = ((d_input.numel() / d_input.size(0)) / d_input.size(1));
 
             // Call RPU function.
             self.finishUpdateCalculations();
@@ -328,7 +328,7 @@ void declare_rpu_tiles_cuda(py::module &m) {
             CHECK_TORCH_CUDA_INPUT(d_input);
 
             int N = d_input.size(0); // batch
-            int d_image_size = d_input.numel()/(d_input.size(0)*d_input.size(1));
+            int d_image_size = d_input.numel() / (d_input.size(0) * d_input.size(1));
 
             // Call RPU function.
             self.finishUpdateCalculations();
diff --git a/src/rpucuda/cuda/rpucuda_difference_device.cu b/src/rpucuda/cuda/rpucuda_difference_device.cu
index 7ce9cbbf..90414cc8 100644
--- a/src/rpucuda/cuda/rpucuda_difference_device.cu
+++ b/src/rpucuda/cuda/rpucuda_difference_device.cu
@@ -100,9 +100,7 @@ void DifferenceRPUDeviceCuda<T>::resetCols(
   VectorRPUDeviceCuda<T>::resetCols(dev_weights, start_col, n_cols, reset_prob);
 }
 
-template <typename T> bool DifferenceRPUDeviceCuda<T>::isInverted() const {
-  return g_plus_ == 0;
-}
+template <typename T> bool DifferenceRPUDeviceCuda<T>::isInverted() const { return g_plus_ == 0; }
 
 template <typename T> void DifferenceRPUDeviceCuda<T>::invert() {
   std::swap(g_plus_, g_minus_);
diff --git a/src/rpucuda/cuda/rpucuda_expstep_device.cu b/src/rpucuda/cuda/rpucuda_expstep_device.cu
index e111c328..5e3420a2 100644
--- a/src/rpucuda/cuda/rpucuda_expstep_device.cu
+++ b/src/rpucuda/cuda/rpucuda_expstep_device.cu
@@ -36,12 +36,12 @@ namespace RPU {
 template <typename T> struct UpdateFunctorExpStep {
 
   __device__ __forceinline__ void operator()(
-      T &w,
+      T &apparent_weight,
       uint32_t n,
       uint32_t negative,
       const float4 par_4,
       const float2 par_2,
-      T &par_1,
+      T &persistent_weight,
       const T *global_pars,
       T noise_std_dw,
       curandState &local_state)
@@ -49,10 +49,13 @@ template <typename T> struct UpdateFunctorExpStep {
   {
     // par_4 order (min_bound, scale_down, max_bound, scale_up )
     // global_pars see below
+    T uw_std = global_pars[6];
     T wmax = par_4.z; //[2];
     T wmin = par_4.x; //[0];
     T b_diff = (wmax - wmin);
 
+    T &w = uw_std > 0 ? persistent_weight : apparent_weight;
+
     if (b_diff > 0) { // only do something when bounds make sense
 
       T A = negative ? global_pars[1] : global_pars[0];        // 1: up, 0: down
@@ -69,12 +72,17 @@ template <typename T> struct UpdateFunctorExpStep {
           UPDATE_ONCE;
         }
       }
+      // add update write noise onto apparent weight
+      if (uw_std > 0) {
+        T stoch_value = curand_normal(&local_state);
+        apparent_weight = persistent_weight + uw_std * stoch_value;
+      }
     }
   }
 };
 #undef UPDATE_ONCE
 
-RPUCUDA_DEVICE_ADD_FUNCTOR_UPDATE_KERNELS(ExpStep, UpdateFunctorExpStep<T>, 6);
+RPUCUDA_DEVICE_ADD_FUNCTOR_UPDATE_KERNELS(ExpStep, UpdateFunctorExpStep<T>, 7);
 
 template class ExpStepRPUDeviceCuda<float>;
 #ifdef RPU_USE_DOUBLE
diff --git a/src/rpucuda/cuda/rpucuda_expstep_device.h b/src/rpucuda/cuda/rpucuda_expstep_device.h
index 69cdc086..5336b47d 100644
--- a/src/rpucuda/cuda/rpucuda_expstep_device.h
+++ b/src/rpucuda/cuda/rpucuda_expstep_device.h
@@ -26,7 +26,7 @@ template <typename T> class ExpStepRPUDeviceCuda : public PulsedRPUDeviceCuda<T>
       ExpStepRPUDeviceCuda,
       ExpStepRPUDevice,
       /*ctor body*/
-      dev_es_par_ = std::unique_ptr<CudaArray<T>>(new CudaArray<T>(this->context_, 6));
+      dev_es_par_ = std::unique_ptr<CudaArray<T>>(new CudaArray<T>(this->context_, 7));
       ,
       /*dtor body*/
       ,
@@ -40,7 +40,7 @@ template <typename T> class ExpStepRPUDeviceCuda : public PulsedRPUDeviceCuda<T>
       swap(a.dev_es_par_, b.dev_es_par_);
       ,
       /*host copy from cpu (rpu_device). Parent device params are copyied automatically*/
-      T es_par_arr[6];
+      T es_par_arr[7];
       auto &par = getPar();
       es_par_arr[0] = par.es_A_down;
       es_par_arr[1] = par.es_A_up;
@@ -48,7 +48,7 @@ template <typename T> class ExpStepRPUDeviceCuda : public PulsedRPUDeviceCuda<T>
       es_par_arr[3] = par.es_gamma_up;
       es_par_arr[4] = par.es_a;
       es_par_arr[5] = par.es_b;
-
+      es_par_arr[6] = getPar().getScaledWriteNoise();
       dev_es_par_->assign(es_par_arr);
       this->context_->synchronize();)
 
@@ -59,6 +59,9 @@ template <typename T> class ExpStepRPUDeviceCuda : public PulsedRPUDeviceCuda<T>
       int use_bo64,
       bool out_trans,
       const PulsedUpdateMetaParameter<T> &up) override;
+  T *get1ParamsData() override {
+    return getPar().usesPersistentWeight() ? this->dev_persistent_weights_->getData() : nullptr;
+  };
 
 private:
   std::unique_ptr<CudaArray<T>> dev_es_par_ = nullptr;
diff --git a/src/rpucuda/cuda/rpucuda_linearstep_device.cu b/src/rpucuda/cuda/rpucuda_linearstep_device.cu
index 0e94b706..c5bf72fc 100644
--- a/src/rpucuda/cuda/rpucuda_linearstep_device.cu
+++ b/src/rpucuda/cuda/rpucuda_linearstep_device.cu
@@ -20,18 +20,21 @@ namespace RPU {
 template <typename T> struct UpdateFunctorLinearStepMult {
 
   __device__ __forceinline__ void operator()(
-      T &w,
+      T &apparent_weight,
       uint32_t n,
       uint32_t negative,
       const float4 par_4,
       const float2 lin_slope,
-      T &par_1,
-      const T *global_par,
+      T &persistent_weight,
+      const T *write_noise_std,
       T noise_std_dw,
       curandState &local_state) {
+    T uw_std = *write_noise_std;
     T lin_dw = (negative > 0) ? (par_4.w) : (-par_4.y);        //[3], [1]
     T lin_a = (negative > 0) ? (lin_slope.y) : (-lin_slope.x); // [1],[0]
 
+    T &w = uw_std > 0 ? persistent_weight : apparent_weight;
+
     // n is larger 0 in any case
     if (n == 1) {
       if (noise_std_dw > 0) {
@@ -59,24 +62,33 @@ template <typename T> struct UpdateFunctorLinearStepMult {
     w = (w > wmax) ? wmax : w;
     T wmin = par_4.x; // [0]
     w = (w < wmin) ? wmin : w;
+
+    // add update write noise onto apparent weight
+    if (uw_std > 0) {
+      T stoch_value = curand_normal(&local_state);
+      apparent_weight = persistent_weight + uw_std * stoch_value;
+    }
   }
 };
 
 template <typename T> struct UpdateFunctorLinearStepAdd {
 
   __device__ __forceinline__ void operator()(
-      T &w,
+      T &apparent_weight,
       uint32_t n,
       uint32_t negative,
       const float4 par_4,
       const float2 lin_slope,
-      T &par_1,
-      const T *global_par,
+      T &persistent_weight,
+      const T *write_noise_std,
       T noise_std_dw,
       curandState &local_state) {
+    T uw_std = *write_noise_std;
     T lin_dw = (negative > 0) ? (par_4.w) : (-par_4.y);        // [3] [1]
     T lin_a = (negative > 0) ? (lin_slope.y) : (-lin_slope.x); //[1],[0]
 
+    T &w = uw_std > 0 ? persistent_weight : apparent_weight;
+
     // n is larger 0 in any case
     if (n == 1) {
       if (noise_std_dw > 0) {
@@ -103,6 +115,12 @@ template <typename T> struct UpdateFunctorLinearStepAdd {
     w = (w > wmax) ? wmax : w;
     T wmin = par_4.x; // [0]
     w = (w < wmin) ? wmin : w;
+
+    // add update write noise onto apparent weight
+    if (uw_std > 0) {
+      T stoch_value = curand_normal(&local_state);
+      apparent_weight = persistent_weight + uw_std * stoch_value;
+    }
   }
 };
 
diff --git a/src/rpucuda/cuda/rpucuda_linearstep_device.h b/src/rpucuda/cuda/rpucuda_linearstep_device.h
index 072f0c2b..2b7f6034 100644
--- a/src/rpucuda/cuda/rpucuda_linearstep_device.h
+++ b/src/rpucuda/cuda/rpucuda_linearstep_device.h
@@ -27,17 +27,21 @@ template <typename T> class LinearStepRPUDeviceCuda : public PulsedRPUDeviceCuda
       LinearStepRPUDevice,
       /*ctor body*/
       dev_slope_ = RPU::make_unique<CudaArray<float>>(this->context_, 2 * this->size_);
+      dev_write_noise_std_ = RPU::make_unique<CudaArray<T>>(this->context_, 1);
       ,
       /*dtor body*/
       ,
       /*copy body*/
       dev_slope_->assign(*other.dev_slope_);
+      dev_write_noise_std_->assign(*other.dev_write_noise_std_);
       ,
       /*move assigment body*/
       dev_slope_ = std::move(other.dev_slope_);
+      dev_write_noise_std_ = std::move(other.dev_write_noise_std_);
       ,
       /*swap body*/
       swap(a.dev_slope_, b.dev_slope_);
+      swap(a.dev_write_noise_std_, b.dev_write_noise_std_);
       ,
       /*host copy from cpu (rpu_device). Parent device params are copyied automatically*/
       int d_size = this->d_size_;
@@ -48,12 +52,13 @@ template <typename T> class LinearStepRPUDeviceCuda : public PulsedRPUDeviceCuda
 
       for (int i = 0; i < d_size; ++i) {
         for (int j = 0; j < x_size; ++j) {
-
           int kk = j * (d_size * 2) + 2 * i;
           tmp_slope[kk] = w_slope_down[i][j];
           tmp_slope[kk + 1] = w_slope_up[i][j];
         }
       } dev_slope_->assign(tmp_slope);
+      dev_write_noise_std_->setConst(getPar().getScaledWriteNoise());
+
       this->context_->synchronize();
       delete[] tmp_slope;);
 
@@ -63,10 +68,15 @@ template <typename T> class LinearStepRPUDeviceCuda : public PulsedRPUDeviceCuda
       int use_bo64,
       bool out_trans,
       const PulsedUpdateMetaParameter<T> &up) override;
+  T *getGlobalParamsData() override { return dev_write_noise_std_->getData(); };
   float *get2ParamsData() override { return dev_slope_->getData(); };
+  T *get1ParamsData() override {
+    return getPar().usesPersistentWeight() ? this->dev_persistent_weights_->getData() : nullptr;
+  };
 
 private:
   std::unique_ptr<CudaArray<float>> dev_slope_ = nullptr;
+  std::unique_ptr<CudaArray<T>> dev_write_noise_std_ = nullptr;
 };
 
 } // namespace RPU
diff --git a/src/rpucuda/cuda/rpucuda_pulsed_device.cu b/src/rpucuda/cuda/rpucuda_pulsed_device.cu
index be040031..f644e83f 100644
--- a/src/rpucuda/cuda/rpucuda_pulsed_device.cu
+++ b/src/rpucuda/cuda/rpucuda_pulsed_device.cu
@@ -32,9 +32,10 @@ PulsedRPUDeviceCuda<T>::PulsedRPUDeviceCuda(CudaContext *c, int x_size, int d_si
 template <typename T> void PulsedRPUDeviceCuda<T>::initialize() {
 
   dev_4params_ = RPU::make_unique<CudaArray<float>>(this->context_, 4 * this->size_);
-  dev_diffusion_rate_ = RPU::make_unique<CudaArray<T>>(this->context_, this->size_);
-  dev_reset_bias_ = RPU::make_unique<CudaArray<T>>(this->context_, this->size_);
   dev_decay_scale_ = RPU::make_unique<CudaArray<T>>(this->context_, this->size_);
+  dev_diffusion_rate_ = nullptr; // on the fly
+  dev_reset_bias_ = nullptr;
+  dev_persistent_weights_ = nullptr;
 
   this->context_->synchronize();
 };
@@ -45,10 +46,21 @@ PulsedRPUDeviceCuda<T>::PulsedRPUDeviceCuda(const PulsedRPUDeviceCuda<T> &other)
   initialize();
 
   dev_4params_->assign(*other.dev_4params_);
-  dev_diffusion_rate_->assign(*other.dev_diffusion_rate_);
-  dev_reset_bias_->assign(*other.dev_reset_bias_);
   dev_decay_scale_->assign(*other.dev_decay_scale_);
 
+  if (other.dev_diffusion_rate_ != nullptr) {
+    dev_diffusion_rate_ = RPU::make_unique<CudaArray<T>>(this->context_, this->size_);
+    dev_diffusion_rate_->assign(*other.dev_diffusion_rate_);
+  }
+  if (other.dev_reset_bias_ != nullptr) {
+    dev_reset_bias_ = RPU::make_unique<CudaArray<T>>(this->context_, this->size_);
+    dev_reset_bias_->assign(*other.dev_reset_bias_);
+  }
+  if (other.dev_persistent_weights_ != nullptr) {
+    dev_persistent_weights_ = RPU::make_unique<CudaArray<T>>(this->context_, this->size_);
+    dev_persistent_weights_->assign(*other.dev_persistent_weights_);
+  }
+
   this->context_->synchronize();
 };
 
@@ -73,6 +85,7 @@ PulsedRPUDeviceCuda<T>::PulsedRPUDeviceCuda(const PulsedRPUDeviceCuda<T> &other)
 //   dev_diffusion_rate_ = std::move(other.dev_diffusion_rate_);
 //   dev_reset_bias_ = std::move(other.dev_reset_bias_);
 //   dev_decay_scale_ = std::move(other.dev_decay_scale_);
+//   dev_persistent_weights_ = std::move(other.dev_persistent_weights_);
 
 //   return *this;
 // };
@@ -100,9 +113,12 @@ void PulsedRPUDeviceCuda<T>::populateFrom(const AbstractRPUDevice<T> &rpu_device
   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 *ds = rpu_device.getDecayScale()[0];
   T *df = rpu_device.getDiffusionRate()[0];
   T *rb = rpu_device.getResetBias()[0];
+  T *pw = rpu_device.getPersistentWeights()[0];
 
   for (int i = 0; i < d_size; ++i) {
     for (int j = 0; j < x_size; ++j) {
@@ -119,25 +135,45 @@ void PulsedRPUDeviceCuda<T>::populateFrom(const AbstractRPUDevice<T> &rpu_device
       tmp_ds[l_t] = ds[l];
       tmp_df[l_t] = df[l];
       tmp_rb[l_t] = rb[l];
+      tmp_pw[l_t] = pw[l];
     }
   }
 
   dev_4params_->assign(tmp);
-  // other parameters
   dev_decay_scale_->assign(tmp_ds);
-  dev_diffusion_rate_->assign(tmp_df);
-  dev_reset_bias_->assign(tmp_rb);
+
+  // other parameters (on the fly)
+  const auto &par = getPar();
+  if (par.diffusion > 0.0) {
+    dev_diffusion_rate_ = RPU::make_unique<CudaArray<T>>(this->context_, size);
+    dev_diffusion_rate_->assign(tmp_df);
+  }
+
+  if (par.needsResetBias()) {
+    dev_reset_bias_ = RPU::make_unique<CudaArray<T>>(this->context_, size);
+    dev_reset_bias_->assign(tmp_rb);
+  }
+
+  if (par.usesPersistentWeight()) {
+    dev_persistent_weights_ = RPU::make_unique<CudaArray<T>>(this->context_, size);
+    dev_persistent_weights_->assign(tmp_pw);
+  }
 
   this->context_->synchronize();
 
   delete[] tmp_ds;
   delete[] tmp_df;
   delete[] tmp_rb;
+  delete[] tmp_pw;
   delete[] tmp;
 }
 
 template <typename T>
 void PulsedRPUDeviceCuda<T>::applyWeightUpdate(T *weights, T *dw_and_current_weight_out) {
+
+  if (getPar().usesPersistentWeight()) {
+    RPU_FATAL("ResetCols is not supported with write_noise_std>0!");
+  }
   RPU::math::elemaddcopysat<T>(
       this->context_, weights, dw_and_current_weight_out, this->size_,
       dev_4params_->getDataConst());
@@ -146,20 +182,34 @@ void PulsedRPUDeviceCuda<T>::applyWeightUpdate(T *weights, T *dw_and_current_wei
 template <typename T>
 void PulsedRPUDeviceCuda<T>::decayWeights(T *weights, T alpha, bool bias_no_decay) {
 
+  T *w = getPar().usesPersistentWeight() ? dev_persistent_weights_->getData() : weights;
+
   RPU::math::elemscalealpha<T>(
-      this->context_, weights, bias_no_decay ? MAX(this->size_ - this->d_size_, 0) : this->size_,
+      this->context_, w, bias_no_decay ? MAX(this->size_ - this->d_size_, 0) : this->size_,
       dev_decay_scale_->getData(), dev_4params_->getData(), alpha);
+
+  applyUpdateWriteNoise(weights);
 }
 
 template <typename T> void PulsedRPUDeviceCuda<T>::decayWeights(T *weights, bool bias_no_decay) {
 
+  T *w = getPar().usesPersistentWeight() ? dev_persistent_weights_->getData() : weights;
+
   RPU::math::elemscale<T>(
-      this->context_, weights, bias_no_decay ? MAX(this->size_ - this->d_size_, 0) : this->size_,
+      this->context_, w, bias_no_decay ? MAX(this->size_ - this->d_size_, 0) : this->size_,
       dev_decay_scale_->getData(), dev_4params_->getData());
+
+  applyUpdateWriteNoise(weights);
 }
 
 template <typename T> void PulsedRPUDeviceCuda<T>::diffuseWeights(T *weights) {
 
+  if (dev_diffusion_rate_ == nullptr) {
+    return; // no diffusion
+  }
+
+  T *w = getPar().usesPersistentWeight() ? dev_persistent_weights_->getData() : weights;
+
   if (this->dev_diffusion_nrnd_ == nullptr) {
     this->initDiffusionRnd();
     this->rnd_context_->randNormal(
@@ -168,20 +218,28 @@ template <typename T> void PulsedRPUDeviceCuda<T>::diffuseWeights(T *weights) {
   this->rnd_context_->synchronize();
 
   RPU::math::elemasb02<T>(
-      this->context_, weights, this->size_, this->dev_diffusion_nrnd_->getData(),
+      this->context_, w, this->size_, this->dev_diffusion_nrnd_->getData(),
       dev_diffusion_rate_->getData(), dev_4params_->getData());
 
   this->rnd_context_->recordWaitEvent(this->context_->getStream());
   this->rnd_context_->randNormal(
       this->dev_diffusion_nrnd_->getData(), this->dev_diffusion_nrnd_->getSize());
+
+  // Note: write noise will use the same rand to save memory. If
+  // diffusion + writenoise is often needed one might want to add an
+  // extra variable for the random numbers
+  applyUpdateWriteNoise(weights);
 }
 
 template <typename T> void PulsedRPUDeviceCuda<T>::clipWeights(T *weights, T clip) {
 
-  RPU::math::elemsat<T>(this->context_, weights, this->size_, dev_4params_->getData());
+  T *w = getPar().usesPersistentWeight() ? dev_persistent_weights_->getData() : weights;
+
+  RPU::math::elemsat<T>(this->context_, w, this->size_, dev_4params_->getData());
   if (clip >= 0) {
-    RPU::math::aclip<T>(this->context_, weights, this->size_, clip);
+    RPU::math::aclip<T>(this->context_, w, this->size_, clip);
   }
+  applyUpdateWriteNoise(weights);
 }
 
 template <typename T> void PulsedRPUDeviceCuda<T>::initResetRnd() {
@@ -197,10 +255,42 @@ template <typename T> void PulsedRPUDeviceCuda<T>::initResetRnd() {
   this->rnd_context_->synchronize();
 }
 
+template <typename T> void PulsedRPUDeviceCuda<T>::applyUpdateWriteNoise(T *dev_weights) {
+
+  const auto &par = getPar();
+
+  if (!par.usesPersistentWeight()) {
+    return;
+  }
+  // re-uses the diffusion rnd
+  if (this->dev_diffusion_nrnd_ == nullptr) {
+    this->initDiffusionRnd();
+    this->rnd_context_->randNormal(
+        this->dev_diffusion_nrnd_->getData(), this->dev_diffusion_nrnd_->getSize());
+  }
+  this->rnd_context_->synchronize();
+
+  RPU::math::elemweightedsum<T>(
+      this->context_, dev_weights, this->size_, dev_persistent_weights_->getData(), (T)1.0,
+      this->dev_diffusion_nrnd_->getData(), par.write_noise_std);
+
+  this->rnd_context_->recordWaitEvent(this->context_->getStream());
+  this->rnd_context_->randNormal(
+      this->dev_diffusion_nrnd_->getData(), this->dev_diffusion_nrnd_->getSize());
+}
+
 template <typename T>
 void PulsedRPUDeviceCuda<T>::resetCols(T *weights, int start_col, int n_cols, T reset_prob) {
   // col-major in CUDA.
 
+  if (dev_reset_bias_ == nullptr) {
+    return; // no reset
+  }
+
+  if (getPar().usesPersistentWeight()) {
+    RPU_FATAL("ResetCols is not supported with write_noise_std>0!");
+  }
+
   if (dev_reset_nrnd_ == nullptr) {
     initResetRnd();
   }
diff --git a/src/rpucuda/cuda/rpucuda_pulsed_device.h b/src/rpucuda/cuda/rpucuda_pulsed_device.h
index 47777634..a32c92af 100644
--- a/src/rpucuda/cuda/rpucuda_pulsed_device.h
+++ b/src/rpucuda/cuda/rpucuda_pulsed_device.h
@@ -107,6 +107,7 @@ template <typename T> class PulsedRPUDeviceCuda : public PulsedRPUDeviceCudaBase
     swap(a.dev_4params_, b.dev_4params_);
     swap(a.dev_reset_nrnd_, b.dev_reset_nrnd_);
     swap(a.dev_reset_flag_, b.dev_reset_flag_);
+    swap(a.dev_persistent_weights_, b.dev_persistent_weights_);
   };
 
   // implement abstract functions
@@ -144,9 +145,12 @@ template <typename T> class PulsedRPUDeviceCuda : public PulsedRPUDeviceCudaBase
   virtual float *get4ParamsData() { return dev_4params_->getData(); }
 
 protected:
+  virtual void applyUpdateWriteNoise(T *dev_weights);
+
   void initResetRnd();
   std::unique_ptr<CudaArray<float>> dev_reset_nrnd_ = nullptr;
   std::unique_ptr<CudaArray<float>> dev_reset_flag_ = nullptr;
+  std::unique_ptr<CudaArray<float>> dev_persistent_weights_ = nullptr;
 
 private:
   void initialize();
diff --git a/src/rpucuda/cuda/rpucuda_transfer_device.cu b/src/rpucuda/cuda/rpucuda_transfer_device.cu
index b8347764..44dcc69b 100644
--- a/src/rpucuda/cuda/rpucuda_transfer_device.cu
+++ b/src/rpucuda/cuda/rpucuda_transfer_device.cu
@@ -122,7 +122,7 @@ void TransferRPUDeviceCuda<T>::populateFrom(const AbstractRPUDevice<T> &rpu_devi
   initialize(); // pwu/iom
 
   current_col_indices_.resize(this->n_devices_ - 1);
-  std::fill(current_col_indices_.begin(), current_col_indices_.end(), (int) 0);
+  std::fill(current_col_indices_.begin(), current_col_indices_.end(), (int)0);
 
   this->current_update_idx_ = 0;
 
diff --git a/src/rpucuda/rpu.h b/src/rpucuda/rpu.h
index 467ecf79..c48f3821 100644
--- a/src/rpucuda/rpu.h
+++ b/src/rpucuda/rpu.h
@@ -598,14 +598,14 @@ template <typename T> class RPUSimple : public RPUAbstract<T> {
 
 private:
   void alpha_warning() {
-    static int alpha_warning_count = 0;
-    alpha_warning_count++;
-    if (alpha_warning_count < 25) {
-      std::cout
-          << "Warning: setting weights with alpha scale. Note that alpha scale is NOT respected "
-             "when getting weights, saving or export, where thus wrong weights are obtained."
-          << std::endl;
-    }
+    DEBUG_CALL(
+        static int alpha_warning_count = 0; alpha_warning_count++; if (alpha_warning_count < 25) {
+          std::cout
+              << "Warning: setting weights with alpha scale. Note that alpha scale is NOT "
+                 "respected "
+                 "when getting weights, saving or export, where thus wrong weights are obtained."
+              << std::endl;
+        })
   }
 
 public:
diff --git a/src/rpucuda/rpu_difference_device.cpp b/src/rpucuda/rpu_difference_device.cpp
index 18a14639..6f6efa98 100644
--- a/src/rpucuda/rpu_difference_device.cpp
+++ b/src/rpucuda/rpu_difference_device.cpp
@@ -95,9 +95,7 @@ void DifferenceRPUDevice<T>::populate(
 /*********************************************************************************/
 /* update */
 
-template <typename T> bool DifferenceRPUDevice<T>::isInverted() const {
-  return g_plus_ == 0;
-}
+template <typename T> bool DifferenceRPUDevice<T>::isInverted() const { return g_plus_ == 0; }
 
 template <typename T> inline void DifferenceRPUDevice<T>::invert() {
   std::swap(g_plus_, g_minus_);
@@ -108,6 +106,9 @@ template <typename T> inline void DifferenceRPUDevice<T>::invert() {
 template <typename T>
 void DifferenceRPUDevice<T>::initUpdateCycle(
     T **weights, const PulsedUpdateMetaParameter<T> &up, T current_lr, int m_batch_info) {
+
+  VectorRPUDevice<T>::initUpdateCycle(weights, up, current_lr, m_batch_info);
+
   if (a_indices_.size() < (size_t)up.desired_BL) {
     a_indices_.resize(up.desired_BL);
     b_indices_.resize(up.desired_BL);
diff --git a/src/rpucuda/rpu_difference_device.h b/src/rpucuda/rpu_difference_device.h
index f388a977..3fda0a11 100644
--- a/src/rpucuda/rpu_difference_device.h
+++ b/src/rpucuda/rpu_difference_device.h
@@ -107,14 +107,8 @@ template <typename T> class DifferenceRPUDevice : public VectorRPUDevice<T> {
 
   void invert();
 
-  // current_update_count not needed:
   void initUpdateCycle(
       T **weights, const PulsedUpdateMetaParameter<T> &up, T current_lr, int m_batch_info) override;
-  void finishUpdateCycle(
-      T **weights,
-      const PulsedUpdateMetaParameter<T> &up,
-      T current_lr,
-      int m_batch_info) override{};
 
   void setHiddenUpdateIdx(int idx) override{};
 
diff --git a/src/rpucuda/rpu_expstep_device.cpp b/src/rpucuda/rpu_expstep_device.cpp
index 84f04d20..16ae8082 100644
--- a/src/rpucuda/rpu_expstep_device.cpp
+++ b/src/rpucuda/rpu_expstep_device.cpp
@@ -23,6 +23,7 @@ void ExpStepRPUDevice<T>::populate(
 template <typename T>
 inline void update_once(
     T &w,
+    T &w_apparent,
     int &sign,
     T &min_bound,
     T &max_bound,
@@ -35,23 +36,28 @@ inline void update_once(
     const T &es_gamma_down,
     const T &es_gamma_up,
     const T &dw_min_std,
+    const T &write_noise_std,
     RNG<T> *rng) {
   T b_diff = (max_bound - min_bound);
   if (b_diff > 0.0) {
     T z = (T)2.0 * w / b_diff * es_a + es_b;
 
     if (sign > 0) {
-      T y_down = MAX((T)1 - es_A_down * exp(es_gamma_down * (-z)), (T)0);
+      T y_down = MAX((T)1 - es_A_down * expf(es_gamma_down * (-z)), (T)0);
       w -= y_down * ((T)1.0 + dw_min_std * rng->sampleGauss()) * scale_down;
 
     } else {
-      T y_up = MAX(1 - es_A_up * exp(es_gamma_up * z), (T)0);
+      T y_up = MAX(1 - es_A_up * expf(es_gamma_up * z), (T)0);
       w += y_up * ((T)1.0 + dw_min_std * rng->sampleGauss()) * scale_up;
     }
 
     // always check both bounds
     w = MIN(w, max_bound);
     w = MAX(w, min_bound);
+
+    if (write_noise_std > (T)0.0) {
+      w_apparent = w + write_noise_std * rng->sampleGauss();
+    }
   }
 }
 
@@ -59,35 +65,41 @@ template <typename T>
 void ExpStepRPUDevice<T>::doSparseUpdate(
     T **weights, int i, const int *x_signed_indices, int x_count, int d_sign, RNG<T> *rng) {
 
+  const auto &par = getPar();
   T *scale_down = this->w_scale_down_[i];
   T *scale_up = this->w_scale_up_[i];
-  T *w = weights[i];
+  T *w = par.usesPersistentWeight() ? this->w_persistent_[i] : weights[i];
+  T *w_apparent = weights[i];
   T *min_bound = this->w_min_bound_[i];
   T *max_bound = this->w_max_bound_[i];
 
-  const auto &par = getPar();
+  T write_noise_std = par.getScaledWriteNoise();
 
   PULSED_UPDATE_W_LOOP(update_once(
-                           w[j], sign, min_bound[j], max_bound[j], scale_down[j], scale_up[j],
-                           par.es_a, par.es_b, par.es_A_down, par.es_A_up, par.es_gamma_down,
-                           par.es_gamma_up, par.dw_min_std, rng););
+                           w[j], w_apparent[j], sign, min_bound[j], max_bound[j], scale_down[j],
+                           scale_up[j], par.es_a, par.es_b, par.es_A_down, par.es_A_up,
+                           par.es_gamma_down, par.es_gamma_up, par.dw_min_std, write_noise_std,
+                           rng););
 };
 
 template <typename T>
 void ExpStepRPUDevice<T>::doDenseUpdate(T **weights, int *coincidences, RNG<T> *rng) {
 
+  const auto &par = getPar();
   T *scale_down = this->w_scale_down_[0];
   T *scale_up = this->w_scale_up_[0];
-  T *w = weights[0];
+  T *w = par.usesPersistentWeight() ? this->w_persistent_[0] : weights[0];
+  T *w_apparent = weights[0];
   T *min_bound = this->w_min_bound_[0];
   T *max_bound = this->w_max_bound_[0];
 
-  const auto &par = getPar();
+  T write_noise_std = par.getScaledWriteNoise();
 
   PULSED_UPDATE_W_LOOP_DENSE(update_once(
-                                 w[j], sign, min_bound[j], max_bound[j], scale_down[j], scale_up[j],
-                                 par.es_a, par.es_b, par.es_A_down, par.es_A_up, par.es_gamma_down,
-                                 par.es_gamma_up, par.dw_min_std, rng););
+                                 w[j], w_apparent[j], sign, min_bound[j], max_bound[j],
+                                 scale_down[j], scale_up[j], par.es_a, par.es_b, par.es_A_down,
+                                 par.es_A_up, par.es_gamma_down, par.es_gamma_up, par.dw_min_std,
+                                 write_noise_std, rng););
 }
 
 template class ExpStepRPUDevice<float>;
diff --git a/src/rpucuda/rpu_expstep_device.h b/src/rpucuda/rpu_expstep_device.h
index 62a6c601..2b9303a7 100644
--- a/src/rpucuda/rpu_expstep_device.h
+++ b/src/rpucuda/rpu_expstep_device.h
@@ -20,26 +20,28 @@ namespace RPU {
 
 template <typename T> class ExpStepRPUDevice;
 
-BUILD_PULSED_DEVICE_META_PARAMETER(ExpStep,
-                                   /*implements*/
-                                   DeviceUpdateType::ExpStep,
-                                   /*parameter def*/
-                                   T es_A_up = (T)0.00081;        /* p_000081 */
-                                   T es_A_down = (T)0.36833;      /* p_036833 */
-                                   T es_gamma_up = (T)12.44625;   /*  p_1244625 */
-                                   T es_gamma_down = (T)12.78785; /* p_1278785 */
-                                   T es_a = (T)0.244;             /* p_0244 */
-                                   T es_b = (T)0.2425;            /*p_02425 */
-                                   ,
-                                   /*print body*/
-                                   ss << "\t es_A_up:\t\t" << es_A_up << std::endl;
-                                   ss << "\t es_A_down:\t\t" << es_A_down << std::endl;
-                                   ss << "\t es_gamma_up:\t\t" << es_gamma_up << std::endl;
-                                   ss << "\t es_gamma_down:\t\t" << es_gamma_down << std::endl;
-                                   ss << "\t es_a:\t\t\t" << es_a << std::endl;
-                                   ss << "\t es_b:\t\t\t" << es_b << std::endl;
-                                   ,
-                                   /*add */
+BUILD_PULSED_DEVICE_META_PARAMETER(
+    ExpStep,
+    /*implements*/
+    DeviceUpdateType::ExpStep,
+    /*parameter def*/
+    T es_A_up = (T)0.00081;        /* p_000081 */
+    T es_A_down = (T)0.36833;      /* p_036833 */
+    T es_gamma_up = (T)12.44625;   /*  p_1244625 */
+    T es_gamma_down = (T)12.78785; /* p_1278785 */
+    T es_a = (T)0.244;             /* p_0244 */
+    T es_b = (T)0.2425;            /*p_02425 */
+    ,
+    /*print body*/
+    ss << "\t es_A_up:\t\t" << es_A_up << std::endl;
+    ss << "\t es_A_down:\t\t" << es_A_down << std::endl;
+    ss << "\t es_gamma_up:\t\t" << es_gamma_up << std::endl;
+    ss << "\t es_gamma_down:\t\t" << es_gamma_down << std::endl;
+    ss << "\t es_a:\t\t\t" << es_a << std::endl;
+    ss << "\t es_b:\t\t\t" << es_b << std::endl;
+    ,
+    /*add */
+    bool implementsWriteNoise() const override { return true; };
 
 );
 
diff --git a/src/rpucuda/rpu_linearstep_device.cpp b/src/rpucuda/rpu_linearstep_device.cpp
index 3b329b72..c61adf6b 100644
--- a/src/rpucuda/rpu_linearstep_device.cpp
+++ b/src/rpucuda/rpu_linearstep_device.cpp
@@ -65,11 +65,14 @@ void LinearStepRPUDevice<T>::populate(
 
 template <typename T> void LinearStepRPUDevice<T>::printDP(int x_count, int d_count) const {
 
-  if (x_count < 0 || x_count > this->x_size_)
+  if (x_count < 0 || x_count > this->x_size_) {
     x_count = this->x_size_;
+  }
 
-  if (d_count < 0 || d_count > this->d_size_)
+  if (d_count < 0 || d_count > this->d_size_) {
     d_count = this->d_size_;
+  }
+  bool persist_if = getPar().usesPersistentWeight();
 
   for (int i = 0; i < d_count; ++i) {
     for (int j = 0; j < x_count; ++j) {
@@ -81,15 +84,24 @@ template <typename T> void LinearStepRPUDevice<T>::printDP(int x_count, int d_co
       std::cout << this->sup_[i][j].scale_down << ">,<";
       std::cout << w_slope_up_[i][j] << ",";
       std::cout << w_slope_down_[i][j] << ">]";
-
-      std::cout << std::endl;
+      std::cout.precision(10);
+      std::cout << this->sup_[i][j].decay_scale << ", ";
+      std::cout.precision(6);
+      std::cout << this->w_diffusion_rate_[i][j] << ", ";
+      std::cout << this->w_reset_bias_[i][j];
+      if (persist_if) {
+        std::cout << ", " << this->w_persistent_[i][j];
+      }
+      std::cout << "]";
     }
+    std::cout << std::endl;
   }
 }
 
 template <typename T>
 inline void update_once_mult(
     T &w,
+    T &w_apparent,
     int &sign,
     T &scale_down,
     T &scale_up,
@@ -98,6 +110,7 @@ inline void update_once_mult(
     T &min_bound,
     T &max_bound,
     const T &dw_min_std,
+    const T &write_noise_std,
     RNG<T> *rng) {
   if (sign > 0) {
     w -= (slope_down * w + scale_down) * ((T)1.0 + dw_min_std * rng->sampleGauss());
@@ -106,11 +119,16 @@ inline void update_once_mult(
   }
   w = MAX(w, min_bound);
   w = MIN(w, max_bound);
+
+  if (write_noise_std > (T)0.0) {
+    w_apparent = w + write_noise_std * rng->sampleGauss();
+  }
 }
 
 template <typename T>
 inline void update_once_add(
     T &w,
+    T &w_apparent,
     int &sign,
     T &scale_down,
     T &scale_up,
@@ -119,6 +137,7 @@ inline void update_once_add(
     T &min_bound,
     T &max_bound,
     const T &dw_min_std,
+    const T &write_noise_std,
     RNG<T> *rng) {
   if (sign > 0) {
     w -= slope_down * w + scale_down * ((T)1.0 + dw_min_std * rng->sampleGauss());
@@ -127,52 +146,66 @@ inline void update_once_add(
   }
   w = MAX(w, min_bound);
   w = MIN(w, max_bound);
+
+  if (write_noise_std > (T)0.0) {
+    w_apparent = w + write_noise_std * rng->sampleGauss();
+  }
 }
 
 template <typename T>
 void LinearStepRPUDevice<T>::doSparseUpdate(
     T **weights, int i, const int *x_signed_indices, int x_count, int d_sign, RNG<T> *rng) {
 
+  const auto &par = getPar();
+
   T *scale_down = this->w_scale_down_[i];
   T *scale_up = this->w_scale_up_[i];
   T *slope_down = w_slope_down_[i];
   T *slope_up = w_slope_up_[i];
-  T *w = weights[i];
+  T *w = par.usesPersistentWeight() ? this->w_persistent_[i] : weights[i];
+  T *w_apparent = weights[i];
   T *min_bound = this->w_min_bound_[i];
   T *max_bound = this->w_max_bound_[i];
-  const auto &par = getPar();
 
+  T write_noise_std = par.getScaledWriteNoise();
   if (par.ls_mult_noise) {
     PULSED_UPDATE_W_LOOP(update_once_mult(
-                             w[j], sign, scale_down[j], scale_up[j], slope_down[j], slope_up[j],
-                             min_bound[j], max_bound[j], par.dw_min_std, rng););
+                             w[j], w_apparent[j], sign, scale_down[j], scale_up[j], slope_down[j],
+                             slope_up[j], min_bound[j], max_bound[j], par.dw_min_std,
+                             write_noise_std, rng););
   } else {
     PULSED_UPDATE_W_LOOP(update_once_add(
-                             w[j], sign, scale_down[j], scale_up[j], slope_down[j], slope_up[j],
-                             min_bound[j], max_bound[j], par.dw_min_std, rng););
+                             w[j], w_apparent[j], sign, scale_down[j], scale_up[j], slope_down[j],
+                             slope_up[j], min_bound[j], max_bound[j], par.dw_min_std,
+                             write_noise_std, rng););
   }
 }
 
 template <typename T>
 void LinearStepRPUDevice<T>::doDenseUpdate(T **weights, int *coincidences, RNG<T> *rng) {
 
+  const auto &par = getPar();
+
   T *scale_down = this->w_scale_down_[0];
   T *scale_up = this->w_scale_up_[0];
   T *slope_down = w_slope_down_[0];
   T *slope_up = w_slope_up_[0];
-  T *w = weights[0];
+  T *w = par.usesPersistentWeight() ? this->w_persistent_[0] : weights[0];
+  T *w_apparent = weights[0];
   T *min_bound = this->w_min_bound_[0];
   T *max_bound = this->w_max_bound_[0];
-  const auto &par = getPar();
+  T write_noise_std = par.getScaledWriteNoise();
 
   if (par.ls_mult_noise) {
     PULSED_UPDATE_W_LOOP_DENSE(update_once_mult(
-                                   w[j], sign, scale_down[j], scale_up[j], slope_down[j],
-                                   slope_up[j], min_bound[j], max_bound[j], par.dw_min_std, rng););
+                                   w[j], w_apparent[j], sign, scale_down[j], scale_up[j],
+                                   slope_down[j], slope_up[j], min_bound[j], max_bound[j],
+                                   par.dw_min_std, write_noise_std, rng););
   } else {
     PULSED_UPDATE_W_LOOP_DENSE(update_once_add(
-                                   w[j], sign, scale_down[j], scale_up[j], slope_down[j],
-                                   slope_up[j], min_bound[j], max_bound[j], par.dw_min_std, rng););
+                                   w[j], w_apparent[j], sign, scale_down[j], scale_up[j],
+                                   slope_down[j], slope_up[j], min_bound[j], max_bound[j],
+                                   par.dw_min_std, write_noise_std, rng););
   }
 }
 
diff --git a/src/rpucuda/rpu_linearstep_device.h b/src/rpucuda/rpu_linearstep_device.h
index a610ecd5..8026ac59 100644
--- a/src/rpucuda/rpu_linearstep_device.h
+++ b/src/rpucuda/rpu_linearstep_device.h
@@ -49,7 +49,7 @@ BUILD_PULSED_DEVICE_META_PARAMETER(
        << std::endl;
     ,
     /*Add*/
-);
+    bool implementsWriteNoise() const override { return true; };);
 
 /*SoftBounds is essentially a specific LinearStep device. It will
   create one. Careful though, one could in principle  set the ls parameters
diff --git a/src/rpucuda/rpu_pulsed_device.cpp b/src/rpucuda/rpu_pulsed_device.cpp
index 47d3b578..f5cb3ad7 100644
--- a/src/rpucuda/rpu_pulsed_device.cpp
+++ b/src/rpucuda/rpu_pulsed_device.cpp
@@ -44,6 +44,10 @@ void PulsedRPUDeviceMetaParameter<T>::printToStream(std::stringstream &ss) const
     ss << "\t resets to:\t\t" << reset << "\t(dtod=" << reset_dtod << ", ctoc=" << reset_std << ")"
        << std::endl;
 
+    if (this->implementsWriteNoise() && write_noise_std > (T)0.0) {
+      ss << "\t write noise std:\t" << write_noise_std << std::endl;
+    }
+
     if (this->lifetime > 0) {
       ss << "\t lifetime [decay]:\t" << this->lifetime << "\t(dtod=" << lifetime_dtod << ")"
          << std::endl;
@@ -90,6 +94,7 @@ template <typename T> void PulsedRPUDevice<T>::allocateContainers() {
   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);
 
   // we better set everything to zero.
   for (int j = 0; j < x_sz; ++j) {
@@ -101,6 +106,7 @@ template <typename T> void PulsedRPUDevice<T>::allocateContainers() {
       w_decay_scale_[i][j] = (T)1.0; // no decay
       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();
@@ -127,6 +133,7 @@ template <typename T> void PulsedRPUDevice<T>::freeContainers() {
     Array_2D_Free<T>(w_decay_scale_);
     Array_2D_Free<T>(w_diffusion_rate_);
     Array_2D_Free<T>(w_reset_bias_);
+    Array_2D_Free<T>(w_persistent_);
 
     containers_allocated_ = false;
   }
@@ -209,6 +216,7 @@ PulsedRPUDevice<T> &PulsedRPUDevice<T>::operator=(PulsedRPUDevice<T> &&other) {
   w_decay_scale_ = other.w_decay_scale_;
   w_diffusion_rate_ = other.w_diffusion_rate_;
   w_reset_bias_ = other.w_reset_bias_;
+  w_persistent_ = other.w_persistent_;
 
   // set pointers to null
   other.sup_ = nullptr;
@@ -222,6 +230,7 @@ PulsedRPUDevice<T> &PulsedRPUDevice<T>::operator=(PulsedRPUDevice<T> &&other) {
   other.w_diffusion_rate_ = nullptr;
 
   other.w_reset_bias_ = nullptr;
+  other.w_persistent_ = nullptr;
 
   return *this;
 }
@@ -238,6 +247,9 @@ template <typename T> void PulsedRPUDevice<T>::getDPNames(std::vector<std::strin
   if (!getPar().legacy_params) {
     names.push_back(std::string("reset_bias"));
   }
+  if (getPar().usesPersistentWeight()) {
+    names.push_back(std::string("hidden_weigths"));
+  }
 }
 
 template <typename T>
@@ -252,14 +264,18 @@ void PulsedRPUDevice<T>::getDeviceParameter(std::vector<T *> &data_ptrs) const {
   }
 
   for (int i = 0; i < this->size_; ++i) {
-    data_ptrs[0][i] = w_max_bound_[0][i];
-    data_ptrs[1][i] = w_min_bound_[0][i];
-    data_ptrs[2][i] = w_scale_up_[0][i];
-    data_ptrs[3][i] = w_scale_down_[0][i];
-    data_ptrs[4][i] = w_decay_scale_[0][i];
-    data_ptrs[5][i] = w_diffusion_rate_[0][i];
+    int n = 0;
+    data_ptrs[n++][i] = w_max_bound_[0][i];
+    data_ptrs[n++][i] = w_min_bound_[0][i];
+    data_ptrs[n++][i] = w_scale_up_[0][i];
+    data_ptrs[n++][i] = w_scale_down_[0][i];
+    data_ptrs[n++][i] = w_decay_scale_[0][i];
+    data_ptrs[n++][i] = w_diffusion_rate_[0][i];
     if (!getPar().legacy_params) {
-      data_ptrs[6][i] = w_reset_bias_[0][i];
+      data_ptrs[n++][i] = w_reset_bias_[0][i];
+    }
+    if (getPar().usesPersistentWeight()) {
+      data_ptrs[n++][i] = w_persistent_[0][i];
     }
   }
 };
@@ -313,8 +329,8 @@ void PulsedRPUDevice<T>::setDeviceParameter(const std::vector<T *> &data_ptrs) {
 template <typename T> void PulsedRPUDevice<T>::decayWeights(T **weights, bool bias_no_decay) {
 
   // maybe a bit overkill to check the bounds...
+  T *w = getPar().usesPersistentWeight() ? w_persistent_[0] : weights[0];
   T *wd = w_decay_scale_[0];
-  T *w = weights[0];
   T *max_bound = w_max_bound_[0];
   T *min_bound = w_min_bound_[0];
 
@@ -334,15 +350,15 @@ template <typename T> void PulsedRPUDevice<T>::decayWeights(T **weights, bool bi
       w[i] = MAX(w[i], min_bound[i]);
     }
   }
+  applyUpdateWriteNoise(weights);
 }
 
 template <typename T>
 void PulsedRPUDevice<T>::decayWeights(T **weights, T alpha, bool bias_no_decay) {
 
   // maybe a bit overkill to check the bounds...
-
+  T *w = getPar().usesPersistentWeight() ? w_persistent_[0] : weights[0];
   T *wd = w_decay_scale_[0];
-  T *w = weights[0];
   T *max_bound = w_max_bound_[0];
   T *min_bound = w_min_bound_[0];
 
@@ -362,12 +378,13 @@ void PulsedRPUDevice<T>::decayWeights(T **weights, T alpha, bool bias_no_decay)
       w[i] = MAX(w[i], min_bound[i]);
     }
   }
+  applyUpdateWriteNoise(weights);
 }
 
 template <typename T> void PulsedRPUDevice<T>::diffuseWeights(T **weights, RNG<T> &rng) {
 
+  T *w = getPar().usesPersistentWeight() ? w_persistent_[0] : weights[0];
   T *diffusion_rate = &(w_diffusion_rate_[0][0]);
-  T *w = &(weights[0][0]);
   T *max_bound = &(w_max_bound_[0][0]);
   T *min_bound = &(w_min_bound_[0][0]);
 
@@ -377,12 +394,12 @@ template <typename T> void PulsedRPUDevice<T>::diffuseWeights(T **weights, RNG<T
     w[i] = MIN(w[i], max_bound[i]);
     w[i] = MAX(w[i], min_bound[i]);
   }
+  applyUpdateWriteNoise(weights);
 }
 
 template <typename T> void PulsedRPUDevice<T>::clipWeights(T **weights, T clip) {
   // apply hard bounds
-
-  T *w = &(weights[0][0]);
+  T *w = getPar().usesPersistentWeight() ? w_persistent_[0] : weights[0];
   T *max_bound = &(w_max_bound_[0][0]);
   T *min_bound = &(w_min_bound_[0][0]);
   if (clip < 0.0) { // only apply bounds
@@ -398,12 +415,17 @@ template <typename T> void PulsedRPUDevice<T>::clipWeights(T **weights, T clip)
       w[i] = MAX(w[i], MAX(min_bound[i], -clip));
     }
   }
+  applyUpdateWriteNoise(weights);
 }
 
 template <typename T>
 void PulsedRPUDevice<T>::resetCols(
     T **weights, int start_col, int n_col, T reset_prob, RealWorldRNG<T> &rng) {
 
+  if (getPar().usesPersistentWeight()) {
+    RPU_FATAL("ResetCols is not supported with write_noise_std>0!");
+  }
+
   T reset_std = getPar().reset_std;
   for (int j = 0; j < this->x_size_; ++j) {
     if ((start_col + n_col <= this->x_size_ && j >= start_col && j < start_col + n_col) ||
@@ -459,9 +481,40 @@ void PulsedRPUDevice<T>::copyInvertDeviceParameter(const PulsedRPUDeviceBase<T>
 }
 
 template <typename T> bool PulsedRPUDevice<T>::onSetWeights(T **weights) {
-  // apply hard bounds
-  this->clipWeights(weights, (T)-1.0);
-  return false; // whether device was changed
+
+  // apply hard bounds to given weights
+  T *w = weights[0];
+  T *max_bound = &(w_max_bound_[0][0]);
+  T *min_bound = &(w_min_bound_[0][0]);
+  PRAGMA_SIMD
+  for (int i = 0; i < this->size_; ++i) {
+    w[i] = MIN(w[i], max_bound[i]);
+    w[i] = MAX(w[i], min_bound[i]);
+  }
+
+  if (getPar().usesPersistentWeight()) {
+    PRAGMA_SIMD
+    for (int i = 0; i < this->size_; i++) {
+      w_persistent_[0][i] = w[i];
+    }
+    applyUpdateWriteNoise(weights);
+    return true; // modified device thus true
+  } else {
+    return false; // whether device was changed
+  }
+}
+
+template <typename T> void PulsedRPUDevice<T>::applyUpdateWriteNoise(T **weights) {
+  // applies new noise to ALL weight values
+  T uw_std = getPar().getScaledWriteNoise();
+
+  if (uw_std <= 0) {
+    return; // nothing to be done, weights assumed to already updated
+  }
+
+  for (int i = 0; i < this->size_; i++) {
+    weights[0][i] = w_persistent_[0][i] + uw_std * write_noise_rng_.sampleGauss();
+  }
 }
 
 /*********************************************************************************/
@@ -545,7 +598,6 @@ void PulsedRPUDevice<T>::populate(const PulsedRPUDeviceMetaParameter<T> &p, Real
 
       //--------------------
       // diffusion
-      par.diffusion = MAX(par.diffusion, (T)0.0);
       {
         T t = fabs(par.diffusion * (1 + par.diffusion_dtod * rng->sampleGauss()));
         sup_[i][j].diffusion_rate = t;
@@ -553,16 +605,13 @@ void PulsedRPUDevice<T>::populate(const PulsedRPUDeviceMetaParameter<T> &p, Real
 
       //--------------------
       // reset
-      par.reset_dtod = MAX(par.reset_dtod, (T)0.0);
-      par.reset_std = MAX(par.reset_std, (T)0.0);
-      {
-        T t = par.reset * ((T)1.0 + par.reset_dtod * rng->sampleGauss());
+      { // additive dtod
+        T t = par.reset + par.reset_dtod * rng->sampleGauss();
         sup_[i][j].reset_bias = t;
       }
 
       //--------------------
       // decay
-      par.lifetime = MAX(par.lifetime, (T)0.0);
       {
         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;
@@ -586,21 +635,31 @@ void PulsedRPUDevice<T>::populate(const PulsedRPUDeviceMetaParameter<T> &p, Real
 template <typename T> void PulsedRPUDevice<T>::printDP(int x_count, int d_count) const {
   int x_count1 = x_count;
   int d_count1 = d_count;
-  if (x_count < 0 || x_count > this->x_size_)
+  if (x_count < 0 || x_count > this->x_size_) {
     x_count1 = this->x_size_;
+  }
 
-  if (d_count < 0 || d_count > this->d_size_)
+  if (d_count < 0 || d_count > this->d_size_) {
     d_count1 = this->d_size_;
+  }
+
+  bool persist_if = getPar().usesPersistentWeight();
 
   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 << "[<" << 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.precision(10);
-      std::cout << sup_[i][j].decay_scale;
+      std::cout << sup_[i][j].decay_scale << ", ";
       std::cout.precision(6);
+      std::cout << w_diffusion_rate_[i][j] << ", ";
+      std::cout << w_reset_bias_[i][j];
+      if (persist_if) {
+        std::cout << ", " << w_persistent_[i][j];
+      }
+      std::cout << "]";
     }
     std::cout << std::endl;
   }
diff --git a/src/rpucuda/rpu_pulsed_device.h b/src/rpucuda/rpu_pulsed_device.h
index f7c7f113..65fd3010 100644
--- a/src/rpucuda/rpu_pulsed_device.h
+++ b/src/rpucuda/rpu_pulsed_device.h
@@ -85,6 +85,8 @@ template <typename T> struct PulsedRPUDeviceMetaParameter : PulsedRPUDeviceMetaP
   T reset_std = (T)0.0;
   T reset_dtod = (T)0.0;
 
+  T write_noise_std = (T)0.0;
+
   void printToStream(std::stringstream &ss) const override;
   using SimpleMetaParameter<T>::print;
   std::string getName() const override { return "PulsedRPUDeviceParameter"; };
@@ -93,6 +95,21 @@ template <typename T> struct PulsedRPUDeviceMetaParameter : PulsedRPUDeviceMetaP
   };
   PulsedRPUDeviceMetaParameter<T> *clone() const override { RPU_FATAL("Needs implementation"); };
   DeviceUpdateType implements() const override { return DeviceUpdateType::Undefined; };
+
+  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();
+    if (!implementsWriteNoise() && usesPersistentWeight()) {
+      RPU_FATAL("Device does not support write noise");
+    }
+    reset_dtod = MAX(reset_dtod, (T)0.0);
+    reset_std = MAX(reset_std, (T)0.0);
+    reset = MAX(reset, (T)0.0);
+  };
 };
 
 template <typename T> class PulsedRPUDeviceBase : public SimpleRPUDevice<T> {
@@ -189,6 +206,7 @@ template <typename T> class PulsedRPUDevice : public PulsedRPUDeviceBase<T> {
     swap(a.w_min_bound_, b.w_min_bound_);
     swap(a.w_decay_scale_, b.w_decay_scale_);
     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_);
 
@@ -202,6 +220,7 @@ template <typename T> class PulsedRPUDevice : public PulsedRPUDeviceBase<T> {
   void setDeviceParameter(const std::vector<T *> &data_ptrs) override;
   void printDP(int x_count, int d_count) const override;
 
+  inline T **getPersistentWeights() const { return w_persistent_; };
   inline T **getMaxBound() const { return w_max_bound_; };
   inline T **getMinBound() const { return w_min_bound_; };
   inline T **getDecayScale() const { return w_decay_scale_; };
@@ -214,19 +233,13 @@ template <typename T> class PulsedRPUDevice : public PulsedRPUDeviceBase<T> {
     return static_cast<PulsedRPUDeviceMetaParameter<T> &>(SimpleRPUDevice<T>::getPar());
   };
 
-  // these are turned off for concrete devices. Only for meta [derived from Base]
-  void initUpdateCycle(
-      T **weights,
-      const PulsedUpdateMetaParameter<T> &up,
-      T current_lr,
-      int m_batch_info) override final{};
-  void finishUpdateCycle(
-      T **weights,
-      const PulsedUpdateMetaParameter<T> &up,
-      T current_lr,
-      int m_batch_info) override final{};
-
   T getDwMin() const override { return this->getPar().dw_min; };
+
+  /* Note: In case of persistent data these weight methods below will
+     affect the perstistent state and ALL apparent weight elements
+     will be re-drawn with noise (even if they were not
+     e.g. clipped) */
+
   void decayWeights(T **weights, bool bias_no_decay) override;
   void decayWeights(T **weights, T alpha, bool bias_no_decay) override;
   void diffuseWeights(T **weights, RNG<T> &rng) override;
@@ -249,6 +262,10 @@ template <typename T> class PulsedRPUDevice : public PulsedRPUDeviceBase<T> {
   T **w_decay_scale_ = nullptr;
   T **w_diffusion_rate_ = nullptr;
   T **w_reset_bias_ = nullptr;
+  T **w_persistent_ = nullptr;
+
+  RealWorldRNG<T> write_noise_rng_{0};
+  virtual void applyUpdateWriteNoise(T **weights);
 
 private:
   void freeContainers();
diff --git a/src/rpucuda/rpu_simple_device.h b/src/rpucuda/rpu_simple_device.h
index 8a429505..a494feb3 100644
--- a/src/rpucuda/rpu_simple_device.h
+++ b/src/rpucuda/rpu_simple_device.h
@@ -64,6 +64,8 @@ template <typename T> struct AbstractRPUDeviceMetaParameter : SimpleMetaParamete
   virtual void initialize() {
     _par_initialized = true;
     _device_parameter_mode_manual = false;
+    this->diffusion = MAX(this->diffusion, (T)0.0);
+    this->lifetime = MAX(this->lifetime, (T)0.0);
   }
 
   friend void
diff --git a/src/rpucuda/rpu_transfer_device.cpp b/src/rpucuda/rpu_transfer_device.cpp
index 09307d79..dd81f7d4 100644
--- a/src/rpucuda/rpu_transfer_device.cpp
+++ b/src/rpucuda/rpu_transfer_device.cpp
@@ -292,7 +292,7 @@ template <typename T> void TransferRPUDevice<T>::setTransferVecs(const T *transf
 
 template <typename T> int TransferRPUDevice<T>::resetCounters(bool force) {
   current_col_indices_.resize(this->n_devices_);
-  std::fill(current_col_indices_.begin(), current_col_indices_.end(), (int) 0);
+  std::fill(current_col_indices_.begin(), current_col_indices_.end(), (int)0);
   return VectorRPUDevice<T>::resetCounters(force);
 }
 /*********************************************************************************/
diff --git a/src/rpucuda/rpu_vector_device.cpp b/src/rpucuda/rpu_vector_device.cpp
index 4ec6a8ba..ba3f32eb 100644
--- a/src/rpucuda/rpu_vector_device.cpp
+++ b/src/rpucuda/rpu_vector_device.cpp
@@ -414,7 +414,15 @@ void VectorRPUDevice<T>::initUpdateCycle(
   }
   default: {
   }
-  };
+  }
+
+  if (par.singleDeviceUpdate()) {
+    rpu_device_vec_[current_device_idx_]->initUpdateCycle(weights, up, current_lr, m_batch_info);
+  } else {
+    for (size_t k = 0; k < rpu_device_vec_.size(); k++) {
+      rpu_device_vec_[k]->initUpdateCycle(weights, up, current_lr, m_batch_info);
+    }
+  }
 }
 
 template <typename T>
@@ -448,6 +456,15 @@ void VectorRPUDevice<T>::doSparseUpdate(
 template <typename T>
 void VectorRPUDevice<T>::finishUpdateCycle(
     T **weights, const PulsedUpdateMetaParameter<T> &up, T current_lr, int m_batch_info) {
+
+  if (getPar().singleDeviceUpdate()) {
+    rpu_device_vec_[current_device_idx_]->finishUpdateCycle(weights, up, current_lr, m_batch_info);
+  } else {
+    for (size_t k = 0; k < rpu_device_vec_.size(); k++) {
+      rpu_device_vec_[k]->finishUpdateCycle(weights, up, current_lr, m_batch_info);
+    }
+  }
+
   // count updates
   current_update_idx_++;
 }


Test patch:
diff --git a/src/rpucuda/cuda/io_iterator_test.cpp b/src/rpucuda/cuda/io_iterator_test.cpp
index 08870e5e..d3d801e5 100644
--- a/src/rpucuda/cuda/io_iterator_test.cpp
+++ b/src/rpucuda/cuda/io_iterator_test.cpp
@@ -19,8 +19,8 @@
 #include <chrono>
 #include <functional>
 #include <memory>
-#include <random>
 #include <numeric>
+#include <random>
 
 #define TOLERANCE 1e-5
 
diff --git a/src/rpucuda/cuda/rpucuda_pulsed_test.cpp b/src/rpucuda/cuda/rpucuda_pulsed_test.cpp
index 236ada92..96c40e83 100644
--- a/src/rpucuda/cuda/rpucuda_pulsed_test.cpp
+++ b/src/rpucuda/cuda/rpucuda_pulsed_test.cpp
@@ -21,8 +21,8 @@
 #include <algorithm>
 #include <chrono>
 #include <memory>
-#include <random>
 #include <numeric>
+#include <random>
 
 #define TOLERANCE 1e-5