Repository: IBM/aihwkit (language: python)

Issue:
Softbound synapses based on Fusi Abott modell
## Description and motivation
This may be caused by my lack of mathematical expertise and a quick solve is obvious. 
The soft-bound synapse model described by Fusi and Abbott [1] has been adopted in models for synaptic resistive devices, e.g. in the detailed study of Frascaroli et al. [2]. In their paper, they give an equation for the weight dependent weight update for potentiation and depression (Eqn. (1) and (2)) as well as respective equations for potentiation and depression in dependence of the pulse number ( Eqn. (3) and (4)). See figures below.

Now to my questions: 
- Is it currently possible to rewrite these equations so they fit in the framework of device models?
- As the model is defined for weights between 0 and 1, one would need to chose the DifferenceUnitCell module to obtain weights between -1 and 1, is that correct? 

I hope this is clear. I appreciate any help!

[1] S. Fusi; L. F. Abbott, Nature Neuroscience, Limits on the memory storage capacity of bounded synapses, 2007 (DOI: 10.1038/nn1859 )
[2] J. Frascaroli; S. Brivio; E. Covi; S. Spiga, Scientific Reports, Evidence of soft bound behaviour in analogue memristive devices for neuromorphic  computing, 2018 (DOI: 10.1038/s41598-018-25376-x )
![Frascaroli_1_2](https://user-images.githubusercontent.com/81355788/112631057-b4585c80-8e36-11eb-8e4b-51fd3fc1a650.PNG)
![Frascaroli_3_4](https://user-images.githubusercontent.com/81355788/112631060-b4f0f300-8e36-11eb-9c3e-bfdfba49ef89.PNG)




Discussion / hints:
Thank you very much for your question! Currently, we do not support this update behavior. We only have a `SoftBoundsDevice` mode, where the update behavior is similiar to (1) and (2) with `gamma=1`, that is the update step size changes linear with weight. 
You won't need `DifferenceUnitCell`, in our models the conductance range is scaled from -|w_min| to w_max, which can be realized with a fixed bias (reference cell) that we usually do not model explicitly (it can be modeled with `ReferenceUnitCell` if needed explicitely).  The `DifferenceUnitCell` is only for devices that are not bi-directional, such as PCM, where one device receives all positive updates and the other the negative updates and a reset and refresh would be needed to avoid saturation.  

That said, implementing your suggested update behavior as a new device model would be relatively straightforward. I'll see whether it can be done relatively quickly.  

Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 8cf64dfa..21205965 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -23,6 +23,9 @@ The format is based on [Keep a Changelog], and this project adheres to
       high precision storage) and then transfers the matrix sequentially to
       the analog device, instead of using the default fully parallel pulsed
       update. (\#159)
+    * A new devicen model class `PowStepDevice` that implements a
+      power-exponent type of non-linearity based on the Fusi & Abott
+      synapse model. (\#192)
     * Option to choose deterministic pulse trains for the rank-1 update of
       analog devices during training. (\#99)
     * More noise types for hardware-aware training for inference
@@ -74,6 +77,7 @@ The format is based on [Keep a Changelog], and this project adheres to
 
 #### Fixed
 
+* Issue of numnber og loop estimations for realistics read (\#192)
 * Fixed small issues that resulted in warnings for windows compilation. (\#99)
 * Faulty backward noise management error message removed for perfect backward
   and CUDA. (\#99)
diff --git a/docs/source/using_simulator.rst b/docs/source/using_simulator.rst
index e1efe507..c45ff1fd 100644
--- a/docs/source/using_simulator.rst
+++ b/docs/source/using_simulator.rst
@@ -108,7 +108,9 @@ Resistive device class                                            Description
 :class:`~aihwkit.simulator.configs.devices.ConstantStepDevice`    pulsed update behavioral model: constant step, where the update step of material is constant throughout the resistive range (up to hard bounds).
 :class:`~aihwkit.simulator.configs.devices.LinearStepDevice`      pulsed update behavioral model: linear step, where the update step response size of the material is linearly dependent with resistance (up to hard bounds).
 :class:`~aihwkit.simulator.configs.devices.SoftBoundsDevice`      pulsed update behavioral model: soft bounds, where the update step response size of the material is linearly dependent and it goes to zero at the bound.
+:class:`~aihwkit.simulator.configs.devices.SoftBoundsPmaxDevice`  same model as in :class:`~aihwkit.simulator.configs.devices.SoftBoundsDevice` but using a more convenient parameterization for easier fits to experimentally measured update response curves.
 :class:`~aihwkit.simulator.configs.devices.ExpStepDevice`         exponential update step or CMOS-like update behavior.
+:class:`~aihwkit.simulator.configs.devices.PowStepDevice`         update step using a power exponent non-linearity.
 ================================================================  ========
 
 Unit cell devices
@@ -118,7 +120,7 @@ Unit cell devices
 Resistive device class                                                Description
 ====================================================================  ========
 :class:`~aihwkit.simulator.configs.devices.VectorUnitCell`            abstract resistive device that combines multiple pulsed resistive devices in a single 'unit cell'.
-:class:`~aihwkit.simulator.configs.devices.DifferenceUnitCell`        abstract device model takes an arbitrary device per crosspoint and implements an explicit plus-minus device pair.
+:class:`~aihwkit.simulator.configs.devices.DifferenceUnitCell`        abstract device model that takes an arbitrary device per crosspoint and implements an one-sided plus-minus device pair, where both devices only receive only positive updates and a refresh or reset needs to be added to avoid saturation.
 :class:`~aihwkit.simulator.configs.devices.ReferenceUnitCell`         abstract device model takes two arbitrary device per cross-point and implements an device with reference pair.
 ====================================================================  ========
 
@@ -129,6 +131,7 @@ Compound devices
 Resistive device class                                                Description
 ====================================================================  ========
 :class:`~aihwkit.simulator.configs.devices.TransferCompound`          abstract device model that takes 2 or more devices per crosspoint and implements a 'transfer' based learning rule such as Tiki-Taka (see `Gokmen & Haensch 2020`_).
+:class:`~aihwkit.simulator.configs.devices.MixedPrecisionCompound`    abstract device model that takes one devices per crosspoint and implements a 'mixed-precision' based learning rule where the rank-update is done in digital instead of using a fully analog parallel write (see `Nandakumar et al. 2020`_).
 ====================================================================  ========
 
 RPU Configurations
@@ -330,6 +333,43 @@ rule instead of plain SGD. Once the configuration is done, the usage
 of this complex analog tile for testing or training from the user
 point of view is however the same as for other tiles.
 
+Mixed Precision Compound
+-------------------------
+
+This abstract device implements an analog SGD optimizer suggested by
+`Nandakumar et al. 2020`_ where the update is not done in analog
+directly, but in digital. Thus is uses a digital rank-update of an
+intermediately stored floating point matrix, which will be used to
+transfer the information to the analog tile that is used in forward
+and backward pass.  This optimizer strategy is in contrast with the
+default mode in the simulator, that uses stochastic pulse trains to
+update in parallel onto the analog tile directly. This will have
+impact on the hardware design as well as expected runtime, as more
+digital computation is needed to be done. For details, see `Nandakumar
+et al. 2020`_.
+
+To enable mixed-precision one defines for example the following ``rpu_config``::
+
+    # Imports from aihwkit.
+    from aihwkit.simulator.configs import DigitalRankUpdateRPUConfig
+    from aihwkit.simulator.configs.devices import SoftBoundsDevice
+
+    rpu_config = DigitalRankUpdateRPUConfig(
+        device=SoftBoundsDevice(),
+
+	# make some adjustments of mixed-precision hyper parameter
+        granularity=0.001,
+	n_x_bins=15,
+	n_d_bins=31,
+    )
+
+    # use tile configuration in model
+    model = AnalogLinear(4, 2, bias=True, rpu_config=rpu_config)
+
+Now this analog tile will use the mixed-precision optimizer with a
+soft bounds device model. 
+
+
 Analog presets
 --------------
 
@@ -354,3 +394,4 @@ instead of manually specifying a ``RPU Configuration``::
 .. _Gokmen & Haensch 2020: https://www.frontiersin.org/articles/10.3389/fnins.2020.00103/full
 .. _Example 7: https://github.com/IBM/aihwkit/blob/master/examples/07_simple_layer_with_other_devices.py
 .. _Example 8: https://github.com/IBM/aihwkit/blob/master/examples/08_simple_layer_with_tiki_taka.py
+.. _Nandakumar et al. 2020: https://www.frontiersin.org/articles/10.3389/fnins.2020.00406/full
diff --git a/src/aihwkit/simulator/configs/devices.py b/src/aihwkit/simulator/configs/devices.py
index 97966820..79ceaa47 100644
--- a/src/aihwkit/simulator/configs/devices.py
+++ b/src/aihwkit/simulator/configs/devices.py
@@ -201,14 +201,14 @@ class PulsedDevice(_PrintableMixin):
     also vary across devices.
 
     :math:`\Delta w_{ij}^d` is set during RPU initialization (for each
-    cross-point `ij`):
+    cross-point :math:`ij`):
 
     .. math::
 
         \Delta w_{ij}^d = d\; \Delta w_\text{min}\, \left(
         1 + d \beta_{ij} + \sigma_\text{d-to-d}\xi\right)
 
-    where \xi is again a standard Gaussian. :math:`\beta_{ij}` is the
+    where :math:`\xi` is again a standard Gaussian. :math:`\beta_{ij}` is the
     directional up `versus` down bias.  At initialization ``up_down_dtod`` and
     ``up_down`` defines this bias term:
 
@@ -217,7 +217,7 @@ class PulsedDevice(_PrintableMixin):
         \beta_{ij} = \beta_\text{up-down} + \xi
         \sigma_\text{up-down-dtod}
 
-    where \xi is again a standard Gaussian number and
+    where :math:`\xi` is again a standard Gaussian number and
     :math:`\beta_\text{up-down}` corresponds to ``up_down``. Note that
     ``up_down_dtod`` is again given in relative units to ``dw_min``.
     """
@@ -665,6 +665,120 @@ class ExpStepDevice(PulsedDevice):
     """
 
 
+@dataclass
+class PowStepDevice(PulsedDevice):
+    r"""Pulsed update behavioral model: power-dependent step.
+
+    Pulsed update behavioral model, where the update step response
+    size of the material has a power-dependent with resistance. This
+    device model implements (a shifted from of) the `Fusi & Abott
+    (2007)`_ synapse model (see also `Frascaroli et al. (2108)`_).
+
+    The model based on :class:`~PulsedDevice` and thus shares most
+    parameters and functionality. However, it implements new `update
+    once` function, where the update step size depends in the
+    following way. If we set :math:`\omega_{ij} =
+    \frac{b_{ij}^\text{max} - w_{ij}}{b_{ij}^\text{max} -
+    b_{ij}^\text{min}}` the relative distance of the current weight to
+    the upper bound, then the update per pulse is for the upwards direction:
+
+    .. math::
+       w_{ij}  \leftarrow  w_{ij} + \Delta w_{ij}^+\,(\omega_{ij})^{\gamma_{ij}^+}
+       \left(1 + \sigma_\text{c-to-c}\,\xi\right)
+
+    and in downwards direction:
+
+    .. math::
+       w_{ij}  \leftarrow  w_{ij} + \Delta w_{ij}^-\,(1 - \omega_{ij})^{\gamma_{ij}^-}
+       \left(1 + \sigma_\text{c-to-c}\,\xi\right)
+
+    Similar to :math:`\Delta w_{ij}^d` the exponent :math:`\gamma_{ij}` can be
+    defined with device-to-device variation and bias in up and down
+    direction:
+
+    .. math::
+
+        \gamma_{ij}^d = d\; \gamma\, \left(1 + d\, \beta_{ij}
+        + \sigma_\text{pow-gamma-d-to-d}\xi\right)
+
+    where :math:`\xi` is again a standard Gaussian. :math:`\beta_{ij}`
+    is the directional up `versus` down bias.  At initialization
+    ``pow_up_down_dtod`` and ``pow_up_down`` defines this bias term:
+
+    .. math::
+
+        \beta_{ij} = \beta_\text{pow-up-down} + \xi\sigma_\text{pow-up-down-dtod}
+
+    where :math:`\xi` is again a standard Gaussian number and
+    :math:`\beta_\text{pow-up-down}` corresponds to ``pow_up_down``.
+
+    Note:
+        The ``pow_gamma_dtod`` and ``pow_up_down_dtod``
+        device-to-device variation parameters are given in relative
+        units to ``pow_gamma``.
+
+    Note:
+        :math:`\Delta w_{ij}^d` is defined as for the
+        :class:`~PulsedDevice`, however, for this device, the update step
+        size will *not* be given by :math:`\Delta w_{ij}` at
+        :math:`w_{ij}=0` as for most other devices models
+
+    ..  _Fusi & Abott (2007): https://www.nature.com/articles/nn1859
+    ..  _Frascaroli et al. (2108): https://www.nature.com/articles/s41598-018-25376-x
+    """
+
+    bindings_class: ClassVar[Type] = devices.PowStepResistiveDeviceParameter
+
+    pow_gamma: float = 1.0
+    r"""The value of :math:`\gamma` as explained above.
+
+    Note:
+        :math:`\gamma` reduces essentially to the
+        :class:`SoftBoundsDevice` (if no device-to-device variation of
+        gamma is used additionally). However, the
+        :class:`SoftBoundsDevice` will be much faster, as it does not
+        need to compute the slow `pow` function.
+    """
+
+    pow_gamma_dtod: float = 0.1
+    r"""Device-to-device variation for ``pow_gamma``.
+
+    i.e. the value of :math:`\gamma_\text{pow-gamma-d-to-d}` given in relative
+    units to ``pow_gamma``.
+    """
+
+    pow_up_down: float = 0.0
+    r"""The up versus down bias of the :math:`\gamma` as described above.
+
+    It is :math:`\gamma^+ = \gamma (1 + \beta_\text{pow-up-down})` and
+    :math:`\gamma^- = \gamma (1 - \beta_\text{pow-up-down})` .
+    """
+
+    pow_up_down_dtod: float = 0.0
+    r"""Device-to-device variation in the up versus down bias of
+    :math:`\gamma` as descibed above.
+
+    In units of ``pow_gamma``.
+    """
+
+    write_noise_std: float = 0.0
+    r"""Whether to use update write noise.
+
+    Whether to use update write noise that is added to the updated
+    devices weight, while the update is done on a hidden persistent weight. The
+    update write noise is then sampled a new when the device is touched
+    again.
+
+    Thus it is:
+
+    .. math::
+       w_\text{apparent}{ij} = w_ij + \sigma_\text{write_noise}\xi
+
+    and the update is done on :math:`w_ij` but the forward sees the
+    :math:`w_\text{apparent}`.
+    """
+
+
 ###############################################################################
 # Specific devices based on ``unit cell``.
 ###############################################################################
diff --git a/src/aihwkit/simulator/rpu_base_src/rpu_base.h b/src/aihwkit/simulator/rpu_base_src/rpu_base.h
index 0e71a9be..8317bd46 100644
--- a/src/aihwkit/simulator/rpu_base_src/rpu_base.h
+++ b/src/aihwkit/simulator/rpu_base_src/rpu_base.h
@@ -17,6 +17,7 @@
 #include "rpu_linearstep_device.h"
 #include "rpu_mixedprec_device.h"
 #include "rpu_mixedprec_device_base.h"
+#include "rpu_powstep_device.h"
 #include "rpu_pulsed.h"
 #include "rpu_simple_device.h"
 #include "rpu_transfer_device.h"
diff --git a/src/aihwkit/simulator/rpu_base_src/rpu_base_devices.cpp b/src/aihwkit/simulator/rpu_base_src/rpu_base_devices.cpp
index 28c3e4d0..8d56a3d5 100644
--- a/src/aihwkit/simulator/rpu_base_src/rpu_base_devices.cpp
+++ b/src/aihwkit/simulator/rpu_base_src/rpu_base_devices.cpp
@@ -26,6 +26,7 @@ void declare_rpu_devices(py::module &m) {
   using DifferenceParam = RPU::DifferenceRPUDeviceMetaParameter<T>;
   using TransferParam = RPU::TransferRPUDeviceMetaParameter<T>;
   using MixedPrecParam = RPU::MixedPrecRPUDeviceMetaParameter<T>;
+  using PowStepParam = RPU::PowStepRPUDeviceMetaParameter<T>;
 
   /*
    * Trampoline classes for allowing inheritance.
@@ -230,6 +231,24 @@ void declare_rpu_devices(py::module &m) {
     }
   };
 
+  class PyPowStepParam : public PowStepParam {
+  public:
+    std::string getName() const override {
+      PYBIND11_OVERLOAD(std::string, PowStepParam, getName, );
+    }
+    PowStepParam *clone() const override {
+      PYBIND11_OVERLOAD(PowStepParam *, PowStepParam, clone, );
+    }
+    RPU::DeviceUpdateType implements() const override {
+      PYBIND11_OVERLOAD(RPU::DeviceUpdateType, PowStepParam, implements, );
+    }
+    RPU::PowStepRPUDevice<T> *
+    createDevice(int x_size, int d_size, RPU::RealWorldRNG<T> *rng) override {
+      PYBIND11_OVERLOAD(
+          RPU::PowStepRPUDevice<T> *, PowStepParam, createDevice, x_size, d_size, rng);
+    }
+  };
+
   /*
    * Python class definitions.
    */
@@ -488,6 +507,19 @@ void declare_rpu_devices(py::module &m) {
         return ss.str();
       });
 
+  py::class_<PowStepParam, PyPowStepParam, PulsedParam>(m, "PowStepResistiveDeviceParameter")
+      .def(py::init<>())
+      .def_readwrite("pow_gamma", &PowStepParam::ps_gamma)
+      .def_readwrite("pow_gamma_dtod", &PowStepParam::ps_gamma_dtod)
+      .def_readwrite("pow_up_down", &PowStepParam::ps_gamma_up_down)
+      .def_readwrite("pow_up_down_dtod", &PowStepParam::ps_gamma_up_down_dtod)
+      .def_readwrite("write_noise_std", &PowStepParam::write_noise_std)
+      .def("__str__", [](PowStepParam &self) {
+        std::stringstream ss;
+        self.printToStream(ss);
+        return ss.str();
+      });
+
   /**
    * Helper enums.
    **/
diff --git a/src/rpucuda/cuda/rpucuda_linearstep_device.cu b/src/rpucuda/cuda/rpucuda_linearstep_device.cu
index c5bf72fc..224e84eb 100644
--- a/src/rpucuda/cuda/rpucuda_linearstep_device.cu
+++ b/src/rpucuda/cuda/rpucuda_linearstep_device.cu
@@ -32,6 +32,8 @@ template <typename T> struct UpdateFunctorLinearStepMult {
     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 wmax = par_4.z;                                          // [2]
+    T wmin = par_4.x;                                          // [0]
 
     T &w = uw_std > 0 ? persistent_weight : apparent_weight;
 
@@ -44,24 +46,27 @@ template <typename T> struct UpdateFunctorLinearStepMult {
       } else {
         w += lin_a * w + lin_dw;
       }
+      w = (w > wmax) ? wmax : w;
+      w = (w < wmin) ? wmin : w;
     } else {
       if (noise_std_dw > 0) {
         for (int i_updates = 0; i_updates < n; i_updates++) {
           T stoch_value = curand_normal(&local_state);
           stoch_value *= noise_std_dw;
           w += (lin_a * w + lin_dw) * (1.0 + stoch_value);
+          // better always check both bounds
+          w = (w > wmax) ? wmax : w;
+          w = (w < wmin) ? wmin : w;
         }
       } else {
         for (int i_updates = 0; i_updates < n; i_updates++) {
           w += lin_a * w + lin_dw;
+          // better always check both bounds
+          w = (w > wmax) ? wmax : w;
+          w = (w < wmin) ? wmin : w;
         }
       }
     }
-    // better always check both bounds
-    T wmax = par_4.z; // [2]
-    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) {
@@ -86,8 +91,9 @@ template <typename T> struct UpdateFunctorLinearStepAdd {
     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;
+    T wmax = par_4.z; // [2]
+    T wmin = par_4.x; // [0]
 
     // n is larger 0 in any case
     if (n == 1) {
@@ -98,24 +104,25 @@ template <typename T> struct UpdateFunctorLinearStepAdd {
       } else {
         w += lin_a * w + lin_dw;
       }
+      w = (w > wmax) ? wmax : w;
+      w = (w < wmin) ? wmin : w;
     } else {
       if (noise_std_dw > 0) {
         for (int i_updates = 0; i_updates < n; i_updates++) {
           T stoch_value = curand_normal(&local_state);
           stoch_value *= noise_std_dw;
           w += lin_a * w + lin_dw * (1.0 + stoch_value);
+          w = (w > wmax) ? wmax : w;
+          w = (w < wmin) ? wmin : w;
         }
       } else {
         for (int i_updates = 0; i_updates < n; i_updates++) {
           w += lin_a * w + lin_dw;
+          w = (w > wmax) ? wmax : w;
+          w = (w < wmin) ? wmin : w;
         }
       }
     }
-    T wmax = par_4.z; // [2]
-    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);
diff --git a/src/rpucuda/cuda/rpucuda_powstep_device.cu b/src/rpucuda/cuda/rpucuda_powstep_device.cu
new file mode 100644
index 00000000..efcaee75
--- /dev/null
+++ b/src/rpucuda/cuda/rpucuda_powstep_device.cu
@@ -0,0 +1,127 @@
+/**
+ * (C) Copyright 2020, 2021 IBM. All Rights Reserved.
+ *
+ * This code is licensed under the Apache License, Version 2.0. You may
+ * obtain a copy of this license in the LICENSE.txt file in the root directory
+ * of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+ *
+ * Any modifications or derivative works of this code must retain this
+ * copyright notice, and modified files need to carry a notice indicating
+ * that they have been altered from the originals.
+ */
+
+#include "pwu_kernel_parameter.h"
+#include "rpu_pulsed_meta_parameter.h"
+#include "rpucuda_powstep_device.h"
+#include <memory>
+
+namespace RPU {
+
+template <typename T> struct UpdateFunctorPowStep {
+
+  __device__ __forceinline__ void operator()(
+      T &apparent_weight,
+      uint32_t n,
+      uint32_t negative,
+      const float4 par_4,         // min_bound, scale_down, max_bound, scale_up,
+      const float2 gamma_down_up, // gamma_down, gamma_up
+      T &persistent_weight,
+      const T *write_noise_std,
+      T noise_std_dw,
+      curandState &local_state) {
+
+    T wmin = par_4.x; // [0]
+    T wmax = par_4.z; // [2]
+    T range = wmax - wmin;
+    if (range == 0) {
+      return;
+    }
+    T uw_std = *write_noise_std;
+    T &w = uw_std > 0 ? persistent_weight : apparent_weight;
+    // negative > 0 means sign < 0 and thus up-direction
+    T scale = (negative > 0) ? (par_4.w) : (-par_4.y);                //[3] (up), [1] (down)
+    T gamma = (negative > 0) ? (gamma_down_up.y) : (gamma_down_up.x); // [1] (up), [0] (down)
+
+    // up direction: ((wmax - w) / range) ^ gamma
+    // down direction:  ((w - wmin) / range) ^ gamma  == (1 - (wmax - w)/range) ^ gamma
+
+    // n is larger 0 in any case
+    if (n == 1) {
+      T x = (wmax - w) / range;
+      T dw = scale * ((negative > 0) ? __powf(x, gamma) : __powf((T)1.0 - x, gamma));
+
+      if (noise_std_dw > 0) {
+        T stoch_value = curand_normal(&local_state);
+        stoch_value *= noise_std_dw;
+        w += dw * ((T)1.0 + stoch_value);
+      } else {
+        w += dw;
+      }
+      w = (w > wmax) ? wmax : w;
+      w = (w < wmin) ? wmin : w;
+
+    } else {
+      if (noise_std_dw > 0) {
+        for (int i_updates = 0; i_updates < n; i_updates++) {
+          T stoch_value = curand_normal(&local_state);
+          stoch_value *= noise_std_dw;
+          T x = (wmax - w) / range;
+          T dw = scale * ((negative > 0) ? __powf(x, gamma) : __powf((T)1.0 - x, gamma));
+          w += dw * ((T)1.0 + stoch_value);
+          // better always check both bounds
+          w = (w > wmax) ? wmax : w;
+          w = (w < wmin) ? wmin : w;
+        }
+      } else {
+        for (int i_updates = 0; i_updates < n; i_updates++) {
+          T x = (wmax - w) / range;
+          w += scale * ((negative > 0) ? __powf(x, gamma) : __powf((T)1.0 - x, gamma));
+          // better always check both bounds
+          w = (w > wmax) ? wmax : w;
+          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;
+    }
+  }
+};
+
+#define ARGS                                                                                       \
+  (this->context_, this->x_size_, this->d_size_, m_batch, nK32, use_bo64, out_trans, up,           \
+   par.getName())
+
+template <typename T>
+pwukpvec_t<T> PowStepRPUDeviceCuda<T>::getUpdateKernels(
+    int m_batch, int nK32, int use_bo64, bool out_trans, const PulsedUpdateMetaParameter<T> &up) {
+
+  pwukpvec_t<T> v;
+  const auto &par = getPar();
+  v.push_back(
+      RPU::make_unique<PWUKernelParameterSingleFunctor<T, UpdateFunctorPowStep<T>, 1>> ARGS);
+  v.push_back(RPU::make_unique<PWUKernelParameterBatchFunctor<T, UpdateFunctorPowStep<T>, 1>> ARGS);
+  v.push_back(
+      RPU::make_unique<PWUKernelParameterBatchSharedFunctor<T, UpdateFunctorPowStep<T>, 1>> ARGS);
+  return v;
+}
+
+#undef ARGS
+
+template <typename T> T PowStepRPUDeviceCuda<T>::getDwMin() const {
+
+  // TODO: this should actually be just be a property from rpu_device. Need to do this in Base.
+  const auto &par = getPar();
+  // just approximately at zero (used for update management)
+  return par.dw_min * pow(0.5, par.ps_gamma);
+};
+
+template class PowStepRPUDeviceCuda<float>;
+#ifdef RPU_USE_DOUBLE
+template class PowStepRPUDeviceCuda<double>;
+#endif
+
+} // namespace RPU
diff --git a/src/rpucuda/cuda/rpucuda_powstep_device.h b/src/rpucuda/cuda/rpucuda_powstep_device.h
new file mode 100644
index 00000000..235e9a61
--- /dev/null
+++ b/src/rpucuda/cuda/rpucuda_powstep_device.h
@@ -0,0 +1,84 @@
+/**
+ * (C) Copyright 2020, 2021 IBM. All Rights Reserved.
+ *
+ * This code is licensed under the Apache License, Version 2.0. You may
+ * obtain a copy of this license in the LICENSE.txt file in the root directory
+ * of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+ *
+ * Any modifications or derivative works of this code must retain this
+ * copyright notice, and modified files need to carry a notice indicating
+ * that they have been altered from the originals.
+ */
+
+#pragma once
+
+#include "pwu_kernel_parameter_base.h"
+#include "rpu_powstep_device.h"
+#include "rpucuda_pulsed_device.h"
+#include <memory>
+
+namespace RPU {
+
+template <typename T> class PowStepRPUDeviceCuda : public PulsedRPUDeviceCuda<T> {
+
+public:
+  BUILD_PULSED_DEVICE_CONSTRUCTORS_CUDA(
+      PowStepRPUDeviceCuda,
+      PowStepRPUDevice,
+      /*ctor body*/
+      dev_gamma_ = 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_gamma_->assign(*other.dev_gamma_);
+      dev_write_noise_std_->assign(*other.dev_write_noise_std_);
+      ,
+      /*move assigment body*/
+      dev_gamma_ = std::move(other.dev_gamma_);
+      dev_write_noise_std_ = std::move(other.dev_write_noise_std_);
+      ,
+      /*swap body*/
+      swap(a.dev_gamma_, b.dev_gamma_);
+      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_;
+      int x_size = this->x_size_;
+      T **w_gamma_up = rpu_device.getGammaUp();
+      T **w_gamma_down = rpu_device.getGammaDown();
+      float *tmp_gamma = new float[2 * this->size_];
+
+      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_gamma[kk] = w_gamma_down[i][j];
+          tmp_gamma[kk + 1] = w_gamma_up[i][j];
+        }
+      } dev_gamma_->assign(tmp_gamma);
+      dev_write_noise_std_->setConst(getPar().getScaledWriteNoise());
+
+      this->context_->synchronize();
+      delete[] tmp_gamma;);
+
+  pwukpvec_t<T> getUpdateKernels(
+      int m_batch,
+      int nK32,
+      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_gamma_->getData(); };
+  T *get1ParamsData() override {
+    return getPar().usesPersistentWeight() ? this->dev_persistent_weights_->getData() : nullptr;
+  };
+
+  T getDwMin() const override;
+
+private:
+  std::unique_ptr<CudaArray<float>> dev_gamma_ = nullptr;
+  std::unique_ptr<CudaArray<T>> dev_write_noise_std_ = nullptr;
+};
+
+} // namespace RPU
diff --git a/src/rpucuda/cuda/rpucuda_pulsed.cu b/src/rpucuda/cuda/rpucuda_pulsed.cu
index dd37259e..82f38245 100644
--- a/src/rpucuda/cuda/rpucuda_pulsed.cu
+++ b/src/rpucuda/cuda/rpucuda_pulsed.cu
@@ -340,7 +340,7 @@ template <typename T> void RPUCudaPulsed<T>::setWeightsReal(const T *weightsptr,
   if (dpar != nullptr) {
     w_min = dpar->w_min;
     w_max = dpar->w_max;
-    static_cast<PulsedRPUDeviceCudaBase<T> &>(*rpucuda_device_).getDwMin();
+    dw_min = static_cast<PulsedRPUDeviceCudaBase<T> &>(*rpucuda_device_).getDwMin();
   }
   T A = 0;
   T B = 0;
diff --git a/src/rpucuda/cuda/rpucuda_simple_device.cu b/src/rpucuda/cuda/rpucuda_simple_device.cu
index bee034e4..f0ea7469 100644
--- a/src/rpucuda/cuda/rpucuda_simple_device.cu
+++ b/src/rpucuda/cuda/rpucuda_simple_device.cu
@@ -20,6 +20,7 @@
 #include "rpucuda_expstep_device.h"
 #include "rpucuda_linearstep_device.h"
 #include "rpucuda_mixedprec_device.h"
+#include "rpucuda_powstep_device.h"
 #include "rpucuda_transfer_device.h"
 #include "rpucuda_vector_device.h"
 
@@ -53,6 +54,8 @@ AbstractRPUDeviceCuda<T>::createFrom(CudaContext *c, const AbstractRPUDevice<T>
     return new SimpleRPUDeviceCuda<T>(c, static_cast<const SimpleRPUDevice<T> &>(rpu_device));
   case DeviceUpdateType::MixedPrec:
     return new MixedPrecRPUDeviceCuda<T>(c, static_cast<const MixedPrecRPUDevice<T> &>(rpu_device));
+  case DeviceUpdateType::PowStep:
+    return new PowStepRPUDeviceCuda<T>(c, static_cast<const PowStepRPUDevice<T> &>(rpu_device));
   default:
     RPU_FATAL("Pulsed device type not implemented in CUDA. Maybe not added to createFrom in "
               "rpucuda_simple_device.cu?");
diff --git a/src/rpucuda/rpu_linearstep_device.cpp b/src/rpucuda/rpu_linearstep_device.cpp
index c61adf6b..2e4030f8 100644
--- a/src/rpucuda/rpu_linearstep_device.cpp
+++ b/src/rpucuda/rpu_linearstep_device.cpp
@@ -98,6 +98,7 @@ template <typename T> void LinearStepRPUDevice<T>::printDP(int x_count, int d_co
   }
 }
 
+namespace {
 template <typename T>
 inline void update_once_mult(
     T &w,
@@ -151,6 +152,7 @@ inline void update_once_add(
     w_apparent = w + write_noise_std * rng->sampleGauss();
   }
 }
+} // namespace
 
 template <typename T>
 void LinearStepRPUDevice<T>::doSparseUpdate(
diff --git a/src/rpucuda/rpu_linearstep_device.h b/src/rpucuda/rpu_linearstep_device.h
index 8026ac59..df149a02 100644
--- a/src/rpucuda/rpu_linearstep_device.h
+++ b/src/rpucuda/rpu_linearstep_device.h
@@ -87,8 +87,8 @@ struct SoftBoundsRPUDeviceMetaParameter : LinearStepRPUDeviceMetaParameter<T> {
   void printToStream(std::stringstream &ss) const override {
     checkSoftBounds();
     PulsedRPUDeviceMetaParameter<T>::printToStream(ss);
-    ss << "   ";
-    ss << "SoftBounds parameter:" << std::endl;
+    // ss << "   ";
+    // ss << getName() << " parameter:" << std::endl;
     ss << "\t ls_mult_noise:\t\t\t" << std::boolalpha << this->ls_mult_noise << std::endl;
   };
 };
diff --git a/src/rpucuda/rpu_powstep_device.cpp b/src/rpucuda/rpu_powstep_device.cpp
new file mode 100644
index 00000000..3da05dba
--- /dev/null
+++ b/src/rpucuda/rpu_powstep_device.cpp
@@ -0,0 +1,183 @@
+/**
+ * (C) Copyright 2020, 2021 IBM. All Rights Reserved.
+ *
+ * This code is licensed under the Apache License, Version 2.0. You may
+ * obtain a copy of this license in the LICENSE.txt file in the root directory
+ * of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+ *
+ * Any modifications or derivative works of this code must retain this
+ * copyright notice, and modified files need to carry a notice indicating
+ * that they have been altered from the originals.
+ */
+
+#include "rpu_powstep_device.h"
+#include "utility_functions.h"
+#include <chrono>
+#include <cmath>
+#include <iostream>
+#include <limits>
+
+namespace RPU {
+
+/********************************************************************************
+ * Linear Step RPU Device
+ *********************************************************************************/
+
+template <typename T>
+void PowStepRPUDevice<T>::populate(
+    const PowStepRPUDeviceMetaParameter<T> &p, RealWorldRNG<T> *rng) {
+
+  PulsedRPUDevice<T>::populate(p, rng); // will clone par
+  auto &par = getPar();
+
+  T gamma = par.ps_gamma;
+  T gain_std = par.ps_gamma_dtod;
+  T up_down_std = par.ps_gamma_up_down_dtod;
+  T up_down = par.ps_gamma_up_down;
+
+  T up_bias = up_down > 0 ? (T)0.0 : up_down;
+  T down_bias = up_down > 0 ? -up_down : (T)0.0;
+
+  for (int i = 0; i < this->d_size_; ++i) {
+    for (int j = 0; j < this->x_size_; ++j) {
+
+      T gain = (T)1.0 + gain_std * rng->sampleGauss();
+      T r = up_down_std * rng->sampleGauss();
+
+      w_gamma_up_[i][j] = (up_bias + gain + r) * gamma;
+      w_gamma_down_[i][j] = (down_bias + gain - r) * gamma;
+
+      if (par.enforce_consistency) {
+        w_gamma_up_[i][j] = fabs(w_gamma_up_[i][j]);
+        w_gamma_down_[i][j] = fabs(w_gamma_down_[i][j]);
+      }
+    }
+  }
+}
+
+template <typename T> void PowStepRPUDevice<T>::printDP(int x_count, int d_count) const {
+
+  if (x_count < 0 || x_count > this->x_size_) {
+    x_count = this->x_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) {
+      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 << 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.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;
+  }
+}
+
+namespace {
+template <typename T>
+inline void update_once(
+    T &w,
+    T &w_apparent,
+    int &sign,
+    T &scale_down,
+    T &scale_up,
+    T &gamma_down,
+    T &gamma_up,
+    T &min_bound,
+    T &max_bound,
+    const T &dw_min_std,
+    const T &write_noise_std,
+    RNG<T> *rng) {
+  T range = max_bound - min_bound;
+  if (range == 0.0) {
+    return;
+  }
+  if (sign > 0) {
+    w -= scale_down * pow((w - min_bound) / range, gamma_down) *
+         ((T)1.0 + dw_min_std * rng->sampleGauss());
+  } else {
+    w += scale_up * pow((max_bound - w) / range, gamma_up) *
+         ((T)1.0 + dw_min_std * rng->sampleGauss());
+  }
+  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();
+  }
+}
+
+} // namespace
+
+template <typename T> T PowStepRPUDevice<T>::getDwMin() const {
+  const auto &par = getPar();
+  // just approximately at zero (used for update management)
+  return par.dw_min * pow(0.5, par.ps_gamma);
+};
+
+template <typename T>
+void PowStepRPUDevice<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 *gamma_down = w_gamma_down_[i];
+  T *gamma_up = w_gamma_up_[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];
+
+  T write_noise_std = par.getScaledWriteNoise();
+  PULSED_UPDATE_W_LOOP(update_once(
+                           w[j], w_apparent[j], sign, scale_down[j], scale_up[j], gamma_down[j],
+                           gamma_up[j], min_bound[j], max_bound[j], par.dw_min_std, write_noise_std,
+                           rng););
+}
+
+template <typename T>
+void PowStepRPUDevice<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 *gamma_down = w_gamma_down_[0];
+  T *gamma_up = w_gamma_up_[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];
+  T write_noise_std = par.getScaledWriteNoise();
+
+  PULSED_UPDATE_W_LOOP_DENSE(update_once(
+                                 w[j], w_apparent[j], sign, scale_down[j], scale_up[j],
+                                 gamma_down[j], gamma_up[j], min_bound[j], max_bound[j],
+                                 par.dw_min_std, write_noise_std, rng););
+}
+
+template class PowStepRPUDevice<float>;
+#ifdef RPU_USE_DOUBLE
+template class PowStepRPUDevice<double>;
+#endif
+
+} // namespace RPU
diff --git a/src/rpucuda/rpu_powstep_device.h b/src/rpucuda/rpu_powstep_device.h
new file mode 100644
index 00000000..fc5a211d
--- /dev/null
+++ b/src/rpucuda/rpu_powstep_device.h
@@ -0,0 +1,128 @@
+/**
+ * (C) Copyright 2020, 2021 IBM. All Rights Reserved.
+ *
+ * This code is licensed under the Apache License, Version 2.0. You may
+ * obtain a copy of this license in the LICENSE.txt file in the root directory
+ * of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+ *
+ * Any modifications or derivative works of this code must retain this
+ * copyright notice, and modified files need to carry a notice indicating
+ * that they have been altered from the originals.
+ */
+
+#pragma once
+
+#include "rng.h"
+#include "rpu_pulsed_device.h"
+#include "utility_functions.h"
+
+namespace RPU {
+
+template <typename T> class PowStepRPUDevice;
+
+BUILD_PULSED_DEVICE_META_PARAMETER(
+    PowStep,
+    /*implements*/
+    DeviceUpdateType::PowStep,
+    /*parameter def*/
+    T ps_gamma = (T)1.0;
+    T ps_gamma_dtod = (T)0.1;
+    T ps_gamma_up_down = (T)0.00;
+    T ps_gamma_up_down_dtod = (T)0.00;
+    ,
+    /*print body*/
+    ss << "\t ps_gamma:\t\t" << ps_gamma << "\t(dtod=" << ps_gamma_dtod << ")" << std::endl;
+    ss << "\t ps_gamma_up_down:\t" << ps_gamma_up_down << "\t(dtod=" << ps_gamma_up_down_dtod << ")"
+       << std::endl;
+    ,
+    /*Add*/
+    bool implementsWriteNoise() const override { return true; };);
+
+template <typename T> class PowStepRPUDevice : public PulsedRPUDevice<T> {
+
+  BUILD_PULSED_DEVICE_CONSTRUCTORS(
+      PowStepRPUDevice,
+      /* ctor*/
+      int x_sz = this->x_size_;
+      int d_sz = this->d_size_;
+
+      w_gamma_down_ = Array_2D_Get<T>(d_sz, x_sz);
+      w_gamma_up_ = Array_2D_Get<T>(d_sz, x_sz);
+
+      for (int j = 0; j < x_sz; ++j) {
+        for (int i = 0; i < d_sz; ++i) {
+          w_gamma_up_[i][j] = (T)0.0;
+          w_gamma_down_[i][j] = (T)0.0;
+        }
+      },
+      /* dtor*/
+      Array_2D_Free<T>(w_gamma_down_);
+      Array_2D_Free<T>(w_gamma_up_);
+      ,
+      /* copy */
+      for (int j = 0; j < other.x_size_; ++j) {
+        for (int i = 0; i < other.d_size_; ++i) {
+          w_gamma_down_[i][j] = other.w_gamma_down_[i][j];
+          w_gamma_up_[i][j] = other.w_gamma_up_[i][j];
+        }
+      },
+      /* move assignment */
+      w_gamma_down_ = other.w_gamma_down_;
+      w_gamma_up_ = other.w_gamma_up_;
+
+      other.w_gamma_down_ = nullptr;
+      other.w_gamma_up_ = nullptr;
+      ,
+      /* swap*/
+      swap(a.w_gamma_up_, b.w_gamma_up_);
+      swap(a.w_gamma_down_, b.w_gamma_down_);
+      ,
+      /* dp names*/
+      names.push_back(std::string("gamma_up"));
+      names.push_back(std::string("gamma_down"));
+      ,
+      /* dp2vec body*/
+      int n_prev = (int)names.size();
+      int size = this->x_size_ * this->d_size_;
+
+      for (int i = 0; i < size; ++i) {
+        data_ptrs[n_prev][i] = w_gamma_up_[0][i];
+        data_ptrs[n_prev + 1][i] = w_gamma_down_[0][i];
+      },
+      /* vec2dp body*/
+      int n_prev = (int)names.size();
+      int size = this->x_size_ * this->d_size_;
+
+      for (int i = 0; i < size; ++i) {
+        w_gamma_up_[0][i] = data_ptrs[n_prev][i];
+        w_gamma_down_[0][i] = data_ptrs[n_prev + 1][i];
+      },
+      /*invert copy DP */
+      T **gamma_down = rpu->getGammaDown();
+      T **gamma_up = rpu->getGammaUp();
+      for (int j = 0; j < this->x_size_; ++j) {
+        for (int i = 0; i < this->d_size_; ++i) {
+          w_gamma_down_[i][j] = gamma_up[i][j];
+          w_gamma_up_[i][j] = gamma_down[i][j];
+        }
+      }
+
+  );
+
+  void printDP(int x_count, int d_count) const override;
+
+  T getDwMin() const override;
+
+  inline T **getGammaUp() const { return w_gamma_up_; };
+  inline T **getGammaDown() const { return w_gamma_down_; };
+
+  void doSparseUpdate(
+      T **weights, int i, const int *x_signed_indices, int x_count, int d_sign, RNG<T> *rng)
+      override;
+  void doDenseUpdate(T **weights, int *coincidences, RNG<T> *rng) override;
+
+private:
+  T **w_gamma_up_ = nullptr;
+  T **w_gamma_down_ = nullptr;
+};
+} // namespace RPU
diff --git a/src/rpucuda/rpu_pulsed.cpp b/src/rpucuda/rpu_pulsed.cpp
index 0f0d1361..7e0abdad 100644
--- a/src/rpucuda/rpu_pulsed.cpp
+++ b/src/rpucuda/rpu_pulsed.cpp
@@ -306,8 +306,8 @@ template <typename T> void RPUPulsed<T>::setWeightsReal(const T *weightsptr, int
   if (dpar != nullptr) {
     w_min = dpar->w_min;
     w_max = dpar->w_max;
-    dynamic_cast<PulsedRPUDeviceBase<T> *>(&*rpu_device_)
-        ->getDwMin(); // this should be safe since we checked the params
+    dw_min = dynamic_cast<PulsedRPUDeviceBase<T> *>(&*rpu_device_)
+                 ->getDwMin(); // this should be safe since we checked the params
   }
   int BL = 0;
   T A = (T)0.0;
diff --git a/src/rpucuda/rpu_pulsed_device.h b/src/rpucuda/rpu_pulsed_device.h
index efefd65f..e2a881b1 100644
--- a/src/rpucuda/rpu_pulsed_device.h
+++ b/src/rpucuda/rpu_pulsed_device.h
@@ -323,7 +323,7 @@ public:
     PulsedRPUDevice<T>::copyInvertDeviceParameter(rpu_device);                                     \
     const auto *rpu = dynamic_cast<const CLASSNAME<T> *>(rpu_device);                              \
     if (rpu == nullptr) {                                                                          \
-      RPU_FATAL("Expect RPU Pulsed device");                                                       \
+      RPU_FATAL("Wrong device class");                                                             \
     };                                                                                             \
     INVERT_COPY_BODY;                                                                              \
   };                                                                                               \
diff --git a/src/rpucuda/rpu_simple_device.h b/src/rpucuda/rpu_simple_device.h
index e76f130a..baad5ea2 100644
--- a/src/rpucuda/rpu_simple_device.h
+++ b/src/rpucuda/rpu_simple_device.h
@@ -40,7 +40,8 @@ enum DeviceUpdateType {
   Vector,
   Difference,
   Transfer,
-  MixedPrec
+  MixedPrec,
+  PowStep,
 };
 
 // inherit from Simple


Test patch:
diff --git a/src/rpucuda/cuda/rpucuda_powstep_test.cpp b/src/rpucuda/cuda/rpucuda_powstep_test.cpp
new file mode 100644
index 00000000..e505d920
--- /dev/null
+++ b/src/rpucuda/cuda/rpucuda_powstep_test.cpp
@@ -0,0 +1,258 @@
+/**
+ * (C) Copyright 2020, 2021 IBM. All Rights Reserved.
+ *
+ * This code is licensed under the Apache License, Version 2.0. You may
+ * obtain a copy of this license in the LICENSE.txt file in the root directory
+ * of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+ *
+ * Any modifications or derivative works of this code must retain this
+ * copyright notice, and modified files need to carry a notice indicating
+ * that they have been altered from the originals.
+ */
+
+#include "cuda.h"
+#include "cuda_util.h"
+#include "io_manager.h"
+#include "rng.h"
+#include "rpucuda_powstep_device.h"
+#include "rpucuda_pulsed.h"
+#include "utility_functions.h"
+#include "gtest/gtest.h"
+#include <chrono>
+#include <memory>
+#include <random>
+
+#define TOLERANCE 1e-5
+
+#ifdef RPU_USE_DOUBLE
+typedef double num_t;
+#else
+typedef float num_t;
+#endif
+
+namespace {
+
+using namespace RPU;
+
+class RPUCudaPowStepTestFixture : public ::testing::TestWithParam<num_t> {
+public:
+  void SetUp() {
+
+    context_container = RPU::make_unique<CudaContext>(-1, false);
+    context = &*context_container;
+
+    x_size = 100;
+    d_size = 110;
+    repeats = 100;
+    K = 10;
+    m_batch = 100; // for the batched versions
+    dim3 = 4;
+
+    num_t bmin = -0.7;
+    num_t bmax = 0.7;
+
+    PulsedMetaParameter<num_t> p;
+    PowStepRPUDeviceMetaParameter<num_t> dp;
+    IOMetaParameter<num_t> p_io;
+
+    p_io.out_noise = 0.0; // no noise in output;
+    dp.dw_min_std = 0.0;
+
+    p.up.desired_BL = K;
+    p.up.pulse_type = PulseType::DeterministicImplicit;
+    p.up.d_res_implicit = 0.01;
+    p.up.x_res_implicit = 0.01;
+
+    dp.w_max = 1;
+    dp.w_min = -1;
+    dp.w_min_dtod = 0.0;
+    dp.w_max_dtod = 0.0;
+
+    dp.dw_min = 0.01;
+
+    // peripheral circuits specs
+    p_io.inp_res = -1;
+    p_io.inp_sto_round = false;
+    p_io.out_res = -1;
+
+    p_io.noise_management = NoiseManagementType::AbsMax;
+    p_io.bound_management = BoundManagementType::None;
+
+    p.f_io = p_io;
+    p.b_io = p_io;
+
+    p.up.update_management = true;
+    p.up.update_bl_management = true;
+
+    dp.lifetime = 0;
+    dp.lifetime_dtod = 0;
+
+    dp.diffusion = 0.00;
+    dp.diffusion_dtod = 0.0;
+
+    // slope
+    dp.ps_gamma = GetParam();
+    dp.ps_gamma_dtod = 0.01;
+    dp.ps_gamma_up_down = 0.1;
+    dp.ps_gamma_up_down_dtod = 0.01;
+
+    dp.print();
+
+    rx.resize(x_size * m_batch);
+    rd.resize(d_size * m_batch);
+
+    num_t lr = 0.4;
+
+    layer_pulsed = RPU::make_unique<RPUPulsed<num_t>>(x_size, d_size);
+    layer_pulsed->populateParameter(&p, &dp);
+    layer_pulsed->setLearningRate(lr);
+    layer_pulsed->setWeightsUniformRandom(bmin, bmax);
+
+    layer_pulsed->disp();
+
+    // culayer
+    culayer_pulsed = RPU::make_unique<RPUCudaPulsed<num_t>>(context, *layer_pulsed);
+    culayer_pulsed->disp();
+
+    unsigned int seed = std::chrono::system_clock::now().time_since_epoch().count();
+    std::default_random_engine generator{seed};
+    std::uniform_real_distribution<num_t> udist(-1., 1.);
+    auto urnd = std::bind(udist, generator);
+
+    // just assign some numbers from the weigt matrix
+    for (int i = 0; i < x_size * m_batch; i++)
+      rx[i] = urnd();
+
+    for (int j = 0; j < d_size * m_batch; j++) {
+      rd[j] = urnd();
+    }
+
+    x_cuvec = RPU::make_unique<CudaArray<num_t>>(context, x_size);
+    x_vec.resize(x_size);
+
+    d_cuvec = RPU::make_unique<CudaArray<num_t>>(context, d_size);
+    d_vec.resize(d_size);
+
+    x_cuvec_batch = RPU::make_unique<CudaArray<num_t>>(context, m_batch * x_size);
+    x_vec_batch.resize(m_batch * x_size);
+
+    d_cuvec_batch = RPU::make_unique<CudaArray<num_t>>(context, m_batch * d_size);
+    d_vec_batch.resize(m_batch * d_size);
+  }
+
+  std::unique_ptr<CudaContext> context_container;
+  CudaContext *context;
+
+  std::unique_ptr<RPUPulsed<num_t>> layer_pulsed;
+  std::unique_ptr<RPUCudaPulsed<num_t>> culayer_pulsed;
+  std::vector<num_t> x_vec, x_vec_batch, d_vec, d_vec_batch, rx, rd;
+  std::unique_ptr<CudaArray<num_t>> x_cuvec, x_cuvec_batch;
+  std::unique_ptr<CudaArray<num_t>> d_cuvec, d_cuvec_batch;
+  int x_size;
+  int d_size;
+  int repeats;
+  int K;
+  int dim3;
+  int m_batch;
+
+  num_t noise_value;
+};
+
+// types
+INSTANTIATE_TEST_CASE_P(GammaTest, RPUCudaPowStepTestFixture, ::testing::Values(1.2, 1.0));
+
+#define RPU_TEST_UPDATE(CUFUN, FUN, NLOOP)                                                         \
+  this->context->synchronizeDevice();                                                              \
+  this->culayer_pulsed->printWeights(3, 1);                                                        \
+  this->layer_pulsed->printWeights(3, 1);                                                          \
+                                                                                                   \
+  std::cout << "RPU Cuda:\n";                                                                      \
+  this->culayer_pulsed->printRPUParameter(1, 1);                                                   \
+  std::cout << "RPU:\n";                                                                           \
+  this->layer_pulsed->printRPUParameter(1, 1);                                                     \
+                                                                                                   \
+  num_t **refweights = Array_2D_Get<num_t>(this->d_size, this->x_size);                            \
+  num_t **w = this->layer_pulsed->getWeights();                                                    \
+  for (int i = 0; i < this->d_size; i++) {                                                         \
+    for (int j = 0; j < this->x_size; j++) {                                                       \
+      refweights[i][j] = w[i][j];                                                                  \
+    }                                                                                              \
+  }                                                                                                \
+                                                                                                   \
+  int nloop = NLOOP;                                                                               \
+  double cudur = 0, dur = 0;                                                                       \
+  int count_diff = 0;                                                                              \
+  auto start_time = std::chrono::high_resolution_clock::now();                                     \
+  for (int loop = 0; loop < nloop; loop++) {                                                       \
+                                                                                                   \
+    start_time = std::chrono::high_resolution_clock::now();                                        \
+    this->culayer_pulsed->CUFUN;                                                                   \
+    this->context->synchronizeDevice();                                                            \
+    auto end_time = std::chrono::high_resolution_clock::now();                                     \
+    cudur += std::chrono::duration_cast<std::chrono::microseconds>(end_time - start_time).count(); \
+                                                                                                   \
+    start_time = std::chrono::high_resolution_clock::now();                                        \
+    this->layer_pulsed->FUN;                                                                       \
+    end_time = std::chrono::high_resolution_clock::now();                                          \
+    dur += std::chrono::duration_cast<std::chrono::microseconds>(end_time - start_time).count();   \
+                                                                                                   \
+    num_t **cuweights = this->culayer_pulsed->getWeights();                                        \
+    num_t **weights = this->layer_pulsed->getWeights();                                            \
+    this->context->synchronizeDevice();                                                            \
+    count_diff = 0;                                                                                \
+    for (int i = 0; i < this->d_size; i++) {                                                       \
+      for (int j = 0; j < this->x_size; j++) {                                                     \
+        ASSERT_NEAR(weights[i][j], cuweights[i][j], TOLERANCE);                                    \
+        count_diff += weights[i][j] != refweights[i][j];                                           \
+      }                                                                                            \
+    }                                                                                              \
+    ASSERT_TRUE(count_diff > 0);                                                                   \
+    if (loop % 3 == 0) {                                                                           \
+      this->culayer_pulsed->setWeights(refweights[0]);                                             \
+      this->context->synchronizeDevice();                                                          \
+      this->layer_pulsed->setWeights(refweights[0]);                                               \
+    }                                                                                              \
+  }                                                                                                \
+  std::cout << "changed weights: " << count_diff << std::endl;                                     \
+  std::cout << BOLD_ON << "\nCUDA Updates done in: " << (num_t)cudur / 1000. / nloop << " msec. "  \
+            << BOLD_OFF << std::endl;                                                              \
+  std::cout << BOLD_ON << "RPU Updates done in: " << (num_t)dur / 1000. / nloop << " msec.\n "     \
+            << BOLD_OFF << std::endl;                                                              \
+                                                                                                   \
+  Array_2D_Free(refweights);
+
+TEST_P(RPUCudaPowStepTestFixture, UpdateVector) {
+
+  this->x_cuvec->assign(this->rx.data());
+  this->x_vec = this->rx;
+
+  this->d_cuvec->assign(this->rd.data());
+  this->d_vec = this->rd;
+
+  this->context->synchronizeDevice();
+
+  RPU_TEST_UPDATE(
+      update(this->x_cuvec->getData(), this->d_cuvec->getData(), false, 1),
+      update(this->x_vec.data(), this->d_vec.data(), false, 1), this->repeats);
+}
+
+TEST_P(RPUCudaPowStepTestFixture, UpdateMatrixBatch) {
+
+  this->x_cuvec_batch->assign(this->rx.data());
+  this->x_vec_batch = this->rx;
+
+  this->d_cuvec_batch->assign(this->rd.data());
+  this->d_vec_batch = this->rd;
+
+  RPU_TEST_UPDATE(
+      update(this->x_cuvec_batch->getData(), this->d_cuvec_batch->getData(), false, this->m_batch),
+      update(this->x_vec_batch.data(), this->d_vec_batch.data(), false, this->m_batch), 10);
+}
+
+} // namespace
+
+int main(int argc, char **argv) {
+  resetCuda();
+  ::testing::InitGoogleTest(&argc, argv);
+  return RUN_ALL_TESTS();
+}
diff --git a/tests/helpers/tiles.py b/tests/helpers/tiles.py
index 098afb40..c843929f 100644
--- a/tests/helpers/tiles.py
+++ b/tests/helpers/tiles.py
@@ -22,6 +22,7 @@
     ExpStepDevice,
     SoftBoundsDevice,
     SoftBoundsPmaxDevice,
+    PowStepDevice,
     IOParameters,
     DifferenceUnitCell,
     VectorUnitCell,
@@ -145,6 +146,21 @@ def get_tile(self, out_size, in_size, rpu_config=None, **kwargs):
         return AnalogTile(out_size, in_size, rpu_config, **kwargs)
 
 
+class PowStep:
+    """AnalogTile with PowStepDevice."""
+
+    simulator_tile_class = tiles.AnalogTile
+    first_hidden_field = 'max_bound'
+    use_cuda = False
+
+    def get_rpu_config(self):
+        return SingleRPUConfig(device=PowStepDevice(w_max_dtod=0, w_min_dtod=0))
+
+    def get_tile(self, out_size, in_size, rpu_config=None, **kwargs):
+        rpu_config = rpu_config or self.get_rpu_config()
+        return AnalogTile(out_size, in_size, rpu_config, **kwargs)
+
+
 class Vector:
     """AnalogTile with VectorUnitCell."""
 
@@ -364,6 +380,21 @@ def get_tile(self, out_size, in_size, rpu_config=None, **kwargs):
         return AnalogTile(out_size, in_size, rpu_config, **kwargs).cuda()
 
 
+class PowStepCuda:
+    """AnalogTile with PowStepDevice."""
+
+    simulator_tile_class = getattr(tiles, 'CudaAnalogTile', None)
+    first_hidden_field = 'max_bound'
+    use_cuda = True
+
+    def get_rpu_config(self):
+        return SingleRPUConfig(device=PowStepDevice(w_max_dtod=0, w_min_dtod=0))
+
+    def get_tile(self, out_size, in_size, rpu_config=None, **kwargs):
+        rpu_config = rpu_config or self.get_rpu_config()
+        return AnalogTile(out_size, in_size, rpu_config, **kwargs).cuda()
+
+
 class VectorCuda:
     """AnalogTile with VectorUnitCell."""
 
diff --git a/tests/test_rpu_configurations.py b/tests/test_rpu_configurations.py
index dd6ec7a4..8b49e523 100644
--- a/tests/test_rpu_configurations.py
+++ b/tests/test_rpu_configurations.py
@@ -20,10 +20,11 @@
 from .helpers.testcases import ParametrizedTestCase
 from .helpers.tiles import (
     FloatingPoint, Ideal, ConstantStep, LinearStep,
-    ExpStep, SoftBounds, SoftBoundsPmax, Vector, Difference, Transfer, MixedPrecision,
+    ExpStep, SoftBounds, SoftBoundsPmax, PowStep, Vector, Difference,
+    Transfer, MixedPrecision,
     FloatingPointCuda, IdealCuda, ConstantStepCuda, LinearStepCuda,
-    ExpStepCuda, SoftBoundsCuda, SoftBoundsPmaxCuda, VectorCuda,
-    DifferenceCuda, TransferCuda, MixedPrecisionCuda,
+    ExpStepCuda, SoftBoundsCuda, SoftBoundsPmaxCuda, PowStepCuda,
+    VectorCuda, DifferenceCuda, TransferCuda, MixedPrecisionCuda
 )
 
 
@@ -64,6 +65,7 @@ def test_config_device_parameters(self):
     ExpStep,
     SoftBounds,
     SoftBoundsPmax,
+    PowStep,
     Vector,
     Difference,
     Transfer,
@@ -74,6 +76,7 @@ def test_config_device_parameters(self):
     ExpStepCuda,
     SoftBoundsCuda,
     SoftBoundsPmaxCuda,
+    PowStepCuda,
     VectorCuda,
     DifferenceCuda,
     TransferCuda,