chore: import upstream snapshot with attribution

.github/workflows/python-app.yml

@@ -0,0 +1,26 @@
+name: Tests
+
+on:
+  push:
+    branches: [ master ]
+  pull_request:
+    branches: [ master ]
+
+jobs:
+  build:
+    timeout-minutes: 5
+    runs-on: ubuntu-latest
+    steps:
+    - uses: actions/checkout@v4
+    - name: Set up Python 3.12
+      uses: actions/setup-python@v4
+      with:
+        python-version: 3.12
+    - name: Install dependencies
+      run: |
+        python -m pip install --upgrade pip
+        pip install flake8 pytest
+        pip install -r requirements.txt
+    - name: Test with pytest
+      run: |
+        pytest

.gitignore

@@ -0,0 +1,7 @@
+*.pyc
+build/
+dist/
+mla.egg-info/
+.cache
+*.swp
+.idea

AUTHORS

@@ -0,0 +1,19 @@
+Artem Golubin <me@rushter.com>
+Anebi Agbo
+Convex Path
+James Chevalier
+Jiancheng
+KaiMin Lai
+Nguyễn Tuấn
+Nicolas Hug
+Xiaochun Ma
+Yiran Sheng
+brady salz
+junwang007
+keineahnung2345
+lucaskolstad
+vincent tang
+xq5he
+LanderTome
+therickli
+Andrew Melnik

Dockerfile

@@ -0,0 +1,10 @@
+FROM python:3
+
+RUN mkdir -p /var/app
+WORKDIR /var/app
+COPY . /var/app
+
+# install scipy & numpy
+# install required packages
+RUN pip install scipy numpy && \
+    pip install .

LICENSE

@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2016-2020 Artem Golubin
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

MANIFEST.in

@@ -0,0 +1 @@
+recursive-include mla/datasets/data *

README.md

@@ -0,0 +1,49 @@
+# Machine learning algorithms
+A collection of minimal and clean implementations of machine learning algorithms.
+
+### Why?
+This project is targeting people who want to learn internals of ml algorithms or implement them from scratch.  
+The code is much easier to follow than the optimized libraries and easier to play with.  
+All algorithms are implemented in Python, using numpy, scipy and autograd.  
+
+### Implemented:
+* [Deep learning (MLP, CNN, RNN, LSTM)](mla/neuralnet)
+* [Linear regression, logistic regression](mla/linear_models.py)
+* [Random Forests](mla/ensemble/random_forest.py)
+* [Support vector machine (SVM) with kernels (Linear, Poly, RBF)](mla/svm)
+* [K-Means](mla/kmeans.py)
+* [Gaussian Mixture Model](mla/gaussian_mixture.py)
+* [K-nearest neighbors](mla/knn.py)
+* [Naive bayes](mla/naive_bayes.py)
+* [Principal component analysis (PCA)](mla/pca.py)
+* [Factorization machines](mla/fm.py)
+* [Restricted Boltzmann machine (RBM)](mla/rbm.py)
+* [t-Distributed Stochastic Neighbor Embedding (t-SNE)](mla/tsne.py)
+* [Gradient Boosting trees (also known as GBDT, GBRT, GBM, XGBoost)](mla/ensemble/gbm.py)
+* [Reinforcement learning (Deep Q learning)](mla/rl)
+
+
+### Installation
+```sh
+        git clone https://github.com/rushter/MLAlgorithms
+        cd MLAlgorithms
+        pip install scipy numpy
+        python setup.py develop
+```
+### How to run examples without installation
+```sh
+        cd MLAlgorithms
+        python -m examples.linear_models
+```
+### How to run examples within Docker
+```sh
+        cd MLAlgorithms
+        docker build -t mlalgorithms .
+        docker run --rm -it mlalgorithms bash
+        python -m examples.linear_models
+```
+### Contributing
+
+Your contributions are always welcome!  
+Feel free to improve existing code, documentation or implement new algorithm.  
+Please open an issue to propose your changes if they are big enough.  

README.wehub.md

@@ -0,0 +1,7 @@
+# WeHub 来源说明
+
+- 原始项目：`rushter/MLAlgorithms`
+- 原始仓库：https://github.com/rushter/MLAlgorithms
+- 导入方式：上游默认分支的最新快照
+- 原作者、版权和许可证信息以原始仓库及本仓库 LICENSE 为准
+- 本文件仅用于记录来源，不代表 WeHub 是原项目作者

examples/__init__.py

@@ -0,0 +1 @@
+# coding: utf-8

examples/gaussian_mixture.py

@@ -0,0 +1,46 @@
+import random
+import numpy as np
+import matplotlib.pyplot as plt
+from sklearn import datasets
+from mla.kmeans import KMeans
+from mla.gaussian_mixture import GaussianMixture
+
+random.seed(1)
+np.random.seed(6)
+
+
+def make_clusters(skew=True, *arg, **kwargs):
+    X, y = datasets.make_blobs(*arg, **kwargs)
+    if skew:
+        nrow = X.shape[1]
+        for i in np.unique(y):
+            X[y == i] = X[y == i].dot(np.random.random((nrow, nrow)) - 0.5)
+    return X, y
+
+
+def KMeans_and_GMM(K):
+    COLOR = "bgrcmyk"
+
+    X, y = make_clusters(skew=True, n_samples=1500, centers=K)
+    _, axes = plt.subplots(1, 3)
+
+    # Ground Truth
+    axes[0].scatter(X[:, 0], X[:, 1], c=[COLOR[int(assignment)] for assignment in y])
+    axes[0].set_title("Ground Truth")
+
+    # KMeans
+    kmeans = KMeans(K=K, init="++")
+    kmeans.fit(X)
+    kmeans.predict()
+    axes[1].set_title("KMeans")
+    kmeans.plot(ax=axes[1], holdon=True)
+
+    # Gaussian Mixture
+    gmm = GaussianMixture(K=K, init="kmeans")
+    gmm.fit(X)
+    axes[2].set_title("Gaussian Mixture")
+    gmm.plot(ax=axes[2])
+
+
+if __name__ == "__main__":
+    KMeans_and_GMM(4)

examples/gbm.py

@@ -0,0 +1,70 @@
+import logging
+
+from sklearn.datasets import make_classification
+from sklearn.datasets import make_regression
+from sklearn.metrics import roc_auc_score
+
+try:
+    from sklearn.model_selection import train_test_split
+except ImportError:
+    from sklearn.cross_validation import train_test_split
+
+from mla.ensemble.gbm import GradientBoostingClassifier, GradientBoostingRegressor
+from mla.metrics.metrics import mean_squared_error
+
+logging.basicConfig(level=logging.DEBUG)
+
+
+def classification():
+    # Generate a random binary classification problem.
+    X, y = make_classification(
+        n_samples=350,
+        n_features=15,
+        n_informative=10,
+        random_state=1111,
+        n_classes=2,
+        class_sep=1.0,
+        n_redundant=0,
+    )
+    X_train, X_test, y_train, y_test = train_test_split(
+        X, y, test_size=0.15, random_state=1111
+    )
+
+    model = GradientBoostingClassifier(
+        n_estimators=50, max_depth=4, max_features=8, learning_rate=0.1
+    )
+    model.fit(X_train, y_train)
+    predictions = model.predict(X_test)
+    print(predictions)
+    print(predictions.min())
+    print(predictions.max())
+    print("classification, roc auc score: %s" % roc_auc_score(y_test, predictions))
+
+
+def regression():
+    # Generate a random regression problem
+    X, y = make_regression(
+        n_samples=500,
+        n_features=5,
+        n_informative=5,
+        n_targets=1,
+        noise=0.05,
+        random_state=1111,
+        bias=0.5,
+    )
+    X_train, X_test, y_train, y_test = train_test_split(
+        X, y, test_size=0.1, random_state=1111
+    )
+
+    model = GradientBoostingRegressor(n_estimators=25, max_depth=5, max_features=3)
+    model.fit(X_train, y_train)
+    predictions = model.predict(X_test)
+    print(
+        "regression, mse: %s"
+        % mean_squared_error(y_test.flatten(), predictions.flatten())
+    )
+
+
+if __name__ == "__main__":
+    classification()
+    # regression()

examples/kmeans.py

@@ -0,0 +1,21 @@
+import numpy as np
+from sklearn.datasets import make_blobs
+
+from mla.kmeans import KMeans
+
+
+def kmeans_example(plot=False):
+    X, y = make_blobs(
+        centers=4, n_samples=500, n_features=2, shuffle=True, random_state=42
+    )
+    clusters = len(np.unique(y))
+    k = KMeans(K=clusters, max_iters=150, init="++")
+    k.fit(X)
+    k.predict()
+
+    if plot:
+        k.plot()
+
+
+if __name__ == "__main__":
+    kmeans_example(plot=True)

examples/linear_models.py

@@ -0,0 +1,60 @@
+import logging
+
+try:
+    from sklearn.model_selection import train_test_split
+except ImportError:
+    from sklearn.cross_validation import train_test_split
+from sklearn.datasets import make_classification
+from sklearn.datasets import make_regression
+
+from mla.linear_models import LinearRegression, LogisticRegression
+from mla.metrics.metrics import mean_squared_error, accuracy
+
+# Change to DEBUG to see convergence
+logging.basicConfig(level=logging.ERROR)
+
+
+def regression():
+    # Generate a random regression problem
+    X, y = make_regression(
+        n_samples=10000,
+        n_features=100,
+        n_informative=75,
+        n_targets=1,
+        noise=0.05,
+        random_state=1111,
+        bias=0.5,
+    )
+    X_train, X_test, y_train, y_test = train_test_split(
+        X, y, test_size=0.25, random_state=1111
+    )
+
+    model = LinearRegression(lr=0.01, max_iters=2000, penalty="l2", C=0.03)
+    model.fit(X_train, y_train)
+    predictions = model.predict(X_test)
+    print("regression mse", mean_squared_error(y_test, predictions))
+
+
+def classification():
+    # Generate a random binary classification problem.
+    X, y = make_classification(
+        n_samples=1000,
+        n_features=100,
+        n_informative=75,
+        random_state=1111,
+        n_classes=2,
+        class_sep=2.5,
+    )
+    X_train, X_test, y_train, y_test = train_test_split(
+        X, y, test_size=0.1, random_state=1111
+    )
+
+    model = LogisticRegression(lr=0.01, max_iters=500, penalty="l1", C=0.01)
+    model.fit(X_train, y_train)
+    predictions = model.predict(X_test)
+    print("classification accuracy", accuracy(y_test, predictions))
+
+
+if __name__ == "__main__":
+    regression()
+    classification()

examples/naive_bayes.py

@@ -0,0 +1,31 @@
+from sklearn.datasets import make_classification
+from sklearn.metrics import roc_auc_score
+from sklearn.model_selection import train_test_split
+
+from mla.naive_bayes import NaiveBayesClassifier
+
+
+def classification():
+    # Generate a random binary classification problem.
+    X, y = make_classification(
+        n_samples=1000,
+        n_features=10,
+        n_informative=10,
+        random_state=1111,
+        n_classes=2,
+        class_sep=2.5,
+        n_redundant=0,
+    )
+    X_train, X_test, y_train, y_test = train_test_split(
+        X, y, test_size=0.1, random_state=1111
+    )
+
+    model = NaiveBayesClassifier()
+    model.fit(X_train, y_train)
+    predictions = model.predict(X_test)[:, 1]
+
+    print("classification accuracy", roc_auc_score(y_test, predictions))
+
+
+if __name__ == "__main__":
+    classification()

examples/nearest_neighbors.py

@@ -0,0 +1,59 @@
+try:
+    from sklearn.model_selection import train_test_split
+except ImportError:
+    from sklearn.cross_validation import train_test_split
+from sklearn.datasets import make_classification
+from sklearn.datasets import make_regression
+from scipy.spatial import distance
+
+from mla import knn
+from mla.metrics.metrics import mean_squared_error, accuracy
+
+
+def regression():
+    # Generate a random regression problem
+    X, y = make_regression(
+        n_samples=500,
+        n_features=5,
+        n_informative=5,
+        n_targets=1,
+        noise=0.05,
+        random_state=1111,
+        bias=0.5,
+    )
+    X_train, X_test, y_train, y_test = train_test_split(
+        X, y, test_size=0.25, random_state=1111
+    )
+
+    model = knn.KNNRegressor(k=5, distance_func=distance.euclidean)
+    model.fit(X_train, y_train)
+    predictions = model.predict(X_test)
+    print("regression mse", mean_squared_error(y_test, predictions))
+
+
+def classification():
+    X, y = make_classification(
+        n_samples=500,
+        n_features=5,
+        n_informative=5,
+        n_redundant=0,
+        n_repeated=0,
+        n_classes=3,
+        random_state=1111,
+        class_sep=1.5,
+    )
+
+    X_train, X_test, y_train, y_test = train_test_split(
+        X, y, test_size=0.1, random_state=1111
+    )
+
+    clf = knn.KNNClassifier(k=5, distance_func=distance.euclidean)
+
+    clf.fit(X_train, y_train)
+    predictions = clf.predict(X_test)
+    print("classification accuracy", accuracy(y_test, predictions))
+
+
+if __name__ == "__main__":
+    regression()
+    classification()

examples/nnet_convnet_mnist.py

@@ -0,0 +1,57 @@
+import logging
+
+from mla.datasets import load_mnist
+from mla.metrics import accuracy
+from mla.neuralnet import NeuralNet
+from mla.neuralnet.layers import (
+    Activation,
+    Convolution,
+    MaxPooling,
+    Flatten,
+    Dropout,
+    Parameters,
+)
+from mla.neuralnet.layers import Dense
+from mla.neuralnet.optimizers import Adadelta
+from mla.utils import one_hot
+
+logging.basicConfig(level=logging.DEBUG)
+
+
+# Load MNIST dataset
+X_train, X_test, y_train, y_test = load_mnist()
+
+# Normalize data
+X_train /= 255.0
+X_test /= 255.0
+
+y_train = one_hot(y_train.flatten())
+y_test = one_hot(y_test.flatten())
+print(X_train.shape, X_test.shape, y_train.shape, y_test.shape)
+
+# Approx. 15-20 min. per epoch
+model = NeuralNet(
+    layers=[
+        Convolution(n_filters=32, filter_shape=(3, 3), padding=(1, 1), stride=(1, 1)),
+        Activation("relu"),
+        Convolution(n_filters=32, filter_shape=(3, 3), padding=(1, 1), stride=(1, 1)),
+        Activation("relu"),
+        MaxPooling(pool_shape=(2, 2), stride=(2, 2)),
+        Dropout(0.5),
+        Flatten(),
+        Dense(128),
+        Activation("relu"),
+        Dropout(0.5),
+        Dense(10),
+        Activation("softmax"),
+    ],
+    loss="categorical_crossentropy",
+    optimizer=Adadelta(),
+    metric="accuracy",
+    batch_size=128,
+    max_epochs=3,
+)
+
+model.fit(X_train, y_train)
+predictions = model.predict(X_test)
+print(accuracy(y_test, predictions))

examples/nnet_mlp.py

@@ -0,0 +1,95 @@
+import logging
+
+try:
+    from sklearn.model_selection import train_test_split
+except ImportError:
+    from sklearn.cross_validation import train_test_split
+from sklearn.datasets import make_classification
+from sklearn.datasets import make_regression
+from sklearn.metrics import roc_auc_score
+
+from mla.metrics.metrics import mean_squared_error
+from mla.neuralnet import NeuralNet
+from mla.neuralnet.constraints import MaxNorm
+from mla.neuralnet.layers import Activation, Dense, Dropout
+from mla.neuralnet.optimizers import Adadelta, Adam
+from mla.neuralnet.parameters import Parameters
+from mla.neuralnet.regularizers import L2
+from mla.utils import one_hot
+
+logging.basicConfig(level=logging.DEBUG)
+
+
+def classification():
+    # Generate a random binary classification problem.
+    X, y = make_classification(
+        n_samples=1000,
+        n_features=100,
+        n_informative=75,
+        random_state=1111,
+        n_classes=2,
+        class_sep=2.5,
+    )
+    y = one_hot(y)
+    X_train, X_test, y_train, y_test = train_test_split(
+        X, y, test_size=0.15, random_state=1111
+    )
+
+    model = NeuralNet(
+        layers=[
+            Dense(256, Parameters(init="uniform", regularizers={"W": L2(0.05)})),
+            Activation("relu"),
+            Dropout(0.5),
+            Dense(128, Parameters(init="normal", constraints={"W": MaxNorm()})),
+            Activation("relu"),
+            Dense(2),
+            Activation("softmax"),
+        ],
+        loss="categorical_crossentropy",
+        optimizer=Adadelta(),
+        metric="accuracy",
+        batch_size=64,
+        max_epochs=25,
+    )
+    model.fit(X_train, y_train)
+    predictions = model.predict(X_test)
+    print("classification accuracy", roc_auc_score(y_test[:, 0], predictions[:, 0]))
+
+
+def regression():
+    # Generate a random regression problem
+    X, y = make_regression(
+        n_samples=5000,
+        n_features=25,
+        n_informative=25,
+        n_targets=1,
+        random_state=100,
+        noise=0.05,
+    )
+    y *= 0.01
+    X_train, X_test, y_train, y_test = train_test_split(
+        X, y, test_size=0.1, random_state=1111
+    )
+
+    model = NeuralNet(
+        layers=[
+            Dense(64, Parameters(init="normal")),
+            Activation("linear"),
+            Dense(32, Parameters(init="normal")),
+            Activation("linear"),
+            Dense(1),
+        ],
+        loss="mse",
+        optimizer=Adam(),
+        metric="mse",
+        batch_size=256,
+        max_epochs=15,
+    )
+    model.fit(X_train, y_train)
+    predictions = model.predict(X_test)
+    print("regression mse", mean_squared_error(y_test, predictions.flatten()))
+
+
+if __name__ == "__main__":
+    classification()
+    regression()

examples/nnet_rnn_binary_add.py

