{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear Algebra\n",
    "\n",
    "Smile provides an MATLAB like environment. In the simplest case, you can use it as a calculator. Besides `java.lang.Math` functions, Smile provides many other important mathematical functions such as `logistic`, `factorial`, `choose`, etc. in `smile.math.MathEx` class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import $ivy.`com.github.haifengl::smile-scala:3.1.0`\n",
    "import $ivy.`org.slf4j:slf4j-simple:2.0.13` \n",
    "\n",
    "import java.lang.Math._\n",
    "import smile.math._\n",
    "import smile.math.matrix._\n",
    "import smile.math.MathEx._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "2+3\n",
    "log(2)\n",
    "logistic(3.0)\n",
    "choose(10, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Special Functions\n",
    "\n",
    "Smile implements `beta`, `erf`, `gamma` and their related functions. Special functions are particular mathematical functions which have more or less established names and notations due to their importance in mathematical analysis, functional analysis, physics, or other applications. Many special functions appear as solutions of differential equations or integrals of elementary functions. For example, the error function `erf` (also called the Gauss error function) is a special function of sigmoid shape which occurs in probability, statistics, materials science, and partial differential equations. The complementary error function, denoted `erfc`, is defined as `erfc(x) = 1 - erf(x)`. The error function and complementary error function are special cases of the incomplete `gamma` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "erf(1.0)\n",
    "digamma(1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Vector\n",
    "\n",
    "Common arithmetic operations on vectors and scalars are similar as in R and Matlab."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "val x = c(1.0, 2.0, 3.0)\n",
    "val y = c(4.0, 3.0, 2.0)\n",
    "\n",
    "x + y\n",
    "\n",
    "1.5 * x - 3.0 * y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that these operations are lazy. The computation is only performed when the results are needed, e.g. when the expression is used where a vector is expected.\n",
    "\n",
    "For a vector, there are multiple functions to calculate its norm such as `norm` (L2 norm), `norm1` (L1 norm), `norm2` (L2 norm), `normInf` (infinity norm), `normFro` (Frobenius norm). We can also `standardize` a vector to mean 0 and variance 1, `unitize` it so that L2 norm be 1, or `unitize1` it so that L1 norm be 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "norm(x)\n",
    "\n",
    "unitize(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a pair of vectors, we can calculate the dot product, distance, divergence, covariance, and correlations with `dot`, `distance`, `kld` (Kullback-Leibler Divergence), `jsd` (Jensen-Shannon Divergence), `cov`, `cor` (Pearson Correlation), `spearman` (Spearman Rank Correlation Coefficient), `kendall` (Kendall Tau Rank Correlation Coefficient)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dot(x, y)\n",
    "\n",
    "cov(x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Matrix\n",
    "\n",
    "Like Matlab, we can use `eye`, `zeros` and `ones` to create identity, zero, or all-ones matrix, respectively. To create a matrix from 2-dimensional array, we can use the constructor `matrix` or the ~ operator. The ~ operator can be applied to 1-dimensional array too, which creates a single column matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "val a = matrix(Array(\n",
    "        Array(0.7220180, 0.07121225, 0.6881997),\n",
    "        Array(-0.2648886, -0.89044952, 0.3700456),\n",
    "        Array(-0.6391588, 0.44947578, 0.6240573))\n",
    "    )\n",
    "\n",
    "val b = matrix(Array(\n",
    "        Array(0.6881997, -0.07121225, 0.7220180),\n",
    "        Array(0.3700456, 0.89044952, -0.2648886),\n",
    "        Array(0.6240573, -0.44947578, -0.6391588))\n",
    "    )\n",
    "\n",
    "val C = Array(\n",
    "        Array(0.9527204, -0.2973347, 0.06257778),\n",
    "        Array(-0.2808735, -0.9403636, -0.19190231),\n",
    "        Array(0.1159052, 0.1652528, -0.97941688)\n",
    "    )\n",
    "\n",
    "val c = ~C // or val c = matrix(C)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Matrix-vector operations are just like in math formula."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a * x + y * 1.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly for matrix-matrix operations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a + b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that `a * b` are element-wise:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a * b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For matrix multiplication, the operator is `%*%`, same as in R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a %*% b %*% c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The method `DenseMatrix.transpose` returns the transpose of matrix, which executes immediately. However, the method `t` is preferred on `MatrixExpression` as it is lazy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a %*% b.t %*% c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Smile has runtime optimization for matrix multiplication chain, which can greatly improve the performance. In the below we generate several random matrices and multiply them together."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "val a = randn( 300,  900)\n",
    "val b = randn( 900,  150)\n",
    "val c = randn( 150, 1800)\n",
    "val d = randn(1800,   30)\n",
    "\n",
    "(a %*% b %*% c %*% d).toMatrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where `randn()` creates a matrix of normally distributed random numbers. Smile tries to load machine optimized BLAS/LAPACK native libraries for most matrix computation. If BLAS/LAPACK is not available, Smile will fall back to pure Java implementation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In linear algebra, a matrix decomposition or matrix factorization is a factorization of a matrix into a product of matrices. There are many different matrix decompositions. In Smile, we provide LU, QR, Cholesky, eigen, and SVD decomposition by functions `lu`, `qr`, `cholesky`, `eigen`, and `svd`, respectively.\n",
    "\n",
    "With these decompositions, many important linear algebra operations can be performed such as calculating matrix rank, determinant, solving linear systems, computing inverse matrix, etc. In fact, Smile has functions `det`, `rank`, `inv` and operator \\ for these common computation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "val a = matrix(Array(\n",
    "  Array(0.7220180, 0.07121225, 0.6881997),\n",
    "  Array(-0.2648886, -0.89044952, 0.3700456),\n",
    "  Array(-0.6391588, 0.44947578, 0.6240573))\n",
    ")\n",
    "\n",
    "val x = Array(1.0, 2.0, 3.0)\n",
    "\n",
    "a \\ x\n",
    "\n",
    "inv(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Scala (2.13)",
   "language": "scala",
   "name": "scala213"
  },
  "language_info": {
   "codemirror_mode": "text/x-scala",
   "file_extension": ".scala",
   "mimetype": "text/x-scala",
   "name": "scala",
   "nbconvert_exporter": "script",
   "version": "2.13.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}