adding new 1/f^alpha generators and switch to smart ptr

core/junction.hpp

@@ -149,7 +149,7 @@ class Layer
 {
 private:
 
-    OneFNoise<T>* ofn;
+    std::shared_ptr<OneFNoise<T>> ofn;
 
     ScalarDriver<T> temperatureDriver;
     ScalarDriver<T> IECDriverTop;
@@ -216,7 +216,7 @@ class Layer
         dWn = CVector<T>(this->distribution);
         dWn.normalize();
         this->cellVolume = this->cellSurface * this->thickness;
-        this->ofn = new OneFNoise<T>(0, 0., 0.);
+        this->ofn = std::shared_ptr<OneFNoise<T>>(new OneFNoise<T>(0, 0., 0.));
     }
 
 public:
@@ -440,7 +440,7 @@ class Layer
     }
 
     void setOneFNoise(unsigned int sources, T bias, T scale) {
-        this->ofn = new OneFNoise<T>(sources, bias, scale);
+        this->ofn = std::shared_ptr<OneFNoise<T>>(new OneFNoise<T>(sources, bias, scale));
         this->pinkNoiseSet = true;
         this->nonStochasticOneFSet = true; // by default turn it on, but in the stochastic sims, we will have to turn it off
     }
@@ -678,6 +678,8 @@ class Layer
         {
             this->I_log = this->currentDriver.getCurrentScalarValue(time);
             // use standard STT formulation
+            // see that literature reports Ms/MAGNETIC_PERMEABILITY
+            // but then the units don't match, we use Ms [T] which works
             const T aJ = HBAR * this->I_log /
                 (ELECTRON_CHARGE * this->Ms * this->thickness);
             // field like

core/noise.hpp

@@ -21,6 +21,10 @@
 #include <numeric>                 // for accumulate
 #include <random>                  // for uniform_real_distribution, geometr...
 #include <vector>                  // for vector
+#include <complex>
+#include <memory>
+#include "cvector.hpp"
+#include "../third_party/kissfft/kissfft.hh"
 
 #define PINK_NOISE_NUM_STAGES 3
 template <typename T>
@@ -118,5 +122,103 @@ class OneFNoise {
     }
 };
 
