diff --git a/mlpp/cost/cost.cpp b/mlpp/cost/cost.cpp
index 9456325..9ae1450 100644
--- a/mlpp/cost/cost.cpp
+++ b/mlpp/cost/cost.cpp
@@ -10,6 +10,582 @@
 #include <cmath>
 #include <iostream>
 
+real_t MLPPCost::msev(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y) {
+	int y_hat_size = y_hat->size();
+
+	ERR_FAIL_COND_V(y_hat_size != y->size(), 0);
+
+	const real_t *y_hat_ptr = y_hat->ptr();
+	const real_t *y_ptr = y->ptr();
+
+	real_t sum = 0;
+	for (int i = 0; i < y_hat_size; ++i) {
+		sum += (y_hat_ptr[i] - y_ptr[i]) * (y_hat_ptr[i] - y_ptr[i]);
+	}
+
+	return sum / 2 * y_hat_size;
+}
+real_t MLPPCost::msem(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y) {
+	int y_hat_data_size = y_hat->data_size();
+
+	ERR_FAIL_COND_V(y_hat_data_size != y->data_size(), 0);
+
+	const real_t *y_hat_ptr = y_hat->ptr();
+	const real_t *y_ptr = y->ptr();
+
+	real_t sum = 0;
+	for (int i = 0; i < y_hat_data_size; ++i) {
+		sum += (y_hat_ptr[i] - y_ptr[i]) * (y_hat_ptr[i] - y_ptr[i]);
+	}
+
+	return sum / 2.0 * static_cast<real_t>(y_hat_data_size);
+}
+
+Ref<MLPPVector> MLPPCost::mse_derivv(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y) {
+	MLPPLinAlg alg;
+	return alg.subtractionnv(y_hat, y);
+}
+
+Ref<MLPPMatrix> MLPPCost::mse_derivm(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y) {
+	MLPPLinAlg alg;
+	return alg.subtractionm(y_hat, y);
+}
+
+real_t MLPPCost::rmsev(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y) {
+	int y_hat_size = y_hat->size();
+
+	ERR_FAIL_COND_V(y_hat_size != y->size(), 0);
+
+	const real_t *y_hat_ptr = y_hat->ptr();
+	const real_t *y_ptr = y->ptr();
+
+	real_t sum = 0;
+	for (int i = 0; i < y_hat_size; ++i) {
+		sum += (y_hat_ptr[i] - y_ptr[i]) * (y_hat_ptr[i] - y_ptr[i]);
+	}
+
+	return Math::sqrt(sum / static_cast<real_t>(y_hat_size));
+}
+real_t MLPPCost::rmsem(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y) {
+	int y_hat_data_size = y_hat->data_size();
+
+	ERR_FAIL_COND_V(y_hat_data_size != y->data_size(), 0);
+
+	const real_t *y_hat_ptr = y_hat->ptr();
+	const real_t *y_ptr = y->ptr();
+
+	real_t sum = 0;
+	for (int i = 0; i < y_hat_data_size; ++i) {
+		sum += (y_hat_ptr[i] - y_ptr[i]) * (y_hat_ptr[i] - y_ptr[i]);
+	}
+
+	return Math::sqrt(sum / static_cast<real_t>(y_hat->size().y));
+}
+
+Ref<MLPPVector> MLPPCost::rmse_derivv(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y) {
+	MLPPLinAlg alg;
+
+	return alg.scalar_multiplynv(1 / (2.0 * Math::sqrt(msev(y_hat, y))), mse_derivv(y_hat, y));
+}
+Ref<MLPPMatrix> MLPPCost::rmse_derivm(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y) {
+	MLPPLinAlg alg;
+
+	return alg.scalar_multiplym(1 / (2.0 / Math::sqrt(msem(y_hat, y))), mse_derivm(y_hat, y));
+}
+
+real_t MLPPCost::maev(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y) {
+	int y_hat_size = y_hat->size();
+
+	ERR_FAIL_COND_V(y_hat_size != y->size(), 0);
+
+	const real_t *y_hat_ptr = y_hat->ptr();
+	const real_t *y_ptr = y->ptr();
+
+	real_t sum = 0;
+	for (int i = 0; i < y_hat_size; i++) {
+		sum += ABS((y_hat_ptr[i] - y_ptr[i]));
+	}
+	return sum / static_cast<real_t>(y_hat_size);
+}
+real_t MLPPCost::maem(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y) {
+	int y_hat_data_size = y_hat->data_size();
+
+	ERR_FAIL_COND_V(y_hat_data_size != y->data_size(), 0);
+
+	const real_t *y_hat_ptr = y_hat->ptr();
+	const real_t *y_ptr = y->ptr();
+
+	real_t sum = 0;
+	for (int i = 0; i < y_hat_data_size; ++i) {
+		sum += ABS((y_hat_ptr[i] - y_ptr[i]));
+	}
+
+	return sum / static_cast<real_t>(y_hat_data_size);
+}
+
+Ref<MLPPVector> MLPPCost::mae_derivv(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y) {
+	int y_hat_size = y_hat->size();
+
+	const real_t *y_hat_ptr = y_hat->ptr();
+
+	Ref<MLPPVector> deriv;
+	deriv.instance();
+	deriv->resize(y_hat_size);
+	real_t *deriv_ptr = deriv->ptrw();
+
+	for (int i = 0; i < y_hat_size; ++i) {
+		int y_hat_ptr_i = y_hat_ptr[i];
+
+		if (y_hat_ptr_i < 0) {
+			deriv_ptr[i] = -1;
+		} else if (y_hat_ptr_i == 0) {
+			deriv_ptr[i] = 0;
+		} else {
+			deriv_ptr[i] = 1;
+		}
+	}
+
+	return deriv;
+}
+Ref<MLPPMatrix> MLPPCost::mae_derivm(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y) {
+	int y_hat_data_size = y_hat->data_size();
+
+	const real_t *y_hat_ptr = y_hat->ptr();
+
+	Ref<MLPPMatrix> deriv;
+	deriv.instance();
+	deriv->resize(y_hat->size());
+	real_t *deriv_ptr = deriv->ptrw();
+
+	for (int i = 0; i < y_hat_data_size; ++i) {
+		int y_hat_ptr_i = y_hat_ptr[i];
+
+		if (y_hat_ptr_i < 0) {
+			deriv_ptr[i] = -1;
+		} else if (y_hat_ptr_i == 0) {
+			deriv_ptr[i] = 0;
+		} else {
+			deriv_ptr[i] = 1;
+		}
+	}
+
+	return deriv;
+}
+
+real_t MLPPCost::mbev(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y) {
+	int y_hat_size = y_hat->size();
+
+	ERR_FAIL_COND_V(y_hat_size != y->size(), 0);
+
+	const real_t *y_hat_ptr = y_hat->ptr();
+	const real_t *y_ptr = y->ptr();
+
+	real_t sum = 0;
+	for (int i = 0; i < y_hat_size; ++i) {
+		sum += (y_hat_ptr[i] - y_ptr[i]);
+	}
+
+	return sum / static_cast<real_t>(y_hat_size);
+}
+real_t MLPPCost::mbem(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y) {
+	int y_hat_data_size = y_hat->data_size();
+
+	ERR_FAIL_COND_V(y_hat_data_size != y->data_size(), 0);
+
+	const real_t *y_hat_ptr = y_hat->ptr();
+	const real_t *y_ptr = y->ptr();
+
+	real_t sum = 0;
+	for (int i = 0; i < y_hat_data_size; ++i) {
+		sum += (y_hat_ptr[i] - y_ptr[i]);
+	}
+
+	return sum / static_cast<real_t>(y_hat_data_size);
+}
+
+Ref<MLPPVector> MLPPCost::mbe_derivv(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y) {
+	MLPPLinAlg alg;
+	return alg.onevecv(y_hat->size());
+}
+Ref<MLPPMatrix> MLPPCost::mbe_derivm(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y) {
+	MLPPLinAlg alg;
+	return alg.onematm(y_hat->size().x, y_hat->size().y);
+}
+
+// Classification Costs
+real_t MLPPCost::log_lossv(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y) {
+	int y_hat_size = y_hat->size();
+
+	ERR_FAIL_COND_V(y_hat_size != y->size(), 0);
+
+	const real_t *y_hat_ptr = y_hat->ptr();
+	const real_t *y_ptr = y->ptr();
+
+	real_t sum = 0;
+	real_t eps = 1e-8;
+	for (int i = 0; i < y_hat_size; ++i) {
+		sum += -(y_ptr[i] * Math::log(y_hat_ptr[i] + eps) + (1 - y_ptr[i]) * Math::log(1 - y_hat_ptr[i] + eps));
+	}
+
+	return sum / static_cast<real_t>(y_hat_size);
+}
+
+real_t MLPPCost::log_lossm(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y) {
+	int y_hat_data_size = y_hat->data_size();
+
+	ERR_FAIL_COND_V(y_hat_data_size != y->data_size(), 0);
+
+	const real_t *y_hat_ptr = y_hat->ptr();
+	const real_t *y_ptr = y->ptr();
+
+	real_t sum = 0;
+	real_t eps = 1e-8;
+	for (int i = 0; i < y_hat_data_size; ++i) {
+		sum += -(y_ptr[i] * Math::log(y_hat_ptr[i] + eps) + (1 - y_ptr[i]) * Math::log(1 - y_hat_ptr[i] + eps));
+	}
+
+	return sum / static_cast<real_t>(y_hat_data_size);
+}
+
+Ref<MLPPVector> MLPPCost::log_loss_derivv(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y) {
+	MLPPLinAlg alg;
+	return alg.additionnv(
+			alg.scalar_multiplynv(-1, alg.element_wise_division(y, y_hat)),
+			alg.element_wise_division(alg.scalar_multiplynv(-1, alg.scalar_addnv(-1, y)), alg.scalar_multiplynv(-1, alg.scalar_addnv(-1, y_hat))));
+}
+
+Ref<MLPPMatrix> MLPPCost::log_loss_derivm(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y) {
+	MLPPLinAlg alg;
+	return alg.additionm(
+			alg.scalar_multiplym(-1, alg.element_wise_divisionm(y, y_hat)),
+			alg.element_wise_divisionm(alg.scalar_multiplym(-1, alg.scalar_addm(-1, y)), alg.scalar_multiplym(-1, alg.scalar_addm(-1, y_hat))));
+}
+
+real_t MLPPCost::cross_entropyv(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y) {
+	int y_hat_size = y_hat->size();
+
+	ERR_FAIL_COND_V(y_hat_size != y->size(), 0);
+
+	const real_t *y_hat_ptr = y_hat->ptr();
+	const real_t *y_ptr = y->ptr();
+
+	real_t sum = 0;
+	for (int i = 0; i < y_hat_size; ++i) {
+		sum += y_ptr[i] * Math::log(y_hat_ptr[i]);
+	}
+
+	return -1 * sum;
+}
+real_t MLPPCost::cross_entropym(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y) {
+	int y_hat_data_size = y_hat->data_size();
+
+	ERR_FAIL_COND_V(y_hat_data_size != y->data_size(), 0);
+
+	const real_t *y_hat_ptr = y_hat->ptr();
+	const real_t *y_ptr = y->ptr();
+
+	real_t sum = 0;
+	for (int i = 0; i < y_hat_data_size; ++i) {
+		sum += y_ptr[i] * Math::log(y_hat_ptr[i]);
+	}
+
+	return -1 * sum;
+}
+
+Ref<MLPPVector> MLPPCost::cross_entropy_derivv(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y) {
+	MLPPLinAlg alg;
+	return alg.scalar_multiplynv(-1, alg.element_wise_division(y, y_hat));
+}
+Ref<MLPPMatrix> MLPPCost::cross_entropy_derivm(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y) {
+	MLPPLinAlg alg;
+	return alg.scalar_multiplym(-1, alg.element_wise_divisionm(y, y_hat));
+}
+
+real_t MLPPCost::huber_lossv(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y, real_t delta) {
+	int y_hat_size = y_hat->size();
+
+	ERR_FAIL_COND_V(y_hat_size != y->size(), 0);
+
+	const real_t *y_hat_ptr = y_hat->ptr();
+	const real_t *y_ptr = y->ptr();
+
+	MLPPLinAlg alg;
+	real_t sum = 0;
+	for (int i = 0; i < y_hat_size; ++i) {
+		if (ABS(y_ptr[i] - y_hat_ptr[i]) <= delta) {
+			sum += (y_ptr[i] - y_hat_ptr[i]) * (y_ptr[i] - y_hat_ptr[i]);
+		} else {
+			sum += 2 * delta * ABS(y_ptr[i] - y_hat_ptr[i]) - delta * delta;
+		}
+	}
+
+	return sum;
+}
+real_t MLPPCost::huber_lossm(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y, real_t delta) {
+	int y_hat_data_size = y_hat->data_size();
+
+	ERR_FAIL_COND_V(y_hat_data_size != y->data_size(), 0);
+
+	const real_t *y_hat_ptr = y_hat->ptr();
+	const real_t *y_ptr = y->ptr();
+
+	MLPPLinAlg alg;
+	real_t sum = 0;
+	for (int i = 0; i < y_hat_data_size; ++i) {
+		if (ABS(y_ptr[i] - y_hat_ptr[i]) <= delta) {
+			sum += (y_ptr[i] - y_hat_ptr[i]) * (y_ptr[i] - y_hat_ptr[i]);
+		} else {
+			sum += 2 * delta * ABS(y_ptr[i] - y_hat_ptr[i]) - delta * delta;
+		}
+	}
+
+	return sum;
+}
+
+Ref<MLPPVector> MLPPCost::huber_loss_derivv(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y, real_t delta) {
+	int y_hat_size = y_hat->size();
+
+	ERR_FAIL_COND_V(y_hat_size != y->size(), 0);
+
+	const real_t *y_hat_ptr = y_hat->ptr();
+	const real_t *y_ptr = y->ptr();
+
+	MLPPLinAlg alg;
+	real_t sum = 0;
+
+	Ref<MLPPVector> deriv;
+	deriv.instance();
+	deriv->resize(y_hat->size());
+
+	real_t *deriv_ptr = deriv->ptrw();
+
+	for (int i = 0; i < y_hat_size; ++i) {
+		if (ABS(y_ptr[i] - y_hat_ptr[i]) <= delta) {
+			deriv_ptr[i] = (-(y_ptr[i] - y_hat_ptr[i]));
+		} else {
+			if (y_hat_ptr[i] > 0 || y_hat_ptr[i] < 0) {
+				deriv_ptr[i] = (2 * delta * (y_hat_ptr[i] / ABS(y_hat_ptr[i])));
+			} else {
+				deriv_ptr[i] = (0);
+			}
+		}
+	}
+
+	return deriv;
+}
+Ref<MLPPMatrix> MLPPCost::huber_loss_derivm(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y, real_t delta) {
+	int y_hat_data_size = y_hat->data_size();
+
+	ERR_FAIL_COND_V(y_hat_data_size != y->data_size(), 0);
+
+	const real_t *y_hat_ptr = y_hat->ptr();
+	const real_t *y_ptr = y->ptr();
+
+	MLPPLinAlg alg;
+	real_t sum = 0;
+
+	Ref<MLPPMatrix> deriv;
+	deriv.instance();
+	deriv->resize(y_hat->size());
+
+	real_t *deriv_ptr = deriv->ptrw();
+
+	for (int i = 0; i < y_hat_data_size; ++i) {
+		if (ABS(y_ptr[i] - y_hat_ptr[i]) <= delta) {
+			deriv_ptr[i] = (-(y_ptr[i] - y_hat_ptr[i]));
+		} else {
+			if (y_hat_ptr[i] > 0 || y_hat_ptr[i] < 0) {
+				deriv_ptr[i] = (2 * delta * (y_hat_ptr[i] / ABS(y_hat_ptr[i])));
+			} else {
+				deriv_ptr[i] = (0);
+			}
+		}
+	}
+
+	return deriv;
+}
+
+real_t MLPPCost::hinge_lossv(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y) {
+	int y_hat_size = y_hat->size();
+
+	ERR_FAIL_COND_V(y_hat_size != y->size(), 0);
+
+	const real_t *y_hat_ptr = y_hat->ptr();
+	const real_t *y_ptr = y->ptr();
+
+	real_t sum = 0;
+	for (int i = 0; i < y_hat_size; ++i) {
+		sum += MAX(0, 1 - y_ptr[i] * y_hat_ptr[i]);
+	}
+
+	return sum / static_cast<real_t>(y_hat_size);
+}
+real_t MLPPCost::hinge_lossm(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y) {
+	int y_hat_data_size = y_hat->data_size();
+
+	ERR_FAIL_COND_V(y_hat_data_size != y->data_size(), 0);
+
+	const real_t *y_hat_ptr = y_hat->ptr();
+	const real_t *y_ptr = y->ptr();
+
+	real_t sum = 0;
+	for (int i = 0; i < y_hat_data_size; ++i) {
+		sum += MAX(0, 1 - y_ptr[i] * y_hat_ptr[i]);
+	}
+
+	return sum / static_cast<real_t>(y_hat_data_size);
+}
+
+Ref<MLPPVector> MLPPCost::hinge_loss_derivv(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y) {
+	int y_hat_size = y_hat->size();
+
+	ERR_FAIL_COND_V(y_hat_size != y->size(), 0);
+
+	const real_t *y_hat_ptr = y_hat->ptr();
+	const real_t *y_ptr = y->ptr();
+
+	Ref<MLPPVector> deriv;
+	deriv.instance();
+	deriv->resize(y_hat->size());
+
+	real_t *deriv_ptr = deriv->ptrw();
+
+	for (int i = 0; i < y_hat_size; ++i) {
+		if (1 - y_ptr[i] * y_hat_ptr[i] > 0) {
+			deriv_ptr[i] = -y_ptr[i];
+		} else {
+			deriv_ptr[i] = 0;
+		}
+	}
+
+	return deriv;
+}
+Ref<MLPPMatrix> MLPPCost::hinge_loss_derivm(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y) {
+	int y_hat_data_size = y_hat->data_size();
+
+	ERR_FAIL_COND_V(y_hat_data_size != y->data_size(), 0);
+
+	const real_t *y_hat_ptr = y_hat->ptr();
+	const real_t *y_ptr = y->ptr();
+
+	Ref<MLPPMatrix> deriv;
+	deriv.instance();
+	deriv->resize(y_hat->size());
+
+	real_t *deriv_ptr = deriv->ptrw();
+
+	for (int i = 0; i < y_hat_data_size; ++i) {
+		if (1 - y_ptr[i] * y_hat_ptr[i] > 0) {
+			deriv_ptr[i] = -y_ptr[i];
+		} else {
+			deriv_ptr[i] = 0;
+		}
+	}
+
+	return deriv;
+}
+
+real_t MLPPCost::hinge_losswv(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y, const Ref<MLPPVector> &weights, real_t C) {
+	MLPPLinAlg alg;
+	MLPPReg regularization;
+
+	return C * hinge_lossv(y_hat, y) + regularization.reg_termv(weights, 1, 0, MLPPReg::REGULARIZATION_TYPE_RIDGE);
+}
+real_t MLPPCost::hinge_losswm(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y, const Ref<MLPPMatrix> &weights, real_t C) {
+	MLPPLinAlg alg;
+	MLPPReg regularization;
+
+	return C * hinge_lossm(y_hat, y) + regularization.reg_termv(weights, 1, 0, MLPPReg::REGULARIZATION_TYPE_RIDGE);
+}
+
+Ref<MLPPVector> MLPPCost::hinge_loss_derivwv(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y, real_t C) {
+	MLPPLinAlg alg;
+	MLPPReg regularization;
+
+	return alg.scalar_multiplynv(C, hinge_loss_derivv(y_hat, y));
+}
+Ref<MLPPMatrix> MLPPCost::hinge_loss_derivwm(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y, real_t C) {
+	MLPPLinAlg alg;
+	MLPPReg regularization;
+
+	return alg.scalar_multiplym(C, hinge_loss_derivm(y_hat, y));
+}
+
+real_t MLPPCost::wasserstein_lossv(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y) {
+	int y_hat_size = y_hat->size();
+
+	ERR_FAIL_COND_V(y_hat_size != y->size(), 0);
+
+	const real_t *y_hat_ptr = y_hat->ptr();
+	const real_t *y_ptr = y->ptr();
+
+	real_t sum = 0;
+	for (int i = 0; i < y_hat_size; ++i) {
+		sum += y_hat_ptr[i] * y_ptr[i];
+	}
+
+	return -sum / static_cast<real_t>(y_hat_size);
+}
+real_t MLPPCost::wasserstein_lossm(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y) {
+	int y_hat_data_size = y_hat->data_size();
+
+	ERR_FAIL_COND_V(y_hat_data_size != y->data_size(), 0);
+
+	const real_t *y_hat_ptr = y_hat->ptr();
+	const real_t *y_ptr = y->ptr();
+
+	real_t sum = 0;
+	for (int i = 0; i < y_hat_data_size; ++i) {
+		sum += y_hat_ptr[i] * y_ptr[i];
+	}
+
+	return -sum / static_cast<real_t>(y_hat_data_size);
+}
+
+Ref<MLPPVector> MLPPCost::wasserstein_loss_derivv(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y) {
+	MLPPLinAlg alg;
+
+	return alg.scalar_multiplynv(-1, y); // Simple.
+}
+Ref<MLPPMatrix> MLPPCost::wasserstein_loss_derivm(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y) {
+	MLPPLinAlg alg;
+
+	return alg.scalar_multiplym(-1, y); // Simple.
+}
+
+real_t MLPPCost::dual_form_svm(const Ref<MLPPVector> &alpha, const Ref<MLPPMatrix> &X, const Ref<MLPPVector> &y) {
+	MLPPLinAlg alg;
+
+	Ref<MLPPMatrix> Y = alg.diagm(y); // Y is a diagnoal matrix. Y[i][j] = y[i] if i = i, else Y[i][j] = 0. Yt = Y.
+	Ref<MLPPMatrix> K = alg.matmultm(X, alg.transposem(X)); // TO DO: DON'T forget to add non-linear kernelizations.
+	Ref<MLPPMatrix> Q = alg.matmultm(alg.matmultm(alg.transposem(Y), K), Y);
+
+	Ref<MLPPMatrix> alpha_m;
+	alpha_m.instance();
+	alpha_m->resize(Size2i(alpha->size(), 1));
+	alpha_m->set_row_mlpp_vector(0, alpha);
+
+	Ref<MLPPMatrix> alpha_m_res = alg.matmultm(alg.matmultm(alpha_m, Q), alg.transposem(alpha_m));
+
+	real_t alphaQ = alpha_m_res->get_element(0, 0);
+	Ref<MLPPVector> one = alg.onevecv(alpha->size());
+
+	return -alg.dotv(one, alpha) + 0.5 * alphaQ;
+}
+
+Ref<MLPPVector> MLPPCost::dual_form_svm_deriv(const Ref<MLPPVector> &alpha, const Ref<MLPPMatrix> &X, const Ref<MLPPVector> &y) {
+	MLPPLinAlg alg;
+
+	Ref<MLPPMatrix> Y = alg.diagm(y); // Y is a diagnoal matrix. Y[i][j] = y[i] if i = i, else Y[i][j] = 0. Yt = Y.
+	Ref<MLPPMatrix> K = alg.matmultm(X, alg.transposem(X)); // TO DO: DON'T forget to add non-linear kernelizations.
+	Ref<MLPPMatrix> Q = alg.matmultm(alg.matmultm(alg.transposem(Y), K), Y);
+	Ref<MLPPVector> alphaQDeriv = alg.mat_vec_multv(Q, alpha);
+	Ref<MLPPVector> one = alg.onevecv(alpha->size());
+
+	return alg.subtractionm(alphaQDeriv, one);
+}
+
+// ======  OLD  ======
 
 real_t MLPPCost::MSE(std::vector<real_t> y_hat, std::vector<real_t> y) {
 	real_t sum = 0;
diff --git a/mlpp/cost/cost.h b/mlpp/cost/cost.h
index dcfd180..aa5ce26 100644
--- a/mlpp/cost/cost.h
+++ b/mlpp/cost/cost.h
@@ -10,11 +10,85 @@
 
 #include "core/math/math_defs.h"
 
+#include "core/object/reference.h"
+
 #include <vector>
 
+#include "../lin_alg/mlpp_matrix.h"
+#include "../lin_alg/mlpp_vector.h"
+
+//void set_weights(const Ref<MLPPMatrix> &val);
+//void set_bias(const Ref<MLPPVector> &val);
+
+class MLPPCost : public Reference {
+	GDCLASS(MLPPCost, Reference);
 
-class MLPPCost {
 public:
+	real_t msev(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y);
+	real_t msem(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y);
+
+	Ref<MLPPVector> mse_derivv(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y);
+	Ref<MLPPMatrix> mse_derivm(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y);
+
+	real_t rmsev(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y);
+	real_t rmsem(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y);
+
+	Ref<MLPPVector> rmse_derivv(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y);
+	Ref<MLPPMatrix> rmse_derivm(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y);
+
+	real_t maev(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y);
+	real_t maem(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y);
+
+	Ref<MLPPVector> mae_derivv(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y);
+	Ref<MLPPMatrix> mae_derivm(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y);
+
+	real_t mbev(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y);
+	real_t mbem(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y);
+
+	Ref<MLPPVector> mbe_derivv(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y);
+	Ref<MLPPMatrix> mbe_derivm(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y);
+
+	// Classification Costs
+	real_t log_lossv(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y);
+	real_t log_lossm(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y);
+
+	Ref<MLPPVector> log_loss_derivv(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y);
+	Ref<MLPPMatrix> log_loss_derivm(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y);
+
+	real_t cross_entropyv(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y);
+	real_t cross_entropym(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y);
+
+	Ref<MLPPVector> cross_entropy_derivv(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y);
+	Ref<MLPPMatrix> cross_entropy_derivm(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y);
+
+	real_t huber_lossv(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y, real_t delta);
+	real_t huber_lossm(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y, real_t delta);
+
+	Ref<MLPPVector> huber_loss_derivv(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y, real_t delta);
+	Ref<MLPPMatrix> huber_loss_derivm(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y, real_t delta);
+
+	real_t hinge_lossv(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y);
+	real_t hinge_lossm(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y);
+
+	Ref<MLPPVector> hinge_loss_derivv(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y);
+	Ref<MLPPMatrix> hinge_loss_derivm(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y);
+
+	real_t hinge_losswv(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y, const Ref<MLPPVector> &weights, real_t C);
+	real_t hinge_losswm(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y, const Ref<MLPPMatrix> &weights, real_t C);
+
+	Ref<MLPPVector> hinge_loss_derivwv(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y, real_t C);
+	Ref<MLPPMatrix> hinge_loss_derivwm(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y, real_t C);
+
+	real_t wasserstein_lossv(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y);
+	real_t wasserstein_lossm(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y);
+
+	Ref<MLPPVector> wasserstein_loss_derivv(const Ref<MLPPVector> &y_hat, const Ref<MLPPVector> &y);
+	Ref<MLPPMatrix> wasserstein_loss_derivm(const Ref<MLPPMatrix> &y_hat, const Ref<MLPPMatrix> &y);
+
+	real_t dual_form_svm(const Ref<MLPPVector> &alpha, const Ref<MLPPMatrix> &X, const Ref<MLPPVector> &y); // TO DO: DON'T forget to add non-linear kernelizations.
+
+	Ref<MLPPVector> dual_form_svm_deriv(const Ref<MLPPVector> &alpha, const Ref<MLPPMatrix> &X, const Ref<MLPPVector> &y);
+
 	// Regression Costs
 	real_t MSE(std::vector<real_t> y_hat, std::vector<real_t> y);
 	real_t MSE(std::vector<std::vector<real_t>> y_hat, std::vector<std::vector<real_t>> y);
@@ -80,9 +154,6 @@ public:
 	real_t dualFormSVM(std::vector<real_t> alpha, std::vector<std::vector<real_t>> X, std::vector<real_t> y); // TO DO: DON'T forget to add non-linear kernelizations.
 
 	std::vector<real_t> dualFormSVMDeriv(std::vector<real_t> alpha, std::vector<std::vector<real_t>> X, std::vector<real_t> y);
-
-private:
 };
 
-
 #endif /* Cost_hpp */
diff --git a/mlpp/lin_alg/lin_alg.cpp b/mlpp/lin_alg/lin_alg.cpp
index da635cf..5832232 100644
--- a/mlpp/lin_alg/lin_alg.cpp
+++ b/mlpp/lin_alg/lin_alg.cpp
@@ -1867,6 +1867,21 @@ real_t MLPPLinAlg::dot(std::vector<real_t> a, std::vector<real_t> b) {
 	return c;
 }
 
+real_t MLPPLinAlg::dotv(const Ref<MLPPVector> &a, const Ref<MLPPVector> &b) {
+	int a_size = a->size();
+
+	ERR_FAIL_COND_V(a_size != b->size(), 0);
+
+	const real_t *a_ptr = a->ptr();
+	const real_t *b_ptr = b->ptr();
+
+	real_t c = 0;
+	for (int i = 0; i < a_size; ++i) {
+		c += a_ptr[i] * b_ptr[i];
+	}
+	return c;
+}
+
 std::vector<real_t> MLPPLinAlg::cross(std::vector<real_t> a, std::vector<real_t> b) {
 	// Cross products exist in R^7 also. Though, I will limit it to R^3 as Wolfram does this.
 	std::vector<std::vector<real_t>> mat = { onevec(3), a, b };
@@ -1905,6 +1920,25 @@ std::vector<std::vector<real_t>> MLPPLinAlg::diag(std::vector<real_t> a) {
 	return B;
 }
 
+Ref<MLPPVector> MLPPLinAlg::diagm(const Ref<MLPPVector> &a) {
+	int a_size = a->size();
+
+	Ref<MLPPMatrix> B;
+	B.instance();
+
+	B->resize(Size2i(a_size, a_size));
+	B->fill(0);
+
+	const real_t *a_ptr = a->ptr();
+	real_t *b_ptr = B->ptrw();
+
+	for (int i = 0; i < a_size; ++i) {
+		b_ptr[B->calculate_index(i, i)] = a_ptr[i];
+	}
+
+	return B;
+}
+
 std::vector<real_t> MLPPLinAlg::full(int n, int k) {
 	std::vector<real_t> full;
 	full.resize(n);
diff --git a/mlpp/lin_alg/lin_alg.h b/mlpp/lin_alg/lin_alg.h
index 1f4acf6..26e02bb 100644
--- a/mlpp/lin_alg/lin_alg.h
+++ b/mlpp/lin_alg/lin_alg.h
@@ -213,6 +213,7 @@ public:
 	Ref<MLPPVector> cbrtv(const Ref<MLPPVector> &a);
 
 	real_t dot(std::vector<real_t> a, std::vector<real_t> b);
+	real_t dotv(const Ref<MLPPVector> &a, const Ref<MLPPVector> &b);
 
 	std::vector<real_t> cross(std::vector<real_t> a, std::vector<real_t> b);
 
@@ -229,6 +230,7 @@ public:
 	Ref<MLPPVector> fullv(int n, int k);
 
 	std::vector<std::vector<real_t>> diag(std::vector<real_t> a);
+	Ref<MLPPVector> diagm(const Ref<MLPPVector> &a);
 
 	std::vector<real_t> sin(std::vector<real_t> a);
 	std::vector<real_t> cos(std::vector<real_t> a);