diff --git a/mlpp/activation/activation.cpp b/mlpp/activation/activation.cpp
index 8d4e490..10fa0f5 100644
--- a/mlpp/activation/activation.cpp
+++ b/mlpp/activation/activation.cpp
@@ -35,33 +35,21 @@ Ref<MLPPMatrix> MLPPActivation::run_activation_deriv_matrix(const ActivationFunc
 	return Ref<MLPPMatrix>();
 }
 
-Ref<MLPPVector> MLPPActivation::activation(const Ref<MLPPVector> &z, real_t (*function)(real_t), const bool deriv) {
-	return Ref<MLPPVector>();
-}
-Ref<MLPPVector> MLPPActivation::activation_norm(const Ref<MLPPVector> &z, real_t (*function)(real_t)) {
-	return Ref<MLPPVector>();
-}
-Ref<MLPPVector> MLPPActivation::activation_deriv(const Ref<MLPPVector> &z, real_t (*function)(real_t)) {
-	return Ref<MLPPVector>();
-}
+Ref<MLPPVector> MLPPActivation::activation(const Ref<MLPPVector> &z, real_t (*function)(real_t)) {
+	Ref<MLPPVector> a;
+	a.instance();
 
-/*
+	int size = z->size();
 
-// TO DO: Implement this template activation
-std::vector<real_t> MLPPActivation::activation(std::vector<real_t> z, bool deriv, real_t (*function)(real_t, bool)) {
-	if (deriv) {
-		std::vector<real_t> deriv;
-		deriv.resize(z.size());
-		for (int i = 0; i < z.size(); i++) {
-			deriv[i] = function(z[i], true);
-		}
-		return deriv;
-	}
-	std::vector<real_t> a;
-	a.resize(z.size());
-	for (int i = 0; i < z.size(); i++) {
-		a[i] = function(z[i], deriv);
+	a->resize(size);
+
+	const real_t *z_ptr = z->ptr();
+	real_t *a_ptr = a->ptrw();
+
+	for (int i = 0; i < size; ++i) {
+		a_ptr[i] = function(z_ptr[i]);
 	}
+
 	return a;
 }
 
@@ -83,11 +71,11 @@ real_t MLPPActivation::linear_deriv(real_t z) {
 }
 Ref<MLPPVector> MLPPActivation::linear_deriv(const Ref<MLPPVector> &z) {
 	MLPPLinAlg alg;
-	return alg.onevec(z.size());
+	return alg.onevecv(z->size());
 }
 Ref<MLPPMatrix> MLPPActivation::linear_deriv(const Ref<MLPPMatrix> &z) {
 	MLPPLinAlg alg;
-	return alg.onemat(z.size(), z[0].size());
+	return alg.onematm(z->size().x, z->size().y);
 }
 
 //SIGMOID
@@ -96,15 +84,14 @@ real_t MLPPActivation::sigmoid_norm(real_t z) {
 }
 Ref<MLPPVector> MLPPActivation::sigmoid_norm(const Ref<MLPPVector> &z) {
 	MLPPLinAlg alg;
-
-	return alg.elementWiseDivision(alg.onevec(z.size()), alg.addition(alg.onevec(z.size()), alg.exp(alg.scalarMultiply(-1, z))));
+	return alg.element_wise_division(alg.onevecv(z->size()), alg.additionm(alg.onevecv(z->size()), alg.expv(alg.scalar_multiplynv(-1, z))));
 }
 Ref<MLPPMatrix> MLPPActivation::sigmoid_norm(const Ref<MLPPMatrix> &z) {
 	MLPPLinAlg alg;
-
-	return alg.elementWiseDivision(alg.onemat(z.size(), z[0].size()), alg.addition(alg.onemat(z.size(), z[0].size()), alg.exp(alg.scalarMultiply(-1, z))));
+	return alg.element_wise_division(alg.onematm(z->size().x, z->size().y), alg.additionm(alg.onematm(z->size().x, z->size().y), alg.expv(alg.scalar_multiplynv(-1, z))));
 }
 
+/*
 real_t MLPPActivation::sigmoid_deriv(real_t z) {
 	return sigmoid_norm(z) * (1 - sigmoid_norm(z));
 }
@@ -1257,9 +1244,7 @@ Ref<MLPPMatrix> MLPPActivation::arcoth_deriv(const Ref<MLPPMatrix> &z) {
 
 	return alg.elementWiseDivision(alg.onemat(z.size(), z[0].size()), alg.subtraction(alg.onemat(z.size(), z[0].size()), alg.hadamard_product(z, z)));
 }
-
 */
-
 //======================== OLD =============================
 
 real_t MLPPActivation::linear(real_t z, bool deriv) {
diff --git a/mlpp/activation/activation.h b/mlpp/activation/activation.h
index d648194..405e0fd 100644
--- a/mlpp/activation/activation.h
+++ b/mlpp/activation/activation.h
@@ -66,9 +66,7 @@ public:
 	Ref<MLPPVector> run_activation_deriv_vector(const ActivationFunction func, const Ref<MLPPVector> &z);
 	Ref<MLPPMatrix> run_activation_deriv_matrix(const ActivationFunction func, const Ref<MLPPMatrix> &z);
 
-	Ref<MLPPVector> activation(const Ref<MLPPVector> &z, real_t (*function)(real_t), const bool deriv = false);
-	Ref<MLPPVector> activation_norm(const Ref<MLPPVector> &z, real_t (*function)(real_t));
-	Ref<MLPPVector> activation_deriv(const Ref<MLPPVector> &z, real_t (*function)(real_t));
+	Ref<MLPPVector> activation(const Ref<MLPPVector> &z, real_t (*function)(real_t));
 
 	//ACTIVATION FUNCTIONS
 
diff --git a/mlpp/lin_alg/lin_alg.cpp b/mlpp/lin_alg/lin_alg.cpp
index d4132fa..2feb86a 100644
--- a/mlpp/lin_alg/lin_alg.cpp
+++ b/mlpp/lin_alg/lin_alg.cpp
@@ -88,6 +88,77 @@ std::vector<std::vector<real_t>> MLPPLinAlg::matmult(std::vector<std::vector<rea
 	return C;
 }
 
+Ref<MLPPMatrix> MLPPLinAlg::additionm(const Ref<MLPPMatrix> &A, const Ref<MLPPMatrix> &B) {
+	ERR_FAIL_COND_V(!A.is_valid() || !B.is_valid(), Ref<MLPPMatrix>());
+	Size2i a_size = A->size();
+	ERR_FAIL_COND_V(a_size != B->size(), Ref<MLPPMatrix>());
+
+	Ref<MLPPMatrix> C;
+	C.instance();
+	C->resize(a_size);
+
+	const real_t *a_ptr = A->ptr();
+	const real_t *b_ptr = B->ptr();
+	real_t *c_ptr = C->ptrw();
+
+	int data_size = A->data_size();
+
+	for (int i = 0; i < data_size; ++i) {
+		c_ptr[i] = a_ptr[i] + b_ptr[i];
+	}
+
+	return C;
+}
+Ref<MLPPMatrix> MLPPLinAlg::subtractionm(const Ref<MLPPMatrix> &A, const Ref<MLPPMatrix> &B) {
+	ERR_FAIL_COND_V(!A.is_valid() || !B.is_valid(), Ref<MLPPMatrix>());
+	Size2i a_size = A->size();
+	ERR_FAIL_COND_V(a_size != B->size(), Ref<MLPPMatrix>());
+
+	Ref<MLPPMatrix> C;
+	C.instance();
+	C->resize(a_size);
+
+	const real_t *a_ptr = A->ptr();
+	const real_t *b_ptr = B->ptr();
+	real_t *c_ptr = C->ptrw();
+
+	int data_size = A->data_size();
+
+	for (int i = 0; i < data_size; ++i) {
+		c_ptr[i] = a_ptr[i] - b_ptr[i];
+	}
+
+	return C;
+}
+Ref<MLPPMatrix> MLPPLinAlg::matmultm(const Ref<MLPPMatrix> &A, const Ref<MLPPMatrix> &B) {
+	ERR_FAIL_COND_V(!A.is_valid() || !B.is_valid(), Ref<MLPPMatrix>());
+	Size2i a_size = A->size();
+	ERR_FAIL_COND_V(a_size != B->size(), Ref<MLPPMatrix>());
+
+	Ref<MLPPMatrix> C;
+	C.instance();
+	C->resize(a_size);
+
+	const real_t *a_ptr = A->ptr();
+	const real_t *b_ptr = B->ptr();
+	real_t *c_ptr = C->ptrw();
+
+	for (int i = 0; i < a_size.y; i++) {
+		for (int k = 0; k < a_size.y; k++) {
+			int ind_i_k = A->calculate_index(i, k);
+
+			for (int j = 0; j < a_size.x; j++) {
+				int ind_i_j = A->calculate_index(i, j);
+				int ind_k_j = A->calculate_index(k, j);
+
+				c_ptr[ind_i_j] += a_ptr[ind_i_k] * b_ptr[ind_k_j];
+			}
+		}
+	}
+
+	return C;
+}
+
 std::vector<std::vector<real_t>> MLPPLinAlg::hadamard_product(std::vector<std::vector<real_t>> A, std::vector<std::vector<real_t>> B) {
 	std::vector<std::vector<real_t>> C;
 	C.resize(A.size());
@@ -402,6 +473,34 @@ std::vector<std::vector<real_t>> MLPPLinAlg::onemat(int n, int m) {
 	return full(n, m, 1);
 }
 
+Ref<MLPPMatrix> MLPPLinAlg::zeromatm(int n, int m) {
+	Ref<MLPPMatrix> mat;
+	mat.instance();
+
+	mat->resize(Size2i(n, m));
+	mat->fill(0);
+
+	return mat;
+}
+Ref<MLPPMatrix> MLPPLinAlg::onematm(int n, int m) {
+	Ref<MLPPMatrix> mat;
+	mat.instance();
+
+	mat->resize(Size2i(n, m));
+	mat->fill(1);
+
+	return mat;
+}
+Ref<MLPPMatrix> MLPPLinAlg::fullm(int n, int m, int k) {
+	Ref<MLPPMatrix> mat;
+	mat.instance();
+
+	mat->resize(Size2i(n, m));
+	mat->fill(k);
+
+	return mat;
+}
+
 std::vector<std::vector<real_t>> MLPPLinAlg::full(int n, int m, int k) {
 	std::vector<std::vector<real_t>> full;
 	full.resize(n);
@@ -995,6 +1094,29 @@ std::vector<real_t> MLPPLinAlg::elementWiseDivision(std::vector<real_t> a, std::
 	return c;
 }
 
+Ref<MLPPVector> MLPPLinAlg::element_wise_division(const Ref<MLPPVector> &a, const Ref<MLPPVector> &b) {
+	ERR_FAIL_COND_V(!a.is_valid() || !b.is_valid(), Ref<MLPPVector>());
+
+	Ref<MLPPVector> out;
+	out.instance();
+
+	int size = a->size();
+
+	ERR_FAIL_COND_V(size != b->size(), Ref<MLPPVector>());
+
+	out->resize(size);
+
+	const real_t *a_ptr = a->ptr();
+	const real_t *b_ptr = b->ptr();
+	real_t *out_ptr = out->ptrw();
+
+	for (int i = 0; i < size; ++i) {
+		out_ptr[i] = a_ptr[i] / b_ptr[i];
+	}
+
+	return out;
+}
+
 std::vector<real_t> MLPPLinAlg::scalarMultiply(real_t scalar, std::vector<real_t> a) {
 	for (int i = 0; i < a.size(); i++) {
 		a[i] *= scalar;
@@ -1210,6 +1332,118 @@ std::vector<real_t> MLPPLinAlg::cbrt(std::vector<real_t> a) {
 	return exponentiate(a, real_t(1) / real_t(3));
 }
 
+Ref<MLPPVector> logv(const Ref<MLPPVector> &a) {
+	ERR_FAIL_COND_V(!a.is_valid(), Ref<MLPPVector>());
+
+	Ref<MLPPVector> out;
+	out.instance();
+
+	int size = a->size();
+	out->resize(size);
+
+	const real_t *a_ptr = a->ptr();
+	real_t *out_ptr = out->ptrw();
+
+	for (int i = 0; i < size; ++i) {
+		out_ptr[i] = Math::log(a_ptr[i]);
+	}
+
+	return out;
+}
+Ref<MLPPVector> log10v(const Ref<MLPPVector> &a) {
+	ERR_FAIL_COND_V(!a.is_valid(), Ref<MLPPVector>());
+
+	Ref<MLPPVector> out;
+	out.instance();
+
+	int size = a->size();
+	out->resize(size);
+
+	const real_t *a_ptr = a->ptr();
+	real_t *out_ptr = out->ptrw();
+
+	for (int i = 0; i < size; ++i) {
+		out_ptr[i] = std::log10(a_ptr[i]);
+	}
+
+	return out;
+}
+Ref<MLPPVector> expv(const Ref<MLPPVector> &a) {
+	ERR_FAIL_COND_V(!a.is_valid(), Ref<MLPPVector>());
+
+	Ref<MLPPVector> out;
+	out.instance();
+
+	int size = a->size();
+	out->resize(size);
+
+	const real_t *a_ptr = a->ptr();
+	real_t *out_ptr = out->ptrw();
+
+	for (int i = 0; i < size; ++i) {
+		out_ptr[i] = Math::exp(a_ptr[i]);
+	}
+
+	return out;
+}
+Ref<MLPPVector> erfv(const Ref<MLPPVector> &a) {
+	ERR_FAIL_COND_V(!a.is_valid(), Ref<MLPPVector>());
+
+	Ref<MLPPVector> out;
+	out.instance();
+
+	int size = a->size();
+	out->resize(size);
+
+	const real_t *a_ptr = a->ptr();
+	real_t *out_ptr = out->ptrw();
+
+	for (int i = 0; i < size; ++i) {
+		out_ptr[i] = std::erf(a_ptr[i]);
+	}
+
+	return out;
+}
+Ref<MLPPVector> exponentiatev(const Ref<MLPPVector> &a, real_t p) {
+	ERR_FAIL_COND_V(!a.is_valid(), Ref<MLPPVector>());
+
+	Ref<MLPPVector> out;
+	out.instance();
+
+	int size = a->size();
+	out->resize(size);
+
+	const real_t *a_ptr = a->ptr();
+	real_t *out_ptr = out->ptrw();
+
+	for (int i = 0; i < size; ++i) {
+		out_ptr[i] = Math::pow(a_ptr[i], p);
+	}
+
+	return out;
+}
+Ref<MLPPVector> sqrtv(const Ref<MLPPVector> &a) {
+	ERR_FAIL_COND_V(!a.is_valid(), Ref<MLPPVector>());
+
+	Ref<MLPPVector> out;
+	out.instance();
+
+	int size = a->size();
+	out->resize(size);
+
+	const real_t *a_ptr = a->ptr();
+	real_t *out_ptr = out->ptrw();
+
+	for (int i = 0; i < size; ++i) {
+		out_ptr[i] = Math::sqrt(a_ptr[i]);
+	}
+
+	return out;
+}
+Ref<MLPPVector> cbrtv(const Ref<MLPPVector> &a) {
+	return exponentiatev(a, static_cast<real_t>(1) / static_cast<real_t>(3));
+}
+
 real_t MLPPLinAlg::dot(std::vector<real_t> a, std::vector<real_t> b) {
 	real_t c = 0;
 	for (int i = 0; i < a.size(); i++) {
@@ -1265,6 +1499,34 @@ std::vector<real_t> MLPPLinAlg::full(int n, int k) {
 	return full;
 }
 
+Ref<MLPPVector> MLPPLinAlg::zerovecv(int n) {
+	Ref<MLPPVector> vec;
+	vec.instance();
+
+	vec->resize(n);
+	vec->fill(0);
+
+	return vec;
+}
+Ref<MLPPVector> MLPPLinAlg::onevecv(int n) {
+	Ref<MLPPVector> vec;
+	vec.instance();
+
+	vec->resize(n);
+	vec->fill(1);
+
+	return vec;
+}
+Ref<MLPPVector> MLPPLinAlg::fullv(int n, int k) {
+	Ref<MLPPVector> vec;
+	vec.instance();
+
+	vec->resize(n);
+	vec->fill(k);
+
+	return vec;
+}
+
 std::vector<real_t> MLPPLinAlg::sin(std::vector<real_t> a) {
 	std::vector<real_t> b;
 	b.resize(a.size());
diff --git a/mlpp/lin_alg/lin_alg.h b/mlpp/lin_alg/lin_alg.h
index f73e105..a8696d9 100644
--- a/mlpp/lin_alg/lin_alg.h
+++ b/mlpp/lin_alg/lin_alg.h
@@ -27,35 +27,27 @@ public:
 	std::vector<std::vector<real_t>> gaussianNoise(int n, int m);
 
 	std::vector<std::vector<real_t>> addition(std::vector<std::vector<real_t>> A, std::vector<std::vector<real_t>> B);
-
 	std::vector<std::vector<real_t>> subtraction(std::vector<std::vector<real_t>> A, std::vector<std::vector<real_t>> B);
-
 	std::vector<std::vector<real_t>> matmult(std::vector<std::vector<real_t>> A, std::vector<std::vector<real_t>> B);
 
+	Ref<MLPPMatrix> additionm(const Ref<MLPPMatrix> &A, const Ref<MLPPMatrix> &B);
+	Ref<MLPPMatrix> subtractionm(const Ref<MLPPMatrix> &A, const Ref<MLPPMatrix> &B);
+	Ref<MLPPMatrix> matmultm(const Ref<MLPPMatrix> &A, const Ref<MLPPMatrix> &B);
+
 	std::vector<std::vector<real_t>> hadamard_product(std::vector<std::vector<real_t>> A, std::vector<std::vector<real_t>> B);
-
 	std::vector<std::vector<real_t>> kronecker_product(std::vector<std::vector<real_t>> A, std::vector<std::vector<real_t>> B);
-
 	std::vector<std::vector<real_t>> elementWiseDivision(std::vector<std::vector<real_t>> A, std::vector<std::vector<real_t>> B);
 
 	std::vector<std::vector<real_t>> transpose(std::vector<std::vector<real_t>> A);
-
 	std::vector<std::vector<real_t>> scalarMultiply(real_t scalar, std::vector<std::vector<real_t>> A);
-
 	std::vector<std::vector<real_t>> scalarAdd(real_t scalar, std::vector<std::vector<real_t>> A);
 
 	std::vector<std::vector<real_t>> log(std::vector<std::vector<real_t>> A);
-
 	std::vector<std::vector<real_t>> log10(std::vector<std::vector<real_t>> A);
-
 	std::vector<std::vector<real_t>> exp(std::vector<std::vector<real_t>> A);
-
 	std::vector<std::vector<real_t>> erf(std::vector<std::vector<real_t>> A);
-
 	std::vector<std::vector<real_t>> exponentiate(std::vector<std::vector<real_t>> A, real_t p);
-
 	std::vector<std::vector<real_t>> sqrt(std::vector<std::vector<real_t>> A);
-
 	std::vector<std::vector<real_t>> cbrt(std::vector<std::vector<real_t>> A);
 
 	std::vector<std::vector<real_t>> matrixPower(std::vector<std::vector<real_t>> A, int n);
@@ -75,13 +67,14 @@ public:
 	std::vector<std::vector<real_t>> pinverse(std::vector<std::vector<real_t>> A);
 
 	std::vector<std::vector<real_t>> zeromat(int n, int m);
-
 	std::vector<std::vector<real_t>> onemat(int n, int m);
-
 	std::vector<std::vector<real_t>> full(int n, int m, int k);
 
-	std::vector<std::vector<real_t>> sin(std::vector<std::vector<real_t>> A);
+	Ref<MLPPMatrix> zeromatm(int n, int m);
+	Ref<MLPPMatrix> onematm(int n, int m);
+	Ref<MLPPMatrix> fullm(int n, int m, int k);
 
+	std::vector<std::vector<real_t>> sin(std::vector<std::vector<real_t>> A);
 	std::vector<std::vector<real_t>> cos(std::vector<std::vector<real_t>> A);
 
 	std::vector<std::vector<real_t>> rotate(std::vector<std::vector<real_t>> A, real_t theta, int axis = -1);
@@ -162,6 +155,7 @@ public:
 	std::vector<real_t> hadamard_product(std::vector<real_t> a, std::vector<real_t> b);
 
 	std::vector<real_t> elementWiseDivision(std::vector<real_t> a, std::vector<real_t> b);
+	Ref<MLPPVector> element_wise_division(const Ref<MLPPVector> &a, const Ref<MLPPVector> &b);
 
 	std::vector<real_t> scalarMultiply(real_t scalar, std::vector<real_t> a);
 	Ref<MLPPVector> scalar_multiplynv(real_t scalar, const Ref<MLPPVector> &a);
@@ -181,19 +175,21 @@ public:
 	std::vector<real_t> subtractMatrixRows(std::vector<real_t> a, std::vector<std::vector<real_t>> B);
 
 	std::vector<real_t> log(std::vector<real_t> a);
-
 	std::vector<real_t> log10(std::vector<real_t> a);
-
 	std::vector<real_t> exp(std::vector<real_t> a);
-
 	std::vector<real_t> erf(std::vector<real_t> a);
-
 	std::vector<real_t> exponentiate(std::vector<real_t> a, real_t p);
-
 	std::vector<real_t> sqrt(std::vector<real_t> a);
-
 	std::vector<real_t> cbrt(std::vector<real_t> a);
 
+	Ref<MLPPVector> logv(const Ref<MLPPVector> &a);
+	Ref<MLPPVector> log10v(const Ref<MLPPVector> &a);
+	Ref<MLPPVector> expv(const Ref<MLPPVector> &a);
+	Ref<MLPPVector> erfv(const Ref<MLPPVector> &a);
+	Ref<MLPPVector> exponentiatev(const Ref<MLPPVector> &a, real_t p);
+	Ref<MLPPVector> sqrtv(const Ref<MLPPVector> &a);
+	Ref<MLPPVector> cbrtv(const Ref<MLPPVector> &a);
+
 	real_t dot(std::vector<real_t> a, std::vector<real_t> b);
 
 	std::vector<real_t> cross(std::vector<real_t> a, std::vector<real_t> b);
@@ -201,13 +197,15 @@ public:
 	std::vector<real_t> abs(std::vector<real_t> a);
 
 	std::vector<real_t> zerovec(int n);
-
 	std::vector<real_t> onevec(int n);
+	std::vector<real_t> full(int n, int k);
+
+	Ref<MLPPVector> zerovecv(int n);
+	Ref<MLPPVector> onevecv(int n);
+	Ref<MLPPVector> fullv(int n, int k);
 
 	std::vector<std::vector<real_t>> diag(std::vector<real_t> a);
 
-	std::vector<real_t> full(int n, int k);
-
 	std::vector<real_t> sin(std::vector<real_t> a);
 
 	std::vector<real_t> cos(std::vector<real_t> a);