+std::mt19937 generator(std::random_device{}());
+template<typename T = double>
+class BufferedAlphaNoise {
+private:
+    std::vector<std::complex<float>> bufferWhite, bufferColoured;
+    std::vector<std::complex<float>> bufferWhiteComplex, bufferColouredComplex;
+    std::vector<float> result;
+    unsigned int bufferSize;
+    std::function <float()> gaussPDF;
+    T alpha = 1.;
+    T scale = 1.;
+    std::shared_ptr<kissfft<float>> fwd, inv; // configs for forward and inverse real fft
+    unsigned int internalCounter = 0;
+    unsigned int refills = 0;
+public:
+
+    /**
+     * @brief Construct a new Buffered Alpha Noise object
+     *
+     * @param bufferSize the size of the buffer
+     * @param alpha the alpha parameter 1/f^alpha
+     * @param std the standard deviation of the gaussian
+     * @param scale the scaling parameter
+     */
+    BufferedAlphaNoise(unsigned int bufferSize, T alpha, T std, T scale) : bufferSize(bufferSize), alpha(alpha), scale(scale) {
+
+        this->bufferColoured.resize(2 * bufferSize);
+        this->bufferWhite.resize(2 * bufferSize);
+        this->result.resize(bufferSize);
+        this->bufferColouredComplex.resize(2 * bufferSize);
+        this->bufferWhiteComplex.resize(2 * bufferSize);
+
+        // these are the filter weights -- we only have to fill it once per alpha and N
+        this->bufferColoured[0] = 1.0;
+        for (unsigned int i = 1; i < this->bufferSize; ++i) {
+            const float weight = (float)(0.5 * alpha + (
+                (float)(i - 1)
+                )) / ((float)i);
+            this->bufferColoured[i] = this->bufferColoured[i - 1] * weight;
+        }
+
+        this->gaussPDF = std::bind(std::normal_distribution<float>(0, std), std::ref(generator));
+
+        this->fwd = std::shared_ptr<kissfft<float>>(new kissfft<float>(2 * this->bufferSize, false));
+        this->inv = std::shared_ptr<kissfft<float>>(new kissfft<float>(2 * this->bufferSize, true));
+    }
+
+    ~BufferedAlphaNoise() {
+    }
+
+    void fillBuffer() {
+        // this is actual generation function
+        // generate random white as a baseline, only N values, rest is 0 padded
+        std::generate(this->bufferWhite.begin(), this->bufferWhite.begin() + this->bufferSize, this->gaussPDF);
+
+        for (unsigned int i = this->bufferSize; i < 2 * this->bufferSize; ++i) {
+            this->bufferColoured[i] = 0;
+            this->bufferWhite[i] = 0;
+        }
+        // perform the fft
+        this->fwd->transform(&this->bufferWhite[0], &this->bufferWhiteComplex[0]);
+        this->fwd->transform(&this->bufferColoured[0], &this->bufferColouredComplex[0]);
+
+        // multiply the two
+        for (unsigned int i = 0; i < this->bufferSize; ++i) {
+            this->bufferColouredComplex[i] = this->bufferColouredComplex[i] * this->bufferWhiteComplex[i];
+        }
+        // invert
+        this->bufferColouredComplex[0] = this->bufferColouredComplex[0] / std::complex<float>(2.0, 0);
+        this->bufferColouredComplex[this->bufferSize - 1] = this->bufferColouredComplex[this->bufferSize - 1] / std::complex<float>(2.0, 0);
+        for (unsigned int i = this->bufferSize; i < 2 * this->bufferSize; ++i) {
+            this->bufferColouredComplex[i] = 0.;
+        }
+        this->inv->transform(&this->bufferColouredComplex[0], &this->bufferWhiteComplex[0]);
+
+        std::transform(
+            this->bufferWhiteComplex.begin(),
+            this->bufferWhiteComplex.begin() + this->bufferSize,
+            this->result.begin(),
+            [&](std::complex<float> x) { return x.real() / (this->bufferSize); }
+        );
+    }
+
+    const std::vector<float>& getFullBuffer() {
+        return this->result;
+    }
+
+    T tick() {
+        // we measure only up to a buffer size, not 2x buffer size
+        if (this->internalCounter == 0) {
+            this->fillBuffer();
+        }
+        const auto ret = this->result[this->internalCounter];
+        this->internalCounter = (this->internalCounter + 1) % this->bufferSize;
+        return this->scale * ret;
+    }
+
+};
 
 #endif

python/cmtj.cpp

@@ -8,6 +8,7 @@
 #include "../core/cvector.hpp"
 #include "../core/drivers.hpp"
 #include "../core/junction.hpp"
+#include "../core/noise.hpp"
 #include <stdio.h>
 #include <vector>
 
@@ -321,4 +322,15 @@ PYBIND11_MODULE(cmtj, m)
         .def("setLayerAnisotropy", &Reservoir::setLayerAnisotropy)
         .def("setLayerExternalField", &Reservoir::setLayerExternalField)
         .def("getMagnetisation", &Reservoir::getMagnetisation);
+
+    // generatormodule
+    py::module generator_module = m.def_submodule("noise", "Submodule with noise generation functions");
+    py::class_<BufferedAlphaNoise<double>>(generator_module, "BufferedAlphaNoise")
+        .def(py::init<unsigned int, double, double, double>(),
+            "bufferSize"_a,
+            "alpha"_a,
+            "std"_a,
+            "scale"_a)
+        .def("fillBuffer", &BufferedAlphaNoise<double>::fillBuffer)
+        .def("tick", &BufferedAlphaNoise<double>::tick);
 }