@@ -0,0 +1,78 @@
+import logging
+from itertools import combinations, islice
+
+import numpy as np
+
+try:
+    from sklearn.model_selection import train_test_split
+except ImportError:
+    from sklearn.cross_validation import train_test_split
+
+from mla.metrics import accuracy
+from mla.neuralnet import NeuralNet
+from mla.neuralnet.layers import Activation, TimeDistributedDense
+from mla.neuralnet.layers.recurrent import LSTM
+from mla.neuralnet.optimizers import Adam
+
+logging.basicConfig(level=logging.DEBUG)
+
+
+def addition_dataset(dim=10, n_samples=10000, batch_size=64):
+    """Generate binary addition dataset.
+    http://devankuleindiren.com/Projects/rnn_arithmetic.php
+    """
+    binary_format = "{:0" + str(dim) + "b}"
+
+    # Generate all possible number combinations
+    combs = list(islice(combinations(range(2 ** (dim - 1)), 2), n_samples))
+
+    # Initialize empty arrays
+    X = np.zeros((len(combs), dim, 2), dtype=np.uint8)
+    y = np.zeros((len(combs), dim, 1), dtype=np.uint8)
+
+    for i, (a, b) in enumerate(combs):
+        # Convert numbers to binary format
+        X[i, :, 0] = list(reversed([int(x) for x in binary_format.format(a)]))
+        X[i, :, 1] = list(reversed([int(x) for x in binary_format.format(b)]))
+
+        # Generate target variable (a+b)
+        y[i, :, 0] = list(reversed([int(x) for x in binary_format.format(a + b)]))
+
+    X_train, X_test, y_train, y_test = train_test_split(
+        X, y, test_size=0.2, random_state=1111
+    )
+
+    # Round number of examples for batch processing
+    train_b = (X_train.shape[0] // batch_size) * batch_size
+    test_b = (X_test.shape[0] // batch_size) * batch_size
+    X_train = X_train[0:train_b]
+    y_train = y_train[0:train_b]
+
+    X_test = X_test[0:test_b]
+    y_test = y_test[0:test_b]
+    return X_train, X_test, y_train, y_test
+
+
+def addition_problem(ReccurentLayer):
+    X_train, X_test, y_train, y_test = addition_dataset(8, 5000)
+
+    print(X_train.shape, X_test.shape)
+    model = NeuralNet(
+        layers=[ReccurentLayer, TimeDistributedDense(1), Activation("sigmoid")],
+        loss="mse",
+        optimizer=Adam(),
+        metric="mse",
+        batch_size=64,
+        max_epochs=15,
+    )
+    model.fit(X_train, y_train)
+    predictions = np.round(model.predict(X_test))
+    predictions = np.packbits(predictions.astype(np.uint8))
+    y_test = np.packbits(y_test.astype(np.int))
+    print(accuracy(y_test, predictions))
+
+
+# RNN
+# addition_problem(RNN(16, parameters=Parameters(constraints={'W': SmallNorm(), 'U': SmallNorm()})))
+# LSTM
+addition_problem(LSTM(16))

examples/nnet_rnn_text_generation.py

@@ -0,0 +1,82 @@
+from __future__ import print_function
+
+import logging
+import random
+
+import numpy as np
+import sys
+
+from mla.datasets import load_nietzsche
+from mla.neuralnet import NeuralNet
+from mla.neuralnet.constraints import SmallNorm
+from mla.neuralnet.layers import Activation, Dense
+from mla.neuralnet.layers.recurrent import LSTM, RNN
+from mla.neuralnet.optimizers import RMSprop
+
+logging.basicConfig(level=logging.DEBUG)
+
+
+# Example taken from: https://github.com/fchollet/keras/blob/master/examples/lstm_text_generation.py
+
+
+def sample(preds, temperature=1.0):
+    # helper function to sample an index from a probability array
+    preds = np.asarray(preds).astype("float64")
+    preds = np.log(preds) / temperature
+    exp_preds = np.exp(preds)
+    preds = exp_preds / np.sum(exp_preds)
+    probas = np.random.multinomial(1, preds, 1)
+    return np.argmax(probas)
+
+
+X, y, text, chars, char_indices, indices_char = load_nietzsche()
+# Round the number of sequences for batch processing
+items_count = X.shape[0] - (X.shape[0] % 64)
+maxlen = X.shape[1]
+X = X[0:items_count]
+y = y[0:items_count]
+
+print(X.shape, y.shape)
+# LSTM OR RNN
+# rnn_layer = RNN(128, return_sequences=False)
+rnn_layer = LSTM(128, return_sequences=False)
+
+model = NeuralNet(
+    layers=[
+        rnn_layer,
+        # Flatten(),
+        # TimeStepSlicer(-1),
+        Dense(X.shape[2]),
+        Activation("softmax"),
+    ],
+    loss="categorical_crossentropy",
+    optimizer=RMSprop(learning_rate=0.01),
+    metric="accuracy",
+    batch_size=64,
+    max_epochs=1,
+    shuffle=False,
+)
+
+for _ in range(25):
+    model.fit(X, y)
+    start_index = random.randint(0, len(text) - maxlen - 1)
+
+    generated = ""
+    sentence = text[start_index : start_index + maxlen]
+    generated += sentence
+    print('----- Generating with seed: "' + sentence + '"')
+    sys.stdout.write(generated)
+    for i in range(100):
+        x = np.zeros((64, maxlen, len(chars)))
+        for t, char in enumerate(sentence):
+            x[0, t, char_indices[char]] = 1.0
+        preds = model.predict(x)[0]
+        next_index = sample(preds, 0.5)
+        next_char = indices_char[next_index]
+
+        generated += next_char
+        sentence = sentence[1:] + next_char
+
+        sys.stdout.write(next_char)
+        sys.stdout.flush()
+    print()

examples/pca.py

@@ -0,0 +1,39 @@
+try:
+    from sklearn.model_selection import train_test_split
+except ImportError:
+    from sklearn.cross_validation import train_test_split
+from sklearn.datasets import make_classification
+
+from mla.linear_models import LogisticRegression
+from mla.metrics import accuracy
+from mla.pca import PCA
+
+# logging.basicConfig(level=logging.DEBUG)
+
+# Generate a random binary classification problem.
+X, y = make_classification(
+    n_samples=1000,
+    n_features=100,
+    n_informative=75,
+    random_state=1111,
+    n_classes=2,
+    class_sep=2.5,
+)
+
+
+X_train, X_test, y_train, y_test = train_test_split(
+    X, y, test_size=0.25, random_state=1111
+)
+
+for s in ["svd", "eigen"]:
+    p = PCA(15, solver=s)
+
+    # fit PCA with training data, not entire dataset
+    p.fit(X_train)
+    X_train_reduced = p.transform(X_train)
+    X_test_reduced = p.transform(X_test)
+
+    model = LogisticRegression(lr=0.001, max_iters=2500)
+    model.fit(X_train_reduced, y_train)
+    predictions = model.predict(X_test_reduced)
+    print("Classification accuracy for %s PCA: %s" % (s, accuracy(y_test, predictions)))

examples/random_forest.py

@@ -0,0 +1,71 @@
+import logging
+
+import numpy as np
+from sklearn.datasets import make_classification
+from sklearn.datasets import make_regression
+from sklearn.metrics import roc_auc_score, accuracy_score
+
+try:
+    from sklearn.model_selection import train_test_split
+except ImportError:
+    from sklearn.cross_validation import train_test_split
+
+from mla.ensemble.random_forest import RandomForestClassifier, RandomForestRegressor
+from mla.metrics.metrics import mean_squared_error
+
+logging.basicConfig(level=logging.DEBUG)
+
+
+def classification():
+    # Generate a random binary classification problem.
+    X, y = make_classification(
+        n_samples=500,
+        n_features=10,
+        n_informative=10,
+        random_state=1111,
+        n_classes=2,
+        class_sep=2.5,
+        n_redundant=0,
+    )
+
+    X_train, X_test, y_train, y_test = train_test_split(
+        X, y, test_size=0.15, random_state=1111
+    )
+
+    model = RandomForestClassifier(n_estimators=10, max_depth=4)
+    model.fit(X_train, y_train)
+
+    predictions_prob = model.predict(X_test)[:, 1]
+    predictions = np.argmax(model.predict(X_test), axis=1)
+    # print(predictions.shape)
+    print("classification, roc auc score: %s" % roc_auc_score(y_test, predictions_prob))
+    print("classification, accuracy score: %s" % accuracy_score(y_test, predictions))
+
+
+def regression():
+    # Generate a random regression problem
+    X, y = make_regression(
+        n_samples=500,
+        n_features=5,
+        n_informative=5,
+        n_targets=1,
+        noise=0.05,
+        random_state=1111,
+        bias=0.5,
+    )
+    X_train, X_test, y_train, y_test = train_test_split(
+        X, y, test_size=0.1, random_state=1111
+    )
+
+    model = RandomForestRegressor(n_estimators=50, max_depth=10, max_features=3)
+    model.fit(X_train, y_train)
+    predictions = model.predict(X_test)
+    print(
+        "regression, mse: %s"
+        % mean_squared_error(y_test.flatten(), predictions.flatten())
+    )
+
+
+if __name__ == "__main__":
+    classification()
+    # regression()

examples/rbm.py

@@ -0,0 +1,25 @@
+import logging
+
+import numpy as np
+
+from mla.rbm import RBM
+
+logging.basicConfig(level=logging.DEBUG)
+
+
+def print_curve(rbm):
+    from matplotlib import pyplot as plt
+
+    def moving_average(a, n=25):
+        ret = np.cumsum(a, dtype=float)
+        ret[n:] = ret[n:] - ret[:-n]
+        return ret[n - 1 :] / n
+
+    plt.plot(moving_average(rbm.errors))
+    plt.show()
+
+
+X = np.random.uniform(0, 1, (1500, 10))
+rbm = RBM(n_hidden=10, max_epochs=200, batch_size=10, learning_rate=0.1)
+rbm.fit(X)
+print_curve(rbm)

examples/rl_deep_q_learning.py

@@ -0,0 +1,37 @@
+import logging
+
+from mla.neuralnet import NeuralNet
+from mla.neuralnet.layers import Activation, Dense
+from mla.neuralnet.optimizers import Adam
+from mla.rl.dqn import DQN
+
+logging.basicConfig(level=logging.CRITICAL)
+
+
+def mlp_model(n_actions, batch_size=64):
+    model = NeuralNet(
+        layers=[Dense(32), Activation("relu"), Dense(n_actions)],
+        loss="mse",
+        optimizer=Adam(),
+        metric="mse",
+        batch_size=batch_size,
+        max_epochs=1,
+        verbose=False,
+    )
+    return model
+
+
+model = DQN(n_episodes=2500, batch_size=64)
+model.init_environment("CartPole-v0")
+model.init_model(mlp_model)
+
+try:
+    # Train the model
+    # It can take from 300 to 2500 episodes to solve CartPole-v0 problem due to randomness of environment.
+    # You can stop training process using Ctrl+C signal
+    # Read more about this problem: https://gym.openai.com/envs/CartPole-v0
+    model.train(render=False)
+except KeyboardInterrupt:
+    pass
+# Render trained model
+model.play(episodes=100)

examples/svm.py

@@ -0,0 +1,42 @@
+import logging
+
+try:
+    from sklearn.model_selection import train_test_split
+except ImportError:
+    from sklearn.cross_validation import train_test_split
+from sklearn.datasets import make_classification
+
+from mla.metrics.metrics import accuracy
+from mla.svm.kernerls import Linear, RBF
+from mla.svm.svm import SVM
+
+logging.basicConfig(level=logging.DEBUG)
+
+
+def classification():
+    # Generate a random binary classification problem.
+    X, y = make_classification(
+        n_samples=1200,
+        n_features=10,
+        n_informative=5,
+        random_state=1111,
+        n_classes=2,
+        class_sep=1.75,
+    )
+    # Convert y to {-1, 1}
+    y = (y * 2) - 1
+    X_train, X_test, y_train, y_test = train_test_split(
+        X, y, test_size=0.2, random_state=1111
+    )
+
+    for kernel in [RBF(gamma=0.1), Linear()]:
+        model = SVM(max_iter=500, kernel=kernel, C=0.6)
+        model.fit(X_train, y_train)
+        predictions = model.predict(X_test)
+        print(
+            "Classification accuracy (%s): %s" % (kernel, accuracy(y_test, predictions))
+        )
+
+
+if __name__ == "__main__":
+    classification()

examples/t-sne.py

@@ -0,0 +1,28 @@
+import logging
+
+import matplotlib.pyplot as plt
+from sklearn.datasets import make_classification
+
+from mla.tsne import TSNE
+
+logging.basicConfig(level=logging.DEBUG)
+
+X, y = make_classification(
+    n_samples=500,
+    n_features=10,
+    n_informative=5,
+    n_redundant=0,
+    random_state=1111,
+    n_classes=2,
+    class_sep=2.5,
+)
+
+p = TSNE(2, max_iter=500)
+X = p.fit_transform(X)
+
+colors = ["red", "green"]
+for t in range(2):
+    t_mask = (y == t).astype(bool)
+    plt.scatter(X[t_mask, 0], X[t_mask, 1], color=colors[t])
+
+plt.show()

mla/__init__.py

@@ -0,0 +1,3 @@
+# coding:utf-8
+# copyright: (c) 2016 by Artem Golubin
+# license: MIT, see LICENSE for more details.

mla/base/__init__.py

@@ -0,0 +1,2 @@
+# coding:utf-8
+from .base import *

mla/base/base.py

@@ -0,0 +1,63 @@
+# coding:utf-8
+import numpy as np
+
+
+class BaseEstimator:
+    y_required = True
+    fit_required = True
+
+    def _setup_input(self, X, y=None):
+        """Ensure inputs to an estimator are in the expected format.
+
+        Ensures X and y are stored as numpy ndarrays by converting from an
+        array-like object if necessary. Enables estimators to define whether
+        they require a set of y target values or not with y_required, e.g.
+        kmeans clustering requires no target labels and is fit against only X.
+
+        Parameters
+        ----------
+        X : array-like
+            Feature dataset.
+        y : array-like
+            Target values. By default is required, but if y_required = false
+            then may be omitted.
+        """
+        if not isinstance(X, np.ndarray):
+            X = np.array(X)
+
+        if X.size == 0:
+            raise ValueError("Got an empty matrix.")
+
+        if X.ndim == 1:
+            self.n_samples, self.n_features = 1, X.shape
+        else:
+            self.n_samples, self.n_features = X.shape[0], np.prod(X.shape[1:])
+
+        self.X = X
+
+        if self.y_required:
+            if y is None:
+                raise ValueError("Missed required argument y")
+
+            if not isinstance(y, np.ndarray):
+                y = np.array(y)
+
+            if y.size == 0:
+                raise ValueError("The targets array must be no-empty.")
+
+        self.y = y
+
+    def fit(self, X, y=None):
+        self._setup_input(X, y)
+
+    def predict(self, X=None):
+        if not isinstance(X, np.ndarray):
+            X = np.array(X)
+
+        if self.X is not None or not self.fit_required:
+            return self._predict(X)
+        else:
+            raise ValueError("You must call `fit` before `predict`")
+
+    def _predict(self, X=None):
+        raise NotImplementedError()

mla/datasets/__init__.py

@@ -0,0 +1,2 @@
+# coding:utf-8
+from mla.datasets.base import *

mla/datasets/base.py

@@ -0,0 +1,78 @@
+# coding:utf-8
+import os
+
+import numpy as np
+
+
+def get_filename(name):
+    return os.path.join(os.path.dirname(__file__), name)
+
+
+def load_mnist():
+    def load(dataset="training", digits=np.arange(10)):
+        import struct
+        from array import array as pyarray
+        from numpy import array, int8, uint8, zeros
+
+        if dataset == "train":
+            fname_img = get_filename("data/mnist/train-images-idx3-ubyte")
+            fname_lbl = get_filename("data/mnist/train-labels-idx1-ubyte")
+        elif dataset == "test":
+            fname_img = get_filename("data/mnist/t10k-images-idx3-ubyte")
+            fname_lbl = get_filename("data/mnist/t10k-labels-idx1-ubyte")
+        else:
+            raise ValueError("Unexpected dataset name: %r" % dataset)
+
+        flbl = open(fname_lbl, "rb")
+        magic_nr, size = struct.unpack(">II", flbl.read(8))
+        lbl = pyarray("b", flbl.read())
+        flbl.close()
+
+        fimg = open(fname_img, "rb")
+        magic_nr, size, rows, cols = struct.unpack(">IIII", fimg.read(16))
+        img = pyarray("B", fimg.read())
+        fimg.close()
+
+        ind = [k for k in range(size) if lbl[k] in digits]
+        N = len(ind)
+
+        images = zeros((N, rows, cols), dtype=uint8)
+        labels = zeros((N, 1), dtype=int8)
+        for i in range(len(ind)):
+            images[i] = array(
+                img[ind[i] * rows * cols : (ind[i] + 1) * rows * cols]
+            ).reshape((rows, cols))
+            labels[i] = lbl[ind[i]]
+
+        return images, labels
+
+    X_train, y_train = load("train")
+    X_test, y_test = load("test")
+
+    X_train = X_train.reshape(X_train.shape[0], 1, 28, 28).astype(np.float32)
+    X_test = X_test.reshape(X_test.shape[0], 1, 28, 28).astype(np.float32)
+
+    return X_train, X_test, y_train, y_test
+
+
+def load_nietzsche():
+    text = open(get_filename("data/nietzsche.txt"), "rt").read().lower()
+    chars = set(list(text))
+    char_indices = {ch: i for i, ch in enumerate(chars)}
+    indices_char = {i: ch for i, ch in enumerate(chars)}
+
+    maxlen = 40
+    step = 3
+    sentences = []
+    next_chars = []
+    for i in range(0, len(text) - maxlen, step):
+        sentences.append(text[i : i + maxlen])
+        next_chars.append(text[i + maxlen])
+
+    X = np.zeros((len(sentences), maxlen, len(chars)), dtype=np.bool)
+    y = np.zeros((len(sentences), len(chars)), dtype=np.bool)
+    for i, sentence in enumerate(sentences):
+        for t, char in enumerate(sentence):
+            X[i, t, char_indices[char]] = 1
+        y[i, char_indices[next_chars[i]]] = 1
+    return X, y, text, chars, char_indices, indices_char

mla/ensemble/__init__.py

@@ -0,0 +1,2 @@
+# coding:utf-8
+from .random_forest import RandomForestClassifier, RandomForestRegressor

mla/ensemble/base.py

@@ -0,0 +1,65 @@
+# coding:utf-8
+import numpy as np
+from scipy import stats
+
+
+def f_entropy(p):
+    # Convert values to probability
+    p = np.bincount(p) / float(p.shape[0])
+
+    ep = stats.entropy(p)
+    if ep == -float("inf"):
+        return 0.0
+    return ep
+
+
+def information_gain(y, splits):
+    splits_entropy = sum(
+        [f_entropy(split) * (float(split.shape[0]) / y.shape[0]) for split in splits]
+    )
+    return f_entropy(y) - splits_entropy
+
+
+def mse_criterion(y, splits):
+    y_mean = np.mean(y)
+    return -sum(
+        [
+            np.sum((split - y_mean) ** 2) * (float(split.shape[0]) / y.shape[0])
+            for split in splits
+        ]
+    )
+
+
+def xgb_criterion(y, left, right, loss):
+    left = loss.gain(left["actual"], left["y_pred"])
+    right = loss.gain(right["actual"], right["y_pred"])
+    initial = loss.gain(y["actual"], y["y_pred"])
+    gain = left + right - initial
+    return gain
+
+
+def get_split_mask(X, column, value):
+    left_mask = X[:, column] < value
+    right_mask = X[:, column] >= value
+    return left_mask, right_mask
+
+
+def split(X, y, value):
+    left_mask = X < value
+    right_mask = X >= value
+    return y[left_mask], y[right_mask]
+
+
+def split_dataset(X, target, column, value, return_X=True):
+    left_mask, right_mask = get_split_mask(X, column, value)
+
+    left, right = {}, {}
+    for key in target.keys():
+        left[key] = target[key][left_mask]
+        right[key] = target[key][right_mask]
+
+    if return_X:
+        left_X, right_X = X[left_mask], X[right_mask]
+        return left_X, right_X, left, right
+    else:
+        return left, right

mla/ensemble/gbm.py

@@ -0,0 +1,152 @@
+# coding:utf-8
+import numpy as np
+
+# logistic function
+from scipy.special import expit
+
+from mla.base import BaseEstimator
+from mla.ensemble.base import mse_criterion
+from mla.ensemble.tree import Tree
+
+"""
+References:
+https://arxiv.org/pdf/1603.02754v3.pdf
+http://www.saedsayad.com/docs/xgboost.pdf
+https://homes.cs.washington.edu/~tqchen/pdf/BoostedTree.pdf
+http://stats.stackexchange.com/questions/202858/loss-function-approximation-with-taylor-expansion
+"""
+
+
+class Loss:
+    """Base class for loss functions."""
+
+    def __init__(self, regularization=1.0):
+        self.regularization = regularization
+
+    def grad(self, actual, predicted):
+        """First order gradient."""
+        raise NotImplementedError()
+
+    def hess(self, actual, predicted):
+        """Second order gradient."""
+        raise NotImplementedError()
+
+    def approximate(self, actual, predicted):
+        """Approximate leaf value."""
+        return self.grad(actual, predicted).sum() / (
+            self.hess(actual, predicted).sum() + self.regularization
+        )
+
+    def transform(self, pred):
+        """Transform predictions values."""
+        return pred
+
+    def gain(self, actual, predicted):
+        """Calculate gain for split search."""
+        nominator = self.grad(actual, predicted).sum() ** 2
+        denominator = self.hess(actual, predicted).sum() + self.regularization
+        return 0.5 * (nominator / denominator)
+
+
+class LeastSquaresLoss(Loss):
+    """Least squares loss"""
+
+    def grad(self, actual, predicted):
+        return actual - predicted
+
+    def hess(self, actual, predicted):
+        return np.ones_like(actual)
+
+
+class LogisticLoss(Loss):
+    """Logistic loss."""
+
+    def grad(self, actual, predicted):
+        return actual * expit(-actual * predicted)
+
+    def hess(self, actual, predicted):
+        expits = expit(predicted)
+        return expits * (1 - expits)
+
+    def transform(self, output):
+        # Apply logistic (sigmoid) function to the output
+        return expit(output)
+
+
+class GradientBoosting(BaseEstimator):
+    """Gradient boosting trees with Taylor's expansion approximation (as in xgboost)."""
+
+    def __init__(
+        self,
+        n_estimators,
+        learning_rate=0.1,
+        max_features=10,
+        max_depth=2,
+        min_samples_split=10,
+    ):
+        self.min_samples_split = min_samples_split
+        self.learning_rate = learning_rate
+        self.max_depth = max_depth
+        self.max_features = max_features
+        self.n_estimators = n_estimators
+        self.trees = []
+        self.loss = None
+
+    def fit(self, X, y=None):
+        self._setup_input(X, y)
+        self.y_mean = np.mean(y)
+        self._train()
+
+    def _train(self):
+        # Initialize model with zeros
+        y_pred = np.zeros(self.n_samples, np.float32)
+        # Or mean
+        # y_pred = np.full(self.n_samples, self.y_mean)
+
+        for n in range(self.n_estimators):
+            residuals = self.loss.grad(self.y, y_pred)
+            tree = Tree(regression=True, criterion=mse_criterion)
+            # Pass multiple target values to the tree learner
+            targets = {
+                # Residual values
+                "y": residuals,
+                # Actual target values
+                "actual": self.y,
+                # Predictions from previous step
+                "y_pred": y_pred,
+            }
+            tree.train(
+                self.X,
+                targets,
+                max_features=self.max_features,
+                min_samples_split=self.min_samples_split,
+                max_depth=self.max_depth,
+                loss=self.loss,
+            )
+            predictions = tree.predict(self.X)
+            y_pred += self.learning_rate * predictions
+            self.trees.append(tree)
+
+    def _predict(self, X=None):
+        y_pred = np.zeros(X.shape[0], np.float32)
+
+        for i, tree in enumerate(self.trees):
+            y_pred += self.learning_rate * tree.predict(X)
+        return y_pred
+
+    def predict(self, X=None):
+        return self.loss.transform(self._predict(X))
+
+
+class GradientBoostingRegressor(GradientBoosting):
+    def fit(self, X, y=None):
+        self.loss = LeastSquaresLoss()
+        super(GradientBoostingRegressor, self).fit(X, y)
+
+
+class GradientBoostingClassifier(GradientBoosting):
+    def fit(self, X, y=None):
+        # Convert labels from {0, 1} to {-1, 1}
+        y = (y * 2) - 1
+        self.loss = LogisticLoss()
+        super(GradientBoostingClassifier, self).fit(X, y)

mla/ensemble/random_forest.py

@@ -0,0 +1,130 @@
+# coding:utf-8
+import numpy as np
+
+from mla.base import BaseEstimator
+from mla.ensemble.base import information_gain, mse_criterion
+from mla.ensemble.tree import Tree
+
+
+class RandomForest(BaseEstimator):
+    def __init__(
+        self,
+        n_estimators=10,
+        max_features=None,
+        min_samples_split=10,
+        max_depth=None,
+        criterion=None,
+    ):
+        """Base class for RandomForest.
+
+        Parameters
+        ----------
+        n_estimators : int
+            The number of decision tree.
+        max_features : int
+            The number of features to consider when looking for the best split.
+        min_samples_split : int
+            The minimum number of samples required to split an internal node.
+        max_depth : int
+            Maximum depth of the tree.
+        criterion : str
+            The function to measure the quality of a split.
+        """
+        self.max_depth = max_depth
+        self.min_samples_split = min_samples_split
+        self.max_features = max_features
+        self.n_estimators = n_estimators
+        self.trees = []
+
+    def fit(self, X, y):
+        self._setup_input(X, y)
+        if self.max_features is None:
+            self.max_features = int(np.sqrt(X.shape[1]))
+        else:
+            assert X.shape[1] > self.max_features
+        self._train()
+
+    def _train(self):
+        for tree in self.trees:
+            tree.train(
+                self.X,
+                self.y,
+                max_features=self.max_features,
+                min_samples_split=self.min_samples_split,
+                max_depth=self.max_depth,
+            )
+
+    def _predict(self, X=None):
+        raise NotImplementedError()
+
+
+class RandomForestClassifier(RandomForest):
+    def __init__(
+        self,
+        n_estimators=10,
+        max_features=None,
+        min_samples_split=10,
+        max_depth=None,
+        criterion="entropy",
+    ):
+        super(RandomForestClassifier, self).__init__(
+            n_estimators=n_estimators,
+            max_features=max_features,
+            min_samples_split=min_samples_split,
+            max_depth=max_depth,
+            criterion=criterion,
+        )
+
+        if criterion == "entropy":
+            self.criterion = information_gain
+        else:
+            raise ValueError()
+
+        # Initialize empty trees
+        for _ in range(self.n_estimators):
+            self.trees.append(Tree(criterion=self.criterion))
+
+    def _predict(self, X=None):
+        y_shape = np.unique(self.y).shape[0]
+        predictions = np.zeros((X.shape[0], y_shape))
+
+        for i in range(X.shape[0]):
+            row_pred = np.zeros(y_shape)
+            for tree in self.trees:
+                row_pred += tree.predict_row(X[i, :])
+
+            row_pred /= self.n_estimators
+            predictions[i, :] = row_pred
+        return predictions
+
+
+class RandomForestRegressor(RandomForest):
+    def __init__(
+        self,
+        n_estimators=10,
+        max_features=None,
+        min_samples_split=10,
+        max_depth=None,
+        criterion="mse",
+    ):
+        super(RandomForestRegressor, self).__init__(
+            n_estimators=n_estimators,
+            max_features=max_features,
+            min_samples_split=min_samples_split,
+            max_depth=max_depth,
+        )
+
+        if criterion == "mse":
+            self.criterion = mse_criterion
+        else:
+            raise ValueError()
+
+        # Initialize empty regression trees
+        for _ in range(self.n_estimators):
+            self.trees.append(Tree(regression=True, criterion=self.criterion))
+
+    def _predict(self, X=None):
+        predictions = np.zeros((X.shape[0], self.n_estimators))
+        for i, tree in enumerate(self.trees):
+            predictions[:, i] = tree.predict(X)
+        return predictions.mean(axis=1)

mla/ensemble/tree.py

@@ -0,0 +1,206 @@
+# coding:utf-8
+import random
+
+import numpy as np
+from scipy import stats
+
+from mla.ensemble.base import split, split_dataset, xgb_criterion
+
+random.seed(111)
+
+
+class Tree(object):
+    """Recursive implementation of decision tree."""
+
+    def __init__(self, regression=False, criterion=None, n_classes=None):
+        self.regression = regression
+        self.impurity = None
+        self.threshold = None
+        self.column_index = None
+        self.outcome = None
+        self.criterion = criterion
+        self.loss = None
+        self.n_classes = n_classes  # Only for classification
+
+        self.left_child = None
+        self.right_child = None
+
+    @property
+    def is_terminal(self):
+        return not bool(self.left_child and self.right_child)
+
+    def _find_splits(self, X):
+        """Find all possible split values."""
+        split_values = set()
+
+        # Get unique values in a sorted order
+        x_unique = list(np.unique(X))
+        for i in range(1, len(x_unique)):
+            # Find a point between two values
+            average = (x_unique[i - 1] + x_unique[i]) / 2.0
+            split_values.add(average)
+
+        return list(split_values)
+
+    def _find_best_split(self, X, target, n_features):
+        """Find best feature and value for a split. Greedy algorithm."""
+
+        # Sample random subset of features
+        subset = random.sample(list(range(0, X.shape[1])), n_features)
+        max_gain, max_col, max_val = None, None, None
+
+        for column in subset:
+            split_values = self._find_splits(X[:, column])
+            for value in split_values:
+                if self.loss is None:
+                    # Random forest
+                    splits = split(X[:, column], target["y"], value)
+                    gain = self.criterion(target["y"], splits)
+                else:
+                    # Gradient boosting
+                    left, right = split_dataset(
+                        X, target, column, value, return_X=False
+                    )
+                    gain = xgb_criterion(target, left, right, self.loss)
+
+                if (max_gain is None) or (gain > max_gain):
+                    max_col, max_val, max_gain = column, value, gain
+        return max_col, max_val, max_gain
+
+    def _train(
+        self,
+        X,
+        target,
+        max_features=None,
+        min_samples_split=10,
+        max_depth=None,
+        minimum_gain=0.01,
+    ):
+        try:
+            # Exit from recursion using assert syntax
+            assert X.shape[0] > min_samples_split
+            assert max_depth > 0
+
+            if max_features is None:
+                max_features = X.shape[1]
+
+            column, value, gain = self._find_best_split(X, target, max_features)
+            assert gain is not None
+            if self.regression:
+                assert gain != 0
+            else:
+                assert gain > minimum_gain
+
+            self.column_index = column
+            self.threshold = value
+            self.impurity = gain
+
+            # Split dataset
+            left_X, right_X, left_target, right_target = split_dataset(
+                X, target, column, value
+            )
+
+            # Grow left and right child
+            self.left_child = Tree(self.regression, self.criterion, self.n_classes)
+            self.left_child._train(
+                left_X,
+                left_target,
+                max_features,
+                min_samples_split,
+                max_depth - 1,
+                minimum_gain,
+            )
+
+            self.right_child = Tree(self.regression, self.criterion, self.n_classes)
+            self.right_child._train(
+                right_X,
+                right_target,
+                max_features,
+                min_samples_split,
+                max_depth - 1,
+                minimum_gain,
+            )
+        except AssertionError:
+            self._calculate_leaf_value(target)
+
+    def train(
+        self,
+        X,
+        target,
+        max_features=None,
+        min_samples_split=10,
+        max_depth=None,
+        minimum_gain=0.01,
+        loss=None,
+    ):
+        """Build a decision tree from training set.
+
+        Parameters
+        ----------
+
+        X : array-like
+            Feature dataset.
+        target : dictionary or array-like
+            Target values.
+        max_features : int or None
+            The number of features to consider when looking for the best split.
+        min_samples_split : int
+            The minimum number of samples required to split an internal node.
+        max_depth : int
+            Maximum depth of the tree.
+        minimum_gain : float, default 0.01
+            Minimum gain required for splitting.
+        loss : function, default None
+            Loss function for gradient boosting.
+        """
+
+        if not isinstance(target, dict):
+            target = {"y": target}
+
+        # Loss for gradient boosting
+        if loss is not None:
+            self.loss = loss
+
+        if not self.regression:
+            self.n_classes = len(np.unique(target["y"]))
+
+        self._train(
+            X,
+            target,
+            max_features=max_features,
+            min_samples_split=min_samples_split,
+            max_depth=max_depth,
+            minimum_gain=minimum_gain,
+        )
+
+    def _calculate_leaf_value(self, targets):
+        """Find optimal value for leaf."""
+        if self.loss is not None:
+            # Gradient boosting
+            self.outcome = self.loss.approximate(targets["actual"], targets["y_pred"])
+        else:
+            # Random Forest
+            if self.regression:
+                # Mean value for regression task
+                self.outcome = np.mean(targets["y"])
+            else:
+                # Probability for classification task
+                self.outcome = (
+                    np.bincount(targets["y"], minlength=self.n_classes)
+                    / targets["y"].shape[0]
+                )
+
+    def predict_row(self, row):
+        """Predict single row."""
+        if not self.is_terminal:
+            if row[self.column_index] < self.threshold:
+                return self.left_child.predict_row(row)
+            else:
+                return self.right_child.predict_row(row)
+        return self.outcome
+
+    def predict(self, X):
+        result = np.zeros(X.shape[0])
+        for i in range(X.shape[0]):
+            result[i] = self.predict_row(X[i, :])
+        return result

mla/fm.py

@@ -0,0 +1,90 @@
+# coding:utf-8
+
+import autograd.numpy as np
+from autograd import elementwise_grad
+
+from mla.base import BaseEstimator
+from mla.metrics import mean_squared_error, binary_crossentropy
+
+
+np.random.seed(9999)
+
+"""
+References:
+Factorization Machines http://www.csie.ntu.edu.tw/~b97053/paper/Rendle2010FM.pdf
+"""
+
+
+class BaseFM(BaseEstimator):
+    def __init__(
+        self,
+        n_components=10,
+        max_iter=100,
+        init_stdev=0.1,
+        learning_rate=0.01,
+        reg_v=0.1,
+        reg_w=0.5,
+        reg_w0=0.0,
+    ):
+        """Simplified factorization machines implementation using SGD optimizer."""
+        self.reg_w0 = reg_w0
+        self.reg_w = reg_w
+        self.reg_v = reg_v
+        self.n_components = n_components
+        self.lr = learning_rate
+        self.init_stdev = init_stdev
+        self.max_iter = max_iter
+        self.loss = None
+        self.loss_grad = None
+
+    def fit(self, X, y=None):
+        self._setup_input(X, y)
+        # bias
+        self.wo = 0.0
+        # Feature weights
+        self.w = np.zeros(self.n_features)
+        # Factor weights
+        self.v = np.random.normal(
+            scale=self.init_stdev, size=(self.n_features, self.n_components)
+        )
+        self._train()
+
+    def _train(self):
+        for epoch in range(self.max_iter):
+            y_pred = self._predict(self.X)
+            loss = self.loss_grad(self.y, y_pred)
+            w_grad = np.dot(loss, self.X) / float(self.n_samples)
+            self.wo -= self.lr * (loss.mean() + 2 * self.reg_w0 * self.wo)
+            self.w -= self.lr * w_grad + (2 * self.reg_w * self.w)
+            self._factor_step(loss)
+
+    def _factor_step(self, loss):
+        for ix, x in enumerate(self.X):
+            for i in range(self.n_features):
+                v_grad = loss[ix] * (x.dot(self.v).dot(x[i])[0] - self.v[i] * x[i] ** 2)
+                self.v[i] -= self.lr * v_grad + (2 * self.reg_v * self.v[i])
+
+    def _predict(self, X=None):
+        linear_output = np.dot(X, self.w)
+        factors_output = (
+            np.sum(np.dot(X, self.v) ** 2 - np.dot(X**2, self.v**2), axis=1) / 2.0
+        )
+        return self.wo + linear_output + factors_output
+
+
+class FMRegressor(BaseFM):
+    def fit(self, X, y=None):
+        super(FMRegressor, self).fit(X, y)
+        self.loss = mean_squared_error
+        self.loss_grad = elementwise_grad(mean_squared_error)
+
+
+class FMClassifier(BaseFM):
+    def fit(self, X, y=None):
+        super(FMClassifier, self).fit(X, y)
+        self.loss = binary_crossentropy
+        self.loss_grad = elementwise_grad(binary_crossentropy)
+
+    def predict(self, X=None):
+        predictions = self._predict(X)
+        return np.sign(predictions)

mla/gaussian_mixture.py

@@ -0,0 +1,185 @@
+# coding:utf-8
+
+import random
+
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.stats import multivariate_normal
+
+from mla.base import BaseEstimator
+from mla.kmeans import KMeans
+
+
+class GaussianMixture(BaseEstimator):
+    """Gaussian Mixture Model: clusters with Gaussian prior.
+
+    Finds clusters by repeatedly performing Expectation–Maximization (EM) algorithm
+    on the dataset. GMM assumes the datasets is distributed in multivariate Gaussian,
+    and tries to find the underlying structure of the Gaussian, i.e. mean and covariance.
+    E-step computes the "responsibility" of the data to each cluster, given the mean
+    and covariance; M-step computes the mean, covariance and weights (prior of each
+    cluster), given the responsibilities. It iterates until the total likelihood
+    changes less than the tolerance.
+
+
+    Parameters
+    ----------
+
+    K : int
+        The number of clusters into which the dataset is partitioned.
+    max_iters: int
+        The maximum iterations of assigning points to the perform EM.
+        Short-circuited by the assignments converging on their own.
+    init: str, default 'random'
+        The name of the method used to initialize the first clustering.
+
+        'random' - Randomly select values from the dataset as the K centroids.
+        'kmeans' - Initialize the centroids, covariances, weights with KMeams's clusters.
+    tolerance: float, default 1e-3
+        The tolerance of difference of the two latest likelihood for convergence.
+    """
+
+    y_required = False
+
+    def __init__(self, K=4, init="random", max_iters=500, tolerance=1e-3):
+        self.K = K
+        self.max_iters = max_iters
+        self.init = init
+        self.assignments = None
+        self.likelihood = []
+        self.tolerance = tolerance
+
+    def fit(self, X, y=None):
+        """Perform Expectation–Maximization (EM) until converged."""
+        self._setup_input(X, y)
+        self._initialize()
+        for _ in range(self.max_iters):
+            self._E_step()
+            self._M_step()
+            if self._is_converged():
+                break
+
+    def _initialize(self):
+        """Set the initial weights, means and covs (with full covariance matrix).
+
+        weights: the prior of the clusters (what percentage of data does a cluster have)
+        means: the mean points of the clusters
+        covs: the covariance matrix of the clusters
+        """
+        self.weights = np.ones(self.K)
+        if self.init == "random":
+            self.means = [
+                self.X[x] for x in random.sample(range(self.n_samples), self.K)
+            ]
+            self.covs = [np.cov(self.X.T) for _ in range(self.K)]
+
+        elif self.init == "kmeans":
+            kmeans = KMeans(K=self.K, max_iters=self.max_iters // 3, init="++")
+            kmeans.fit(self.X)
+            self.assignments = kmeans.predict()
+            self.means = kmeans.centroids
+            self.covs = []
+            for i in np.unique(self.assignments):
+                self.weights[int(i)] = (self.assignments == i).sum()
+                self.covs.append(np.cov(self.X[self.assignments == i].T))
+        else:
+            raise ValueError("Unknown type of init parameter")
+        self.weights /= self.weights.sum()
+
+    def _E_step(self):
+        """Expectation(E-step) for Gaussian Mixture."""
+        likelihoods = self._get_likelihood(self.X)
+        self.likelihood.append(likelihoods.sum())
+        weighted_likelihoods = self._get_weighted_likelihood(likelihoods)
+        self.assignments = weighted_likelihoods.argmax(axis=1)
+        weighted_likelihoods /= weighted_likelihoods.sum(axis=1)[:, np.newaxis]
+        self.responsibilities = weighted_likelihoods
+
+    def _M_step(self):
+        """Maximization (M-step) for Gaussian Mixture."""
+        weights = self.responsibilities.sum(axis=0)
+        for assignment in range(self.K):
+            resp = self.responsibilities[:, assignment][:, np.newaxis]
+            self.means[assignment] = (resp * self.X).sum(axis=0) / resp.sum()
+            self.covs[assignment] = (self.X - self.means[assignment]).T.dot(
+                (self.X - self.means[assignment]) * resp
+            ) / weights[assignment]
+        self.weights = weights / weights.sum()
+
+    def _is_converged(self):
+        """Check if the difference of the latest two likelihood is less than the tolerance."""
+        if (len(self.likelihood) > 1) and (
+            self.likelihood[-1] - self.likelihood[-2] <= self.tolerance
+        ):
+            return True
+        return False
+
+    def _predict(self, X):
+        """Get the assignments for X with GMM clusters."""
+        if not X.shape:
+            return self.assignments
+        likelihoods = self._get_likelihood(X)
+        weighted_likelihoods = self._get_weighted_likelihood(likelihoods)
+        assignments = weighted_likelihoods.argmax(axis=1)
+        return assignments
+
+    def _get_likelihood(self, data):
+        n_data = data.shape[0]
+        likelihoods = np.zeros([n_data, self.K])
+        for c in range(self.K):
+            likelihoods[:, c] = multivariate_normal.pdf(
+                data, self.means[c], self.covs[c]
+            )
+        return likelihoods
+
+    def _get_weighted_likelihood(self, likelihood):
+        return self.weights * likelihood
+
+    def plot(self, data=None, ax=None, holdon=False):
+        """Plot contour for 2D data."""
+        if not (len(self.X.shape) == 2 and self.X.shape[1] == 2):
+            raise AttributeError("Only support for visualizing 2D data.")
+
+        if ax is None:
+            _, ax = plt.subplots()
+
+        if data is None:
+            data = self.X
+            assignments = self.assignments
+        else:
+            assignments = self.predict(data)
+
+        COLOR = "bgrcmyk"
+        cmap = lambda assignment: COLOR[int(assignment) % len(COLOR)]
+
+        # generate grid
+        delta = 0.025
+        margin = 0.2
+        xmax, ymax = self.X.max(axis=0) + margin
+        xmin, ymin = self.X.min(axis=0) - margin
+        axis_X, axis_Y = np.meshgrid(
+            np.arange(xmin, xmax, delta), np.arange(ymin, ymax, delta)
+        )
+
+        def grid_gaussian_pdf(mean, cov):
+            grid_array = np.array(list(zip(axis_X.flatten(), axis_Y.flatten())))
+            return multivariate_normal.pdf(grid_array, mean, cov).reshape(axis_X.shape)
+
+        # plot scatters
+        if assignments is None:
+            c = None
+        else:
+            c = [cmap(assignment) for assignment in assignments]
+        ax.scatter(data[:, 0], data[:, 1], c=c)
+
+        # plot contours
+        for assignment in range(self.K):
+            ax.contour(
+                axis_X,
+                axis_Y,
+                grid_gaussian_pdf(self.means[assignment], self.covs[assignment]),
+                colors=cmap(assignment),
+            )
+
+        if not holdon:
+            plt.show()

mla/kmeans.py

@@ -0,0 +1,158 @@
+# coding:utf-8
+
+import random
+
+import matplotlib.pyplot as plt
+import numpy as np
+import seaborn as sns
+
+from mla.base import BaseEstimator
+from mla.metrics.distance import euclidean_distance
+
+random.seed(1111)
+
+
+class KMeans(BaseEstimator):
+    """Partition a dataset into K clusters.
+
+    Finds clusters by repeatedly assigning each data point to the cluster with
+    the nearest centroid and iterating until the assignments converge (meaning
+    they don't change during an iteration) or the maximum number of iterations
+    is reached.
+
+    Parameters
+    ----------
+
+    K : int
+        The number of clusters into which the dataset is partitioned.
+    max_iters: int
+        The maximum iterations of assigning points to the nearest cluster.
+        Short-circuited by the assignments converging on their own.
+    init: str, default 'random'
+        The name of the method used to initialize the first clustering.
+
+        'random' - Randomly select values from the dataset as the K centroids.
+        '++' - Select a random first centroid from the dataset, then select
+               K - 1 more centroids by choosing values from the dataset with a
+               probability distribution proportional to the squared distance
+               from each point's closest existing cluster. Attempts to create
+               larger distances between initial clusters to improve convergence
+               rates and avoid degenerate cases.
+    """
+
+    y_required = False
+
+    def __init__(self, K=5, max_iters=100, init="random"):
+        self.K = K
+        self.max_iters = max_iters
+        self.clusters = [[] for _ in range(self.K)]
+        self.centroids = []
+        self.init = init
+
+    def _initialize_centroids(self, init):
+        """Set the initial centroids."""
+
+        if init == "random":
+            self.centroids = [
+                self.X[x] for x in random.sample(range(self.n_samples), self.K)
+            ]
+        elif init == "++":
+            self.centroids = [random.choice(self.X)]
+            while len(self.centroids) < self.K:
+                self.centroids.append(self._choose_next_center())
+        else:
+            raise ValueError("Unknown type of init parameter")
+
+    def _predict(self, X=None):
+        """Perform clustering on the dataset."""
+        self._initialize_centroids(self.init)
+        centroids = self.centroids
+
+        # Optimize clusters
+        for _ in range(self.max_iters):
+            self._assign(centroids)
+            centroids_old = centroids
+            centroids = [self._get_centroid(cluster) for cluster in self.clusters]
+
+            if self._is_converged(centroids_old, centroids):
+                break
+
+        self.centroids = centroids
+
+        return self._get_predictions()
+
+    def _get_predictions(self):
+        predictions = np.empty(self.n_samples)
+
+        for i, cluster in enumerate(self.clusters):
+            for index in cluster:
+                predictions[index] = i
+        return predictions
+
+    def _assign(self, centroids):
+        for row in range(self.n_samples):
+            for i, cluster in enumerate(self.clusters):
+                if row in cluster:
+                    self.clusters[i].remove(row)
+                    break
+
+            closest = self._closest(row, centroids)
+            self.clusters[closest].append(row)
+
+    def _closest(self, fpoint, centroids):
+        """Find the closest centroid for a point."""
+        closest_index = None
+        closest_distance = None
+        for i, point in enumerate(centroids):
+            dist = euclidean_distance(self.X[fpoint], point)
+            if closest_index is None or dist < closest_distance:
+                closest_index = i
+                closest_distance = dist
+        return closest_index
+
+    def _get_centroid(self, cluster):
+        """Get values by indices and take the mean."""
+        return [np.mean(np.take(self.X[:, i], cluster)) for i in range(self.n_features)]
+
+    def _dist_from_centers(self):
+        """Calculate distance from centers."""
+        return np.array(
+            [min([euclidean_distance(x, c) for c in self.centroids]) for x in self.X]
+        )
+
+    def _choose_next_center(self):
+        distances = self._dist_from_centers()
+        squared_distances = distances**2
+        probs = squared_distances / squared_distances.sum()
+        ind = np.random.choice(self.X.shape[0], 1, p=probs)[0]
+        return self.X[ind]
+
+    def _is_converged(self, centroids_old, centroids):
+        """Check if the distance between old and new centroids is zero."""
+        distance = 0
+        for i in range(self.K):
+            distance += euclidean_distance(centroids_old[i], centroids[i])
+        return distance == 0
+
+    def plot(self, ax=None, holdon=False):
+        sns.set(style="white")
+        palette = sns.color_palette("hls", self.K + 1)
+        data = self.X
+
+        if ax is None:
+            _, ax = plt.subplots()
+
+        for i, index in enumerate(self.clusters):
+            point = np.array(data[index]).T
+            ax.scatter(
+                *point,
+                c=[
+                    palette[i],
+                ],
+            )
+
+        for point in self.centroids:
+            ax.scatter(*point, marker="x", linewidths=10)
+
+        if not holdon:
+            plt.show()

mla/knn.py

@@ -0,0 +1,74 @@
+# coding:utf-8
+
+from collections import Counter
+
+import numpy as np
+from scipy.spatial.distance import euclidean
+
+from mla.base import BaseEstimator
+
+
+class KNNBase(BaseEstimator):
+    def __init__(self, k=5, distance_func=euclidean):
+        """Base class for Nearest neighbors classifier and regressor.
+
+        Parameters
+        ----------
+        k : int, default 5
+            The number of neighbors to take into account. If 0, all the
+            training examples are used.
+        distance_func : function, default euclidean distance
+            A distance function taking two arguments. Any function from
+            scipy.spatial.distance will do.
+        """
+
+        self.k = None if k == 0 else k  # l[:None] returns the whole list
+        self.distance_func = distance_func
+
+    def aggregate(self, neighbors_targets):
+        raise NotImplementedError()
+
+    def _predict(self, X=None):
+        predictions = [self._predict_x(x) for x in X]
+
+        return np.array(predictions)
+
+    def _predict_x(self, x):
+        """Predict the label of a single instance x."""
+
+        # compute distances between x and all examples in the training set.
+        distances = (self.distance_func(x, example) for example in self.X)
+
+        # Sort all examples by their distance to x and keep their target value.
+        neighbors = sorted(
+            ((dist, target) for (dist, target) in zip(distances, self.y)),
+            key=lambda x: x[0],
+        )
+
+        # Get targets of the k-nn and aggregate them (most common one or
+        # average).
+        neighbors_targets = [target for (_, target) in neighbors[: self.k]]
+
+        return self.aggregate(neighbors_targets)
+
+
+class KNNClassifier(KNNBase):
+    """Nearest neighbors classifier.
+
+    Note: if there is a tie for the most common label among the neighbors, then
+    the predicted label is arbitrary."""
+
+    def aggregate(self, neighbors_targets):
+        """Return the most common target label."""
+
+        most_common_label = Counter(neighbors_targets).most_common(1)[0][0]
+        return most_common_label
+
+
+class KNNRegressor(KNNBase):
+    """Nearest neighbors regressor."""
+
+    def aggregate(self, neighbors_targets):
+        """Return the mean of all targets."""
+
+        return np.mean(neighbors_targets)

mla/linear_models.py

@@ -0,0 +1,138 @@
+# coding:utf-8
+
+import logging
+
+import autograd.numpy as np
+from autograd import grad
+
+from mla.base import BaseEstimator
+from mla.metrics.metrics import mean_squared_error, binary_crossentropy
+
+np.random.seed(1000)
+
+
+class BasicRegression(BaseEstimator):
+    def __init__(
+        self, lr=0.001, penalty="None", C=0.01, tolerance=0.0001, max_iters=1000
+    ):
+        """Basic class for implementing continuous regression estimators which
+        are trained with gradient descent optimization on their particular loss
+        function.
+
+        Parameters
+        ----------
+        lr : float, default 0.001
+            Learning rate.
+        penalty : str, {'l1', 'l2', None'}, default None
+            Regularization function name.
+        C : float, default 0.01
+            The regularization coefficient.
+        tolerance : float, default 0.0001
+            If the gradient descent updates are smaller than `tolerance`, then
+            stop optimization process.
+        max_iters : int, default 10000
+            The maximum number of iterations.
+        """
+        self.C = C
+        self.penalty = penalty
+        self.tolerance = tolerance
+        self.lr = lr
+        self.max_iters = max_iters
+        self.errors = []
+        self.theta = []
+        self.n_samples, self.n_features = None, None
+        self.cost_func = None
+
+    def _loss(self, w):
+        raise NotImplementedError()
+
+    def init_cost(self):
+        raise NotImplementedError()
+
+    def _add_penalty(self, loss, w):
+        """Apply regularization to the loss."""
+        if self.penalty == "l1":
+            loss += self.C * np.abs(w[1:]).sum()
+        elif self.penalty == "l2":
+            loss += (0.5 * self.C) * (w[1:] ** 2).sum()
+        return loss
+
+    def _cost(self, X, y, theta):
+        prediction = X.dot(theta)
+        error = self.cost_func(y, prediction)
+        return error
+
+    def fit(self, X, y=None):
+        self._setup_input(X, y)
+        self.init_cost()
+        self.n_samples, self.n_features = X.shape
+
+        # Initialize weights + bias term
+        self.theta = np.random.normal(size=(self.n_features + 1), scale=0.5)
+
+        # Add an intercept column
+        self.X = self._add_intercept(self.X)
+
+        self._train()
+
+    @staticmethod
+    def _add_intercept(X):
+        b = np.ones([X.shape[0], 1])
+        return np.concatenate([b, X], axis=1)
+
+    def _train(self):
+        self.theta, self.errors = self._gradient_descent()
+        logging.info(" Theta: %s" % self.theta.flatten())
+
+    def _predict(self, X=None):
+        X = self._add_intercept(X)
+        return X.dot(self.theta)
+
+    def _gradient_descent(self):
+        theta = self.theta
+        errors = [self._cost(self.X, self.y, theta)]
+        # Get derivative of the loss function
+        cost_d = grad(self._loss)
+        for i in range(1, self.max_iters + 1):
+            # Calculate gradient and update theta
+            delta = cost_d(theta)
+            theta -= self.lr * delta
+
+            errors.append(self._cost(self.X, self.y, theta))
+            logging.info("Iteration %s, error %s" % (i, errors[i]))
+
+            error_diff = np.linalg.norm(errors[i - 1] - errors[i])
+            if error_diff < self.tolerance:
+                logging.info("Convergence has reached.")
+                break
+        return theta, errors
+
+
+class LinearRegression(BasicRegression):
+    """Linear regression with gradient descent optimizer."""
+
+    def _loss(self, w):
+        loss = self.cost_func(self.y, np.dot(self.X, w))
+        return self._add_penalty(loss, w)
+
+    def init_cost(self):
+        self.cost_func = mean_squared_error
+
+
+class LogisticRegression(BasicRegression):
+    """Binary logistic regression with gradient descent optimizer."""
+
+    def init_cost(self):
+        self.cost_func = binary_crossentropy
+
+    def _loss(self, w):
+        loss = self.cost_func(self.y, self.sigmoid(np.dot(self.X, w)))
+        return self._add_penalty(loss, w)
+
+    @staticmethod
+    def sigmoid(x):
+        return 0.5 * (np.tanh(0.5 * x) + 1)
+
+    def _predict(self, X=None):
+        X = self._add_intercept(X)
+        return self.sigmoid(X.dot(self.theta))

mla/metrics/__init__.py

@@ -0,0 +1,2 @@
+# coding:utf-8
+from .metrics import *

mla/metrics/base.py

@@ -0,0 +1,25 @@
+# coding:utf-8
+import numpy as np
+
+
+def check_data(a, b):
+    if not isinstance(a, np.ndarray):
+        a = np.array(a)
+
+    if not isinstance(b, np.ndarray):
+        b = np.array(b)
+
+    if type(a) != type(b):
+        raise ValueError("Type mismatch: %s and %s" % (type(a), type(b)))
+
+    if a.size != b.size:
+        raise ValueError("Arrays must be equal in length.")
+    return a, b
+
+
+def validate_input(function):
+    def wrapper(a, b):
+        a, b = check_data(a, b)
+        return function(a, b)
+
+    return wrapper

mla/metrics/distance.py

@@ -0,0 +1,17 @@
+# coding:utf-8
+import math
+
+import numpy as np
+
+
+def euclidean_distance(a, b):
+    if isinstance(a, list) and isinstance(b, list):
+        a = np.array(a)
+        b = np.array(b)
+
+    return math.sqrt(sum((a - b) ** 2))
+
+
+def l2_distance(X):
+    sum_X = np.sum(X * X, axis=1)
+    return (-2 * np.dot(X, X.T) + sum_X).T + sum_X

mla/metrics/metrics.py

@@ -0,0 +1,90 @@
+# coding:utf-8
+import autograd.numpy as np
+
+EPS = 1e-15
+
+
+def unhot(function):
+    """Convert one-hot representation into one column."""
+
+    def wrapper(actual, predicted):
+        if len(actual.shape) > 1 and actual.shape[1] > 1:
+            actual = actual.argmax(axis=1)
+        if len(predicted.shape) > 1 and predicted.shape[1] > 1:
+            predicted = predicted.argmax(axis=1)
+        return function(actual, predicted)
+
+    return wrapper
+
+
+def absolute_error(actual, predicted):
+    return np.abs(actual - predicted)
+
+
+@unhot
+def classification_error(actual, predicted):
+    return (actual != predicted).sum() / float(actual.shape[0])
+
+
+@unhot
+def accuracy(actual, predicted):
+    return 1.0 - classification_error(actual, predicted)
+
+
+def mean_absolute_error(actual, predicted):
+    return np.mean(absolute_error(actual, predicted))
+
+
+def squared_error(actual, predicted):
+    return (actual - predicted) ** 2
+
+
+def squared_log_error(actual, predicted):
+    return (np.log(np.array(actual) + 1) - np.log(np.array(predicted) + 1)) ** 2
+
+
+def mean_squared_log_error(actual, predicted):
+    return np.mean(squared_log_error(actual, predicted))
+
+
+def mean_squared_error(actual, predicted):
+    return np.mean(squared_error(actual, predicted))
+
+
+def root_mean_squared_error(actual, predicted):
+    return np.sqrt(mean_squared_error(actual, predicted))
+
+
+def root_mean_squared_log_error(actual, predicted):
+    return np.sqrt(mean_squared_log_error(actual, predicted))
+
+
+def logloss(actual, predicted):
+    predicted = np.clip(predicted, EPS, 1 - EPS)
+    loss = -np.sum(actual * np.log(predicted))
+    return loss / float(actual.shape[0])
+
+
+def hinge(actual, predicted):
+    return np.mean(np.max(1.0 - actual * predicted, 0.0))
+
+
+def binary_crossentropy(actual, predicted):
+    predicted = np.clip(predicted, EPS, 1 - EPS)
+    return np.mean(
+        -np.sum(actual * np.log(predicted) + (1 - actual) * np.log(1 - predicted))
+    )
+
+
+# aliases
+mse = mean_squared_error
+rmse = root_mean_squared_error
+mae = mean_absolute_error
+
+
+def get_metric(name):
+    """Return metric function by name"""
+    try:
+        return globals()[name]
+    except Exception:
+        raise ValueError("Invalid metric function.")

mla/metrics/tests/__init__.py

@@ -0,0 +1 @@
+# coding:utf-8

mla/metrics/tests/test_metrics.py

@@ -0,0 +1,88 @@
+from __future__ import division
+
+import numpy as np
+import pytest
+from numpy.testing import assert_almost_equal
+
+from mla.metrics.base import check_data, validate_input
+from mla.metrics.metrics import get_metric
+
+
+def test_data_validation():
+    with pytest.raises(ValueError):
+        check_data([], 1)
+
+    with pytest.raises(ValueError):
+        check_data([1, 2, 3], [3, 2])
+
+    a, b = check_data([1, 2, 3], [3, 2, 1])
+
+    assert np.all(a == np.array([1, 2, 3]))
+    assert np.all(b == np.array([3, 2, 1]))
+
+
+def metric(name):
+    return validate_input(get_metric(name))
+
+
+def test_classification_error():
+    f = metric("classification_error")
+    assert f([1, 2, 3, 4], [1, 2, 3, 4]) == 0
+    assert f([1, 2, 3, 4], [1, 2, 3, 5]) == 0.25
+    assert f([1, 1, 1, 0, 0, 0], [1, 1, 1, 1, 0, 0]) == (1.0 / 6)
+
+
+def test_absolute_error():
+    f = metric("absolute_error")
+    assert f([3], [5]) == [2]
+    assert f([-1], [-4]) == [3]
+
+
+def test_mean_absolute_error():
+    f = metric("mean_absolute_error")
+    assert f([1, 2, 3], [1, 2, 3]) == 0
+    assert f([1, 2, 3], [3, 2, 1]) == 4 / 3
+
+
+def test_squared_error():
+    f = metric("squared_error")
+    assert f([1], [1]) == [0]
+    assert f([3], [1]) == [4]
+
+
+def test_squared_log_error():
+    f = metric("squared_log_error")
+    assert f([1], [1]) == [0]
+    assert f([3], [1]) == [np.log(2) ** 2]
+    assert f([np.exp(2) - 1], [np.exp(1) - 1]) == [1.0]
+
+
+def test_mean_squared_log_error():
+    f = metric("mean_squared_log_error")
+    assert f([1, 2, 3], [1, 2, 3]) == 0
+    assert f([1, 2, 3, np.exp(1) - 1], [1, 2, 3, np.exp(2) - 1]) == 0.25
+
+
+def test_root_mean_squared_log_error():
+    f = metric("root_mean_squared_log_error")
+    assert f([1, 2, 3], [1, 2, 3]) == 0
+    assert f([1, 2, 3, np.exp(1) - 1], [1, 2, 3, np.exp(2) - 1]) == 0.5
+
+
+def test_mean_squared_error():
+    f = metric("mean_squared_error")
+    assert f([1, 2, 3], [1, 2, 3]) == 0
+    assert f(range(1, 5), [1, 2, 3, 6]) == 1
+
+
+def test_root_mean_squared_error():
+    f = metric("root_mean_squared_error")
+    assert f([1, 2, 3], [1, 2, 3]) == 0
+    assert f(range(1, 5), [1, 2, 3, 5]) == 0.5
+
+
+def test_multiclass_logloss():
+    f = metric("logloss")
+    assert_almost_equal(f([1], [1]), 0)
+    assert_almost_equal(f([1, 1], [1, 1]), 0)
+    assert_almost_equal(f([1], [0.5]), -np.log(0.5))

mla/naive_bayes.py

@@ -0,0 +1,61 @@
+# coding:utf-8
+
+import numpy as np
+
+from mla.base import BaseEstimator
+from mla.neuralnet.activations import softmax
+
+
+class NaiveBayesClassifier(BaseEstimator):
+    """Gaussian Naive Bayes."""
+
+    # Binary problem.
+    n_classes = 2
+
+    def fit(self, X, y=None):
+        self._setup_input(X, y)
+        # Check target labels
+        assert list(np.unique(y)) == [0, 1]
+
+        # Mean and variance for each class and feature combination
+        self._mean = np.zeros((self.n_classes, self.n_features), dtype=np.float64)
+        self._var = np.zeros((self.n_classes, self.n_features), dtype=np.float64)
+
+        self._priors = np.zeros(self.n_classes, dtype=np.float64)
+
+        for c in range(self.n_classes):
+            # Filter features by class
+            X_c = X[y == c]
+
+            # Calculate mean, variance, prior for each class
+            self._mean[c, :] = X_c.mean(axis=0)
+            self._var[c, :] = X_c.var(axis=0)
+            self._priors[c] = X_c.shape[0] / float(X.shape[0])
+
+    def _predict(self, X=None):
+        # Apply _predict_proba for each row
+        predictions = np.apply_along_axis(self._predict_row, 1, X)
+
+        # Normalize probabilities so that each row will sum up to 1.0
+        return softmax(predictions)
+
+    def _predict_row(self, x):
+        """Predict log likelihood for given row."""
+        output = []
+        for y in range(self.n_classes):
+            prior = np.log(self._priors[y])
+            posterior = np.log(self._pdf(y, x)).sum()
+            prediction = prior + posterior
+
+            output.append(prediction)
+        return output
+
+    def _pdf(self, n_class, x):
+        """Calculate Gaussian PDF for each feature."""
+
+        mean = self._mean[n_class]
+        var = self._var[n_class]
+
+        numerator = np.exp(-((x - mean) ** 2) / (2 * var))
+        denominator = np.sqrt(2 * np.pi * var)
+        return numerator / denominator

mla/neuralnet/__init__.py

@@ -0,0 +1 @@
+from .nnet import NeuralNet

mla/neuralnet/activations.py

@@ -0,0 +1,64 @@
+import autograd.numpy as np
+
+"""
+References:
+https://en.wikipedia.org/wiki/Activation_function
+"""
+
+
+def sigmoid(z):
+    return 1.0 / (1.0 + np.exp(-z))
+
+
+def softmax(z):
+    # Avoid numerical overflow by removing max
+    e = np.exp(z - np.amax(z, axis=1, keepdims=True))
+    return e / np.sum(e, axis=1, keepdims=True)
+
+
+def linear(z):
+    return z
+
+
+def softplus(z):
+    """Smooth relu."""
+    # Avoid numerical overflow, see:
+    # https://docs.scipy.org/doc/numpy/reference/generated/numpy.logaddexp.html
+    return np.logaddexp(0.0, z)
+
+
+def softsign(z):
+    return z / (1 + np.abs(z))
+
+
+def tanh(z):
+    return np.tanh(z)
+
+
+def relu(z):
+    return np.maximum(0, z)
+
+
+def leakyrelu(z, a=0.01):
+    return np.maximum(z * a, z)
+
+
+def gelu(z):
+    """
+    Gaussian Error Linear Unit (GELU)
+    """
+    # mainly used in transformers smoother version of relu
+    
+    return 0.5 * z * (
+        1.0 + np.tanh(
+            np.sqrt(2.0 / np.pi) * (z + 0.044715 * np.power(z, 3))
+        )
+    )
+
+
+def get_activation(name):
+    """Return activation function by name"""
+    try:
+        return globals()[name]
+    except Exception:
+        raise ValueError("Invalid activation function.")

mla/neuralnet/constraints.py

@@ -0,0 +1,40 @@
+# coding:utf-8
+import numpy as np
+
+EPSILON = 10e-8
+
+
+class Constraint(object):
+    def clip(self, p):
+        return p
+
+
+class MaxNorm(object):
+    def __init__(self, m=2, axis=0):
+        self.axis = axis
+        self.m = m
+
+    def clip(self, p):
+        norms = np.sqrt(np.sum(p**2, axis=self.axis))
+        desired = np.clip(norms, 0, self.m)
+        p = p * (desired / (EPSILON + norms))
+        return p
+
+
+class NonNeg(object):
+    def clip(self, p):
+        p[p < 0.0] = 0.0
+        return p
+
+
+class SmallNorm(object):
+    def clip(self, p):
+        return np.clip(p, -5, 5)
+
+
+class UnitNorm(Constraint):
+    def __init__(self, axis=0):
+        self.axis = axis
+
+    def clip(self, p):
+        return p / (EPSILON + np.sqrt(np.sum(p**2, axis=self.axis)))

mla/neuralnet/initializations.py

@@ -0,0 +1,75 @@
+import numpy as np
+
+"""
+References:
+http://jmlr.org/proceedings/papers/v9/glorot10a/glorot10a.pdf
+
+"""
+
+
+def normal(shape, scale=0.5):
+    return np.random.normal(size=shape, scale=scale)
+
+
+def uniform(shape, scale=0.5):
+    return np.random.uniform(size=shape, low=-scale, high=scale)
+
+
+def zero(shape, **kwargs):
+    return np.zeros(shape)
+
+
+def one(shape, **kwargs):
+    return np.ones(shape)
+
+
+def orthogonal(shape, scale=0.5):
+    flat_shape = (shape[0], np.prod(shape[1:]))
+    array = np.random.normal(size=flat_shape)
+    u, _, v = np.linalg.svd(array, full_matrices=False)
+    array = u if u.shape == flat_shape else v
+    return np.reshape(array * scale, shape)
+
+
+def _glorot_fan(shape):
+    assert len(shape) >= 2
+
+    if len(shape) == 4:
+        receptive_field_size = np.prod(shape[2:])
+        fan_in = shape[1] * receptive_field_size
+        fan_out = shape[0] * receptive_field_size
+    else:
+        fan_in, fan_out = shape[:2]
+    return float(fan_in), float(fan_out)
+
+
+def glorot_normal(shape, **kwargs):
+    fan_in, fan_out = _glorot_fan(shape)
+    s = np.sqrt(2.0 / (fan_in + fan_out))
+    return normal(shape, s)
+
+
+def glorot_uniform(shape, **kwargs):
+    fan_in, fan_out = _glorot_fan(shape)
+    s = np.sqrt(6.0 / (fan_in + fan_out))
+    return uniform(shape, s)
+
+
+def he_normal(shape, **kwargs):
+    fan_in, fan_out = _glorot_fan(shape)
+    s = np.sqrt(2.0 / fan_in)
+    return normal(shape, s)
+
+
+def he_uniform(shape, **kwargs):
+    fan_in, fan_out = _glorot_fan(shape)
+    s = np.sqrt(6.0 / fan_in)
+    return uniform(shape, s)
+
+
+def get_initializer(name):
+    """Returns initialization function by the name."""
+    try:
+        return globals()[name]
+    except Exception:
+        raise ValueError("Invalid initialization function.")

mla/neuralnet/layers/__init__.py

@@ -0,0 +1,4 @@
+# coding:utf-8
+from .basic import *
+from .convnet import *
+from .normalization import *

mla/neuralnet/layers/basic.py

@@ -0,0 +1,183 @@
+# coding:utf-8
+import autograd.numpy as np
+from autograd import elementwise_grad
+
+from mla.neuralnet.activations import get_activation
+from mla.neuralnet.parameters import Parameters
+
+np.random.seed(9999)
+
+
+class Layer(object):
+    def setup(self, X_shape):
+        """Allocates initial weights."""
+        pass
+
+    def forward_pass(self, x):
+        raise NotImplementedError()
+
+    def backward_pass(self, delta):
+        raise NotImplementedError()
+
+    def shape(self, x_shape):
+        """Returns shape of the current layer."""
+        raise NotImplementedError()
+
+
+class ParamMixin(object):
+    @property
+    def parameters(self):
+        return self._params
+
+
+class PhaseMixin(object):
+    _train = False
+
+    @property
+    def is_training(self):
+        return self._train
+
+    @is_training.setter
+    def is_training(self, is_train=True):
+        self._train = is_train
+
+    @property
+    def is_testing(self):
+        return not self._train
+
+    @is_testing.setter
+    def is_testing(self, is_test=True):
+        self._train = not is_test
+
+
+class Dense(Layer, ParamMixin):
+    def __init__(self, output_dim, parameters=None):
+        """A fully connected layer.
+
+        Parameters
+        ----------
+        output_dim : int
+        """
+        self._params = parameters
+        self.output_dim = output_dim
+        self.last_input = None
+
+        if parameters is None:
+            self._params = Parameters()
+
+    def setup(self, x_shape):
+        self._params.setup_weights((x_shape[1], self.output_dim))
+
+    def forward_pass(self, X):
+        self.last_input = X
+        return self.weight(X)
+
+    def weight(self, X):
+        W = np.dot(X, self._params["W"])
+        return W + self._params["b"]
+
+    def backward_pass(self, delta):
+        dW = np.dot(self.last_input.T, delta)
+        db = np.sum(delta, axis=0)
+
+        # Update gradient values
+        self._params.update_grad("W", dW)
+        self._params.update_grad("b", db)
+        return np.dot(delta, self._params["W"].T)
+
+    def shape(self, x_shape):
+        return x_shape[0], self.output_dim
+
+
+class Activation(Layer):
+    def __init__(self, name):
+        self.last_input = None
+        self.activation = get_activation(name)
+        # Derivative of activation function
+        self.activation_d = elementwise_grad(self.activation)
+
+    def forward_pass(self, X):
+        self.last_input = X
+        return self.activation(X)
+
+    def backward_pass(self, delta):
+        return self.activation_d(self.last_input) * delta
+
+    def shape(self, x_shape):
+        return x_shape
+
+
+class Dropout(Layer, PhaseMixin):
+    """Randomly set a fraction of `p` inputs to 0 at each training update."""
+
+    def __init__(self, p=0.1):
+        self.p = p
+        self._mask = None
+
+    def forward_pass(self, X):
+        assert self.p > 0
+        if self.is_training:
+            self._mask = np.random.uniform(size=X.shape) > self.p
+            y = X * self._mask
+        else:
+            y = X * (1.0 - self.p)
+
+        return y
+
+    def backward_pass(self, delta):
+        return delta * self._mask
+
+    def shape(self, x_shape):
+        return x_shape
+
+
+class TimeStepSlicer(Layer):
+    """Take a specific time step from 3D tensor."""
+
+    def __init__(self, step=-1):
+        self.step = step
+
+    def forward_pass(self, x):
+        return x[:, self.step, :]
+
+    def backward_pass(self, delta):
+        return np.repeat(delta[:, np.newaxis, :], 2, 1)
+
+    def shape(self, x_shape):
+        return x_shape[0], x_shape[2]
+
+
+class TimeDistributedDense(Layer):
+    """Apply regular Dense layer to every timestep."""
+
+    def __init__(self, output_dim):
+        self.output_dim = output_dim
+        self.n_timesteps = None
+        self.dense = None
+        self.input_dim = None
+
+    def setup(self, X_shape):
+        self.dense = Dense(self.output_dim)
+        self.dense.setup((X_shape[0], X_shape[2]))
+        self.input_dim = X_shape[2]
+
+    def forward_pass(self, X):
+        n_timesteps = X.shape[1]
+        X = X.reshape(-1, X.shape[-1])
+        y = self.dense.forward_pass(X)
+        y = y.reshape((-1, n_timesteps, self.output_dim))
+        return y
+
+    def backward_pass(self, delta):
+        n_timesteps = delta.shape[1]
+        X = delta.reshape(-1, delta.shape[-1])
+        y = self.dense.backward_pass(X)
+        y = y.reshape((-1, n_timesteps, self.input_dim))
+        return y
+
+    @property
+    def parameters(self):
+        return self.dense._params
+
+    def shape(self, x_shape):
+        return x_shape[0], x_shape[1], self.output_dim

mla/neuralnet/layers/convnet.py

@@ -0,0 +1,231 @@
+# coding:utf-8
+import autograd.numpy as np
+
+from mla.neuralnet.layers import Layer, ParamMixin
+from mla.neuralnet.parameters import Parameters
+
+
+class Convolution(Layer, ParamMixin):
+    def __init__(
+        self,
+        n_filters=8,
+        filter_shape=(3, 3),
+        padding=(0, 0),
+        stride=(1, 1),
+        parameters=None,
+    ):
+        """A 2D convolutional layer.
+        Input shape: (n_images, n_channels, height, width)
+
+        Parameters
+        ----------
+        n_filters : int, default 8
+            The number of filters (kernels).
+        filter_shape : tuple(int, int), default (3, 3)
+            The shape of the filters. (height, width)
+        parameters : Parameters instance, default None
+        stride : tuple(int, int), default (1, 1)
+            The step of the convolution. (height, width).
+        padding : tuple(int, int), default (0, 0)
+            The number of pixel to add to each side of the input. (height, weight)
+
+        """
+        self.padding = padding
+        self._params = parameters
+        self.stride = stride
+        self.filter_shape = filter_shape
+        self.n_filters = n_filters
+        if self._params is None:
+            self._params = Parameters()
+
+    def setup(self, X_shape):
+        n_channels, self.height, self.width = X_shape[1:]
+
+        W_shape = (self.n_filters, n_channels) + self.filter_shape
+        b_shape = self.n_filters
+        self._params.setup_weights(W_shape, b_shape)
+
+    def forward_pass(self, X):
+        n_images, n_channels, height, width = self.shape(X.shape)
+        self.last_input = X
+        self.col = image_to_column(X, self.filter_shape, self.stride, self.padding)
+        self.col_W = self._params["W"].reshape(self.n_filters, -1).T
+
+        out = np.dot(self.col, self.col_W) + self._params["b"]
+        out = out.reshape(n_images, height, width, -1).transpose(0, 3, 1, 2)
+        return out
+
+    def backward_pass(self, delta):
+        delta = delta.transpose(0, 2, 3, 1).reshape(-1, self.n_filters)
+
+        d_W = np.dot(self.col.T, delta).transpose(1, 0).reshape(self._params["W"].shape)
+        d_b = np.sum(delta, axis=0)
+        self._params.update_grad("b", d_b)
+        self._params.update_grad("W", d_W)
+
+        d_c = np.dot(delta, self.col_W.T)
+        return column_to_image(
+            d_c, self.last_input.shape, self.filter_shape, self.stride, self.padding
+        )
+
+    def shape(self, x_shape):
+        height, width = convoltuion_shape(
+            self.height, self.width, self.filter_shape, self.stride, self.padding
+        )
+        return x_shape[0], self.n_filters, height, width
+
+
+class MaxPooling(Layer):
+    def __init__(self, pool_shape=(2, 2), stride=(1, 1), padding=(0, 0)):
+        """Max pooling layer.
+        Input shape: (n_images, n_channels, height, width)
+
+        Parameters
+        ----------
+        pool_shape : tuple(int, int), default (2, 2)
+        stride : tuple(int, int), default (1,1)
+        padding : tuple(int, int), default (0,0)
+        """
+        self.pool_shape = pool_shape
+        self.stride = stride
+        self.padding = padding
+
+    def forward_pass(self, X):
+        self.last_input = X
+
+        out_height, out_width = pooling_shape(self.pool_shape, X.shape, self.stride)
+        n_images, n_channels, _, _ = X.shape
+
+        col = image_to_column(X, self.pool_shape, self.stride, self.padding)
+        col = col.reshape(-1, self.pool_shape[0] * self.pool_shape[1])
+
+        arg_max = np.argmax(col, axis=1)
+        out = np.max(col, axis=1)
+        self.arg_max = arg_max
+        return out.reshape(n_images, out_height, out_width, n_channels).transpose(
+            0, 3, 1, 2
+        )
+
+    def backward_pass(self, delta):
+        delta = delta.transpose(0, 2, 3, 1)
+
+        pool_size = self.pool_shape[0] * self.pool_shape[1]
+        y_max = np.zeros((delta.size, pool_size))
+        y_max[np.arange(self.arg_max.size), self.arg_max.flatten()] = delta.flatten()
+        y_max = y_max.reshape(delta.shape + (pool_size,))
+
+        dcol = y_max.reshape(y_max.shape[0] * y_max.shape[1] * y_max.shape[2], -1)
+        return column_to_image(
+            dcol, self.last_input.shape, self.pool_shape, self.stride, self.padding
+        )
+
+    def shape(self, x_shape):
+        h, w = convoltuion_shape(
+            x_shape[2], x_shape[3], self.pool_shape, self.stride, self.padding
+        )
+        return x_shape[0], x_shape[1], h, w
+
+
+class Flatten(Layer):
+    """Flattens multidimensional input into 2D matrix."""
+
+    def forward_pass(self, X):
+        self.last_input_shape = X.shape
+        return X.reshape((X.shape[0], -1))
+
+    def backward_pass(self, delta):
+        return delta.reshape(self.last_input_shape)
+
+    def shape(self, x_shape):
+        return x_shape[0], np.prod(x_shape[1:])
+
+
+def image_to_column(images, filter_shape, stride, padding):
+    """Rearrange image blocks into columns.
+
+    Parameters
+    ----------
+
+    filter_shape : tuple(height, width)
+    images : np.array, shape (n_images, n_channels, height, width)
+    padding: tuple(height, width)
+    stride : tuple (height, width)
+
+    """
+    n_images, n_channels, height, width = images.shape
+    f_height, f_width = filter_shape
+    out_height, out_width = convoltuion_shape(
+        height, width, (f_height, f_width), stride, padding
+    )
+    images = np.pad(images, ((0, 0), (0, 0), padding, padding), mode="constant")
+
+    col = np.zeros((n_images, n_channels, f_height, f_width, out_height, out_width))
+    for y in range(f_height):
+        y_bound = y + stride[0] * out_height
+        for x in range(f_width):
+            x_bound = x + stride[1] * out_width
+            col[:, :, y, x, :, :] = images[
+                :, :, y : y_bound : stride[0], x : x_bound : stride[1]
+            ]
+
+    col = col.transpose(0, 4, 5, 1, 2, 3).reshape(n_images * out_height * out_width, -1)
+    return col
+
+
+def column_to_image(columns, images_shape, filter_shape, stride, padding):
+    """Rearrange columns into image blocks.
+
+    Parameters
+    ----------
+    columns
+    images_shape : tuple(n_images, n_channels, height, width)
+    filter_shape : tuple(height, _width)
+    stride : tuple(height, width)
+    padding : tuple(height, width)
+    """
+    n_images, n_channels, height, width = images_shape
+    f_height, f_width = filter_shape
+
+    out_height, out_width = convoltuion_shape(
+        height, width, (f_height, f_width), stride, padding
+    )
+    columns = columns.reshape(
+        n_images, out_height, out_width, n_channels, f_height, f_width
+    ).transpose(0, 3, 4, 5, 1, 2)
+
+    img_h = height + 2 * padding[0] + stride[0] - 1
+    img_w = width + 2 * padding[1] + stride[1] - 1
+    img = np.zeros((n_images, n_channels, img_h, img_w))
+    for y in range(f_height):
+        y_bound = y + stride[0] * out_height
+        for x in range(f_width):
+            x_bound = x + stride[1] * out_width
+            img[:, :, y : y_bound : stride[0], x : x_bound : stride[1]] += columns[
+                :, :, y, x, :, :
+            ]
+
+    return img[:, :, padding[0] : height + padding[0], padding[1] : width + padding[1]]
+
+
+def convoltuion_shape(img_height, img_width, filter_shape, stride, padding):
+    """Calculate output shape for convolution layer."""
+    height = (img_height + 2 * padding[0] - filter_shape[0]) / float(stride[0]) + 1
+    width = (img_width + 2 * padding[1] - filter_shape[1]) / float(stride[1]) + 1
+
+    assert height % 1 == 0
+    assert width % 1 == 0
+
+    return int(height), int(width)
+
+
+def pooling_shape(pool_shape, image_shape, stride):
+    """Calculate output shape for pooling layer."""
+    n_images, n_channels, height, width = image_shape
+
+    height = (height - pool_shape[0]) / float(stride[0]) + 1
+    width = (width - pool_shape[1]) / float(stride[1]) + 1
+
+    assert height % 1 == 0
+    assert width % 1 == 0
+
+    return int(height), int(width)

mla/neuralnet/layers/normalization.py

@@ -0,0 +1,158 @@
+# coding:utf-8
+import numpy as np
+
+from mla.neuralnet.layers import Layer, PhaseMixin, ParamMixin
+from mla.neuralnet.parameters import Parameters
+
+"""
+References:
+https://kratzert.github.io/2016/02/12/understanding-the-gradient-flow-through-the-batch-normalization-layer.html
+"""
+
+
+class BatchNormalization(Layer, ParamMixin, PhaseMixin):
+    def __init__(self, momentum=0.9, eps=1e-5, parameters=None):
+        super().__init__()
+        self._params = parameters
+        if self._params is None:
+            self._params = Parameters()
+        self.momentum = momentum
+        self.eps = eps
+        self.ema_mean = None
+        self.ema_var = None
+
+    def setup(self, x_shape):
+        self._params.setup_weights((1, x_shape[1]))
+
+    def _forward_pass(self, X):
+        gamma = self._params["W"]
+        beta = self._params["b"]
+
+        if self.is_testing:
+            mu = self.ema_mean
+            xmu = X - mu
+            var = self.ema_var
+            sqrtvar = np.sqrt(var + self.eps)
+            ivar = 1.0 / sqrtvar
+            xhat = xmu * ivar
+            gammax = gamma * xhat
+            return gammax + beta
+
+        N, D = X.shape
+
+        # step1: calculate mean
+        mu = 1.0 / N * np.sum(X, axis=0)
+
+        # step2: subtract mean vector of every trainings example
+        xmu = X - mu
+
+        # step3: following the lower branch - calculation denominator
+        sq = xmu**2
+
+        # step4: calculate variance
+        var = 1.0 / N * np.sum(sq, axis=0)
+
+        # step5: add eps for numerical stability, then sqrt
+        sqrtvar = np.sqrt(var + self.eps)
+
+        # step6: invert sqrtwar
+        ivar = 1.0 / sqrtvar
+
+        # step7: execute normalization
+        xhat = xmu * ivar
+
+        # step8: Nor the two transformation steps
+        gammax = gamma * xhat
+
+        # step9
+        out = gammax + beta
+
+        # store running averages of mean and variance during training for use during testing
+        if self.ema_mean is None or self.ema_var is None:
+            self.ema_mean = mu
+            self.ema_var = var
+        else:
+            self.ema_mean = self.momentum * self.ema_mean + (1 - self.momentum) * mu
+            self.ema_var = self.momentum * self.ema_var + (1 - self.momentum) * var
+        # store intermediate
+        self.cache = (xhat, gamma, xmu, ivar, sqrtvar, var)
+
+        return out
+
+    def forward_pass(self, X):
+        if len(X.shape) == 2:
+            # input is a regular layer
+            return self._forward_pass(X)
+        elif len(X.shape) == 4:
+            # input is a convolution layer
+            N, C, H, W = X.shape
+            x_flat = X.transpose(0, 2, 3, 1).reshape(-1, C)
+            out_flat = self._forward_pass(x_flat)
+            return out_flat.reshape(N, H, W, C).transpose(0, 3, 1, 2)
+        else:
+            raise NotImplementedError(
+                "Unknown model with dimensions = {}".format(len(X.shape))
+            )
+
+    def _backward_pass(self, delta):
+        # unfold the variables stored in cache
+        xhat, gamma, xmu, ivar, sqrtvar, var = self.cache
+
+        # get the dimensions of the input/output
+        N, D = delta.shape
+
+        # step9
+        dbeta = np.sum(delta, axis=0)
+        dgammax = delta  # not necessary, but more understandable
+
+        # step8
+        dgamma = np.sum(dgammax * xhat, axis=0)
+        dxhat = dgammax * gamma
+
+        # step7
+        divar = np.sum(dxhat * xmu, axis=0)
+        dxmu1 = dxhat * ivar
+
+        # step6
+        dsqrtvar = -1.0 / (sqrtvar**2) * divar
+
+        # step5
+        dvar = 0.5 * 1.0 / np.sqrt(var + self.eps) * dsqrtvar
+
+        # step4
+        dsq = 1.0 / N * np.ones((N, D)) * dvar
+
+        # step3
+        dxmu2 = 2 * xmu * dsq
+
+        # step2
+        dx1 = dxmu1 + dxmu2
+        dmu = -1 * np.sum(dxmu1 + dxmu2, axis=0)
+
+        # step1
+        dx2 = 1.0 / N * np.ones((N, D)) * dmu
+
+        # step0
+        dx = dx1 + dx2
+
+        # Update gradient values
+        self._params.update_grad("W", dgamma)
+        self._params.update_grad("b", dbeta)
+
+        return dx
+
+    def backward_pass(self, X):
+        if len(X.shape) == 2:
+            # input is a regular layer
+            return self._backward_pass(X)
+        elif len(X.shape) == 4:
+            # input is a convolution layer
+            N, C, H, W = X.shape
+            x_flat = X.transpose(0, 2, 3, 1).reshape(-1, C)
+            out_flat = self._backward_pass(x_flat)
+            return out_flat.reshape(N, H, W, C).transpose(0, 3, 1, 2)
+        else:
+            raise NotImplementedError("Unknown model shape: {}".format(X.shape))
+
+    def shape(self, x_shape):
+        return x_shape

mla/neuralnet/layers/recurrent/__init__.py

@@ -0,0 +1,3 @@
+# coding:utf-8
+from .lstm import *
+from .rnn import *

mla/neuralnet/layers/recurrent/lstm.py

@@ -0,0 +1,195 @@
+# coding:utf-8
+import autograd.numpy as np
+from autograd import elementwise_grad
+
+from mla.neuralnet.activations import sigmoid
+from mla.neuralnet.initializations import get_initializer
+from mla.neuralnet.layers import Layer, get_activation, ParamMixin
+from mla.neuralnet.parameters import Parameters
+
+"""
+References:
+Understanding LSTM Networks http://colah.github.io/posts/2015-08-Understanding-LSTMs/
+A Critical Review of Recurrent Neural Networks for Sequence Learning http://arxiv.org/pdf/1506.00019v4.pdf
+"""
+
+
+class LSTM(Layer, ParamMixin):
+    def __init__(
+        self,
+        hidden_dim,
+        activation="tanh",
+        inner_init="orthogonal",
+        parameters=None,
+        return_sequences=True,
+    ):
+        self.return_sequences = return_sequences
+        self.hidden_dim = hidden_dim
+        self.inner_init = get_initializer(inner_init)
+        self.activation = get_activation(activation)
+        self.activation_d = elementwise_grad(self.activation)
+        self.sigmoid_d = elementwise_grad(sigmoid)
+
+        if parameters is None:
+            self._params = Parameters()
+        else:
+            self._params = parameters
+
+        self.last_input = None
+        self.states = None
+        self.outputs = None
+        self.gates = None
+        self.hprev = None
+        self.input_dim = None
+        self.W = None
+        self.U = None
+
+    def setup(self, x_shape):
+        """
+        Naming convention:
+        i : input gate
+        f : forget gate
+        c : cell
+        o : output gate
+
+        Parameters
+        ----------
+        x_shape : np.array(batch size, time steps, input shape)
+        """
+        self.input_dim = x_shape[2]
+        # Input -> Hidden
+        W_params = ["W_i", "W_f", "W_o", "W_c"]
+        # Hidden -> Hidden
+        U_params = ["U_i", "U_f", "U_o", "U_c"]
+        # Bias terms
+        b_params = ["b_i", "b_f", "b_o", "b_c"]
+
+        # Initialize params
+        for param in W_params:
+            self._params[param] = self._params.init((self.input_dim, self.hidden_dim))
+
+        for param in U_params:
+            self._params[param] = self.inner_init((self.hidden_dim, self.hidden_dim))
+
+        for param in b_params:
+            self._params[param] = np.full((self.hidden_dim,), self._params.initial_bias)
+
+        # Combine weights for simplicity
+        self.W = [self._params[param] for param in W_params]
+        self.U = [self._params[param] for param in U_params]
+
+        # Init gradient arrays for all weights
+        self._params.init_grad()
+
+        self.hprev = np.zeros((x_shape[0], self.hidden_dim))
+        self.oprev = np.zeros((x_shape[0], self.hidden_dim))
+
+    def forward_pass(self, X):
+        n_samples, n_timesteps, input_shape = X.shape
+        p = self._params
+        self.last_input = X
+
+        self.states = np.zeros((n_samples, n_timesteps + 1, self.hidden_dim))
+        self.outputs = np.zeros((n_samples, n_timesteps + 1, self.hidden_dim))
+        self.gates = {
+            k: np.zeros((n_samples, n_timesteps, self.hidden_dim))
+            for k in ["i", "f", "o", "c"]
+        }
+
+        self.states[:, -1, :] = self.hprev
+        self.outputs[:, -1, :] = self.oprev
+
+        for i in range(n_timesteps):
+            t_gates = np.dot(X[:, i, :], self.W) + np.dot(
+                self.outputs[:, i - 1, :], self.U
+            )
+
+            # Input
+            self.gates["i"][:, i, :] = sigmoid(t_gates[:, 0, :] + p["b_i"])
+            # Forget
+            self.gates["f"][:, i, :] = sigmoid(t_gates[:, 1, :] + p["b_f"])
+            # Output
+            self.gates["o"][:, i, :] = sigmoid(t_gates[:, 2, :] + p["b_o"])
+            # Cell
+            self.gates["c"][:, i, :] = self.activation(t_gates[:, 3, :] + p["b_c"])
+
+            # (previous state * forget) + input + cell
+            self.states[:, i, :] = (
+                self.states[:, i - 1, :] * self.gates["f"][:, i, :]
+                + self.gates["i"][:, i, :] * self.gates["c"][:, i, :]
+            )
+            self.outputs[:, i, :] = self.gates["o"][:, i, :] * self.activation(
+                self.states[:, i, :]
+            )
+
+        self.hprev = self.states[:, n_timesteps - 1, :].copy()
+        self.oprev = self.outputs[:, n_timesteps - 1, :].copy()
+
+        if self.return_sequences:
+            return self.outputs[:, 0:-1, :]
+        else:
+            return self.outputs[:, -2, :]
+
+    def backward_pass(self, delta):
+        if len(delta.shape) == 2:
+            delta = delta[:, np.newaxis, :]
+
+        n_samples, n_timesteps, input_shape = delta.shape
+
+        # Temporal gradient arrays
+        grad = {k: np.zeros_like(self._params[k]) for k in self._params.keys()}
+
+        dh_next = np.zeros((n_samples, input_shape))
+        output = np.zeros((n_samples, n_timesteps, self.input_dim))
+
+        # Backpropagation through time
+        for i in reversed(range(n_timesteps)):
+            dhi = (
+                delta[:, i, :]
+                * self.gates["o"][:, i, :]
+                * self.activation_d(self.states[:, i, :])
+                + dh_next
+            )
+
+            og = delta[:, i, :] * self.activation(self.states[:, i, :])
+            de_o = og * self.sigmoid_d(self.gates["o"][:, i, :])
+
+            grad["W_o"] += np.dot(self.last_input[:, i, :].T, de_o)
+            grad["U_o"] += np.dot(self.outputs[:, i - 1, :].T, de_o)
+            grad["b_o"] += de_o.sum(axis=0)
+
+            de_f = (dhi * self.states[:, i - 1, :]) * self.sigmoid_d(
+                self.gates["f"][:, i, :]
+            )
+            grad["W_f"] += np.dot(self.last_input[:, i, :].T, de_f)
+            grad["U_f"] += np.dot(self.outputs[:, i - 1, :].T, de_f)
+            grad["b_f"] += de_f.sum(axis=0)
+
+            de_i = (dhi * self.gates["c"][:, i, :]) * self.sigmoid_d(
+                self.gates["i"][:, i, :]
+            )
+            grad["W_i"] += np.dot(self.last_input[:, i, :].T, de_i)
+            grad["U_i"] += np.dot(self.outputs[:, i - 1, :].T, de_i)
+            grad["b_i"] += de_i.sum(axis=0)
+
+            de_c = (dhi * self.gates["i"][:, i, :]) * self.activation_d(
+                self.gates["c"][:, i, :]
+            )
+            grad["W_c"] += np.dot(self.last_input[:, i, :].T, de_c)
+            grad["U_c"] += np.dot(self.outputs[:, i - 1, :].T, de_c)
+            grad["b_c"] += de_c.sum(axis=0)
+
+            dh_next = dhi * self.gates["f"][:, i, :]
+
+        # TODO: propagate error to the next layer
+
+        # Change actual gradient arrays
+        for k in grad.keys():
+            self._params.update_grad(k, grad[k])
+        return output
+
+    def shape(self, x_shape):
+        if self.return_sequences:
+            return x_shape[0], x_shape[1], self.hidden_dim
+        else:
+            return x_shape[0], self.hidden_dim

mla/neuralnet/layers/recurrent/rnn.py

@@ -0,0 +1,110 @@
+# coding:utf-8
+import autograd.numpy as np
+from autograd import elementwise_grad
+
+from mla.neuralnet.initializations import get_initializer
+from mla.neuralnet.layers import Layer, get_activation, ParamMixin
+from mla.neuralnet.parameters import Parameters
+
+
+class RNN(Layer, ParamMixin):
+    """Vanilla RNN."""
+
+    def __init__(
+        self,
+        hidden_dim,
+        activation="tanh",
+        inner_init="orthogonal",
+        parameters=None,
+        return_sequences=True,
+    ):
+        self.return_sequences = return_sequences
+        self.hidden_dim = hidden_dim
+        self.inner_init = get_initializer(inner_init)
+        self.activation = get_activation(activation)
+        self.activation_d = elementwise_grad(self.activation)
+        if parameters is None:
+            self._params = Parameters()
+        else:
+            self._params = parameters
+        self.last_input = None
+        self.states = None
+        self.hprev = None
+        self.input_dim = None
+
+    def setup(self, x_shape):
+        """
+        Parameters
+        ----------
+        x_shape : np.array(batch size, time steps, input shape)
+        """
+        self.input_dim = x_shape[2]
+
+        # Input -> Hidden
+        self._params["W"] = self._params.init((self.input_dim, self.hidden_dim))
+        # Bias
+        self._params["b"] = np.full((self.hidden_dim,), self._params.initial_bias)
+        # Hidden -> Hidden layer
+        self._params["U"] = self.inner_init((self.hidden_dim, self.hidden_dim))
+
+        # Init gradient arrays
+        self._params.init_grad()
+
+        self.hprev = np.zeros((x_shape[0], self.hidden_dim))
+
+    def forward_pass(self, X):
+        self.last_input = X
+        n_samples, n_timesteps, input_shape = X.shape
+        states = np.zeros((n_samples, n_timesteps + 1, self.hidden_dim))
+        states[:, -1, :] = self.hprev.copy()
+        p = self._params
+
+        for i in range(n_timesteps):
+            states[:, i, :] = np.tanh(
+                np.dot(X[:, i, :], p["W"])
+                + np.dot(states[:, i - 1, :], p["U"])
+                + p["b"]
+            )
+
+        self.states = states
+        self.hprev = states[:, n_timesteps - 1, :].copy()
+        if self.return_sequences:
+            return states[:, 0:-1, :]
+        else:
+            return states[:, -2, :]
+
+    def backward_pass(self, delta):
+        if len(delta.shape) == 2:
+            delta = delta[:, np.newaxis, :]
+        n_samples, n_timesteps, input_shape = delta.shape
+        p = self._params
+
+        # Temporal gradient arrays
+        grad = {k: np.zeros_like(p[k]) for k in p.keys()}
+
+        dh_next = np.zeros((n_samples, input_shape))
+        output = np.zeros((n_samples, n_timesteps, self.input_dim))
+
+        # Backpropagation through time
+        for i in reversed(range(n_timesteps)):
+            dhi = self.activation_d(self.states[:, i, :]) * (delta[:, i, :] + dh_next)
+
+            grad["W"] += np.dot(self.last_input[:, i, :].T, dhi)
+            grad["b"] += delta[:, i, :].sum(axis=0)
+            grad["U"] += np.dot(self.states[:, i - 1, :].T, dhi)
+
+            dh_next = np.dot(dhi, p["U"].T)
+
+            d = np.dot(delta[:, i, :], p["U"].T)
+            output[:, i, :] = np.dot(d, p["W"].T)
+
+        # Change actual gradient arrays
+        for k in grad.keys():
+            self._params.update_grad(k, grad[k])
+        return output
+
+    def shape(self, x_shape):
+        if self.return_sequences:
+            return x_shape[0], x_shape[1], self.hidden_dim
+        else:
+            return x_shape[0], self.hidden_dim

mla/neuralnet/loss.py

@@ -0,0 +1,11 @@
+from ..metrics import mse, logloss, mae, hinge, binary_crossentropy
+
+categorical_crossentropy = logloss
+
+
+def get_loss(name):
+    """Returns loss function by the name."""
+    try:
+        return globals()[name]
+    except KeyError:
+        raise ValueError("Invalid metric function.")

mla/neuralnet/nnet.py

@@ -0,0 +1,181 @@
+import logging
+
+import numpy as np
+from autograd import elementwise_grad
+
+from mla.base import BaseEstimator
+from mla.metrics.metrics import get_metric
+from mla.neuralnet.layers import PhaseMixin
+from mla.neuralnet.loss import get_loss
+from mla.utils import batch_iterator
+
+np.random.seed(9999)
+
+"""
+Architecture inspired from:
+
+    https://github.com/fchollet/keras
+    https://github.com/andersbll/deeppy
+"""
+
+
+class NeuralNet(BaseEstimator):
+    fit_required = False
+
+    def __init__(
+        self,
+        layers,
+        optimizer,
+        loss,
+        max_epochs=10,
+        batch_size=64,
+        metric="mse",
+        shuffle=False,
+        verbose=True,
+    ):
+        self.verbose = verbose
+        self.shuffle = shuffle
+        self.optimizer = optimizer
+
+        self.loss = get_loss(loss)
+
+        # TODO: fix
+        if loss == "categorical_crossentropy":
+            self.loss_grad = lambda actual, predicted: -(actual - predicted)
+        else:
+            self.loss_grad = elementwise_grad(self.loss, 1)
+        self.metric = get_metric(metric)
+        self.layers = layers
+        self.batch_size = batch_size
+        self.max_epochs = max_epochs
+        self._n_layers = 0
+        self.log_metric = True if loss != metric else False
+        self.metric_name = metric
+        self.bprop_entry = self._find_bprop_entry()
+        self.training = False
+        self._initialized = False
+
+    def _setup_layers(self, x_shape):
+        """Initialize model's layers."""
+        x_shape = list(x_shape)
+        x_shape[0] = self.batch_size
+
+        for layer in self.layers:
+            layer.setup(x_shape)
+            x_shape = layer.shape(x_shape)
+
+        self._n_layers = len(self.layers)
+        # Setup optimizer
+        self.optimizer.setup(self)
+        self._initialized = True
+        logging.info("Total parameters: %s" % self.n_params)
+
+    def _find_bprop_entry(self):
+        """Find entry layer for back propagation."""
+
+        if len(self.layers) > 0 and not hasattr(self.layers[-1], "parameters"):
+            return -1
+        return len(self.layers)
+
+    def fit(self, X, y=None):
+        if not self._initialized:
+            self._setup_layers(X.shape)
+
+        if y.ndim == 1:
+            # Reshape vector to matrix
+            y = y[:, np.newaxis]
+        self._setup_input(X, y)
+
+        self.is_training = True
+        # Pass neural network instance to an optimizer
+        self.optimizer.optimize(self)
+        self.is_training = False
+
+    def update(self, X, y):
+        # Forward pass
+        y_pred = self.fprop(X)
+
+        # Backward pass
+        grad = self.loss_grad(y, y_pred)
+        for layer in reversed(self.layers[: self.bprop_entry]):
+            grad = layer.backward_pass(grad)
+        return self.loss(y, y_pred)
+
+    def fprop(self, X):
+        """Forward propagation."""
+        for layer in self.layers:
+            X = layer.forward_pass(X)
+        return X
+
+    def _predict(self, X=None):
+        if not self._initialized:
+            self._setup_layers(X.shape)
+
+        y = []
+        X_batch = batch_iterator(X, self.batch_size)
+        for Xb in X_batch:
+            y.append(self.fprop(Xb))
+        return np.concatenate(y)
+
+    @property
+    def parametric_layers(self):
+        for layer in self.layers:
+            if hasattr(layer, "parameters"):
+                yield layer
+
+    @property
+    def parameters(self):
+        """Returns a list of all parameters."""
+        params = []
+        for layer in self.parametric_layers:
+            params.append(layer.parameters)
+        return params
+
+    def error(self, X=None, y=None):
+        """Calculate an error for given examples."""
+        training_phase = self.is_training
+        if training_phase:
+            # Temporally disable training.
+            # Some layers work differently while training (e.g. Dropout).
+            self.is_training = False
+        if X is None and y is None:
+            y_pred = self._predict(self.X)
+            score = self.metric(self.y, y_pred)
+        else:
+            y_pred = self._predict(X)
+            score = self.metric(y, y_pred)
+        if training_phase:
+            self.is_training = True
+        return score
+
+    @property
+    def is_training(self):
+        return self.training
+
+    @is_training.setter
+    def is_training(self, train):
+        self.training = train
+        for layer in self.layers:
+            if isinstance(layer, PhaseMixin):
+                layer.is_training = train
+
+    def shuffle_dataset(self):
+        """Shuffle rows in the dataset."""
+        n_samples = self.X.shape[0]
+        indices = np.arange(n_samples)
+        np.random.shuffle(indices)
+        self.X = self.X.take(indices, axis=0)
+        self.y = self.y.take(indices, axis=0)
+
+    @property
+    def n_layers(self):
+        """Returns the number of layers."""
+        return self._n_layers
+
+    @property
+    def n_params(self):
+        """Return the number of trainable parameters."""
+        return sum([layer.parameters.n_params for layer in self.parametric_layers])
+
+    def reset(self):
+        self._initialized = False

mla/neuralnet/optimizers.py

@@ -0,0 +1,251 @@
+import logging
+import time
+from collections import defaultdict
+
+import numpy as np
+from tqdm import tqdm
+
+from mla.utils import batch_iterator
+
+"""
+References:
+
+    Gradient descent optimization algorithms  https://ruder.io/optimizing-gradient-descent/
+"""
+
+
+class Optimizer(object):
+    def optimize(self, network):
+        loss_history = []
+        for i in range(network.max_epochs):
+            if network.shuffle:
+                network.shuffle_dataset()
+
+            start_time = time.time()
+            loss = self.train_epoch(network)
+            loss_history.append(loss)
+            if network.verbose:
+                msg = "Epoch:%s, train loss: %s" % (i, loss)
+                if network.log_metric:
+                    msg += ", train %s: %s" % (network.metric_name, network.error())
+                msg += ", elapsed: %s sec." % (time.time() - start_time)
+                logging.info(msg)
+        return loss_history
+
+    def update(self, network):
+        """Performs an update of parameters."""
+        raise NotImplementedError
+
+    def train_epoch(self, network):
+        losses = []
+
+        # Create batch iterator
+        X_batch = batch_iterator(network.X, network.batch_size)
+        y_batch = batch_iterator(network.y, network.batch_size)
+
+        batch = zip(X_batch, y_batch)
+        if network.verbose:
+            batch = tqdm(
+                batch, total=int(np.ceil(network.n_samples / network.batch_size))
+            )
+
+        for X, y in batch:
+            loss = np.mean(network.update(X, y))
+            self.update(network)
+            losses.append(loss)
+
+        epoch_loss = np.mean(losses)
+        return epoch_loss
+
+    def train_batch(self, network, X, y):
+        loss = np.mean(network.update(X, y))
+        self.update(network)
+        return loss
+
+    def setup(self, network):
+        """Creates additional variables.
+        Note: Must be called before optimization process."""
+        raise NotImplementedError
+
+
+class SGD(Optimizer):
+    def __init__(self, learning_rate=0.01, momentum=0.9, decay=0.0, nesterov=False):
+        self.nesterov = nesterov
+        self.decay = decay
+        self.momentum = momentum
+        self.lr = learning_rate
+        self.iteration = 0
+        self.velocity = None
+
+    def update(self, network):
+        lr = self.lr * (1.0 / (1.0 + self.decay * self.iteration))
+
+        for i, layer in enumerate(network.parametric_layers):
+            for n in layer.parameters.keys():
+                # Get gradient values
+                grad = layer.parameters.grad[n]
+                update = self.momentum * self.velocity[i][n] - lr * grad
+                self.velocity[i][n] = update
+                if self.nesterov:
+                    # Adjust using updated velocity
+                    update = self.momentum * self.velocity[i][n] - lr * grad
+                layer.parameters.step(n, update)
+        self.iteration += 1
+
+    def setup(self, network):
+        self.velocity = defaultdict(dict)
+        for i, layer in enumerate(network.parametric_layers):
+            for n in layer.parameters.keys():
+                self.velocity[i][n] = np.zeros_like(layer.parameters[n])
+
+
+class Adagrad(Optimizer):
+    def __init__(self, learning_rate=0.01, epsilon=1e-8):
+        self.eps = epsilon
+        self.lr = learning_rate
+
+    def update(self, network):
+        for i, layer in enumerate(network.parametric_layers):
+            for n in layer.parameters.keys():
+                grad = layer.parameters.grad[n]
+                self.accu[i][n] += grad**2
+                step = self.lr * grad / (np.sqrt(self.accu[i][n]) + self.eps)
+                layer.parameters.step(n, -step)
+
+    def setup(self, network):
+        # Accumulators
+        self.accu = defaultdict(dict)
+        for i, layer in enumerate(network.parametric_layers):
+            for n in layer.parameters.keys():
+                self.accu[i][n] = np.zeros_like(layer.parameters[n])
+
+
+class Adadelta(Optimizer):
+    def __init__(self, learning_rate=1.0, rho=0.95, epsilon=1e-8):
+        self.rho = rho
+        self.eps = epsilon
+        self.lr = learning_rate
+
+    def update(self, network):
+        for i, layer in enumerate(network.parametric_layers):
+            for n in layer.parameters.keys():
+                grad = layer.parameters.grad[n]
+                self.accu[i][n] = (
+                    self.rho * self.accu[i][n] + (1.0 - self.rho) * grad**2
+                )
+                step = (
+                    grad
+                    * np.sqrt(self.d_accu[i][n] + self.eps)
+                    / np.sqrt(self.accu[i][n] + self.eps)
+                )
+
+                layer.parameters.step(n, -step * self.lr)
+                # Update delta accumulator
+                self.d_accu[i][n] = (
+                    self.rho * self.d_accu[i][n] + (1.0 - self.rho) * step**2
+                )
+
+    def setup(self, network):
+        # Accumulators
+        self.accu = defaultdict(dict)
+        self.d_accu = defaultdict(dict)
+        for i, layer in enumerate(network.parametric_layers):
+            for n in layer.parameters.keys():
+                self.accu[i][n] = np.zeros_like(layer.parameters[n])
+                self.d_accu[i][n] = np.zeros_like(layer.parameters[n])
+
+
+class RMSprop(Optimizer):
+    def __init__(self, learning_rate=0.001, rho=0.9, epsilon=1e-8):
+        self.eps = epsilon
+        self.rho = rho
+        self.lr = learning_rate
+
+    def update(self, network):
+        for i, layer in enumerate(network.parametric_layers):
+            for n in layer.parameters.keys():
+                grad = layer.parameters.grad[n]
+                self.accu[i][n] = (self.rho * self.accu[i][n]) + (1.0 - self.rho) * (
+                    grad**2
+                )
+                step = self.lr * grad / (np.sqrt(self.accu[i][n]) + self.eps)
+                layer.parameters.step(n, -step)
+
+    def setup(self, network):
+        # Accumulators
+        self.accu = defaultdict(dict)
+        for i, layer in enumerate(network.parametric_layers):
+            for n in layer.parameters.keys():
+                self.accu[i][n] = np.zeros_like(layer.parameters[n])
+
+
+class Adam(Optimizer):
+    def __init__(self, learning_rate=0.001, beta_1=0.9, beta_2=0.999, epsilon=1e-8):
+        self.epsilon = epsilon
+        self.beta_2 = beta_2
+        self.beta_1 = beta_1
+        self.lr = learning_rate
+        self.iterations = 0
+        self.t = 1
+
+    def update(self, network):
+        for i, layer in enumerate(network.parametric_layers):
+            for n in layer.parameters.keys():
+                grad = layer.parameters.grad[n]
+                self.ms[i][n] = (self.beta_1 * self.ms[i][n]) + (
+                    1.0 - self.beta_1
+                ) * grad
+                self.vs[i][n] = (self.beta_2 * self.vs[i][n]) + (
+                    1.0 - self.beta_2
+                ) * grad**2
+                lr = (
+                    self.lr
+                    * np.sqrt(1.0 - self.beta_2**self.t)
+                    / (1.0 - self.beta_1**self.t)
+                )
+
+                step = lr * self.ms[i][n] / (np.sqrt(self.vs[i][n]) + self.epsilon)
+                layer.parameters.step(n, -step)
+        self.t += 1
+
+    def setup(self, network):
+        # Accumulators
+        self.ms = defaultdict(dict)
+        self.vs = defaultdict(dict)
+        for i, layer in enumerate(network.parametric_layers):
+            for n in layer.parameters.keys():
+                self.ms[i][n] = np.zeros_like(layer.parameters[n])
+                self.vs[i][n] = np.zeros_like(layer.parameters[n])
+
+
+class Adamax(Optimizer):
+    def __init__(self, learning_rate=0.002, beta_1=0.9, beta_2=0.999, epsilon=1e-8):
+        self.epsilon = epsilon
+        self.beta_2 = beta_2
+        self.beta_1 = beta_1
+        self.lr = learning_rate
+        self.t = 1
+
+    def update(self, network):
+        for i, layer in enumerate(network.parametric_layers):
+            for n in layer.parameters.keys():
+                grad = layer.parameters.grad[n]
+                self.ms[i][n] = self.beta_1 * self.ms[i][n] + (1.0 - self.beta_1) * grad
+                self.us[i][n] = np.maximum(self.beta_2 * self.us[i][n], np.abs(grad))
+
+                step = (
+                    self.lr
+                    / (1 - self.beta_1**self.t)
+                    * self.ms[i][n]
+                    / (self.us[i][n] + self.epsilon)
+                )
+                layer.parameters.step(n, -step)
+        self.t += 1
+
+    def setup(self, network):
+        self.ms = defaultdict(dict)
+        self.us = defaultdict(dict)
+        for i, layer in enumerate(network.parametric_layers):
+            for n in layer.parameters.keys():
+                self.ms[i][n] = np.zeros_like(layer.parameters[n])
+                self.us[i][n] = np.zeros_like(layer.parameters[n])

mla/neuralnet/parameters.py

@@ -0,0 +1,98 @@
+# coding:utf-8
+import numpy as np
+
+from mla.neuralnet.initializations import get_initializer
+
+
+class Parameters(object):
+    def __init__(
+        self,
+        init="glorot_uniform",
+        scale=0.5,
+        bias=1.0,
+        regularizers=None,
+        constraints=None,
+    ):
+        """A container for layer's parameters.
+
+        Parameters
+        ----------
+        init : str, default 'glorot_uniform'.
+            The name of the weight initialization function.
+        scale : float, default 0.5
+        bias : float, default 1.0
+            Initial values for bias.
+        regularizers : dict
+            Weight regularizers.
+            >>> {'W' : L2()}
+        constraints : dict
+            Weight constraints.
+            >>> {'b' : MaxNorm()}
+        """
+        if constraints is None:
+            self.constraints = {}
+        else:
+            self.constraints = constraints
+
+        if regularizers is None:
+            self.regularizers = {}
+        else:
+            self.regularizers = regularizers
+
+        self.initial_bias = bias
+        self.scale = scale
+        self.init = get_initializer(init)
+
+        self._params = {}
+        self._grads = {}
+
+    def setup_weights(self, W_shape, b_shape=None):
+        if "W" not in self._params:
+            self._params["W"] = self.init(shape=W_shape, scale=self.scale)
+            if b_shape is None:
+                self._params["b"] = np.full(W_shape[1], self.initial_bias)
+            else:
+                self._params["b"] = np.full(b_shape, self.initial_bias)
+        self.init_grad()
+
+    def init_grad(self):
+        """Init gradient arrays corresponding to each weight array."""
+        for key in self._params.keys():
+            if key not in self._grads:
+                self._grads[key] = np.zeros_like(self._params[key])
+
+    def step(self, name, step):
+        """Increase specific weight by amount of the step parameter."""
+        self._params[name] += step
+
+        if name in self.constraints:
+            self._params[name] = self.constraints[name].clip(self._params[name])
+
+    def update_grad(self, name, value):
+        """Update gradient values."""
+        self._grads[name] = value
+
+        if name in self.regularizers:
+            self._grads[name] += self.regularizers[name](self._params[name])
+
+    @property
+    def n_params(self):
+        """Count the number of parameters in this layer."""
+        return sum([np.prod(self._params[x].shape) for x in self._params.keys()])
+
+    def keys(self):
+        return self._params.keys()
+
+    @property
+    def grad(self):
+        return self._grads
+
+    # Allow access to the fields using dict syntax, e.g. parameters['W']
+    def __getitem__(self, item):
+        if item in self._params:
+            return self._params[item]
+        else:
+            raise ValueError
+
+    def __setitem__(self, key, value):
+        self._params[key] = value

mla/neuralnet/regularizers.py

@@ -0,0 +1,35 @@
+# coding:utf-8
+import numpy as np
+from autograd import elementwise_grad
+
+
+class Regularizer(object):
+    def __init__(self, C=0.01):
+        self.C = C
+        self._grad = elementwise_grad(self._penalty)
+
+    def _penalty(self, weights):
+        raise NotImplementedError()
+
+    def grad(self, weights):
+        return self._grad(weights)
+
+    def __call__(self, weights):
+        return self.grad(weights)
+
+
+class L1(Regularizer):
+    def _penalty(self, weights):
+        return self.C * np.abs(weights)
+
+
+class L2(Regularizer):
+    def _penalty(self, weights):
+        return self.C * weights**2
+
+
+class ElasticNet(Regularizer):
+    """Linear combination of L1 and L2 penalties."""
+
+    def _penalty(self, weights):
+        return 0.5 * self.C * weights**2 + (1.0 - self.C) * np.abs(weights)

mla/neuralnet/tests/test_activations.py

@@ -0,0 +1,21 @@
+import sys
+
+import numpy as np
+
+from mla.neuralnet.activations import *
+
+
+def test_softplus():
+    # np.exp(z_max) will overflow
+    z_max = np.log(sys.float_info.max) + 1.0e10
+    # 1.0 / np.exp(z_min) will overflow
+    z_min = np.log(sys.float_info.min) - 1.0e10
+    inputs = np.array([0.0, 1.0, -1.0, z_min, z_max])
+    # naive implementation of np.log(1 + np.exp(z_max)) will overflow
+    # naive implementation of z + np.log(1 + 1 / np.exp(z_min)) will
+    # throw ZeroDivisionError
+    outputs = np.array(
+        [np.log(2.0), np.log1p(np.exp(1.0)), np.log1p(np.exp(-1.0)), 0.0, z_max]
+    )
+
+    assert np.allclose(outputs, softplus(inputs))

mla/neuralnet/tests/test_optimizers.py

@@ -0,0 +1,72 @@
+from sklearn.datasets import make_classification
+from sklearn.metrics import roc_auc_score
+from sklearn.model_selection import train_test_split
+
+from mla.neuralnet import NeuralNet
+from mla.neuralnet.layers import Dense, Activation, Dropout, Parameters
+from mla.neuralnet.optimizers import *
+from mla.utils import one_hot
+
+
+def clasifier(optimizer):
+    X, y = make_classification(
+        n_samples=1000,
+        n_features=100,
+        n_informative=75,
+        random_state=1111,
+        n_classes=2,
+        class_sep=2.5,
+    )
+    y = one_hot(y)
+
+    X -= np.mean(X, axis=0)
+    X /= np.std(X, axis=0)
+    X_train, X_test, y_train, y_test = train_test_split(
+        X, y, test_size=0.15, random_state=1111
+    )
+
+    model = NeuralNet(
+        layers=[
+            Dense(128, Parameters(init="uniform")),
+            Activation("relu"),
+            Dropout(0.5),
+            Dense(64, Parameters(init="normal")),
+            Activation("relu"),
+            Dense(2),
+            Activation("softmax"),
+        ],
+        loss="categorical_crossentropy",
+        optimizer=optimizer,
+        metric="accuracy",
+        batch_size=64,
+        max_epochs=10,
+    )
+    model.fit(X_train, y_train)
+    predictions = model.predict(X_test)
+    return roc_auc_score(y_test[:, 0], predictions[:, 0])
+
+
+def test_adadelta():
+    assert clasifier(Adadelta()) > 0.9
+
+
+def test_adam():
+    assert clasifier(Adam()) > 0.9
+
+
+def test_adamax():
+    assert clasifier(Adamax()) > 0.9
+
+
+def test_rmsprop():
+    assert clasifier(RMSprop()) > 0.9
+
+
+def test_adagrad():
+    assert clasifier(Adagrad()) > 0.9
+
+
+def test_sgd():
+    assert clasifier(SGD(learning_rate=0.0001)) > 0.9
+    assert clasifier(SGD(learning_rate=0.0001, nesterov=True, momentum=0.9)) > 0.9
+    assert clasifier(SGD(learning_rate=0.0001, nesterov=False, momentum=0.0)) > 0.9

mla/pca.py

@@ -0,0 +1,63 @@
+# coding:utf-8
+import logging
+
+import numpy as np
+from scipy.linalg import svd
+
+from mla.base import BaseEstimator
+
+np.random.seed(1000)
+
+
+class PCA(BaseEstimator):
+    y_required = False
+
+    def __init__(self, n_components, solver="svd"):
+        """Principal component analysis (PCA) implementation.
+
+        Transforms a dataset of possibly correlated values into n linearly
+        uncorrelated components. The components are ordered such that the first
+        has the largest possible variance and each following component as the
+        largest possible variance given the previous components. This causes
+        the early components to contain most of the variability in the dataset.
+
+        Parameters
+        ----------
+        n_components : int
+        solver : str, default 'svd'
+            {'svd', 'eigen'}
+        """
+        self.solver = solver
+        self.n_components = n_components
+        self.components = None
+        self.mean = None
+
+    def fit(self, X, y=None):
+        self.mean = np.mean(X, axis=0)
+        self._decompose(X)
+
+    def _decompose(self, X):
+        # Mean centering
+        X = X.copy()
+        X -= self.mean
+
+        if self.solver == "svd":
+            _, s, Vh = svd(X, full_matrices=True)
+        elif self.solver == "eigen":
+            s, Vh = np.linalg.eig(np.cov(X.T))
+            Vh = Vh.T
+
+        s_squared = s**2
+        variance_ratio = s_squared / s_squared.sum()
+        logging.info(
+            "Explained variance ratio: %s" % (variance_ratio[0 : self.n_components])
+        )
+        self.components = Vh[0 : self.n_components]
+
+    def transform(self, X):
+        X = X.copy()
+        X -= self.mean
+        return np.dot(X, self.components.T)
+
+    def _predict(self, X=None):
+        return self.transform(X)

mla/rbm.py

@@ -0,0 +1,101 @@
+# coding:utf-8
+import logging
+
+import numpy as np
+from scipy.special import expit
+
+from mla.base import BaseEstimator
+from mla.utils import batch_iterator
+
+np.random.seed(9999)
+sigmoid = expit
+
+"""
+References:
+A Practical Guide to Training Restricted Boltzmann Machines https://www.cs.toronto.edu/~hinton/absps/guideTR.pdf
+"""
+
+
+class RBM(BaseEstimator):
+    y_required = False
+
+    def __init__(self, n_hidden=128, learning_rate=0.1, batch_size=10, max_epochs=100):
+        """Bernoulli Restricted Boltzmann Machine (RBM)
+
+        Parameters
+        ----------
+
+        n_hidden : int, default 128
+            The number of hidden units.
+        learning_rate : float, default 0.1
+        batch_size : int, default 10
+        max_epochs : int, default 100
+        """
+        self.max_epochs = max_epochs
+        self.batch_size = batch_size
+        self.lr = learning_rate
+        self.n_hidden = n_hidden
+
+    def fit(self, X, y=None):
+        self.n_visible = X.shape[1]
+        self._init_weights()
+        self._setup_input(X, y)
+        self._train()
+
+    def _init_weights(self):
+        self.W = np.random.randn(self.n_visible, self.n_hidden) * 0.1
+
+        # Bias for visible and hidden units
+        self.bias_v = np.zeros(self.n_visible, dtype=np.float32)
+        self.bias_h = np.zeros(self.n_hidden, dtype=np.float32)
+
+        self.errors = []
+
+    def _train(self):
+        """Use CD-1 training procedure, basically an exact inference for `positive_associations`,
+        followed by a "non burn-in" block Gibbs Sampling for the `negative_associations`."""
+
+        for i in range(self.max_epochs):
+            error = 0
+            for batch in batch_iterator(self.X, batch_size=self.batch_size):
+                positive_hidden = sigmoid(np.dot(batch, self.W) + self.bias_h)
+                hidden_states = self._sample(positive_hidden)  # sample hidden state h1
+                positive_associations = np.dot(batch.T, positive_hidden)
+
+                negative_visible = sigmoid(
+                    np.dot(hidden_states, self.W.T) + self.bias_v
+                )
+                negative_visible = self._sample(
+                    negative_visible
+                )  # use the sampled hidden state h1 to sample v1
+                negative_hidden = sigmoid(
+                    np.dot(negative_visible, self.W) + self.bias_h
+                )
+                negative_associations = np.dot(negative_visible.T, negative_hidden)
+
+                lr = self.lr / float(batch.shape[0])
+                self.W += lr * (
+                    (positive_associations - negative_associations)
+                    / float(self.batch_size)
+                )
+                self.bias_h += lr * (
+                    negative_hidden.sum(axis=0) - negative_associations.sum(axis=0)
+                )
+                self.bias_v += lr * (
+                    np.asarray(batch.sum(axis=0)).squeeze()
+                    - negative_visible.sum(axis=0)
+                )
+
+                error += np.sum((batch - negative_visible) ** 2)
+
+            self.errors.append(error)
+            logging.info("Iteration %s, error %s" % (i, error))
+        logging.debug("Weights: %s" % self.W)
+        logging.debug("Hidden bias: %s" % self.bias_h)
+        logging.debug("Visible bias: %s" % self.bias_v)
+
+    def _sample(self, X):
+        return X > np.random.random_sample(size=X.shape)
+
+    def _predict(self, X=None):
+        return sigmoid(np.dot(X, self.W) + self.bias_h)

mla/rl/__init__.py

@@ -0,0 +1 @@
+# coding:utf-8

mla/rl/dqn.py

@@ -0,0 +1,158 @@
+# coding:utf-8
+import logging
+import random
+
+import gym
+import numpy as np
+from gym import wrappers
+
+np.random.seed(9999)
+
+logger = logging.getLogger()
+logger.setLevel(logging.INFO)
+
+"""
+References:
+    Sutton, Barto (2017). Reinforcement Learning: An Introduction. MIT Press, Cambridge, MA.
+"""
+
+
+class DQN(object):
+    def __init__(
+        self,
+        n_episodes=500,
+        gamma=0.99,
+        batch_size=32,
+        epsilon=1.0,
+        decay=0.005,
+        min_epsilon=0.1,
+        memory_limit=500,
+    ):
+        """Deep Q learning implementation.
+
+        Parameters
+        ----------
+
+        min_epsilon : float
+            Minimal value for epsilon.
+        epsilon : float
+            ε-greedy value.
+        decay : float
+            Epsilon decay rate.
+        memory_limit : int
+            Limit of experience replay memory.
+
+        """
+
+        self.memory_limit = memory_limit
+        self.min_epsilon = min_epsilon
+        self.gamma = gamma
+        self.epsilon = epsilon
+        self.n_episodes = n_episodes
+        self.batch_size = batch_size
+        self.decay = decay
+
+    def init_environment(self, name="CartPole-v0", monitor=False):
+        self.env = gym.make(name)
+        if monitor:
+            self.env = wrappers.Monitor(
+                self.env, name, force=True, video_callable=False
+            )
+
+        self.n_states = self.env.observation_space.shape[0]
+        self.n_actions = self.env.action_space.n
+
+        # Experience replay
+        self.replay = []
+
+    def init_model(self, model):
+        self.model = model(self.n_actions, self.batch_size)
+
+    def train(self, render=False):
+        max_reward = 0
+
+        for ep in range(self.n_episodes):
+            state = self.env.reset()
+
+            total_reward = 0
+
+            while True:
+                if render:
+                    self.env.render()
+
+                if np.random.rand() <= self.epsilon:
+                    # Exploration
+                    action = np.random.randint(self.n_actions)
+                else:
+                    # Exploitation
+                    action = np.argmax(self.model.predict(state[np.newaxis, :])[0])
+
+                # Run one timestep of the environment
+                new_state, reward, done, _ = self.env.step(action)
+                self.replay.append([state, action, reward, new_state, done])
+
+                # Sample batch from experience replay
+                batch_size = min(len(self.replay), self.batch_size)
+                batch = random.sample(self.replay, batch_size)
+
+                X = np.zeros((batch_size, self.n_states))
+                y = np.zeros((batch_size, self.n_actions))
+
+                states = np.array([b[0] for b in batch])
+                new_states = np.array([b[3] for b in batch])
+
+                Q = self.model.predict(states)
+                new_Q = self.model.predict(new_states)
+
+                # Construct training data
+                for i in range(batch_size):
+                    state_r, action_r, reward_r, new_state_r, done_r = batch[i]
+                    target = Q[i]
+
+                    if done_r:
+                        target[action_r] = reward_r
+                    else:
+                        target[action_r] = reward_r + self.gamma * np.amax(new_Q[i])
+
+                    X[i, :] = state_r
+                    y[i, :] = target
+
+                # Train deep learning model
+                self.model.fit(X, y)
+
+                total_reward += reward
+                state = new_state
+
+                if done:
+                    # Exit from current episode
+                    break
+
+            # Remove old entries from replay memory
+            while len(self.replay) > self.memory_limit:
+                self.replay.pop(0)
+
+            self.epsilon = self.min_epsilon + (1.0 - self.min_epsilon) * np.exp(
+                -self.decay * ep
+            )
+
+            max_reward = max(max_reward, total_reward)
+            logger.info(
+                "Episode: %s, reward %s,  epsilon %s, max reward %s"
+                % (ep, total_reward, self.epsilon, max_reward)
+            )
+        logging.info("Training finished.")
+
+    def play(self, episodes):
+        for i in range(episodes):
+            state = self.env.reset()
+            total_reward = 0
+
+            while True:
+                self.env.render()
+                action = np.argmax(self.model.predict(state[np.newaxis, :])[0])
+                state, reward, done, _ = self.env.step(action)
+                total_reward += reward
+                if done:
+                    break
+            logger.info("Episode: %s, reward %s" % (i, total_reward))
+        self.env.close()

mla/svm/__init__.py

@@ -0,0 +1 @@
+# coding:utf-8

mla/svm/kernerls.py

@@ -0,0 +1,35 @@
+# coding:utf-8
+import numpy as np
+import scipy.spatial.distance as dist
+
+
+class Linear(object):
+    def __call__(self, x, y):
+        return np.dot(x, y.T)
+
+    def __repr__(self):
+        return "Linear kernel"
+
+
+class Poly(object):
+    def __init__(self, degree=2):
+        self.degree = degree
+
+    def __call__(self, x, y):
+        return np.dot(x, y.T) ** self.degree
+
+    def __repr__(self):
+        return "Poly kernel"
+
+
+class RBF(object):
+    def __init__(self, gamma=0.1):
+        self.gamma = gamma
+
+    def __call__(self, x, y):
+        x = np.atleast_2d(x)
+        y = np.atleast_2d(y)
+        return np.exp(-self.gamma * dist.cdist(x, y) ** 2).flatten()
+
+    def __repr__(self):
+        return "RBF kernel"

mla/svm/svm.py

@@ -0,0 +1,145 @@
+# coding:utf-8
+import logging
+
+import numpy as np
+
+from mla.base import BaseEstimator
+from mla.svm.kernerls import Linear
+
+np.random.seed(9999)
+
+"""
+References:
+The Simplified SMO Algorithm http://cs229.stanford.edu/materials/smo.pdf
+"""
+
+
+class SVM(BaseEstimator):
+    def __init__(self, C=1.0, kernel=None, tol=1e-3, max_iter=100):
+        """Support vector machines implementation using simplified SMO optimization.
+
+        Parameters
+        ----------
+        C : float, default 1.0
+        kernel : Kernel object
+        tol : float , default 1e-3
+        max_iter : int, default 100
+        """
+        self.C = C
+        self.tol = tol
+        self.max_iter = max_iter
+        if kernel is None:
+            self.kernel = Linear()
+        else:
+            self.kernel = kernel
+
+        self.b = 0
+        self.alpha = None
+        self.K = None
+
+    def fit(self, X, y=None):
+        self._setup_input(X, y)
+        self.K = np.zeros((self.n_samples, self.n_samples))
+        for i in range(self.n_samples):
+            self.K[:, i] = self.kernel(self.X, self.X[i, :])
+        self.alpha = np.zeros(self.n_samples)
+        self.sv_idx = np.arange(0, self.n_samples)
+        return self._train()
+
+    def _train(self):
+        iters = 0
+        while iters < self.max_iter:
+            iters += 1
+            alpha_prev = np.copy(self.alpha)
+
+            for j in range(self.n_samples):
+                # Pick random i
+                i = self.random_index(j)
+
+                eta = 2.0 * self.K[i, j] - self.K[i, i] - self.K[j, j]
+                if eta >= 0:
+                    continue
+                L, H = self._find_bounds(i, j)
+
+                # Error for current examples
+                e_i, e_j = self._error(i), self._error(j)
+
+                # Save old alphas
+                alpha_io, alpha_jo = self.alpha[i], self.alpha[j]
+
+                # Update alpha
+                self.alpha[j] -= (self.y[j] * (e_i - e_j)) / eta
+                self.alpha[j] = self.clip(self.alpha[j], H, L)
+
+                self.alpha[i] = self.alpha[i] + self.y[i] * self.y[j] * (
+                    alpha_jo - self.alpha[j]
+                )
+
+                # Find intercept
+                b1 = (
+                    self.b
+                    - e_i
+                    - self.y[i] * (self.alpha[i] - alpha_io) * self.K[i, i]
+                    - self.y[j] * (self.alpha[j] - alpha_jo) * self.K[i, j]
+                )
+                b2 = (
+                    self.b
+                    - e_j
+                    - self.y[j] * (self.alpha[j] - alpha_jo) * self.K[j, j]
+                    - self.y[i] * (self.alpha[i] - alpha_io) * self.K[i, j]
+                )
+                if 0 < self.alpha[i] < self.C:
+                    self.b = b1
+                elif 0 < self.alpha[j] < self.C:
+                    self.b = b2
+                else:
+                    self.b = 0.5 * (b1 + b2)
+
+            # Check convergence
+            diff = np.linalg.norm(self.alpha - alpha_prev)
+            if diff < self.tol:
+                break
+        logging.info("Convergence has reached after %s." % iters)
+
+        # Save support vectors index
+        self.sv_idx = np.where(self.alpha > 0)[0]
+
+    def _predict(self, X=None):
+        n = X.shape[0]
+        result = np.zeros(n)
+        for i in range(n):
+            result[i] = np.sign(self._predict_row(X[i, :]))
+        return result
+
+    def _predict_row(self, X):
+        k_v = self.kernel(self.X[self.sv_idx], X)
+        return np.dot((self.alpha[self.sv_idx] * self.y[self.sv_idx]).T, k_v.T) + self.b
+
+    def clip(self, alpha, H, L):
+        if alpha > H:
+            alpha = H
+        if alpha < L:
+            alpha = L
+        return alpha
+
+    def _error(self, i):
+        """Error for single example."""
+        return self._predict_row(self.X[i]) - self.y[i]
+
+    def _find_bounds(self, i, j):
+        """Find L and H such that L <= alpha <= H.
+        Also, alpha must satisfy the constraint 0 <= αlpha <= C.
+        """
+        if self.y[i] != self.y[j]:
+            L = max(0, self.alpha[j] - self.alpha[i])
+            H = min(self.C, self.C - self.alpha[i] + self.alpha[j])
+        else:
+            L = max(0, self.alpha[i] + self.alpha[j] - self.C)
+            H = min(self.C, self.alpha[i] + self.alpha[j])
+        return L, H
+
+    def random_index(self, z):
+        i = z
+        while i == z:
+            i = np.random.randint(0, self.n_samples - 1)
+        return i

mla/tests/test_classification_accuracy.py

@@ -0,0 +1,114 @@
+from sklearn.metrics import roc_auc_score
+
+from mla.ensemble import RandomForestClassifier
+from mla.ensemble.gbm import GradientBoostingClassifier
+from mla.knn import KNNClassifier
+from mla.linear_models import LogisticRegression
+from mla.metrics import accuracy
+from mla.naive_bayes import NaiveBayesClassifier
+from mla.neuralnet import NeuralNet
+from mla.neuralnet.constraints import MaxNorm
+from mla.neuralnet.layers import Activation, Dense, Dropout
+from mla.neuralnet.optimizers import Adadelta
+from mla.neuralnet.parameters import Parameters
+from mla.neuralnet.regularizers import L2
+from mla.svm.kernerls import RBF, Linear
+from mla.svm.svm import SVM
+from mla.utils import one_hot
+
+try:
+    from sklearn.model_selection import train_test_split
+except ImportError:
+    from sklearn.cross_validation import train_test_split
+from sklearn.datasets import make_classification
+
+# Generate a random regression problem
+X, y = make_classification(
+    n_samples=750,
+    n_features=10,
+    n_informative=8,
+    random_state=1111,
+    n_classes=2,
+    class_sep=2.5,
+    n_redundant=0,
+)
+X_train, X_test, y_train, y_test = train_test_split(
+    X, y, test_size=0.12, random_state=1111
+)
+
+
+# All classifiers except convnet, RNN, LSTM.
+
+
+def test_linear_model():
+    model = LogisticRegression(lr=0.01, max_iters=500, penalty="l1", C=0.01)
+    model.fit(X_train, y_train)
+    predictions = model.predict(X_test)
+    assert roc_auc_score(y_test, predictions) >= 0.95
+
+
+def test_random_forest():
+    model = RandomForestClassifier(n_estimators=10, max_depth=4)
+    model.fit(X_train, y_train)
+    predictions = model.predict(X_test)[:, 1]
+    assert roc_auc_score(y_test, predictions) >= 0.95
+
+
+def test_svm_classification():
+    y_signed_train = (y_train * 2) - 1
+    y_signed_test = (y_test * 2) - 1
+
+    for kernel in [RBF(gamma=0.05), Linear()]:
+        model = SVM(max_iter=500, kernel=kernel)
+        model.fit(X_train, y_signed_train)
+        predictions = model.predict(X_test)
+        assert accuracy(y_signed_test, predictions) >= 0.8
+
+
+def test_mlp():
+    y_train_onehot = one_hot(y_train)
+    y_test_onehot = one_hot(y_test)
+
+    model = NeuralNet(
+        layers=[
+            Dense(256, Parameters(init="uniform", regularizers={"W": L2(0.05)})),
+            Activation("relu"),
+            Dropout(0.5),
+            Dense(128, Parameters(init="normal", constraints={"W": MaxNorm()})),
+            Activation("relu"),
+            Dense(2),
+            Activation("softmax"),
+        ],
+        loss="categorical_crossentropy",
+        optimizer=Adadelta(),
+        metric="accuracy",
+        batch_size=64,
+        max_epochs=25,
+    )
+    model.fit(X_train, y_train_onehot)
+    predictions = model.predict(X_test)
+    assert roc_auc_score(y_test_onehot[:, 0], predictions[:, 0]) >= 0.95
+
+
+def test_gbm():
+    model = GradientBoostingClassifier(
+        n_estimators=25, max_depth=3, max_features=5, learning_rate=0.1
+    )
+    model.fit(X_train, y_train)
+    predictions = model.predict(X_test)
+    assert roc_auc_score(y_test, predictions) >= 0.95
+
+
+def test_naive_bayes():
+    model = NaiveBayesClassifier()
+    model.fit(X_train, y_train)
+    predictions = model.predict(X_test)[:, 1]
+    assert roc_auc_score(y_test, predictions) >= 0.95
+
+
+def test_knn():
+    clf = KNNClassifier(k=5)
+
+    clf.fit(X_train, y_train)
+    predictions = clf.predict(X_test)
+    assert accuracy(y_test, predictions) >= 0.95

mla/tests/test_reduction.py

@@ -0,0 +1,46 @@
+# coding=utf-8
+import pytest
+from sklearn.datasets import make_classification
+from sklearn.metrics import roc_auc_score
+
+try:
+    from sklearn.model_selection import train_test_split
+except ImportError:
+    from sklearn.cross_validation import train_test_split
+
+from mla.ensemble import RandomForestClassifier
+from mla.pca import PCA
+
+
+@pytest.fixture
+def dataset():
+    # Generate a random binary classification problem.
+    return make_classification(
+        n_samples=1000,
+        n_features=100,
+        n_informative=75,
+        random_state=1111,
+        n_classes=2,
+        class_sep=2.5,
+    )
+
+
+# TODO: fix
+@pytest.mark.skip()
+def test_PCA(dataset):
+    X, y = dataset
+    X_train, X_test, y_train, y_test = train_test_split(
+        X, y, test_size=0.25, random_state=1111
+    )
+    p = PCA(50, solver="eigen")
+
+    # fit PCA with training set, not the entire dataset
+    p.fit(X_train)
+    X_train_reduced = p.transform(X_train)
+    X_test_reduced = p.transform(X_test)
+
+    model = RandomForestClassifier(n_estimators=25, max_depth=5)
+    model.fit(X_train_reduced, y_train)
+    predictions = model.predict(X_test_reduced)[:, 1]
+    score = roc_auc_score(y_test, predictions)
+    assert score >= 0.75

mla/tests/test_regression_accuracy.py

@@ -0,0 +1,61 @@
+try:
+    from sklearn.model_selection import train_test_split
+except ImportError:
+    from sklearn.cross_validation import train_test_split
+from sklearn.datasets import make_regression
+
+from mla.knn import KNNRegressor
+from mla.linear_models import LinearRegression
+from mla.metrics.metrics import mean_squared_error
+from mla.neuralnet import NeuralNet
+from mla.neuralnet.layers import Activation, Dense
+from mla.neuralnet.optimizers import Adam
+from mla.neuralnet.parameters import Parameters
+
+# Generate a random regression problem
+X, y = make_regression(
+    n_samples=1000,
+    n_features=10,
+    n_informative=10,
+    n_targets=1,
+    noise=0.05,
+    random_state=1111,
+    bias=0.5,
+)
+X_train, X_test, y_train, y_test = train_test_split(
+    X, y, test_size=0.25, random_state=1111
+)
+
+
+def test_linear():
+    model = LinearRegression(lr=0.01, max_iters=2000, penalty="l2", C=0.003)
+    model.fit(X_train, y_train)
+    predictions = model.predict(X_test)
+    assert mean_squared_error(y_test, predictions) < 0.25
+
+
+def test_mlp():
+    model = NeuralNet(
+        layers=[
+            Dense(16, Parameters(init="normal")),
+            Activation("linear"),
+            Dense(8, Parameters(init="normal")),
+            Activation("linear"),
+            Dense(1),
+        ],
+        loss="mse",
+        optimizer=Adam(),
+        metric="mse",
+        batch_size=64,
+        max_epochs=150,
+    )
+    model.fit(X_train, y_train)
+    predictions = model.predict(X_test)
+    assert mean_squared_error(y_test, predictions.flatten()) < 1.0
+
+
+def test_knn():
+    model = KNNRegressor(k=5)
+    model.fit(X_train, y_train)
+    predictions = model.predict(X_test)
+    assert mean_squared_error(y_test, predictions) < 10000

mla/tsne.py

@@ -0,0 +1,136 @@
+# coding:utf-8
+import logging
+
+import numpy as np
+
+from mla.base import BaseEstimator
+from mla.metrics.distance import l2_distance
+
+np.random.seed(999)
+
+"""
+References:
+https://lvdmaaten.github.io/tsne/
+Based on:
+https://lvdmaaten.github.io/tsne/code/tsne_python.zip
+"""
+
+
+class TSNE(BaseEstimator):
+    y_required = False
+
+    def __init__(
+        self, n_components=2, perplexity=30.0, max_iter=200, learning_rate=500
+    ):
+        """A t-Distributed Stochastic Neighbor Embedding implementation.
+
+        Parameters
+        ----------
+        max_iter : int, default 200
+        perplexity : float, default 30.0
+        n_components : int, default 2
+        """
+        self.max_iter = max_iter
+        self.perplexity = perplexity
+        self.n_components = n_components
+        self.initial_momentum = 0.5
+        self.final_momentum = 0.8
+        self.min_gain = 0.01
+        self.lr = learning_rate
+        self.tol = 1e-5
+        self.perplexity_tries = 50
+
+    def fit_transform(self, X, y=None):
+        self._setup_input(X, y)
+
+        Y = np.random.randn(self.n_samples, self.n_components)
+        velocity = np.zeros_like(Y)
+        gains = np.ones_like(Y)
+
+        P = self._get_pairwise_affinities(X)
+
+        iter_num = 0
+        while iter_num < self.max_iter:
+            iter_num += 1
+
+            D = l2_distance(Y)
+            Q = self._q_distribution(D)
+
+            # Normalizer q distribution
+            Q_n = Q / np.sum(Q)
+
+            # Early exaggeration & momentum
+            pmul = 4.0 if iter_num < 100 else 1.0
+            momentum = 0.5 if iter_num < 20 else 0.8
+
+            # Perform gradient step
+            grads = np.zeros(Y.shape)
+            for i in range(self.n_samples):
+                grad = 4 * np.dot((pmul * P[i] - Q_n[i]) * Q[i], Y[i] - Y)
+                grads[i] = grad
+
+            gains = (gains + 0.2) * ((grads > 0) != (velocity > 0)) + (gains * 0.8) * (
+                (grads > 0) == (velocity > 0)
+            )
+            gains = gains.clip(min=self.min_gain)
+
+            velocity = momentum * velocity - self.lr * (gains * grads)
+            Y += velocity
+            Y = Y - np.mean(Y, 0)
+
+            error = np.sum(P * np.log(P / Q_n))
+            logging.info("Iteration %s, error %s" % (iter_num, error))
+        return Y
+
+    def _get_pairwise_affinities(self, X):
+        """Computes pairwise affinities."""
+        affines = np.zeros((self.n_samples, self.n_samples), dtype=np.float32)
+        target_entropy = np.log(self.perplexity)
+        distances = l2_distance(X)
+
+        for i in range(self.n_samples):
+            affines[i, :] = self._binary_search(distances[i], target_entropy)
+
+        # Fill diagonal with near zero value
+        np.fill_diagonal(affines, 1.0e-12)
+
+        affines = affines.clip(min=1e-100)
+        affines = (affines + affines.T) / (2 * self.n_samples)
+        return affines
+
+    def _binary_search(self, dist, target_entropy):
+        """Performs binary search to find suitable precision."""
+        precision_min = 0
+        precision_max = 1.0e15
+        precision = 1.0e5
+
+        for _ in range(self.perplexity_tries):
+            denom = np.sum(np.exp(-dist[dist > 0.0] / precision))
+            beta = np.exp(-dist / precision) / denom
+
+            # Exclude zeros
+            g_beta = beta[beta > 0.0]
+            entropy = -np.sum(g_beta * np.log2(g_beta))
+
+            error = entropy - target_entropy
+
+            if error > 0:
+                # Decrease precision
+                precision_max = precision
+                precision = (precision + precision_min) / 2.0
+            else:
+                # Increase precision
+                precision_min = precision
+                precision = (precision + precision_max) / 2.0
+
+            if np.abs(error) < self.tol:
+                break
+
+        return beta
+
+    def _q_distribution(self, D):
+        """Computes Student t-distribution."""
+        Q = 1.0 / (1.0 + D)
+        np.fill_diagonal(Q, 0.0)
+        Q = Q.clip(min=1e-100)
+        return Q

mla/utils/__init__.py

@@ -0,0 +1,3 @@
+# coding:utf-8
+
+from .main import *

mla/utils/main.py

@@ -0,0 +1,25 @@
+# coding:utf-8
+import numpy as np
+
+
+def one_hot(y):
+    n_values = np.max(y) + 1
+    return np.eye(n_values)[y]
+
+
+def batch_iterator(X, batch_size=64):
+    """Splits X into equal sized chunks."""
+    n_samples = X.shape[0]
+    n_batches = n_samples // batch_size
+    batch_end = 0
+
+    for b in range(n_batches):
+        batch_begin = b * batch_size
+        batch_end = batch_begin + batch_size
+
+        X_batch = X[batch_begin:batch_end]
+
+        yield X_batch
+
+    if n_batches * batch_size < n_samples:
+        yield X[batch_end:]

requirements.txt

@@ -0,0 +1,8 @@
+tqdm
+matplotlib>=1.5.1
+numpy>=1.11.1
+scikit-learn>=0.18
+scipy>=0.18.0
+seaborn>=0.7.1
+autograd>=1.1.7
+gym

setup.cfg

@@ -0,0 +1,8 @@
+[bdist_wheel]
+universal=1
+
+[metadata]
+description-file=README.md
+
+[flake8]
+max-line-length = 120

setup.py

@@ -0,0 +1,35 @@
+from setuptools import setup, find_packages
+from codecs import open
+from os import path
+
+__version__ = '0.0.1'
+
+here = path.abspath(path.dirname(__file__))
+
+# Get the long description from the README file
+with open(path.join(here, 'README.md'), encoding='utf-8') as f:
+    long_description = f.read()
+
+# get the dependencies and installs
+with open(path.join(here, 'requirements.txt'), encoding='utf-8') as f:
+    all_reqs = f.read().split('\n')
+
+install_requires = [x.strip() for x in all_reqs if 'git+' not in x]
+dependency_links = [x.strip().replace('git+', '') for x in all_reqs if x.startswith('git+')]
+
+setup(
+    name='mla',
+    version=__version__,
+    description='A collection of minimal and clean implementations of machine learning algorithms.',
+    long_description=long_description,
+    url='https://github.com/rushter/mla',
+    download_url='https://github.com/rushter/mla/tarball/' + __version__,
+    license='MIT',
+    packages=find_packages(exclude=['docs', 'tests*']),
+    include_package_data=True,
+    author='Artem Golubin',
+    install_requires=install_requires,
+    setup_requires=['numpy>=1.10', 'scipy>=0.17'],
+    dependency_links=dependency_links,
+    author_email='gh@rushter.com'
+)