Upload files to "ex1G_dense_matrices_access"

ex1G_dense_matrices_access/DenseMatrix.h

@@ -0,0 +1,71 @@
+#pragma once
+#include "sigmoid.h"
+#include <iostream>
+#include <vector>
+using namespace std;
+
+class DenseMatrix
+{
+    private:
+        vector<double> M;
+        size_t n;
+        size_t m;
+
+
+    public:
+        vector<double> Mult(const vector<double> &x) const
+        {
+            vector<double> y(n,0);
+            for(size_t i = 0; i < n; ++i)   // iterate row
+            {
+                for(size_t j = 0; j < m; ++j)   // iterate column
+                {
+                    y[i] += M[i*m + j]*x[j];
+                }
+            }
+            return y;
+        }
+
+        vector<double> MultT(const vector<double> &y) const
+        {
+            vector<double> x(m,0);
+            for(size_t j = 0; j < m; ++j)   // iterate column
+            {
+                for(size_t i = 0; i < n; ++i)   // iterate row
+                {
+                    x[j] += M[i*m + j]*y[i];
+                }
+            }
+            return x;
+        }
+
+        void Print() const
+        {
+            for(size_t i = 0; i < n; ++i)   // iterate row
+            {
+                for(size_t j = 0; j < m; ++j)   // iterate column
+                {
+                    cout << M[i*m + j] << "  ";
+                }
+                cout << endl;
+            }
+            cout << endl;
+        }
+
+
+        DenseMatrix(size_t n, size_t m)
+        {
+            this->n = n;
+            this->m = m;
+            M = vector<double>(n*m);
+            size_t nm = max(n,m);
+
+            for(size_t i = 0; i < n; ++i)   // iterate row
+            {
+                for(size_t j = 0; j < m; ++j)   // iterate column
+                {
+                    M[i*m + j] = sigmoid(x_entry(i,nm))*sigmoid(x_entry(j,nm));
+                }
+            }
+        }
+};

ex1G_dense_matrices_access/ProductMatrix.h

@@ -0,0 +1,52 @@
+#pragma once
+#include <cmath>
+#include <iostream>
+#include <vector>
+using namespace std;
+
+class ProductMatrix
+{
+    private:
+        vector<double> u;
+        vector<double> v;
+        size_t n;
+        size_t m;
+
+    public:
+        vector<double> Mult(const vector<double> &x) const
+        {
+            vector<double> y(n,0);
+            for(int i = 0; i < n; ++i)
+            {
+                for(int j = 0; j < m; ++j)
+                {
+                    y[i] += v[j]*x[j];
+                }
+                y[i] *= u[i];
+            }
+            return y;
+        }
+
+        vector<double> MultT(const vector<double> &y) const
+        {
+            vector<double> x(m,0);
+            for(int j = 0; j < m; ++j)
+            {
+                for(int i = 0; i < n; ++i)
+                {
+                    x[j] += y[i]*u[i];
+                }
+                x[j] *= v[j];
+            }
+            return x;
+        }
+
+        ProductMatrix(const vector<double> &u, const vector<double> &v)
+        {
+            n = u.size();
+            m = v.size();
+            this->u = u;
+            this->v = v;
+            
+        }
+};

ex1G_dense_matrices_access/main.cpp

@@ -0,0 +1,93 @@
+#include "../utils/timing.h"
+#include "DenseMatrix.h"
+#include "ProductMatrix.h"
+#include <algorithm>
+
+int main()
+{
+    // b) ------------------------------------------------------------------------------------------------------
+    DenseMatrix const M(5,3);
+    vector<double> const u{{1, 2, 3}};
+    vector<double> f1 = M.Mult(u);
+    vector<double> const v{{-1, 2, -3, 4, -5}};
+    vector<double> f2 = M.MultT(v);
+
+    M.Print();
+
+    for(size_t i = 0; i < f1.size(); ++i)
+        cout << f1[i] << endl;
+    cout << endl;
+
+
+    for(size_t j = 0; j < f2.size(); ++j)
+        cout << f2[j] << "  ";
+    cout << endl << "-------------------------------------------------" << endl;
+
+    // c) ------------------------------------------------------------------------------------------------------
+    size_t n = pow(10,3);
+    DenseMatrix const M_1(n,n);
+    vector<double> x(n, 1.0);
+
+    size_t n_loops = 100;
+    vector<double> y_1;
+    vector<double> y_2;
+
+    double time_1 = 0;
+    double time_2 = 0;
+
+    tic();
+    for(int l = 0; l < n_loops; ++l)
+        y_1 = M_1.Mult(x);
+    time_1 += toc();
+
+    tic();
+    for(int l = 0; l < n_loops; ++l)
+        y_2 = M_1.MultT(x);
+    time_2 += toc();
+
+    vector<double> error_vec(n,0);
+    for(int i = 0; i < n; ++i)
+        error_vec[i] = abs(y_1[i] - y_2[i]);
+    double sup_error = *max_element(error_vec.begin(), error_vec.end());
+
+
+    cout << "n = " << n << endl;
+    cout << "Average duration for Mult: " << time_1/n_loops << endl;
+    cout << "Average duration for MultT: " << time_2/n_loops << endl;
+    cout << "sup-error: " << sup_error << endl;
+    cout << "-------------------------------------------------" << endl;
+
+    // d) ------------------------------------------------------------------------------------------------------
+    vector<double> u_M(n,0);
+    for(int i = 0; i < n; ++i)
+        u_M[i] = sigmoid(x_entry(i, n));
+
+    ProductMatrix const M_2(u_M, u_M);
+
+    time_1 = 0;
+    time_2 = 0;
+
+    tic();
+    for(int l = 0; l < n_loops; ++l)
+        y_1 = M_2.Mult(x);
+    time_1 += toc();
+
+    tic();
+    for(int l = 0; l < n_loops; ++l)
+        y_2 = M_2.MultT(x);
+    time_2 += toc();
+
+    for(int i = 0; i < n; ++i)
+        error_vec[i] = abs(y_1[i] - y_2[i]);
+    sup_error = *max_element(error_vec.begin(), error_vec.end());
+
+
+    cout << "n = " << n << endl;
+    cout << "Average duration for Mult: " << time_1/n_loops << endl;
+    cout << "Average duration for MultT: " << time_2/n_loops << endl;
+    cout << "sup-error: " << sup_error << endl;
+    cout << "-------------------------------------------------" << endl;
+
+
+    return 0;
+}

ex1G_dense_matrices_access/sigmoid.h

@@ -0,0 +1,12 @@
+#pragma once
+#include <cmath>
+
+double sigmoid(double x)
+{
+    return 1./(1. + exp(-x));
+}
+
+double x_entry(size_t k, size_t nm)
+{
+    return (10.*k)/(nm - 1) - 5.;
+}