add backends, filters, learning models, optimizers

README.md

@@ -1,9 +1,8 @@
 # PyRoboLearn
 
 This repository contains the code for the *PyRoboLearn* (PRL) framework: a Python framework for Robot Learning.
-This framework revolves mainly around 7 axes: simulators, worlds, robots, interfaces, learning tasks (= environment 
-+ policy), learning models, and learning algorithms. 
+This framework revolves mainly around 7 axes: simulators, worlds, robots, interfaces, learning tasks (= environment and policy), learning models, and learning algorithms. 
 
 ## Requirements
 
-The framework has been tested with Python 2.7 and Ubuntu 16.04.
+The framework has been tested with Python 2.7 and Ubuntu 16.04.

pyrobolearn/__init__.py

@@ -52,6 +52,9 @@ import tasks
 
 # import metrics
 
+# import optimizers
+# import optim
+
 # import algos
 import algos
 
@@ -62,7 +65,7 @@ import algos
 __author__ = "Brian Delhaisse"
 __copyright__ = "Copyright 2018, PyRoboLearn"
 __credits__ = ["Brian Delhaisse"]
-__license__ = "(c) Brian Delhaisse"
+__license__ = "MIT"
 __version__ = "1.0.0"
 __maintainer__ = "Brian Delhaisse"
 __email__ = "briandelhaisse@gmail.com"

pyrobolearn/actions/action.py

@@ -15,7 +15,7 @@ from pyrobolearn.utils.data_structures.orderedset import OrderedSet
 __author__ = "Brian Delhaisse"
 __copyright__ = "Copyright 2018, PyRoboLearn"
 __credits__ = ["Brian Delhaisse"]
-__license__ = "(c) Brian Delhaisse"
+__license__ = "MIT"
 __version__ = "1.0.0"
 __maintainer__ = "Brian Delhaisse"
 __email__ = "briandelhaisse@gmail.com"

pyrobolearn/actions/gym_actions.py

@@ -13,7 +13,7 @@ from action import Action
 __author__ = "Brian Delhaisse"
 __copyright__ = "Copyright 2018, PyRoboLearn"
 __credits__ = ["Brian Delhaisse"]
-__license__ = "(c) Brian Delhaisse"
+__license__ = "MIT"
 __version__ = "1.0.0"
 __maintainer__ = "Brian Delhaisse"
 __email__ = "briandelhaisse@gmail.com"

pyrobolearn/actions/robot_actions/joint_actions.py

@@ -12,7 +12,7 @@ from robot_actions import RobotAction
 __author__ = "Brian Delhaisse"
 __copyright__ = "Copyright 2018, PyRoboLearn"
 __credits__ = ["Brian Delhaisse"]
-__license__ = "(c) Brian Delhaisse"
+__license__ = "MIT"
 __version__ = "1.0.0"
 __maintainer__ = "Brian Delhaisse"
 __email__ = "briandelhaisse@gmail.com"

pyrobolearn/actions/robot_actions/link_actions.py

@@ -11,7 +11,7 @@ from robot_actions import RobotAction
 __author__ = "Brian Delhaisse"
 __copyright__ = "Copyright 2018, PyRoboLearn"
 __credits__ = ["Brian Delhaisse"]
-__license__ = "(c) Brian Delhaisse"
+__license__ = "MIT"
 __version__ = "1.0.0"
 __maintainer__ = "Brian Delhaisse"
 __email__ = "briandelhaisse@gmail.com"

pyrobolearn/actions/robot_actions/robot_actions.py

@@ -16,7 +16,7 @@ from pyrobolearn.robots import Robot
 __author__ = "Brian Delhaisse"
 __copyright__ = "Copyright 2018, PyRoboLearn"
 __credits__ = ["Brian Delhaisse"]
-__license__ = "(c) Brian Delhaisse"
+__license__ = "MIT"
 __version__ = "1.0.0"
 __maintainer__ = "Brian Delhaisse"
 __email__ = "briandelhaisse@gmail.com"

pyrobolearn/backends/__init__.py

@@ -0,0 +1,27 @@
+# The various backends allows to specify which framework we want to use
+
+import os
+
+# check if the 'PYROBOLEARN_BACKEND' environment variable has been defined
+if 'PYROBOLEARN_BACKEND' not in os.environ:
+    os.environ['PYROBOLEARN_BACKEND'] = 'torch'
+
+# provide the various backends
+if os.environ['PYROBOLEARN_BACKEND'] == 'torch':
+    import torch_backend as backend
+elif os.environ['PYROBOLEARN_BACKEND'] == 'numpy':
+    import numpy_backend as backend
+else:
+    pass
+
+# alias
+b = backend
+
+
+# Tests #
+
+# create array
+a = b.array([1, 2, 3])
+print(type(a))
+print(a)
+print(a.reshape((3, 1)))

pyrobolearn/backends/numpy_backend.py

@@ -0,0 +1,20 @@
+
+from autograd import *
+from autograd.numpy import *
+import autograd.numpy as np
+import numpy
+
+
+def array(data, dtype=None, copy=True, device=None, requires_grad=False, ndmin=0):
+    return np.array(data, dtype=dtype, copy=copy, ndmin=ndmin)
+
+
+def inv(data, out=None):
+    # TODO: inverse not in autograd.numpy
+    if out is None:
+        return numpy.linalg.inv(data)
+
+    # inplace
+    output = numpy.linalg.inv(data)
+    for i in range(len(out)):
+        out[i] = output[i]

pyrobolearn/backends/torch_backend.py

@@ -0,0 +1,18 @@
+
+import torch
+from torch import *
+
+
+def array(data, dtype=None, copy=True, device=None, requires_grad=False, ndmin=0):
+    return torch.tensor(data, dtype=dtype, device=device, requires_grad=requires_grad)
+
+
+def inv(data, out=None):
+    return torch.inverse(data, out=out)
+
+
+def concatenate(data, axis=0, out=None):
+    return torch.cat(data, axis, out)
+
+# np.size vs torch.nelement() vs torch.size()
+

pyrobolearn/filters/__init__.py

@@ -0,0 +1,18 @@
+
+# import basic filter
+from filter import *
+
+# import histogram filter
+from histogram_filter import HistogramFilter
+
+# import kalman filter
+from kalman_filter import KalmanFilter
+
+# import extended kalman filter
+from extended_kalman_filter import EKF
+
+# import unscented kalman filter
+from unscented_kalman_filter import UKF
+
+# import particle filter
+from particle_filter import ParticleFilter

pyrobolearn/filters/description.txt

@@ -0,0 +1 @@
+filters --> state estimators

pyrobolearn/filters/extended_kalman_filter.py

@@ -0,0 +1,142 @@
+#!/usr/bin/env python
+"""Define the Extented Kalman Filter.
+"""
+
+import numpy as np
+
+__author__ = "Brian Delhaisse"
+__copyright__ = "Copyright 2018, PyRoboLearn"
+__credits__ = ["Brian Delhaisse"]
+__license__ = "MIT"
+__version__ = "1.0.0"
+__maintainer__ = "Brian Delhaisse"
+__email__ = "briandelhaisse@gmail.com"
+__status__ = "Development"
+
+
+class EKF(object):
+    r"""Extended Kalman Filter (EKF)
+
+    Type: Gaussian filter (=Bayes filter for continuous spaces)
+
+    This is an extension and generalization to the standard Kalman filter which also works with nonlinear systems.
+    This is more practical as state transitions and measurements are seldomly linear. The EKF calculates a Gaussian
+    distribution to the true belief, and represents thus an approximate of the belief (i.e. posterior).
+
+    The dynamic model :math:`p(x_t | x_{t-1}, u_t)` and the measurement model :math:`p(z_t | x_t)` are given by:
+
+    .. math::
+
+        x_t = f(x_{t-1}, u_t) + \epsilon_t \qquad \mbox{where} \qquad \epsilon_t \sim \mathcal{N}(0, R_t)
+        z_t = h(x_t) + \delta_t \qquad \mbox{where} \qquad \delta_t \sim \mathcal{N}(0, Q_t)
+
+    where :math:`f` and :math:`h` are nonlinear functions.
+
+    Because EKF cannot compute the true statistics in closed form, it uses linearization (via Taylor expansion).
+    That is, it approximates :math:`f` and :math:`h` by a linear function. Note that another way would be to compute
+    a Monte-carlo estimate of the Gaussian but would require a higher time-complexity (dependent on the number of
+    samples).
+
+    In a mathematical language, the linearization via Taylor expansion is given by:
+
+    .. math::
+
+        f(x_{t-1}, u_t) = f(\mu_{t-1}, u_t) + f'(\mu_{t-1}, u_t) (x_{t-1} - \mu_{t-1})
+        h(x_t) = h(\mu_t) + h'(\mu_t) (x_t - \mu_t)
+
+    where :math:`f'(\mu_{t-1}, u_t)` and :math:` h'(\mu_t)` represent the Jacobians evaluated at the means.
+
+    Notes: EKF requires to compute the Jacobians for the dynamical and measurement models.
+
+    Complexity: :math:`O(d^3 + n^2)` because of the matrix inversion when computing the kalman gain. :math:`d` is the
+    dimension of the measurement vector, and :math:`n` is the dimension of the state space.
+
+    References:
+        [1] "Probabilistic Robotics", Thrun et al., 2006 (sec 3.3)
+    """
+
+    def __init__(self, mean, covariance, dynamic_noise_cov, measurement_noise_cov):
+        """
+        Initialize the Extended Kalman Filter.
+
+        Args:
+            mean (float[N]): initial mean of the state belief
+            covariance (float[N,N]): initial covariance matrix of the state belief
+            dynamic_noise_cov (float[N,N]): dynamic noise covariance matrix
+            measurement_noise_cov (float[K,K]): measurement noise covariance matrix
+        """
+        self.mu = mean
+        self.Sigma = covariance
+        self.R = dynamic_noise_cov
+        self.Q = measurement_noise_cov
+        self.I = np.identity(self.Sigma.shape[0])
+
+    def predict(self, f, u):
+        """
+        Predict (a priori) the next state (without taking into account measurements) using the nonlinear dynamic model
+        with added Gaussian noise.
+
+        Args:
+            f (callable class): (nonlinear) dynamical function. This should have a method `jacobian`.
+            u (np.array): control array
+
+        Returns:
+            float[N]: a priori mean of the Gaussian belief
+            float[N,N]: a priori covariance of the Gaussian belief
+        """
+        F = f.jacobian(self.mu, u)
+
+        # predict a priori the next state (by only using the nonlinear dynamic model)
+        self.mu = f(self.mu, u)
+        self.Sigma = F.dot(self.Sigma).dot(F.T) + self.R
+
+        # return the a priori Gaussian belief on the next state
+        return self.mu, self.Sigma
+
+    def measurement_update(self, h, z):
+        """
+        Predict (a posteriori) the next state by incorporating the measurement :math:`z_t`.
+
+        Args:
+            h (callable class): (nonlinear) measurement function. This should have a method `jacobian`.
+            z (np.array): measurement array
+
+        Returns:
+            float[N]: a posteriori mean of the Gaussian belief
+            float[N,N]: a posteriori covariance matrix of the Gaussian belief
+        """
+        H = h.jacobian(self.mu)
+
+        # deviation error between the measurement and the a priori predicted state (aka innovation)
+        y = z - h(self.mu)
+
+        # Innovation covariance matrix
+        S = H.dot(self.Sigma).dot(H.T) + self.Q
+
+        # Kalman gain (which specifies the degree to which the measurement is incorporated into the new state)
+        K = self.Sigma.dot(H.T).dot(np.linalg.inv(S))
+
+        # predict a posteriori the next state
+        self.mu = self.mu + K.dot(y)
+        self.Sigma = (self.I - K.dot(H)).dot(self.Sigma)
+
+        # return the a posteriori Gaussian belief on the next state
+        return self.mu, self.Sigma
+
+    def compute(self, f, u, h, z):
+        """
+        Perform a prediction and measurement update step.
+
+        Args:
+            f (callable class): (nonlinear) dynamical function. This should have a method `jacobian`.
+            u (np.array): control array
+            h (callable class): (nonlinear) measurement function. This should have a method `jacobian`.
+            z (np.array): measurement array
+
+        Returns:
+            float[N]: a posteriori mean vector of the Gaussian belief
+            float[N,N]: a posteriori covariance matrix of the Gaussian belief
+        """
+        self.predict(f, u)
+        self.measurement_update(h, z)
+        return self.mu, self.Sigma

pyrobolearn/filters/filter.py

@@ -0,0 +1,8 @@
+# This file provides some common filters used in signal processing
+
+import scipy
+
+
+class Filter(object):
+    r"""Filter abstract class"""
+    pass

pyrobolearn/filters/histogram_filter.py

@@ -0,0 +1,100 @@
+#!/usr/bin/env python
+"""Define the Histogram Filter.
+"""
+
+import numpy as np
+
+__author__ = "Brian Delhaisse"
+__copyright__ = "Copyright 2018, PyRoboLearn"
+__credits__ = ["Brian Delhaisse"]
+__license__ = "MIT"
+__version__ = "1.0.0"
+__maintainer__ = "Brian Delhaisse"
+__email__ = "briandelhaisse@gmail.com"
+__status__ = "Development"
+
+
+class HistogramFilter(object):
+    r"""Histogram Filter (HF)
+
+    Type: Non-parametric filter
+
+    The histogram filter decomposes the continuous state space into finite regions, and represents the posterior
+    by a histogram. [1]
+
+    Notes:
+        - this filter is well-suited to represent complex multimodal belief [1]
+        - the complexity depends on the number of parameters
+
+    Complexity: :math:`O(M^N)` where :math:`M` is the number of regions/bins, and :math:`N` is the dimensionality
+        of the state.
+
+    References:
+        [1] "Probabilistic Robotics", Thrun et al., 2006 (sec 4.1)
+    """
+
+    def __init__(self, state_dim=1, num_bins_per_dim=10):
+        self.num_bins = num_bins_per_dim
+        self.state_dim = state_dim
+
+        # total number of parameters (=regions)
+        self.num_bins = num_bins_per_dim**state_dim
+
+        # initial belief: uniform distribution
+        p = 1./self.num_bins
+        self.p = np.full([num_bins_per_dim]*state_dim, p)
+
+    def predict(self, f, u):
+        """
+        Predict (a priori) the next state (without taking into account measurements) using the probabilistic
+        nonlinear dynamic model.
+
+        Args:
+            f (callable class/function): probabilistic (nonlinear) dynamical function that stacks on the axis 0 the
+                probability distribution for each region. That is the size of axis 0 should be
+                :math:`num(bins)^{dim(state)}`.
+            u (np.array): control array
+        """
+        # probability to arrive in region k given control input u from any region i (convolution operation)
+        s = ''.join([chr(i) for i in range(97, 97 + self.state_dim + 1)])
+        self.p = np.einsum(s+','+s[1:]+'->'+s[1:], f(self.p.shape, u), self.p)
+
+        # return belief
+        return self.p
+
+    def measurement_update(self, h, z):
+        """
+        Predict (a posteriori) the next state by incorporating the measurement :math:`z_t`.
+
+        Args:
+            h (callable class/function): probabilistic (nonlinear) measurement function to see measurement from
+                region k. This function should return an array of the same shape as the one provided in argument
+                where each cell contains the probability to see the measurement `z` from that region (i.e. index).
+            z (np.array): measurement array
+        """
+        # probability to see measurement z from region k (multiplication operation)
+        self.p = h(z, self.p.shape) * self.p
+
+        # normalize to get a proper probability distribution
+        self.p /= np.sum(self.p)
+
+        # return belief
+        return self.p
+
+    def compute(self, f, u, h, z):
+        """
+        Perform a prediction and measurement update step.
+
+        Args:
+            f (callable class/function): probabilistic (nonlinear) dynamical function that stacks on the axis 0 the
+                probability distribution for each region. That is the size of axis 0 should be
+                :math:`num(bins)^{dim(state)}`.
+            u (np.array): control array
+            h (callable class/function):  probabilistic (nonlinear) measurement function to see measurement from
+                region k. This function should return an array of the same shape as the one provided in argument
+                where each cell contains the probability to see the measurement `z` from that region (i.e. index).
+            z (np.array): measurement array
+        """
+        self.predict(f, u)
+        self.measurement_update(h, z)
+        return self.p

pyrobolearn/filters/kalman_filter.py

@@ -0,0 +1,141 @@
+#!/usr/bin/env python
+"""Define the Kalman Filter state estimator.
+"""
+
+import numpy as np
+
+__author__ = "Brian Delhaisse"
+__copyright__ = "Copyright 2018, PyRoboLearn"
+__credits__ = ["Brian Delhaisse"]
+__license__ = "MIT"
+__version__ = "1.0.0"
+__maintainer__ = "Brian Delhaisse"
+__email__ = "briandelhaisse@gmail.com"
+__status__ = "Development"
+
+
+class KalmanFilter(object):
+    r"""Kalman Filter (KF)
+
+    Type: Gaussian filter (=Bayes filter for continuous spaces)
+
+    The Kalman filter (also known as the linear quadratic estimator (LQE)) is a state estimator which "uses a series
+    of measurements observed over time, containing statistical noise and other inaccuracies, and produces estimates
+    of unknown variables that tend to be more accurate than those based on a single measurement alone, by estimating
+    a joint probability distribution over the variables for each timeframe." (Wikipedia)
+
+    Assumptions:
+        * Markov assumption
+        * linear dynamic model with added Gaussian noise
+        * linear measurement model with added Gaussian noise
+        * initial belief on the state is normally distributed
+
+    The Kalman filter consists of two steps:
+        * Prediction step: produces the a priori estimate of the current state variables, along with their
+            corresponding uncertainties, based on the dynamical model.
+        * Measurement update step: produces the a posteriori estimate of the state along with its uncertainty,
+            by taking into account the received measurements, comparing and incorporating them with the above
+            prediction.
+
+    The dynamic model :math:`p(x_t | x_{t-1}, u_t)` and the measurement model :math:`p(z_t | x_t)` are given by:
+
+    .. math::
+
+            x_t = A_t x_{t-1} + B_t u_t + \epsilon_t \qquad \mbox{where} \qquad \epsilon_t \sim \mathcal{N}(0, R_t)
+            z_t = C_t x_t + \delta_t \qquad \mbox{where} \qquad \delta_t \sim \mathcal{N}(0, Q_t)
+
+    Notes: The use of the KF " does not assume that the errors are Gaussian. However, the filter yields the exact
+    conditional probability estimate in the special case that all errors are Gaussian." (Wikipedia)
+
+    Complexity: :math:`O(d^3 + n^2)` because of the matrix inversion when computing the kalman gain. :math:`d` is the
+    dimension of the measurement vector, and :math:`n` is the dimension of the state space.
+
+    References:
+        [1] "A New Approach to Linear Filtering and Prediction Problems", Kalman, 1960
+        [2] "Probabilistic Robotics", Thrun et al., 2006 (sec 3.2)
+    """
+
+    def __init__(self, mean, covariance, dynamic_noise_cov, measurement_noise_cov):
+        """
+        Initialize the Kalman filter.
+
+        Args:
+            mean (float[N]): initial mean of the state belief
+            covariance (float[N,N]): initial covariance matrix of the state belief
+            dynamic_noise_cov (float[N,N]): dynamic noise covariance matrix
+            measurement_noise_cov (float[K,K]): measurement noise covariance matrix
+        """
+        # self.belief = Gaussian(...)
+        self.mu = mean
+        self.Sigma = covariance
+        self.R = dynamic_noise_cov
+        self.Q = measurement_noise_cov
+        self.I = np.identity(self.Sigma.shape[0])
+
+    def predict(self, A, B, u):
+        """
+        Predict (a priori) the next state (without taking into account measurements) using the linear dynamic model
+        with added Gaussian noise.
+
+        Args:
+            A (float[N,N]): linear state transformation matrix
+            B (float[N,M]): linear control transformation matrix
+            u (float[M]): control vector
+
+        Returns:
+            float[N]: a priori mean of the Gaussian belief
+            float[N,N]: a priori covariance of the Gaussian belief
+        """
+        # predict a priori the next state (by only using the linear dynamic model)
+        self.mu = A.dot(self.mu) + B.dot(u)
+        self.Sigma = A.dot(self.Sigma).dot(A.T) + self.R
+
+        # return the a priori Gaussian belief on the next state
+        return self.mu, self.Sigma
+
+    def measurement_update(self, C, z):
+        """
+        Predict (a posteriori) the next state by incorporating the measurement :math:`z_t`.
+
+        Args:
+            C (float[K,N]): linear state-measurement transformation matrix
+            z (float[K]): measurement vector
+
+        Returns:
+            float[N]: a posteriori mean vector of the Gaussian belief
+            float[N,N]: a posteriori covariance matrix of the Gaussian belief
+        """
+        # deviation error between the measurement and the a priori predicted state (aka innovation)
+        y = z - C.dot(self.mu)
+
+        # Innovation covariance matrix
+        S = C.dot(self.Sigma).dot(C.T) + self.Q
+
+        # Kalman gain (which specifies the degree to which the measurement is incorporated into the new state)
+        K = self.Sigma.dot(C.T).dot(np.linalg.inv(S))
+
+        # predict a posteriori the next state
+        self.mu = self.mu + K.dot(y)
+        self.Sigma = (self.I - K.dot(C)).dot(self.Sigma)
+
+        # return the a posteriori Gaussian belief on the next state
+        return self.mu, self.Sigma
+
+    def compute(self, A, B, u, C, z):
+        """
+        Perform a prediction and measurement update step.
+
+        Args:
+            A (float[N,N]): linear state transformation matrix
+            B (float[N,M]): linear control transformation matrix
+            u (float[M]): control vector
+            C (float[K,N]): linear state-measurement transformation matrix
+            z (float[K]): measurement vector
+
+        Returns:
+            float[N]: a posteriori mean vector of the Gaussian belief
+            float[N,N]: a posteriori covariance matrix of the Gaussian belief
+        """
+        self.predict(A, B, u)
+        self.measurement_update(C, z)
+        return self.mu, self.Sigma

pyrobolearn/filters/particle_filter.py

@@ -0,0 +1,96 @@
+#!/usr/bin/env python
+"""Define the Particle Filter.
+"""
+
+import numpy as np
+import bisect
+
+__author__ = "Brian Delhaisse"
+__copyright__ = "Copyright 2018, PyRoboLearn"
+__credits__ = ["Brian Delhaisse"]
+__license__ = "MIT"
+__version__ = "1.0.0"
+__maintainer__ = "Brian Delhaisse"
+__email__ = "briandelhaisse@gmail.com"
+__status__ = "Development"
+
+
+class Particle(object):
+    r"""Particle
+    """
+    def __init__(self, state):
+        self.state = state
+
+
+class ParticleFilter(object):
+    r"""Particle Filter
+
+    Type: Non-parametric filter
+
+    "Particle filtering uses a genetic mutation-selection sampling approach, with a set of particles
+    (also called samples) to represent the posterior distribution of some stochastic process given noisy
+    and/or partial observations." (Wikipedia)
+    In other words, "a particle is a hypothesis as to what the true world state may be at time t." [1]
+
+    Complexity: O(exp(n))
+
+    References:
+        [1] "Probabilistic Robotics", Thrun et al., 2006 (sec 4.3)
+    """
+
+    def __init__(self, state_min, state_max, num_particles=100):
+        # initialize particles randomly between the 2 given bounds (i.e. state_min and state_max)
+        self.particles = np.random.uniform(state_min, state_max, size=(num_particles, state_min.size))
+        self.weights = np.ones(num_particles)
+
+    @property
+    def num_particles(self):
+        return len(self.particles)
+
+    def predict(self, f, u):
+        # for each particle, predict next state of the particle given the control input
+        self.particles = [f(state, u) for state in self.particles]
+        return self.particles
+
+    def measurement_update(self, h, z):
+        # for each particle, compute the probability to see measurement z from the particle state
+        self.weights = [h(z, state) for state in self.particles]
+        return self.weights
+
+    def sampling(self):
+        # importance sampling (survival of the fittest)
+        particles = []
+
+        # resampling
+
+        # Wheel algorithm (from Udacity)
+        # idx = np.random.randint(0, self.num_particles)
+        # b, wmax = 0., max(self.weights)
+        # for i in range(self.num_particles):
+        #     b += np.random.random() * 2 * wmax
+        #     while self.weights[idx] < b:
+        #         b -= self.weights[idx]
+        #         idx = (idx+1) % self.num_particles
+        #     particles.append(self.particles[idx])
+
+        # resampling O(N*log(N)) algorithm: I observe that this one was better than the Wheel algorithm
+        # compute cumulative probability
+        sumProb = sum(self.weights)
+        cumulativeProb = [w/sumProb for w in self.weights]      # normalize
+        for i in range(1, self.num_particles):
+            cumulativeProb[i] += cumulativeProb[i-1]
+
+        # resample
+        for i in range(self.num_particles):
+            idx = bisect.bisect(cumulativeProb, np.random.uniform(0,1))
+            particles.append(self.particles[idx])
+
+        # return particles (=states)
+        self.particles = particles
+        return self.particles
+
+    def compute(self, f, u, h, z):
+        self.predict(f, u)
+        self.measurement_update(h, z)
+        self.sampling()
+        return self.particles

pyrobolearn/filters/unscented_kalman_filter.py

@@ -0,0 +1,177 @@
+#!/usr/bin/env python
+"""Define the Unscented Kalman Filter.
+"""
+
+import numpy as np
+
+__author__ = "Brian Delhaisse"
+__copyright__ = "Copyright 2018, PyRoboLearn"
+__credits__ = ["Brian Delhaisse"]
+__license__ = "MIT"
+__version__ = "1.0.0"
+__maintainer__ = "Brian Delhaisse"
+__email__ = "briandelhaisse@gmail.com"
+__status__ = "Development"
+
+
+class UKF(object):
+    r"""Unscented Kalman Filter (UKF)
+
+    Type: Derivative-free Gaussian filter (=Bayes filter for continuous spaces)
+
+    This is an extension and generalization to the standard Kalman filter which also works with nonlinear systems.
+    This filter is an alternative to deal with non-linear process/measurement models, using :math:`\sigma`-points
+    to approximate the probability distribution. Specifically, as EKF, it linearizes the transformation of a Gaussian,
+    but instead of doing it via a Taylor expansion, it "performs a stochastic linearization through the use of
+    a weighted statistical linear regression process". [1]
+
+    The dynamic model :math:`p(x_t | x_{t-1}, u_t)` and the measurement model :math:`p(z_t | x_t)` are given by:
+
+    .. math::
+
+        x_t = f(x_{t-1}, u_t) + \epsilon_t \qquad \mbox{where} \qquad \epsilon_t \sim \mathcal{N}(0, R_t)
+        z_t = h(x_t) + \delta_t \qquad \mbox{where} \qquad \delta_t \sim \mathcal{N}(0, Q_t)
+
+    where :math:`f` and :math:`h` are nonlinear functions.
+
+    The UKF works by generating :math:`2n+1` :math:`sigma`-points, and pass them through the corresponding dynamical
+    and measurement functions, and compute from them a weighted empirical mean and covariance.
+
+    Notes:
+        - Compared to EKF, this filter does not require to compute the Jacobians for the dynamical and measurement
+        models, and is thus a derivative-free filter. "If the belief is highly non-Gaussian, then the UKF
+        representation is too restrictive and the filter can perform arbitrarily poorly" [1]
+        - An extension to multimodal posteriors is known as multi-hypothesis Kalman filter which uses a mixture of
+        Gaussians. See also, the histogram and particle filters that are well-suited to represent complex multimodal
+        beliefs. [1]
+
+    Complexity: "The asympototic complexity of the UKF algorithm is the same as for the EKF. In practice, the EKF is
+        often slightly faster than the UKF. Nevertheless, the UKF is still highly efficient" [1]
+
+    References:
+        [1] "Probabilistic Robotics", Thrun et al., 2006 (sec 3.4)
+    """
+
+    def __init__(self, mean, covariance, dynamic_noise_cov, measurement_noise_cov, alpha=1., kappa=0, beta=2):
+        """
+        Initialize the Unscented Kalman Filter.
+
+        Args:
+            mean (float[N]): initial mean of the state belief
+            covariance (float[N,N]): initial covariance matrix of the state belief
+            dynamic_noise_cov (float[N,N]): dynamic noise covariance matrix
+            measurement_noise_cov (float[K,K]): measurement noise covariance matrix
+            alpha (float): control the spread of the sigma points (recommended value: 1e-3 or 1)
+            kappa (float): control the spread of the sigma points
+            beta (int): this parameter can be chosen to encode additional (higher order) knowledge about the
+                distribution underlying the Gaussian representation.
+        """
+        self.mu = mean
+        self.Sigma = covariance
+        self.n = self.mu.size
+        self.tau = alpha**2 * (self.n + kappa) - self.n
+        self.weight_mean_0 = self.tau / (self.n + self.tau)
+        self.weight_cov_0 = self.tau / (self.n + self.tau) + (1 - alpha**2 + beta)
+        self.weight = 1./(2.*(self.n + self.tau))
+        self.gamma = np.sqrt(self.n + self.tau)
+
+    def generate_sigma_points(self):
+        """
+        Generate the sigma points for the UKF.
+
+        Returns:
+            float[2N+1, N]: sigma points
+        """
+        term = self.gamma * np.sqrt(self.Sigma.T)
+        sigma_points = np.vstack((self.mu, self.mu + term, self.mu - term))
+        return sigma_points
+
+    def predict(self, f, u):
+        """
+        Predict (a priori) the next state (without taking into account measurements) using the nonlinear dynamic model
+        with added Gaussian noise.
+
+        Args:
+            f (callable class/function): (nonlinear) dynamical function.
+            u (np.array): control array
+
+        Returns:
+            float[N]: a priori mean of the Gaussian belief
+            float[N,N]: a priori covariance of the Gaussian belief
+        """
+        # generate sigma points
+        sigma_points = self.generate_sigma_points()
+
+        # feed each sigma point to the dynamical function
+        X = np.array([f(sigma, u) for sigma in sigma_points])
+
+        # compute the weighted empirical mean and covariance
+        self.mu = self.weight_mean_0 * X[0] + self.weight * X[1:].sum(axis=0)
+        diff0 = (X[0] - self.mu).reshape(-1,1)
+        diff = (X[1:] - self.mu).T
+        self.Sigma = self.weight_cov_0 * diff0.dot(diff0.T) + self.weight * diff.dot(diff.T) + self.R
+
+        # return the a priori Gaussian belief on the next state
+        return self.mu, self.Sigma
+
+    def measurement_update(self, h, z):
+        """
+        Predict (a posteriori) the next state by incorporating the measurement :math:`z_t`.
+
+        Args:
+            h (callable class/function): (nonlinear) measurement function.
+            z (np.array): measurement array
+
+        Returns:
+            float[N]: a posteriori mean of the Gaussian belief
+            float[N,N]: a posteriori covariance matrix of the Gaussian belief
+        """
+        # generate sigma points
+        sigma_points = self.generate_sigma_points()
+
+        # feed each sigma point to the measurement function
+        Z = np.array([h(sigma) for sigma in sigma_points])
+
+        # weighted empirical mean for the measurement
+        z_mean = self.weight_mean_0 * Z[0] + self.weight * Z[1:].sum(axis=0)
+
+        # deviation error between the measurement and the a priori predicted state (aka innovation)
+        y = z - z_mean
+
+        # compute the weighted cross-covariance
+        term0 = ((sigma_points[0] - self.mu).reshape(-1,1)).dot( (Z[0] - z_mean).reshape(1, -1) )
+        term = ((sigma_points[1:] - self.mu).T).dot(Z[1:] - z_mean)
+        C = self.weight_cov_0 * term0 + self.weight * term
+
+        # Innovation covariance matrix
+        diff0 = (Z[0] - z_mean).reshape(-1, 1)
+        diff = (Z[1:] - z_mean).T
+        S = self.weight_cov_0 * diff0.dot(diff0.T) + self.weight * diff.dot(diff.T) + self.Q
+
+        # Kalman gain (which specifies the degree to which the measurement is incorporated into the new state)
+        K = C.dot(np.linalg.inv(S))
+
+        # compute the weighted empirical mean and covariance
+        self.mu = self.mu + K.dot(y)
+        self.Sigma = self.Sigma - K.dot(S.dot(K.T))
+
+        # return the a posteriori Gaussian belief on the next state
+        return self.mu, self.Sigma
+
+    def compute(self, f, u, h, z):
+        """
+        Perform a prediction and measurement update step.
+
+        Args:
+            f (callable class/function): (nonlinear) dynamical function.
+            u (np.array): control array
+            h (callable class/function): (nonlinear) measurement function.
+            z (np.array): measurement array
+
+        Returns:
+            float[N]: a posteriori mean vector of the Gaussian belief
+            float[N,N]: a posteriori covariance matrix of the Gaussian belief
+        """
+        self.predict(f, u)
+        self.measurement_update(h, z)
+        return self.mu, self.Sigma

pyrobolearn/models/README.md

@@ -0,0 +1,3 @@
+## Learning models
+
+In this folder, we provide the various learning models. These can be categorized into two categories: movement primitives and general function approximators.

pyrobolearn/models/__init__.py

@@ -0,0 +1,44 @@
+
+# General model
+from model import Model
+
+# General Learning models #
+
+# Linear
+from linear import Linear
+
+# PCA
+from pca import PCA
+
+# Polynomial
+from polynomial import Polynomial, PolynomialFunction
+
+# Gaussian
+from gaussian import Gaussian, MVN      # MVN is an alias
+
+# GMM/GMR
+from gmm import GMM
+
+# GP
+from gp import GPR
+
+# HMM
+from hmm import *
+
+# DNN
+from nn import *
+
+
+# Learning model for trajectories (movement primitives) #
+
+# CPG
+from cpg import *
+
+# DMP
+from dmp import *
+
+# ProMP
+from promp import *
+
+# KMP
+from kmp import *

pyrobolearn/models/cpg.py

@@ -0,0 +1,716 @@
+#!/usr/bin/env python
+"""Define the Central Pattern Generator Node and Network classes
+
+Central Pattern Generators (CPGs) allows to model rhythmic movement primitives, and are often used in locomotion.
+This file provides the CPG Node and Network, where the network is composed of CPG nodes connected / phased together.
+"""
+
+import numpy as np
+import torch
+import matplotlib.pyplot as plt
+from matplotlib.widgets import Slider, CheckButtons
+from matplotlib.animation import FuncAnimation
+
+
+__author__ = "Brian Delhaisse"
+__copyright__ = "Copyright 2018, PyRoboLearn"
+__credits__ = ["Brian Delhaisse"]
+__license__ = "MIT"
+__version__ = "1.0.0"
+__maintainer__ = "Brian Delhaisse"
+__email__ = "briandelhaisse@gmail.com"
+__status__ = "Development"
+
+
+class CPGNode(object):
+    r"""Central Pattern Generator Node
+
+    Each CPG node :math:`i` is driven by the following differential equations:
+    * phase: :math:`\dot{\phi}_i &= \omega_i + \sum_j a_j w_{ij} \sin(\phi_j - \phi_i - \varphi_{ij})`
+    * amplitude: :math:`\ddot{a}_i &= K_a (A_i - a_i) - D_a \dot{a}_i`
+    * offset: :math:`\ddot{o}_i &= K_o (O_i - o_i) - D_o \dot{o}_i`
+    * target angle: :math:`\theta_i &= o_i + a_i \cos(\phi_i)`
+
+    References:
+        [1] "Central pattern generators for locomotion control in animals and robots: a review", Ijspeert, 2008
+    """
+
+    def __init__(self, id, phi=0, offset=0, amplitude=1., timesteps=100, freq=None):
+        self.id = id
+        self.timesteps = timesteps
+        self.dt = 1./self.timesteps
+        self.t = 0.
+
+        # gains
+        self.D_amp, self.D_offset = 20., 20.
+        self.K_amp, self.K_offset = self.D_amp**2/4., self.D_offset**2/4.
+
+        # state parameters
+        self.phi, self.curr_phi = phi, phi
+        self.amp, self.damp, self.curr_amp = amplitude, 0, amplitude
+        self.offset, self.doffset, self.curr_offset = offset, 0, offset
+
+        # control parameters
+        self.des_offset = offset
+        self.des_amp = amplitude
+        self.des_freq = 1. if freq is None else freq  # 1./self.timesteps if freq is None else freq
+        self.des_omega = 2 * np.pi * self.des_freq
+
+        # init values
+        self.init_phi = self.phi
+        # self.init_offset = self.offset
+        # self.init_amp = self.amp
+
+        # angle
+        self.theta = self.offset + self.amp * np.cos(self.phi)
+
+        # dict containing the CPG nodes connected to this node with their coupling weight and bias
+        self.nodes = {}
+        self.sliderNodes = {}
+
+        # plot
+        self.frame = None
+        self.fig = None
+        self.do_plot_offset, self.do_plot_amp, self.do_plot_phi, self.do_plot_theta = False, False, False, True
+        self.updated = True
+
+        # keep in memory the previous `timesteps` values of offset, amplitude, phi, and theta
+        self.offsets = [self.offset] * self.timesteps
+        self.amps = [self.amp] * self.timesteps
+        self.phis = [self.phi] * self.timesteps
+        self.thetas = [self.theta] * self.timesteps
+
+    ##############
+    # Properties #
+    ##############
+
+    @property
+    def des_offset(self):
+        """Return the desired offset."""
+        return self._des_offset
+
+    @des_offset.setter
+    def des_offset(self, offset):
+        """Set the desired offset."""
+        if offset is None:
+            offset = 0.
+        if not isinstance(offset, (int, float)):
+            raise TypeError("Expecting the desired offset to be an integer or float, got instead: "
+                            "{}".format(type(offset)))
+        self._des_offset = offset
+
+    @property
+    def des_amp(self):
+        """Return the desired amplitude."""
+        return self._des_amp
+
+    @des_amp.setter
+    def des_amp(self, amplitude):
+        """Set the desired amplitude."""
+        if amplitude is None:
+            amplitude = 0.
+        if not isinstance(amplitude, (int, float)):
+            raise TypeError("Expecting the desired offset to be an integer or float, got instead: "
+                            "{}".format(type(amplitude)))
+        self._des_amp = amplitude
+
+    @property
+    def des_freq(self):
+        """Return the desired frequency."""
+        return self._des_freq
+
+    @des_freq.setter
+    def des_freq(self, freq):
+        """Set the desired frequency."""
+        if freq is None:
+            freq = 1.
+        if not isinstance(freq, (int, float)):
+            raise TypeError("Expecting the desired frequency to be an integer or float, got instead: "
+                            "{}".format(type(freq)))
+        self._des_freq = freq
+        self.des_omega = 2 * np.pi * self._des_freq
+
+    @property
+    def num_coupling_nodes(self):
+        """Return the number of coupling nodes."""
+        return len(self.nodes)
+
+    # @property
+    # def phi(self):
+    #     """Return the phase."""
+    #     return self._phi
+    #
+    # @phi.setter
+    # def phi(self, phi):
+    #     """Set the phase."""
+    #     self.curr_phi, self.phi = phi, phi
+    #     self.theta = self.offset + self.amp * np.cos(self.phi)
+
+    def set_phi(self, phi):
+        self.curr_phi, self.phi = phi, phi
+        self.theta = self.offset + self.amp * np.cos(self.phi)
+
+    @property
+    def num_parameters(self):
+        return 3 + 2 * len(self.nodes)
+
+    ###########
+    # Methods #
+    ###########
+
+    def add_node(self, cpg_node, weight=0., bias=0.):
+        """Add a coupling node."""
+        self.nodes[cpg_node] = {'weight': weight, 'bias': bias}
+
+    def remove_node(self, cpg_node):
+        """Remove a coupling node"""
+        if cpg_node in self.nodes:
+            self.nodes.pop(cpg_node)
+
+    def reset(self):
+        """Reset the phase of the CPG node; this can be useful for phase resetting."""
+        # Re-initialize the CPG
+        self.curr_phi, self.phi = self.init_phi, self.init_phi
+        self.theta = self.offset + self.amp * np.cos(self.phi)
+
+    def parameters(self):
+        """Returns an iterator over the model parameters."""
+        # proper node parameters
+        yield self.des_amp
+        yield self.des_offset
+        yield self.des_freq
+
+        # coupling parameters
+        for node in self.nodes:
+            yield self.nodes[node]['weight']
+            yield self.nodes[node]['bias']
+
+    def named_parameters(self):
+        """Returns an iterator over the model parameters, yielding both the name and the parameter itself"""
+        # proper node parameters
+        yield str(self) + ':amplitude', self.des_amp
+        yield str(self) + ':offset', self.des_offset
+        yield str(self) + ':frequency', self.des_freq
+
+        # coupling parameters
+        for node in self.nodes:
+            yield str(self) + ':weight:' + str(node), self.nodes[node]['weight']
+            yield str(self) + ':bias:' + str(node), self.nodes[node]['bias']
+
+    def list_parameters(self):
+        """Return a list of parameters"""
+        return list(self.parameters())
+
+    def get_vectorized_parameters(self, to_numpy=True):
+        """Return a vectorized form of the parameters (weights, biases, offsets)."""
+        if to_numpy:
+            return np.array(self.list_parameters())
+        else:
+            return torch.from_numpy(np.array(self.list_parameters()))
+
+    def set_vectorized_parameters(self, vector):
+        """Set the vector parameters."""
+        # set the parameters from the vectorized one
+        if len(vector) != self.num_parameters:
+            raise ValueError("Expecting the size of the vectorized parameters to match the number of parameters "
+                             "of this node. Instead of having {}, I got {}.".format(self.num_parameters, len(vector)))
+        self.des_amp = vector[0]
+        self.des_offset = vector[1]
+        self.des_freq = vector[2]
+
+        # coupling parameters
+        for idx, node in enumerate(self.nodes):
+            self.nodes[node]['weight'] = vector[2*idx + 3]
+            self.nodes[node]['bias'] = vector[2*idx + 4]
+
+    def step(self):
+        """
+        Perform a step by integrating (i.e. Euler integration) the differential equations governing the CPG.
+
+        The CPG equations for node :math:`i` are:
+
+        .. math::
+
+            \dot{\phi}_i &= \omega_i + \sum_j a_j w_{ij} \sin(\phi_j - \phi_i - \varphi_{ij}) \\
+            \ddot{a}_i &= K_a (A_i - a_i) - D_a \dot{a}_i \\
+            \ddot{o}_i &= K_o (O_i - o_i) - D_o \dot{o}_i \\
+            \theta_i &= o_i + a_i \cos(\phi_i) \\
+
+        where
+        :math:`\phi` is the phase,
+        :math:`\omega` is the desired angular velocity (desired frequency),
+        :math:`A` and :math:`a` are the desired and current amplitude,
+        :math:`O` and :math:`o` are the desired and current offset,
+        :math:`K` and :math:`D` are the stiffness and damping gains (which are normally related such that the system
+        is critically damped),
+        :math:`w_{ij}` and :math:`\varphi_{ij}` are the coupling weights and phase biases, and finally,
+        :math:`\theta` is the resulting (joint) angle (to be sent to the controller).
+
+        .. note:: for each node, update() has to be called only after step() has been called for all the nodes!
+        """
+
+        # offset
+        ddoffset = self.K_offset * (self.des_offset - self.offset) - self.D_offset * self.doffset
+        self.doffset += ddoffset * self.dt
+        self.curr_offset += self.doffset * self.dt
+
+        # amplitude
+        ddamp = self.K_amp * (self.des_amp - self.amp) - self.D_amp * self.damp
+        self.damp += ddamp * self.dt
+        self.curr_amp += self.damp * self.dt
+
+        # phase
+        dphi = self.des_omega
+        for node, coupling in self.nodes.items():
+            dphi += node.amp * coupling['weight'] * np.sin(node.phi - self.phi - coupling['bias'])
+        self.curr_phi = (self.curr_phi + dphi * self.dt) % (2*np.pi)
+
+        # self.t += self.dt
+
+    def update(self):
+        curr_t = self.t + self.dt
+        curr_theta = self.offset + self.amp * np.cos(self.phi)
+
+        # update plot
+        if self.fig is not None:
+            if self.do_plot_theta:
+                self.ax_plot.plot([self.t, curr_t], [self.theta, curr_theta], 'b')  # , animated=True)
+            if self.do_plot_amp:
+                self.ax_plot.plot([self.t, curr_t], [self.amp, self.curr_amp], 'g')  # , animated=True)
+            if self.do_plot_offset:
+                self.ax_plot.plot([self.t, curr_t], [self.offset, self.curr_offset], 'r')  # , animated=True)
+            if self.do_plot_phi:
+                self.ax_plot.plot([self.t, curr_t], [self.phi, self.curr_phi], 'm')  # , animated=True)
+
+            # move axis to get a feeling of online plotting. Warning: this slows down rendering
+            # self.ax_plot.set_xlim(curr_t - 0.5, curr_t + 0.5)
+
+        # update state params
+        self.phi = self.curr_phi
+        self.amp = self.curr_amp
+        self.offset = self.curr_offset
+        self.theta = curr_theta
+        self.t = curr_t
+
+        # update data
+        # if len(self.phis) >= self.timesteps: self.phis = self.phis[1:]
+        # if len(self.amps) >= self.timesteps: self.amps = self.amps[1:]
+        # if len(self.offsets) >= self.timesteps: self.offsets = self.offsets[1:]
+        # if len(self.thetas) >= self.timesteps: self.thetas = self.thetas[1:]
+        # self.phis.append(self.phi)
+        # self.amps.append(self.amp)
+        # self.offsets.append(self.offset)
+        # self.thetas.append(self.theta)
+        # self.updated = True
+
+        # update sliders
+        if self.fig is not None:
+            self.slider_offset.set_val(self.offset)
+            self.slider_amp.set_val(self.amp)
+            self.slider_phi.set_val(self.phi)
+
+        # update canvas
+        # self.fig.canvas.draw()
+        # plt.show(block=False)
+
+    def plot(self):
+        if self.fig is None:
+            self.fig = plt.figure('CPG '+str(self.id))
+
+        # Plot signals #
+        h = 0.65
+        self.ax_plot = self.fig.add_axes([0.1, h, 0.8, 0.3])
+        self.ax_plot.set_ylim(-np.pi, np.pi)
+        self.ax_plot.set_xlim(0, 1)
+        ax = self.fig.add_axes([0.91, h+0.06, 0.08, 0.18])
+        # ax.axis('off')
+        self.check_offset = CheckButtons(ax, ['o', 'a', 'phi', 'theta'],
+                                         [self.do_plot_offset, self.do_plot_amp, self.do_plot_phi, self.do_plot_theta])
+
+        def check_buttons(val):
+            if val == 'o':
+                self.do_plot_offset = not self.do_plot_offset
+                # self.offset_line.set_visible(self.do_plot_offset)
+            elif val == 'a':
+                self.do_plot_amp = not self.do_plot_amp
+                # self.amp_line.set_visible(self.do_plot_amp)
+            elif val == 'phi':
+                self.do_plot_phi = not self.do_plot_phi
+                # self.phi_line.set_visible(self.do_plot_phi)
+            elif val == 'theta':
+                self.do_plot_theta = not self.do_plot_theta
+                # self.theta_line.set_visible(self.do_plot_theta)
+        self.check_offset.on_clicked(check_buttons)
+
+        # x = np.linspace(0., 1., self.timesteps)
+        # self.theta_line = self.ax_plot.plot(x, self.thetas)[0]
+        # self.amp_line = self.ax_plot.plot(x, self.amps)[0]
+        # self.offset_line = self.ax_plot.plot(x, self.offsets)[0]
+        # self.phi_line = self.ax_plot.plot(x, self.phis)[0]
+        #
+        # self.theta_line.set_visible(self.do_plot_theta)
+        # self.amp_line.set_visible(self.do_plot_amp)
+        # self.offset_line.set_visible(self.do_plot_offset)
+        # self.phi_line.set_visible(self.do_plot_phi)
+        #
+        # def update_plot(_):
+        #     lines = []
+        #     if self.updated:
+        #         if self.do_plot_theta:
+        #             self.theta_line.set_ydata(self.thetas)
+        #             lines.append(self.theta_line)
+        #         if self.do_plot_amp:
+        #             self.amp_line.set_ydata(self.amps)
+        #             lines.append(self.amp_line)
+        #         if self.do_plot_offset:
+        #             self.offset_line.set_ydata(self.offsets)
+        #             lines.append(self.offset_line)
+        #         if self.do_plot_phi:
+        #             self.phi_line.set_ydata(self.phis)
+        #             lines.append(self.phi_line)
+        #         self.updated = False
+        #     return lines
+        #
+        # self.anim = FuncAnimation(self.fig, update_plot, interval=100, blit=True)
+
+        # State parameters #
+        # offset
+        h -= 0.1
+        # ax = self.fig.add_axes([0.2, h, 0.6, 0.03])
+        ax = self.fig.add_axes([0.1, h, 0.35, 0.03])
+        self.slider_offset = Slider(ax, 'o', -np.pi, np.pi, valinit=self.offset, dragging=False)
+
+        def set_offset(val):
+            self.offset = self.slider_offset.val
+
+        self.slider_offset.on_changed(set_offset)
+
+        # amplitude
+        # h -= 0.05
+        # ax = self.fig.add_axes([0.2, h, 0.6, 0.03])
+        ax = self.fig.add_axes([0.55, h, 0.35, 0.03])
+        self.slider_amp = Slider(ax, 'a', 0, np.pi, valinit=self.amp, dragging=False)
+
+        def set_amp(val):
+            self.amp = self.slider_amp.val
+
+        self.slider_amp.on_changed(set_amp)
+
+        # phi
+        h -= 0.05
+        # ax = self.fig.add_axes([0.2, h, 0.6, 0.03])
+        ax = self.fig.add_axes([0.1, h, 0.35, 0.03])
+        self.slider_phi = Slider(ax, 'phi', 0, 2*np.pi, valinit=self.phi, dragging=False)
+
+        def set_phi(val):
+            self.phi = self.slider_phi.val
+
+        self.slider_phi.on_changed(set_phi)
+
+        # Control parameters #
+        # desired offset
+        h -= 0.05
+        # ax = self.fig.add_axes([0.2, h, 0.6, 0.03])
+        ax = self.fig.add_axes([0.1, h, 0.35, 0.03])
+        self.slider_des_offset = Slider(ax, 'O', -np.pi, np.pi, valinit=self.des_offset)
+
+        def set_des_offset(val):
+            self.des_offset = self.slider_des_offset.val
+
+        self.slider_des_offset.on_changed(set_des_offset)
+
+        # desired amplitude
+        # h -= 0.05
+        # ax = self.fig.add_axes([0.2, h, 0.6, 0.03])
+        ax = self.fig.add_axes([0.55, h, 0.35, 0.03])
+        self.slider_des_amp = Slider(ax, 'A', 0, np.pi, valinit=self.des_amp)
+
+        def set_des_amp(val):
+            self.des_amp = self.slider_des_amp.val
+
+        self.slider_des_amp.on_changed(set_des_amp)
+
+        # desired frequency
+        h -= 0.05
+        # ax = self.fig.add_axes([0.2, h, 0.6, 0.03])
+        ax = self.fig.add_axes([0.1, h, 0.35, 0.03])
+        self.slider_des_freq = Slider(ax, 'f', 0, 50, valinit=self.des_freq, valfmt='%0.0f')
+
+        def set_des_freq(val):
+            self.des_freq = self.slider_des_freq.val
+            self.des_omega = 2 * np.pi * self.des_freq
+
+        self.slider_des_freq.on_changed(set_des_freq)
+
+        # Coupling parameters #
+        class Coupling(object):
+            """
+            Class to be used with sliders for the coupling parameters (weights and biases)
+            """
+            def __init__(self, node, nodes, sliderNodes):
+                self.node = node
+                self.nodes = nodes
+                self.sliderNodes = sliderNodes
+
+            def set_weight(self, val):
+                self.nodes[self.node]['weight'] = self.sliderNodes[self.node]['weight'].val
+
+            def set_bias(self, val):
+                self.nodes[self.node]['bias'] = self.sliderNodes[self.node]['bias'].val
+
+        for node, coupling in self.nodes.items():
+            h -= 0.05
+            # slider for coupling weight
+            ax = self.fig.add_axes([0.1, h, 0.35, 0.03])
+            slider = Slider(ax, 'w'+str(self.id)+str(node.id), -5, 5, valinit=coupling['weight'], valfmt='%0.1f')
+            c = Coupling(node, self.nodes, self.sliderNodes)
+            slider.on_changed(c.set_weight)
+            self.sliderNodes[node] = {'weight': slider}
+            # slider for phase bias
+            ax = self.fig.add_axes([0.55, h, 0.35, 0.03])
+            slider = Slider(ax, 'b' + str(self.id) + str(node.id), -np.pi, np.pi, valinit=coupling['bias'], valfmt='%0.2f')
+            slider.on_changed(c.set_bias)
+            self.sliderNodes[node]['bias'] = slider
+
+        plt.show(block=False)
+        # plt.draw()
+        # self.fig.show()
+        # self.fig.canvas.draw() # Problem: not reactive to keyboard/mouse when moving sliders...
+
+    #############
+    # Operators #
+    #############
+
+    def __str__(self):
+        """Return a string describing the class."""
+        return self.__class__.__name__ + "(" + str(self.id) + ")"
+
+
+class CPGNetwork(object):
+    r"""Central Pattern Generator Network
+
+    """
+
+    def __init__(self, nodes, timesteps=100):
+        nodes = nodes if isinstance(nodes, dict) else self.fully_connected_network(nodes)
+
+        # create each node
+        self.nodes, init_params = {}, {'phi', 'offset', 'amplitude', 'freq'}
+        for node_id in nodes.keys():
+            d = {key: val for key, val in nodes[node_id].items() if key in init_params}
+            self.nodes[node_id] = CPGNode(node_id, timesteps=timesteps, **d)
+
+        # couple the nodes
+        for node_id, node in self.nodes.items():
+            if 'nodes' in nodes[node_id]:  # check if 'nodes' in dict
+                for coupling_node in nodes[node_id]['nodes']:  # coupling_node = {'id':..., 'weight':..., 'bias': ...}
+                    w = coupling_node['weight'] if 'weight' in coupling_node else 0.
+                    b = coupling_node['bias'] if 'bias' in coupling_node else 0.
+                    node.add_node(self.nodes[coupling_node['id']], weight=w, bias=b)
+
+        self.fig = None
+        self.nodes_id = self.nodes.keys()
+        self.nodes_id.sort()
+
+    @property
+    def num_parameters(self):
+        """Return the total number of parameters in this CPG network."""
+        return sum([node.num_parameters for node in self.nodes.values()])
+
+    def parameters(self):
+        """Returns an iterator over the model parameters."""
+        for node in self.nodes.values():
+            yield np.array(list(node.parameters()))
+
+    def named_parameters(self):
+        """Returns an iterator over the model parameters, yielding both the name and the parameter itself"""
+        for node in self.nodes.values():
+            yield str(node), np.array(list(node.parameters()))
+
+    def list_parameters(self):
+        """Return a list of parameters"""
+        return list(self.parameters())
+
+    def get_vectorized_parameters(self, to_numpy=True):
+        """Return a vectorized form of the parameters (amplitudes, offsets, frequencies, weights, biases)."""
+        return np.concatenate([node.get_vectorized_parameters(to_numpy=to_numpy) for node in self.nodes.values()])
+
+    def set_vectorized_parameters(self, vector):
+        """Set the vector parameters."""
+        # set the parameters from the vectorized one
+        if len(vector) != self.num_parameters:
+            raise ValueError("Expecting the size of the vectorized parameters to match the number of parameters "
+                             "of this node. Instead of having {}, I got {}.".format(self.num_parameters, len(vector)))
+
+        # convert from torch tensor to numpy array if necessary
+        if isinstance(vector, torch.Tensor):
+            if vector.requires_grad:
+                vector = vector.detach().numpy()
+            else:
+                vector = vector.numpy()
+
+        # set the parameters from the vectorized one
+        idx = 0
+        for node in self.nodes.values():
+            size = node.num_parameters
+            node.set_vectorized_parameters(vector[idx:idx+size])
+            idx += size
+
+    def add_node(self, node):
+        """Add a node to the CPG network."""
+        # TODO: update self.nodes_id
+        self.nodes[node.id] = node
+
+    def reset(self):
+        """Reset the phase of each CPG node in the network; this can be useful for phase resetting."""
+        for node in self.nodes.values():
+            node.reset()
+
+    def step(self):
+        """Perform a step with the CPG network."""
+        # perform one step
+        for node in self.nodes.values():
+            node.step()
+        # perform update
+        for node in self.nodes.values():
+            node.update()
+        return np.array([self.nodes[idx].theta for idx in self.nodes_id])
+
+    def plot(self):
+        """Plot each node."""
+        for node in self.nodes.values():
+            node.plot()
+
+    @staticmethod
+    def fully_connected_network(num_nodes, init_phi=0, offset=0, amplitude=1., weight=0, bias=0, freq=1.):
+        """
+        Create an initial dictionary describing a fully connected network of CPG nodes.
+
+        Args:
+            num_nodes (int): the total number of nodes in the network
+            init_phi (float): initial phase of CPG node
+            offset (float): desired and initial offset of CPG node
+            amplitude (float): desired and initial amplitude of CPG node
+            weight (float): weight between 2 nodes
+            bias (float): phase bias between 2 nodes
+
+        Returns:
+            dict: dictionary describing how the nodes are connected
+        """
+        node_ids = range(1, num_nodes+1)
+        d = {i: {'phi': init_phi,
+                 'offset': offset,
+                 'amplitude': amplitude,
+                 'freq': freq,
+                 'nodes': [{'id': n, 'weight': weight, 'bias': bias} for n in (node_ids[:i-1]+node_ids[i:])]}
+             for i in node_ids}
+        return d
+
+
+# Tests
+if __name__ == "__main__":
+    # Define and parse command line arguments
+    import argparse
+    parser = argparse.ArgumentParser()
+    parser.add_argument('-t', '--test', help='What to test', type=str,
+                        choices=['1_node', '2_nodes', 'coman_gazebo', 'minitaur_pybullet'],
+                        default='minitaur_pybullet')
+    parser.add_argument('-T', '--nb_steps', help='Total time number of steps to run the simulation',
+                        type=int, default=100)
+    parser.add_argument('-dt', '--dt', help='Time to wait before next step', type=float, default=0.01)
+    args = parser.parse_args()
+
+    dt = args.dt
+    T = args.nb_steps
+
+    # Check a single node #
+    if args.test == '1_node':
+        # create CPG node
+        node = CPGNode(1)
+        node.plot()
+
+        # Run for T steps
+        for _ in range(T):
+            node.step()
+            node.update()
+            # time.sleep(0.1) # Don't use time, instead use plt.pause!
+            plt.pause(dt)
+        plt.show()
+
+    # Check 2 nodes #
+    elif args.test == '2_nodes':
+        # create CPG network
+        nodes = {1: {'phi': -np.pi / 2, 'nodes': [{'id': 2, 'weight': 0, 'bias': 0}]},
+                 2: {'phi': np.pi / 2, 'nodes': [{'id': 1, 'weight': 0, 'bias': 0}]}}
+        network = CPGNetwork(nodes)
+        network.plot()
+
+        # Run for T steps
+        for _ in range(T):
+            network.step()
+            plt.pause(dt)
+        plt.show()
+
+    # Check coman robot in gazebo #
+    # Warning: Don't forget to launch gazebo with coman before running this code!
+    # $ roslaunch coman_gazebo coman_world.launch
+    elif args.test == 'coman_gazebo':
+        from pyrobolearn.robots.ros.coman.comanpublisher import ComanPublisher
+        import rospy
+
+        # create CPG network of 2 nodes
+        pub = ComanPublisher(joints=['RHipSag', 'LHipSag'])
+        nodes = {1: {'phi': -np.pi/2, 'nodes': [{'id': 2, 'weight': 1, 'bias': 0}]},
+                 2: {'phi': np.pi/2, 'nodes': [{'id': 1, 'weight': 0, 'bias': 0}]}}
+        network = CPGNetwork(nodes)
+        # network.plot()
+
+        # Run for T steps
+        # rate = rospy.Rate(30)
+        # while not rospy.is_shutdown():
+        for _ in range(T):
+            network.step()
+            pub.send({'RHipSag': network.nodes[1].theta,
+                     'LHipSag': network.nodes[2].theta})
+            plt.pause(dt)
+            # rate.sleep()
+
+    # Check CPG with Minitaur #
+    elif args.test == 'minitaur_pybullet':
+        from pybullet_envs.bullet.bullet_client import BulletClient
+        from pybullet_envs.bullet.minitaur import Minitaur
+        import pybullet
+        import pybullet_data
+        import time
+
+        # create and configure pybullet simulator
+        client = BulletClient(connection_mode=pybullet.GUI)
+        client.setAdditionalSearchPath(pybullet_data.getDataPath())
+        floor = client.loadURDF('plane.urdf')
+        client.setGravity(0, 0, -9.81)
+        minitaur = Minitaur(client, urdf_root=pybullet_data.getDataPath())
+
+        # create CPG network
+        phis = [0, 0, 0, 0]
+        nodes = {1: {'phi': phis[0], 'offset': -np.pi/2, 'amplitude': np.pi/4, 'nodes': []},  # front left leg
+                 2: {'phi': phis[1], 'offset': -np.pi/2, 'amplitude': np.pi/4, 'nodes': []},  # back left leg
+                 3: {'phi': phis[2], 'offset': np.pi/2, 'amplitude': np.pi/4, 'nodes': []},  # front right leg
+                 4: {'phi': phis[3], 'offset': np.pi/2, 'amplitude': np.pi/4, 'nodes': []},  # back right leg
+                 }
+        network = CPGNetwork(nodes)
+
+        # Run simulation
+        for _ in range(10*T):
+            # Get angles from CPG network and set them to the minitaur
+            act = network.step()
+            for i in range(len(act)//2): # Left
+                minitaur._SetDesiredMotorAngleById(minitaur._motor_id_list[2*i], act[i])
+                minitaur._SetDesiredMotorAngleById(minitaur._motor_id_list[2*i+1], -np.pi-act[i])
+            for i in range(len(act)//2, len(act)): # Right
+                minitaur._SetDesiredMotorAngleById(minitaur._motor_id_list[2*i], act[i])
+                minitaur._SetDesiredMotorAngleById(minitaur._motor_id_list[2*i+1], np.pi - act[i])
+
+            # run one-step forward the simulation
+            client.stepSimulation()
+            time.sleep(dt)
+            # print(client.getEulerFromQuaternion(minitaur.GetBaseOrientation()))

pyrobolearn/models/gp.py

@@ -0,0 +1,649 @@
+#!/usr/bin/env python
+"""Define the Gaussian Process (GP) learning model.
+
+This file provides the Gaussian Process (GP) model; a non-parametric, discriminative, and probabilistic model.
+
+As for neural networks, several frameworks can be used, such as `GPy` (which uses `numpy`), `GPyTorch` (which uses
+`pytorch`), or `GPFlow` (which uses `tensorflow`). We decided to use `GPyTorch` because of the `pytorch` framework
+popularity in the research community field, its flexibility, its similarity with numpy (but with automatic
+differentiation: autograd), GPU capabilities, and more Pythonic approach.
+"""
+
+import copy
+
+try:
+    import cPickle as pickle
+except ImportError as e:
+    import pickle
+
+import numpy as np
+import torch
+import gpytorch
+# import GPy
+# from model import Model
+
+__author__ = "Brian Delhaisse"
+__copyright__ = "Copyright 2018, PyRoboLearn"
+__credits__ = ["Brian Delhaisse"]
+__license__ = "MIT"
+__version__ = "1.0.0"
+__maintainer__ = "Brian Delhaisse"
+__email__ = "briandelhaisse@gmail.com"
+__status__ = "Development"
+
+
+class GP(object):
+    r"""Gaussian Process model
+
+    This is a wrapper around the GPyTorch gaussian process models. The Gaussian process is a generalization of
+    the multivariate Gaussian distribution. It is a non-parametric, probabilistic, and discriminative model.
+
+    This works by putting a prior distribution on the function:
+
+    .. math:: f ~ GP(0, K(X,X))
+
+    where :math:`K(.,.)` is the kernel matrix where each entry contains :math:`k(x_i, x_j)`, i.e. the kernel function
+    evaluated at the corresponding points. As it can be seen the kernel matrix grows with the number of samples.
+
+    See Also:
+        - `GPy` (which uses numpy) [4]
+        - `GPFlow` (which uses TensorFlow) [5]
+
+    References:
+        [1] "Gaussian Processes for Machine Learning", Rasmussen and Williams, 2006
+        [2] GPyTorch: https://github.com/cornellius-gp/gpytorch
+        [3] GPyTorch examples: https://github.com/cornellius-gp/gpytorch/tree/master/examples
+        [4] GPy: https://gpy.readthedocs.io/en/deploy/
+        [5] GPFlow: http://gpflow.readthedocs.io/en/latest/intro.html
+    """
+    pass
+
+
+class GPC(GP):
+    r"""Gaussian Process Classification
+
+    References:
+        [1] "Gaussian Processes for Machine Learning", Rasmussen and Williams, 2006
+        [2] GPyTorch: https://github.com/cornellius-gp/gpytorch
+        [3] GPyTorch examples: https://github.com/cornellius-gp/gpytorch/tree/master/examples
+    """
+    pass
+
+
+class ExactGPModel(gpytorch.models.ExactGP):
+    r"""Create GP prior model.
+    """
+
+    def __init__(self, mean, kernel, likelihood=None, x=None, y=None):
+        # create likelihood if not already set.
+        if likelihood is None:
+            likelihood = gpytorch.likelihoods.GaussianLikelihood()
+        super(ExactGPModel, self).__init__(train_inputs=x, train_targets=y, likelihood=likelihood)
+
+        # set mean and kernel covariance function
+        self.mean = mean
+        self.kernel = kernel
+
+    @property
+    def mean(self):
+        r"""Return the GP prior mean; that is :math:`\mu(x)` from :math:`p(f|x) = N(\mu(x), K(x,x))`."""
+        return self._mean
+
+    @mean.setter
+    def mean(self, mean):
+        r"""Set the GP prior mean; that is :math:`\mu(x)` from :math:`p(f|x) = \mathcal{N}(\mu(x), K(x,x))`."""
+        if mean is None:
+            mean = gpytorch.means.ConstantMean()
+        if not isinstance(mean, gpytorch.means.Mean):
+            raise TypeError("Expecting the mean to be an instance of `gpytorch.means.Mean`, got instead "
+                            "{}".format(type(mean)))
+        self._mean = mean
+
+    @property
+    def kernel(self):
+        r"""Return the kernel function :math:`K(.,.)` (=prior covariance of the GP)."""
+        return self._kernel
+
+    @kernel.setter
+    def kernel(self, kernel):
+        r"""Set the kernel function :math:`K(.,.)` (=prior covariance of the GP)."""
+        if kernel is None:
+            kernel = gpytorch.kernels.RBFKernel()  # + gpytorch.kernels.WhiteNoiseKernel()
+            kernel = gpytorch.kernels.ScaleKernel(kernel)
+        if not isinstance(kernel, gpytorch.kernels.Kernel):
+            raise TypeError("Expecting the kernel to be an instance of `gpytorch.kernels.Kernel`, got instead "
+                            "{}".format(type(kernel)))
+        self._kernel = kernel
+
+    def forward(self, x):
+        r"""Return the prior probability density function :math:`p(f|x) = \mathcal{N}(. | \mu(x), K(x,x))`."""
+        mean_x = self.mean(x)
+        covar_x = self.kernel(x)
+        # return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)
+        return gpytorch.random_variables.GaussianRandomVariable(mean_x, covar_x)
+
+
+class GPR(GP):
+    r"""Gaussian Process Regression
+
+    The Gaussian process is a generalization of the multivariate Gaussian distribution. It is a non-parametric,
+    probabilistic, and discriminative model.
+
+    This works by putting a prior distribution on the function:
+
+    .. math:: f|X ~ GP(0, K(X,X))
+
+    where :math:`K(.,.)` is the kernel matrix where each entry contains :math:`k(x_i, x_j)`, i.e. the kernel function
+    evaluated at the corresponding points. As it can be seen the kernel matrix grows with the number of samples.
+
+    The log likelihood is given by:
+
+    .. math:: \log p(y | f) = \mathcal{N}(y | 0, K(X,X) + \sigma I)
+
+    Learning the hyperparameters of the kernel are carried out by maximizing the marginal log likelihood, which is
+    given by:
+
+    .. math:: \log p(y | X) = \int p(y | f) p(f | X) df
+
+    The predictive distribution is carried out by assuming that the observed target values :math:`y` and
+    the function values :math:`f^*` at the test locations :math:`X^*` are from the same joint Gaussian distribution.
+    By conditioning this distribution with respect to the old dataset :math:`X, y` and the test locations :math:`X^*`,
+    we can derive :math:`p(f^* | X, y, X^*)` which is the predictive output distribution given the new data points
+    :math:`X^*`.
+
+    Notes:
+        * GP takes into account correlations in the input domain but not in the output space
+        * The time complexity to learn a GP is :math:`O(N^3)` because of the matrix inversion during training.
+        * GMM vs GP:
+            * Both are probabilistic models.
+            * GMM is a generative semi-parametric model while GP is a discriminative non-parametric model.
+            * GMM captures the correlation between the inputs and outputs, while GP only captures correlation in the
+                input space.
+            * GMR models the variability/correlation between the predicted outputs while the GP provides uncertainty
+                on the predicted outputs. The predicted outputs in a GP are independent unless using a heteroscedastic
+                GP or a generalized Wishart process is used.
+
+    GPyTorch::
+
+    For most GP regression models, you will need to construct the following GPyTorch objects:
+    1. A GP Model (`gpytorch.models.ExactGP`) - This handles most of the inference.
+    2. A Likelihood (`gpytorch.likelihoods.GaussianLikelihood`) - This is the most common likelihood used for GP
+        regression.
+    3. A Mean - This defines the prior mean of the GP. If you don't know which mean to use, a
+        `gpytorch.means.ConstantMean` is a good place to start.
+    4. A Kernel - This defines the prior covariance of the GP. If you don't know which kernel to use, a
+        `gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel())` is a good place to start.
+    5. A MultivariateNormal Distribution (`gpytorch.distributions.MultivariateNormal`) - This is the object used to
+        represent multivariate normal distributions.
+
+
+    References:
+        [1] "Gaussian Processes for Machine Learning", Rasmussen and Williams, 2006
+        [2] GPy: https://gpy.readthedocs.io/en/deploy/
+        [3] GPyTorch: https://github.com/cornellius-gp/gpytorch
+        [4] GPFlow: http://gpflow.readthedocs.io/en/latest/intro.html
+    """
+
+    def __init__(self,  mean=None, kernel=None, model=None, likelihood=None):
+        """
+        Initialize the GPR.
+
+        Args:
+            mean (None, gpytorch.means.Mean): mean prior. If None, it will be set to `gpytorch.means.ConstantMean()`.
+            kernel (None, gpytorch.kernels.Kernel): kernel prior. If None it will be set to
+                `gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel() + gpytorch.kernels.WhiteNoiseKernel())`
+            model (None, gpytorch.module.Module): the prior GP model. If None, it will create `ExactGPModel()`, a GP
+                model using the provided mean, kernel, and likelihood.
+            likelihood (None, gpytorch.likelihoods.Likelihood): the likelihood pdf. If None, it will use the
+                `gpytorch.likelihoods.GaussianLikelihood()`
+        """
+        # check model
+        if model is None:
+            self.model = ExactGPModel(mean, kernel, likelihood)
+        else:
+            self.model = model
+
+        # set model into evaluation mode
+        self.eval()
+
+        # set the marginal log likelihood pdf: p(y|x) = \int p(y|f,x) p(f|x) df
+        self.mll = gpytorch.mlls.ExactMarginalLogLikelihood(self.likelihood_prob, self.prior)
+
+    ##############
+    # Properties #
+    ##############
+
+    @property
+    def model(self):
+        """Return the GP model; i.e. the prior probability density function p(f|x)."""
+        return self._model
+
+    @model.setter
+    def model(self, model):
+        """Set the GP model; i.e. the prior probability density function p(f|x)."""
+        if not isinstance(model, gpytorch.module.Module):
+            raise TypeError("Expecting the GP model to be an instance of `gpytorch.module.Module`, got instead "
+                            "{}".format(type(model)))
+        self._model = model
+
+    # alias
+    prior = model
+
+    @property
+    def likelihood_prob(self):
+        """Return the likelihood probability density function p(y|f,x)."""
+        return self.model.likelihood
+
+    @likelihood_prob.setter
+    def likelihood_prob(self, likelihood):
+        r"""Set the likelihood probability density function p(y|f,x)."""
+        if likelihood is None:
+            likelihood = gpytorch.likelihoods.GaussianLikelihood()
+        if not isinstance(likelihood, gpytorch.likelihoods.Likelihood):
+            raise TypeError("Expecting the likelihood to be an instance of `gpytorch.likelihood.Likelihood`, got "
+                            "instead {}".format(type(likelihood)))
+        self.model.likelihood = likelihood
+
+    @property
+    def log_marginal_likelihood_prob(self):
+        r"""Return the log marginal log likelihood pdf: log p(y|x)."""
+        return self.mll
+
+    @property
+    def mean(self):
+        r"""Return the GP prior mean; that is :math:`\mu(x)` from :math:`p(f|x) = N(\mu(x), K(x,x))`."""
+        return self.model.mean
+
+    @mean.setter
+    def mean(self, mean):
+        r"""Set the GP prior mean; that is :math:`\mu(x)` from :math:`p(f|x) = \mathcal{N}(\mu(x), K(x,x))`."""
+        if mean is None:
+            mean = gpytorch.means.ConstantMean()
+        if not isinstance(mean, gpytorch.means.Mean):
+            raise TypeError("Expecting the mean to be an instance of `gpytorch.means.Mean`, got instead "
+                            "{}".format(type(mean)))
+        self.model.mean = mean
+
+    @property
+    def kernel(self):
+        r"""Return the kernel function :math:`K(.,.)` (=prior covariance of the GP)."""
+        return self.model.kernel
+
+    @kernel.setter
+    def kernel(self, kernel):
+        r"""Set the kernel function :math:`K(.,.)` (=prior covariance of the GP)."""
+        if kernel is None:
+            kernel = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel() + gpytorch.kernels.WhiteNoiseKernel())
+        if not isinstance(kernel, gpytorch.kernels.Kernel):
+            raise TypeError("Expecting the kernel to be an instance of `gpytorch.kernels.Kernel`, got instead "
+                            "{}".format(type(kernel)))
+        self.model.kernel = kernel
+
+    @property
+    def dim(self):
+        """Return the dimension of the kernel"""
+        return self.kernel.dim
+
+    ##################
+    # Static Methods #
+    ##################
+
+    @staticmethod
+    def copy(other, deep=True):
+        """copy the other GP"""
+        if not isinstance(other, GPR):
+            raise TypeError("Expecting a GPR model.")
+        if deep:
+            return copy.deepcopy(other)
+        return copy.copy(other)
+
+    @staticmethod
+    def is_parametric():
+        """The Gaussian process is a non parametric model."""
+        return False
+
+    @staticmethod
+    def is_linear():
+        """The Gaussian process does not have any parameters and thus no linear parameters"""
+        return False
+
+    @staticmethod
+    def is_recurrent():
+        """The Gaussian process is not a recurrent model"""
+        return False
+
+    @staticmethod
+    def is_probabilistic():
+        """The Gaussian process is a probabilistic model"""
+        return True
+
+    @staticmethod
+    def is_discriminative():
+        """The Gaussian process is a discriminative model which predicts :math:`p(y|x)`"""
+        return True
+
+    @staticmethod
+    def is_generative():
+        """The Gaussian process is not a generative model, and thus we can not sample from it"""
+        # TODO: actually we can sample a function from it given the initial data (in kernel matrix)
+        return False
+
+    @staticmethod
+    def load(filename):
+        """
+        Load a model from memory.
+
+        Args:
+            filename (str): file that contains the model.
+        """
+        return pickle.load(filename)
+
+    @staticmethod
+    def _convert_to_torch(x):
+        if isinstance(x, np.ndarray):
+            return torch.from_numpy(x).float()
+        return x
+
+    @staticmethod
+    def _convert_to_numpy(x):
+        if isinstance(x, torch.Tensor):
+            if x.requires_grad:
+                return x.detach().numpy()
+            return x.numpy()
+        return x
+
+    @staticmethod
+    def _convert(x, to_numpy=True):
+        if to_numpy:
+            return GPR._convert_to_numpy(x)
+        return x
+
+    ###########
+    # Methods #
+    ###########
+
+    def train(self, mode=True):
+        """Set the model into training mode."""
+        if mode:
+            self.model.train(True)
+            self.likelihood_prob.train(True)
+        else:
+            self.model.train(False)
+            self.likelihood_prob.train(False)
+
+    def eval(self):
+        """Set the mode into evaluation mode."""
+        self.model.eval()
+        self.likelihood_prob.eval()
+
+    def hyperparameters(self):
+        """Return an iterator over the hyperparameters"""
+        return self.model.hyperparameters()
+
+    def named_hyperparameters(self):
+        """Return an iterator over the model parameters, yielding both the name and the parameter itself."""
+        return self.model.named_hyperparameters()
+
+    def likelihood(self, x, y, to_numpy=False):
+        r"""Evaluate the likelihood p(y|f,x)."""
+        likelihood = torch.exp(self.log_likelihood(x, y, to_numpy=False))
+        return self._convert(likelihood, to_numpy=to_numpy)
+
+    def log_likelihood(self, x, y, to_numpy=False):
+        r"""Evaluate the log likelihood: log p(y|f,x)."""
+        x = self._convert_to_torch(x)
+        y = self._convert_to_torch(y)
+        f = self.model(x)
+        log_likelihood = self.likelihood_prob.log_probability(f, y)
+        return self._convert(log_likelihood, to_numpy=to_numpy)
+
+    def marginal_likelihood(self, x, y, to_numpy=False):
+        r"""Evaluate the marginal likelihood: p(y|x)."""
+        ml = torch.exp(self.log_marginal_likelihood(x, y, to_numpy=False))
+        return self._convert(ml, to_numpy=to_numpy)
+
+    def log_marginal_likelihood(self, x, y, to_numpy=False):
+        r"""Evaluate the log marginal likelihood: log p(y|x)."""
+        x = self._convert_to_torch(x)
+        y = self._convert_to_torch(y)
+        f = self.model(x)
+        mll = self.mll(f, y)
+        return self._convert(mll[0], to_numpy=to_numpy)
+
+    def fit(self, x, y, num_iters=100, tolerance=1e-5, optimizer=None, verbose=False):
+        r"""Fit the input and output data; find optimal model hyperparameters."""
+        # check input and output data
+        x = self._convert_to_torch(x)
+        y = self._convert_to_torch(y)
+
+        # set training data
+        # self.model.set_train_data(x, y)
+        self.model = ExactGPModel(self.mean, self.kernel, self.likelihood_prob, x, y)
+
+        # set into training mode
+        self.train(mode=True)
+
+        # define optimizer
+        need_closure = False
+        if optimizer is None or (isinstance(optimizer, str) and optimizer.lower() == 'lbfgs'):
+            optimizer = torch.optim.LBFGS([{'params': self.model.parameters()}], max_iter=num_iters,
+                                          tolerance_change=tolerance)
+            need_closure = True
+
+            def closure():
+                optimizer.zero_grad()
+                output = self.model(x)
+                loss = -self.mll(output, y)
+                loss.backward()
+                return loss
+
+        elif isinstance(optimizer, str) and optimizer.lower() == 'adam':
+            optimizer = torch.optim.Adam([{'params': self.model.parameters()}], lr=0.1)
+
+        if not isinstance(optimizer, torch.optim.Optimizer):
+            raise TypeError("Expecting the optimizer to be an instance of `torch.optim.Optimizer`.")
+
+        # optimize
+        for i in range(num_iters):
+            # zero gradients from previous iteration
+            optimizer.zero_grad()
+
+            # predict output from model p(f|x)
+            output = self.model(x)
+
+            # compute loss
+            loss = - self.mll(output, y)
+
+            # call backward on the loss to fill the gradients
+            loss.backward()
+
+            # print info if specified
+            if verbose:
+                # print('Iter %d/%d - Loss: %.3f   lengthscale: %.3f   noise: %.3f' % (i + 1, num_iters, loss.item(),
+                #       self.model.kernel.base_kernel.lengthscale.item(), self.model.likelihood.noise.item()))
+                print('Iter %d/%d - Loss: %.3f' % (i + 1, num_iters, loss.item()))
+
+            # perform a step with the optimizer
+            if need_closure:
+                optimizer.step(closure)
+            else:
+                optimizer.step()
+
+        # set into evaluation mode (=predictive posterior mode)
+        self.eval()
+
+    def predict(self, x, to_numpy=True):
+        r"""
+        Predict mean output array :math:`\mu(y)` given input array :math:`x`. That is, it returns the mean of the
+        predictive posterior distribution :math:`E[p(y|x)]`.
+
+        Args:
+            x (np.ndarray, torch.Tensor): input array
+            to_numpy (bool): if True, return a np.array
+
+        Returns:
+            np.ndarray, torch.Tensor: output mean array
+        """
+        x = self._convert_to_torch(x)
+
+        # compute p(f|x)
+        f = self.model(x)
+
+        # compute p(y|f,x)
+        y = self.likelihood_prob(f)
+
+        # return mean
+        if to_numpy:
+            return self._convert_to_numpy(y.mean())
+        return y.mean()
+
+    def predict_prob(self, x, to_numpy=True):
+        r"""
+        Predict p(y|x) by returning the mean and the covariance arrays.
+
+        Args:
+            x (np.ndarray, torch.Tensor): input array
+            to_numpy (bool): if True, return a np.array
+
+        Returns:
+            np.ndarray, torch.Tensor: output mean array
+            np.ndarray, torch.Tensor: output covariance array
+        """
+        x = self._convert_to_torch(x)
+
+        # compute p(f|x)
+        f = self.model(x)
+
+        # compute p(y|f,x)
+        y = self.likelihood_prob(f)
+
+        # return mean and covariance
+        if to_numpy:
+            return self._convert_to_numpy(y.mean()), self._convert_to_numpy(y.var())  # y.covar())
+        return y.mean(), y.var()  # y.covar()
+
+    def forward(self, x):
+        r"""
+        Return the predictive distribution p(y|x).
+
+        Args:
+            x (np.ndarray, torch.Tensor): input array
+
+        Returns:
+            gpytorch.random_variables.GaussianRandomVariable: multivariate normal (Gaussian) distribution
+        """
+        x = self._convert_to_torch(x)
+
+        # compute p(f|x)
+        f = self.model(x)
+
+        # return p(y|f,x)
+        return self.likelihood_prob(f)
+
+    def sample(self, x, num_samples=1, to_numpy=True):
+        """Sample the function vector from the GP; i.e. f ~ p(f|x)."""
+        x = self._convert_to_torch(x)
+        f = self.model(x)
+        if to_numpy:
+            return self._convert_to_numpy(f.sample(num_samples))
+        return f.sample(num_samples)
+
+    #############
+    # Operators #
+    #############
+
+    def __str__(self):
+        """Return name of class"""
+        return self.__class__.__name__
+
+
+# class GPRTorch(GPR):
+#     r"""Gaussian Process Regression using GPyTorch
+#
+#     This provides a wrapper around the GPyTorch module.
+#
+#     References:
+#         [1] "GPyTorch: Blackbox Matrix-Matrix Gaussian Process Inference with GPU Acceleration", Gardner et al., 2018
+#         [2] GPyTorch: https://github.com/cornellius-gp/gpytorch
+#         [3] GPyTorch Examples: https://github.com/cornellius-gp/gpytorch/tree/master/examples
+#     """
+#
+#     def __init__(self, model):
+#         super(GPRTorch, self).__init__(model)
+#
+#     def _predict(self, x=None):
+#         return self.model(x)
+
+
+# class GPRy(GPR):
+#     r"""Gaussian Process Regression using GPy.
+#
+#     References:
+#         [1] "GPy: A Gaussian process framework in python", Sheffield, 2014
+#         [2] GPy: https://github.com/SheffieldML/GPy
+#     """
+#
+#     def __init__(self, model):
+#         super(GPRy, self).__init__(model)
+#
+#     def _predict(self, x=None, full_cov=False):
+#         return self.model.predict(x, full_cov)[0]
+
+
+# TESTS
+if __name__ == '__main__':
+    import matplotlib.pyplot as plt
+    from utils.converter import torch_to_numpy
+
+    # create input and output data
+    x = torch.linspace(0, 1, 100)
+    y = torch.sin(x * (2 * np.pi)) + torch.randn(x.size()) * 0.2
+
+    x, y = x.numpy(), y.numpy()
+
+    # plot true data
+    plt.plot(x, y, 'x')
+
+    # create GPR
+    model = GPR()
+
+    # plot prior possible functions
+    f = model.sample(x, num_samples=10, to_numpy=True)
+    plt.plot(x, f.T)
+
+    # compute log likelihoods
+    print("\nBefore training:")
+    print("Log likelihood: {}".format(model.log_likelihood(x, y, to_numpy=True)))
+    print("Log marginal likelihood: {}".format(model.log_marginal_likelihood(x, y, to_numpy=True)))
+
+    # fit the data
+    optimizer = 'adam'  # 'lbfgs'
+    model.fit(x, y, num_iters=100, optimizer=optimizer, verbose=True)
+
+    # compute log likelihoods
+    print("\nAfter training:")
+    print("Log likelihood: {}".format(model.log_likelihood(x, y, to_numpy=True)))
+    print("Log marginal likelihood: {}".format(model.log_marginal_likelihood(x, y, to_numpy=True)))
+
+    # sample function and plot it
+    f = model.sample(x, num_samples=1, to_numpy=True)
+    plt.plot(x, f.T, 'k', linewidth=2.)
+
+    # predict prob
+    x_test = torch.linspace(0, 1, 51).numpy()
+    mean_y, var_y = model.predict_prob(x_test, to_numpy=True)
+    std_y = np.sqrt(var_y)
+    plt.plot(x_test, mean_y, 'b')
+    plt.fill_between(x_test, mean_y-2*std_y, mean_y+2*std_y, facecolor='green', alpha=0.5)
+    plt.ylim([-3, 3])
+    plt.show()
+
+    # Another way to predict
+    pred = model.forward(x_test)
+    lower, upper = pred.confidence_region()
+
+    plt.plot(x, y, 'k*')
+    plt.plot(x_test, torch_to_numpy(pred.mean()), 'b')
+    plt.fill_between(x_test, torch_to_numpy(lower), torch_to_numpy(upper), alpha=0.5)
+    plt.ylim([-3, 3])
+    plt.show()

pyrobolearn/models/hmm.py

@@ -0,0 +1,130 @@
+# This file describes the Hidden Markov Model
+
+from gaussian import Gaussian
+from model import Model
+from hmmlearn.hmm import GaussianHMM
+
+
+class HMM(object):
+    r"""Hidden Markov Models
+
+    Description: emission probabilities, transition probabilities,...
+
+    References:
+        [1] "Pattern Recognition and Machine Learning" (chap 13), Bishop, 2006
+
+    The code was inspired by the following codes:
+    * `hmmlearn`: https://github.com/hmmlearn/hmmlearn
+    * `ghmm`: http://ghmm.sourceforge.net/
+    * `pbdlib`: https://gitlab.idiap.ch/rli/pbdlib-python/tree/master/pbdlib
+    """
+
+    def __init__(self, emission_prob=None):
+        if emission_prob is None or (isinstance(emission_prob, str) and emission_prob.lower() == 'gaussian'):
+            emission_prob = Gaussian()
+
+    ##############
+    # Properties #
+    ##############
+
+
+    ##################
+    # Static Methods #
+    ##################
+
+    @staticmethod
+    def copy(other):
+        if not isinstance(other, HMM):
+            raise TypeError("Trying to copy an object which is not a HMM")
+
+    @staticmethod
+    def isParametric():
+        """The HMM is a parametric model"""
+        return True
+
+    @staticmethod
+    def isLinear():
+        """The HMM is a non-linear model"""
+        return True
+
+    @staticmethod
+    def isRecurrent():
+        """The HMM is recurrent; current outputs depends on previous inputs and states"""
+        return True
+
+    @staticmethod
+    def isProbabilistic():
+        """The HMM a probabilistic model"""
+        return False
+
+    @staticmethod
+    def isDiscriminative():
+        """The HMM is a discriminative model"""
+        return False
+
+    @staticmethod
+    def isGenerative():
+        """The HMM is a generative model which models the joint distributions on states and outputs.
+        This means we can sample from it."""
+        return True
+
+    ###########
+    # Methods #
+    ###########
+
+    def likelihood(self):
+        pass
+
+    # alias
+    pdf = likelihood
+
+    def joint_pdf(self, X, Z):
+        pass
+
+    def sample(self, size=None, seed=None):
+        """
+        Sample from the HMM.
+
+        Args:
+            size:
+            seed:
+
+        Returns:
+
+        """
+        pass
+
+    def expectation_step(self):
+        """
+        Expectation step in the expectation-maximization algorithm.
+
+        Returns:
+
+        """
+        pass
+
+    def maximization_step(self):
+        pass
+
+    def expectation_maximization(self, X):
+        """Expectation-Maximization (EM) algorithm"""
+        pass
+
+    def forward_backward(self):
+        """Forward backward algorithm"""
+        pass
+
+    def sum_product(self):
+        """Sum-product algorithm"""
+        pass
+
+    def viterbi(self):
+        """Viterbi algorithm"""
+        pass
+
+
+class HSMM(HMM):
+    r"""Hidden semi-Markov Models
+
+    """
+    pass

pyrobolearn/models/linear.py

@@ -0,0 +1,254 @@
+#!/usr/bin/env python
+"""Provide the Linear Model.
+
+The linear model is a discriminative deterministic model given by: :math:`y = f(x) = w^T x`.
+"""
+
+import copy
+
+try:
+    import cPickle as pickle
+except ImportError as e:
+    import pickle
+
+import numpy as np
+import torch
+# from model import Model
+
+__author__ = "Brian Delhaisse"
+__copyright__ = "Copyright 2018, PyRoboLearn"
+__credits__ = ["Brian Delhaisse"]
+__license__ = "MIT"
+__version__ = "1.0.0"
+__maintainer__ = "Brian Delhaisse"
+__email__ = "briandelhaisse@gmail.com"
+__status__ = "Development"
+
+
+class Linear(object):
+    r"""Linear Model
+
+    This class describes the linear parametric model: :math:`y = W x + b` where :math:`x` and :math:`y` are
+    respectively the input and output vectors, :math:`W` is the weight matrix, and :math:`b` is the bias/intercept.
+    """
+
+    def __init__(self, num_inputs, num_outputs, add_bias=True):
+        r"""
+        Initialize the affine/linear model described mathematically by:
+
+        .. math:: y = W x + b
+
+        Args:
+            num_inputs (int): dimension of the input
+            num_outputs (int):  dimension of the output
+            add_bias (bool): if True, it will add a bias to the prediction
+        """
+        super(Linear, self).__init__()
+        self.model = torch.nn.Linear(num_inputs, num_outputs, bias=add_bias)
+        self._num_parameters = len(self.get_vectorized_parameters())
+
+    ##############
+    # Properties #
+    ##############
+
+    @property
+    def input_dims(self):
+        """Return the input dimension of the model"""
+        return self.model.weight.shape[1]
+
+    @property
+    def output_dims(self):
+        """Return the output dimension of the model"""
+        return self.model.weight.shape[0]
+
+    @property
+    def input_shape(self):
+        """Return the input shape of the model"""
+        return tuple([self.input_dims])
+
+    @property
+    def output_shape(self):
+        """Return the output shape of the model"""
+        return tuple([self.output_dims])
+
+    @property
+    def num_parameters(self):
+        """Return the total number of parameters"""
+        return self._num_parameters
+
+    ##################
+    # Static Methods #
+    ##################
+
+    @staticmethod
+    def copy(other, deep=True):
+        """Return another copy of the linear model"""
+        if not isinstance(other, Linear):
+            raise TypeError("Trying to copy an object which is not a Linear model")
+        if deep:
+            return copy.deepcopy(other)
+        return copy.copy(other)
+
+    @staticmethod
+    def is_parametric():
+        """The linear model is a parametric model"""
+        return True
+
+    @staticmethod
+    def is_linear():
+        """By definition, a linear model is linear"""
+        return True
+
+    @staticmethod
+    def is_recurrent():
+        """The linear model is not recurrent; current outputs do not depend on previous inputs"""
+        return False
+
+    @staticmethod
+    def is_deterministic():
+        """The linear model is a deterministic model"""
+        return True
+
+    @staticmethod
+    def is_probabilistic():
+        """The linear model is not a probabilistic model; it is a deterministic one"""
+        return False
+
+    @staticmethod
+    def is_discriminative():
+        """The linear model is a discriminative model which predicts :math:`y = Wx + b` where :math:`x` is the input,
+        and :math:`y` is the output"""
+        return True
+
+    @staticmethod
+    def is_generative():
+        """The linear model is not a generative model, and thus we can not sample from it"""
+        return False
+
+    @staticmethod
+    def load(filename):
+        """
+        Load a model from memory.
+
+        Args:
+            filename (str): file that contains the model.
+        """
+        return pickle.load(filename)
+
+    ###########
+    # Methods #
+    ###########
+
+    def parameters(self):
+        """Returns an iterator over the model parameters."""
+        return self.model.parameters()
+
+    def named_parameters(self):
+        """Returns an iterator over the model parameters, yielding both the name and the parameter itself"""
+        return self.model.named_parameters()
+
+    def list_parameters(self):
+        """Return a list of parameters"""
+        return list(self.parameters())
+
+    def get_vectorized_parameters(self, to_numpy=True):
+        """Return a vectorized form (1 dimensional array) of the parameters."""
+        parameters = self.parameters()
+        vector = torch.cat([parameter.reshape(-1) for parameter in parameters])  # np.concatenate = torch.cat
+        if to_numpy:
+            return vector.detach().numpy()
+        return vector
+
+    def set_vectorized_parameters(self, vector):
+        """Set the vector parameters."""
+        # convert the vector to torch array
+        if isinstance(vector, np.ndarray):
+            vector = torch.from_numpy(vector).float()
+
+        # set the parameters from the vectorized one
+        idx = 0
+        for parameter in self.parameters():
+            size = parameter.nelement()
+            parameter.data = vector[idx:idx+size].reshape(parameter.shape)
+            idx += size
+
+    def predict(self, x, to_numpy=True):
+        r"""
+        Predict output vector :math:`y` given input vector :math:`x`, using the formula: :math:`y = W x + b`
+
+        Args:
+            x (np.ndarray, torch.Tensor): input vector
+            to_numpy (bool): if True, return a np.array
+
+        Returns:
+            np.ndarray, torch.Tensor: output vector
+        """
+        # convert from numpy to pytorch if necessary
+        if isinstance(x, np.ndarray):
+            x = torch.from_numpy(x).float()
+
+        # predict
+        y = self.model(x)
+
+        # return the output and convert it if necessary
+        if to_numpy:
+            if y.requires_grad:
+                return y.detach().numpy()
+            return y.numpy()
+        return y
+
+    def fit(self, X, Y):
+        r"""Train the linear model using Linear Regression.
+
+        This is performed by minimizing the L2 loss with respect to the parameters:
+
+        .. math:: \min_W || Y - XW ||^2
+
+        where :math:`X \in \mathcal{R}^{N \times (D_x+1)}` is the augmented input data matrix, and
+        :math:`Y \in \mathcal{R}^{N \times (D_y)}` is the output data matrix.
+
+        The best set of weights can be obtained by solving in closed-loop the above optimization process.
+        The optimal solution is given by the pseudo-inverse: :math:`W^* = (X^\top X)^{-1} X^T Y`.
+
+        Args:
+            X (np.array[N,Dx], torch.Tensor[N,Dx]): input data matrix.
+            Y (np.array[N,Dy], torch.Tensor[N,Dy]): output data matrix.
+        """
+        pass
+
+    def reset(self):
+        pass
+
+    def __call__(self, x, to_numpy=True):
+        return self.predict(x, to_numpy=to_numpy)
+
+    # def concatenate(self, other):
+    #     """
+    #     Concatenate a linear model with another one.
+    #
+    #     Args:
+    #         other (Linear): the other linear model
+    #     """
+    #     if not isinstance(other, Linear):
+    #         raise TypeError("Expecting the other model to be also linear")
+    #     # self.model.add_module()
+
+
+# Tests
+if __name__ == '__main__':
+    # test with numpy
+    x = np.array(range(3))
+    model = Linear(num_inputs=len(x), num_outputs=2, add_bias=True)
+    print("Linear model's input size: {}".format(model.input_dims))
+    print("Linear model's output size: {}".format(model.output_dims))
+    y = model(x)
+    print("Linear input: {}".format(x))
+    print("Linear output: {}".format(y))
+
+    # test with pytorch
+    x = torch.from_numpy(x).float()
+    y = model(x, to_numpy=False)
+    print("Linear input: {}".format(x))
+    print("Linear torch output: {}".format(y))
+    y = model(x, to_numpy=True)
+    print("Linear numpy output: {}".format(y))

pyrobolearn/models/model.py

@@ -0,0 +1,373 @@
+#!/usr/bin/env python
+"""Define the Model abstract class from which all learning models inherit from.
+
+Dependencies: None
+"""
+
+
+from abc import ABCMeta, abstractmethod
+
+__author__ = "Brian Delhaisse"
+__copyright__ = "Copyright 2018, PyRoboLearn"
+__credits__ = ["Brian Delhaisse"]
+__license__ = "MIT"
+__version__ = "1.0.0"
+__maintainer__ = "Brian Delhaisse"
+__email__ = "briandelhaisse@gmail.com"
+__status__ = "Development"
+
+
+class Model(object):
+    r"""(Learning Base) Model (aka parametrized model)
+
+    This abstract class must be inherited by any learning models, and provides a common API for all these models.
+    It is often a simple wrapper over a specific learning model implemented in a certain library.
+
+    The learning model has no direct knowledge about the inputs and outputs (such as states / actions), and can be
+    used out of the box. In our framework, we refer these learning models as the 'inner models', and the class that
+    connects the inputs/outputs (such as states and actions) with the inner model as the 'outer model'.
+    Outer models are thus learning models that maps states / actions to states / actions, while inner models often
+    maps arrays to arrays (where an array can represent a scalar, a vector, a matrix, a tensor,...). In other words,
+    the outer model is a wrapper around the inner model but knows how to deal with state / action inputs and outputs.
+
+    For instance, neural networks are a popular kind of inner models which map arrays to arrays. These inner models
+    can be used to represent outer models such as rl (which map states to actions), dynamic models (which map
+    states and actions to the next states), value estimators (which, for example, map states to a real number),
+    transformation mappings (which map states to states, or actions to actions). Transformation mappings can
+    for instance be used to map a human kinematic state to a robot kinematic state.
+
+    .. seealso:
+        * https://www.codecademy.com/en/forum_questions/512cd091ffeb9e603b005713
+    """
+    __metaclass__ = ABCMeta
+
+    def __init__(self):
+        """
+        Initialize the learning model.
+        """
+        self._models = []  # TODO: should be a directed graph
+
+        self._input_shape = None
+        self._output_shape = None
+
+    ##############
+    # Properties #
+    ##############
+
+    @property
+    def models(self):
+        return self._models
+
+    @property
+    def input_shape(self):
+        return self._input_shape
+
+    @property
+    def output_shape(self):
+        return self._output_shape
+
+    ##################
+    # Static Methods #
+    ##################
+
+    @staticmethod
+    def is_parametric():
+        """
+        Return True if the model is parametric.
+
+        Returns:
+            bool: True if the model is parametric.
+        """
+        raise NotImplementedError
+
+    @staticmethod
+    def is_linear():
+        """
+        Return True if the model is linear (wrt the parameters). This can be for instance useful for some learning
+        algorithms (some only works on linear models).
+
+        Returns:
+            bool: True if it is a linear model
+        """
+        raise NotImplementedError
+
+    @staticmethod
+    def is_recurrent():
+        """
+        Return True if the model is recurrent. This can be for instance useful for some learning algorithms which
+        change their behavior when they deal with recurrent learning models.
+
+        Returns:
+            bool: True if it is a recurrent model.
+        """
+        raise NotImplementedError
+
+    @staticmethod
+    def is_deterministic():
+        """
+        Return True if the model is deterministic.
+
+        Returns:
+            bool: True if the model is deterministic.
+        """
+        return not Model.is_probabilistic()
+
+    @staticmethod
+    def is_probabilistic():  # is_stochastic
+        """
+        Return True if the model is probabilistic.
+
+        Returns:
+            bool: True if the model is probabilistic.
+        """
+        raise NotImplementedError
+
+    @staticmethod
+    def is_discriminative():
+        """
+        Return True if the model is discriminative, that is, if the model estimates the conditional probability
+        :math:`p(y|x)`.
+
+        Returns:
+            bool: True if the model is discriminative.
+        """
+        raise NotImplementedError
+
+    @staticmethod
+    def is_generative():
+        """
+        Return True if the model is generative, that is, if the model estimates the joint distribution of the input
+        and output :math:`p(x,y)`. A generative model allows to sample from it.
+
+        Returns:
+            bool: True if the model is generative.
+        """
+        raise NotImplementedError
+
+    # TODO: isClassifier, isRegressive, isSequential
+
+    ###########
+    # Methods #
+    ###########
+
+    def has_models(self):
+        return len(self._models) > 0
+
+    def add_model(self, model):
+        if not isinstance(model, Model):
+            raise TypeError("Expecting the model to be an instance of Model.")
+        if self.has_models():
+            # check that the output size of the last model is equal to the input size of the new model
+            last_model = self._models[-1]
+            if last_model.output_dims() != model.input_dims():
+                # TODO
+                pass
+        self._models.append(model)
+
+    @abstractmethod
+    def _predict(self, x=None):
+        """
+        Given a possible input, predict the output. The input doesn't always have to be given. For instance,
+        generative models generate data without any inputs.
+        """
+        raise NotImplementedError
+
+    def predict(self, x):
+        if self.has_models():
+            return [model.predict(x) for model in self.models]
+        return self._predict(x)
+
+    @abstractmethod
+    def parameters(self):
+        """
+        Return an iterator over the parameters of the model.
+        """
+        raise NotImplementedError
+
+    @abstractmethod
+    def named_parameters(self):
+        """
+        Return an iterator over the model parameters, yielding both the name and the parameter itself.
+        """
+        raise NotImplementedError
+
+    @abstractmethod
+    def get_params(self):
+        """
+        Return the parameters in the form of a list or dictionary.
+        """
+        raise NotImplementedError
+
+    @abstractmethod
+    def hyperparameters(self):
+        """
+        Return an iterator over the hyperparameters.
+        """
+        raise NotImplementedError
+
+    @abstractmethod
+    def get_hyperparams(self):
+        """
+        Return the hyperparameters in the form of a list or dictionary.
+        """
+        raise NotImplementedError
+
+    @abstractmethod
+    def get_input_dims(self):
+        raise NotImplementedError
+
+    # alias
+    # input_dims = getInputDims
+
+    @abstractmethod
+    def get_output_dims(self):
+        raise NotImplementedError
+
+    # alias
+    # output_dims = getOutputDims
+
+    @abstractmethod
+    def save(self, filename):
+        """
+        Save the model in memory.
+
+        Args:
+            filename (str): file to save the model in.
+
+        Returns:
+            None
+        """
+        raise NotImplementedError
+
+    @abstractmethod
+    def load(self, filename):
+        """
+        Load a model from memory.
+
+        Args:
+            filename (str): file that contains the model.
+
+        Returns:
+            None
+        """
+        raise NotImplementedError
+
+    def latex(self):
+        """
+        Returns the latex equations that describe the learning model.
+        This function can also be called __repr__ and/or __str__.
+        """
+        raise NotImplementedError
+
+    def bibtex(self):
+        """
+        Returns the references of a learning model in the bibtex format.
+        """
+        raise NotImplementedError
+
+    def concatenate(self, model):
+        """
+        Concatenate sequentially two models.
+        """
+        raise NotImplementedError
+
+    #############
+    # Operators #
+    #############
+
+    def __repr__(self):
+        if self.has_models():
+            lst = [self.__class__.__name__ + '(']
+            for model in self._models:
+                lst.append('\t' + model.__repr__() + ',')
+            lst.append(')')
+            return '\n'.join(lst)
+        else:
+            return self.__class__.__name__
+
+    def __len__(self):
+        if self.has_models():
+            return len(self._models)
+
+    def __call__(self, *args, **kwargs):
+        return self.predict(*args, **kwargs)
+
+    def __add__(self, other):
+        """
+        Define how to combine two models together.
+        There are two ways to combine a model, sequentially (one after another in time), or in parallel.
+        """
+        pass
+
+    def __lshift__(self, other):
+        """
+        Define how to concatenate/sequence two models inline.
+
+        Examples:
+            nn = Model()
+            nn1 = MLP()
+            nn2 = MLP()
+            nn << nn1 << nn2  # same as nn << (nn1 >> nn2)
+        """
+        pass
+
+    def __rshift__(self, other):
+        """
+        Define how to concatenate/sequence two models. The input of the first model is given to the output
+        of the second model.
+        """
+        # if same model check `rshift()` fct in the corresponding class
+        # if different models, it is defined here
+        if isinstance(self, NN):
+            if isinstance(self, NN):
+                pass  # TODO: call rshift()
+            elif isinstance(other, CPG):
+                pass
+            elif isinstance(other, DMP):
+                pass
+            else:
+                raise NotImplementedError("Do not know how to concatenate {} and {}.".format(type(self).__name__,
+                                                                                             type(other).__name__))
+        elif isinstance(self, GP):
+            if isinstance(other, GP):
+                pass  # TODO: call rshift
+            elif isinstance(other, DMP):
+                pass
+            else:
+                raise NotImplementedError("Do not know how to concatenate {} and {}.".format(type(self).__name__,
+                                                                                             type(other).__name__))
+
+        elif isinstance(self, GMM):
+            if isinstance(other, GMM):
+                pass
+            elif isinstance(other, DMP):
+                pass
+            else:
+                raise NotImplementedError("Do not know how to concatenate {} and {}.".format(type(self).__name__,
+                                                                                             type(other).__name__))
+
+        else:
+            raise NotImplementedError("Do not know how to concatenate {} and {}.".format(type(self).__name__,
+                                                                                         type(other).__name__))
+
+    def __floordiv__(self, other):
+        """
+        Define how to recursively sequence two models.
+        """
+        pass
+
+    def __mod__(self, other):
+        """
+        Define how much to wait when producing the output.
+        Basically, the output of a model is put into a FIFO queue of size :math:`s` which is specified by the given
+        argument. Thus at any time steps :math:`t`, the output actually produced is :math:`o_{t-s}`.
+        """
+        if not isinstance(other, int):
+            raise TypeError("Expecting an integer")
+
+    def __mul__(self, other):
+        pass
+
+    def __getitem__(self, key):
+        pass
+

pyrobolearn/models/nn/__init__.py

@@ -0,0 +1,24 @@
+
+# import general neural network model
+from dnn import NN
+
+# import multilayer perceptron model
+from mlp import *
+
+# import NEAT model
+from neat_model import NEATModel
+
+# import convolutional neural network
+# from cnn import *
+
+# import recurrent neural network
+# from rnn import *
+
+# import auto-encoder
+# from ae import *
+
+# import variational auto-encoder
+# from vae import *
+
+# import generative adversarial networks
+# from gan import *

pyrobolearn/models/nn/ae.py

@@ -0,0 +1,111 @@
+#!/usr/bin/env python
+"""Define the Auto-Encoder (AE) learning model.
+
+This file provides the AE model; a parametric, generally non-linear, non-recurrent, discriminative,
+and deterministic model. This model is a latent variable model which projects the input data into a latent lower
+dimensional space through an encoder, and re-projects it to the original data space through the use of the decoder.
+
+We decided to use the `pytorch` framework because of its popularity in the research community field, flexibility,
+similarity with numpy (but with automatic differentiation: autograd), GPU capabilities, and more Pythonic approach.
+While we thought about using other frameworks (such as Keras, Tensorflow, and others) as well, it would have
+unnecessarily complexify the whole framework, as these frameworks would not only influence the learning models,
+but also the losses, optimizers, and other modules. While we could have written some interfaces that makes the bridge
+between these various frameworks and ours, we came to the conclusion that this would take a considerable amount of
+efforts and time that we do not have for the moment.
+
+References:
+    [1] "Deep Learning" (http://www.deeplearningbook.org/), Goodfellow et al., 2016
+    [2] PyTorch: https://pytorch.org/
+"""
+
+import copy
+import inspect
+import numpy as np
+import torch
+
+from dnn import NN
+
+__author__ = "Brian Delhaisse"
+__copyright__ = "Copyright 2018, PyRoboLearn"
+__credits__ = ["Brian Delhaisse"]
+__license__ = "(c) Brian Delhaisse"
+__version__ = "1.0.0"
+__maintainer__ = "Brian Delhaisse"
+__email__ = "briandelhaisse@gmail.com"
+__status__ = "Development"
+
+
+class AE(NN):
+    r"""Auto-Encoder
+
+    References:
+        [1] "Deep Learning" (http://www.deeplearningbook.org/), Goodfellow et al., 2016
+    """
+
+    def __init__(self):
+        pass
+
+
+class _AETorch(torch.nn.Module):
+    r"""Auto-Encoder written in Pytorch (that inherits from `torch.nn.Module`)
+    """
+
+    def __init__(self, layer_sizes=[], activation_fct=None, dropout=None, encoder=None, decoder=None):
+        super(_AETorch, self).__init__()
+
+        if encoder is None and decoder is None:
+            # nb of layers (the input layer doesn't count)
+            self.num_layers = len(layer_sizes) - 1
+            layers = [torch.nn.Linear(layer_sizes[i], layer_sizes[i+1]) for i in range(self.num_layers)]
+
+            # check activation function and insert it after each linear layer
+            if activation_fct is not None:
+                if isinstance(activation_fct, str):
+                    activation_fct = getattr(torch.nn, activation_fct)()
+                elif activation_fct.__module__ == 'torch.nn.modules.activation':
+                    if inspect.isclass(activation_fct):
+                        activation_fct = activation_fct()
+                else:
+                    raise ValueError("activation_fct should be a string or belong to torch.nn.modules.activation")
+
+                # add activation layer
+                for i in range(self.num_layers-1, 0, -1):
+                    layers.insert(activation_fct)
+
+            # check dropout
+            if dropout is not None:
+                if isinstance(dropout, float):
+                    dropout = torch.nn.Dropout(dropout)
+                elif dropout.__module__ == 'torch.nn.modules.dropout':
+                    raise ValueError("Dropout should be a float or belong to torch.nn.modules.dropout")
+
+                # add dropout layer
+                for i in range(self.num_layers-1, 0, -2):
+                    layers.insert(dropout)
+
+            # Encoder
+            encoder = torch.nn.Sequential(*layers[:len(layers)//2])
+            # Decoder
+            decoder = torch.nn.Sequential(*layers[len(layers)//2:])
+
+        self.encoder = encoder
+        self.decoder = decoder
+
+    def forward(self, x):
+        x = self.encoder(x)
+        x = self.decoder(x)
+        return x
+
+    def encode(self, x):
+        x = self.encoder(x)
+        return x
+
+    def decode(self, x):
+        x = self.decoder(x)
+        return x
+
+
+class AETorch(NNTorch):
+    r"""Auto-Encoder in Pytorch
+    """
+    pass

pyrobolearn/models/nn/cnn.py

@@ -0,0 +1,49 @@
+#!/usr/bin/env python
+"""Define the Convolutional Neural Network (CNN) learning model.
+
+This file provides the CNN model; a parametric, generally non-linear, non-recurrent, discriminative,
+and deterministic model. This model is convenient for data arrays/tensors that have cells that have a spatial
+relationship between them. For instance, pictures are 2D or 3D arrays where each pixel is related with its neighbors.
+
+We decided to use the `pytorch` framework because of its popularity in the research community field, flexibility,
+similarity with numpy (but with automatic differentiation: autograd), GPU capabilities, and more Pythonic approach.
+While we thought about using other frameworks (such as Keras, Tensorflow, and others) as well, it would have
+unnecessarily complexify the whole framework, as these frameworks would not only influence the learning models,
+but also the losses, optimizers, and other modules. While we could have written some interfaces that makes the bridge
+between these various frameworks and ours, we came to the conclusion that this would take a considerable amount of
+efforts and time that we do not have for the moment.
+
+References:
+    [1] "Deep Learning" (http://www.deeplearningbook.org/), Goodfellow et al., 2016
+    [2] PyTorch: https://pytorch.org/
+"""
+
+import copy
+import inspect
+import numpy as np
+import torch
+
+from dnn import NN
+
+__author__ = "Brian Delhaisse"
+__copyright__ = "Copyright 2018, PyRoboLearn"
+__credits__ = ["Brian Delhaisse"]
+__license__ = "(c) Brian Delhaisse"
+__version__ = "1.0.0"
+__maintainer__ = "Brian Delhaisse"
+__email__ = "briandelhaisse@gmail.com"
+__status__ = "Development"
+
+
+class CNN(NN):
+    r"""Convolutional Neural Network
+
+    Feed-forward CNN.
+    """
+    pass
+
+
+class CNNTorch(NNTorch):
+    r"""Convolutional Neural Network in PyTorch
+    """
+    pass

pyrobolearn/models/nn/dnn.py

@@ -0,0 +1,257 @@
+#!/usr/bin/env python
+"""Define the Deep Neural Network (DNN) learning model.
+
+This file provides the DNN model; a parametric, generally non-linear, possibly recurrent, discriminative/generative,
+and deterministic/probabilistic model.
+
+We decided to use the `pytorch` framework because of its popularity in the research community field, flexibility,
+similarity with numpy (but with automatic differentiation: autograd), GPU capabilities, and more Pythonic approach.
+While we thought about using other frameworks (such as Keras, Tensorflow, and others) as well, it would have
+unnecessarily complexify the whole framework, as these frameworks would not only influence the learning models,
+but also the losses, optimizers, and other modules. While we could have written some interfaces that makes the bridge
+between these various frameworks and ours, we came to the conclusion that this would take a considerable amount of
+efforts and time that we do not have for the moment.
+
+References:
+    [1] "Deep Learning" (http://www.deeplearningbook.org/), Goodfellow et al., 2016
+    [2] PyTorch: https://pytorch.org/
+"""
+
+import copy
+import inspect
+import numpy as np
+import torch
+
+# from pyrobolearn.models import Model
+
+__author__ = "Brian Delhaisse"
+__copyright__ = "Copyright 2018, PyRoboLearn"
+__credits__ = ["Brian Delhaisse"]
+__license__ = "(c) Brian Delhaisse"
+__version__ = "1.0.0"
+__maintainer__ = "Brian Delhaisse"
+__email__ = "briandelhaisse@gmail.com"
+__status__ = "Development"
+
+
+class NN(object):  # Model
+    r"""Neural Network class
+
+    This class describes the neural network model. It is basically a wrapper around deep learning frameworks such
+    as pytorch, Keras, tensorflow and others. This class is inherited by any other neural network classes, such as
+    convolution neural networks, recurrent neural networks, and so on.
+
+    Note that we currently mainly focus on `PyTorch`.
+
+    * PyTorch (https://pytorch.org/)
+        * PyTorch is ... PyTorch allows for dynamic ...
+        * torch.nn.Module: it represents the base class for all the neural networks / layers
+        * torch.nn.modules.loss(torch.nn.Module): it contains the definition of some popular losses
+        * torch.optim.Optimizer: it is the base class for all the optimizers
+
+    Examples::
+
+        import torch.nn as nn
+
+        model = nn.Sequential(nn.Conv2d(in_channels=3, out_channels=10, kernel_size=3), Flatten(), nn.Linear(320, 10))
+        model = NN(model, input_dims=..., output_dims=...)
+
+    References:
+        [1] "Deep Learning" (http://www.deeplearningbook.org/), Goodfellow et al., 2016
+        [2] PyTorch: https://pytorch.org/
+        nn.Module: https://pytorch.org/docs/master/_modules/torch/nn/modules/module.html#Module
+# - nn.Sequential: https://pytorch.org/docs/master/_modules/torch/nn/modules/container.html#Sequential
+    """
+
+    def __init__(self, model, input_dims, output_dims, framework=None):
+        r"""Initialize the NN model.
+
+        Args:
+            model (torch.nn.Module or keras.base_layer.Layer): pytorch/keras model
+            input_dims (int, tuple/list of int): dimensions of the input
+            output_dims (int, tuple/list of int): dimensions of the output
+        """
+        super(NN, self).__init__()
+
+        # check if given model is valid
+        # if model is not None:
+        #     if isinstance(model, torch.nn.Module):
+        #         self.framework = 'pytorch'
+        #     elif isinstance(model, keras.models.Model):
+        #         self.framework = 'keras'
+        #     else:
+        #         raise TypeError("Model should be an instance of torch.nn.Module")
+
+        # set model (written in the specified framework)
+        self.model = model
+        self.input_dims = input_dims
+        self.output_dims = output_dims
+
+        # TODO: infer the framework based on the model
+        self.framework = framework
+
+    ##############
+    # Properties #
+    ##############
+
+    @property
+    def model(self):
+        return self._model
+
+    @model.setter
+    def model(self, model):
+        if model is not None:
+            if not (isinstance(model, torch.nn.Module) or isinstance(model, keras.models.Model)):
+                raise TypeError("The model should be an instance of torch.nn.Module or keras.models.Model")
+        self._model = model
+
+    @property
+    def input_shape(self):  # TODO
+        return self.input_dims
+
+    @property
+    def output_shape(self):  # TODO
+        return self.output_dims
+
+    ##################
+    # Static Methods #
+    ##################
+
+    @staticmethod
+    def is_parametric():
+        """A neural network is a parametric model"""
+        return True
+
+    @staticmethod
+    def is_linear():  # unless all layers are linear
+        """A neural network is in general non-linear, where non-linear activation functions are applied on each
+         layer output. If all the activation layers are linear, then the NN is linear."""
+        return False
+
+    @staticmethod
+    def is_recurrent():  # unless RNN
+        """Unless the neural network is a RNN, it is not recurrent."""
+        return False
+
+    @staticmethod
+    def is_probabilistic():  # unless variational methods (dropouts,...)
+        """The neural network is not a probabilistic model per se. However, it can be simulated by using dropouts,
+        or using a probabilistic distribution on the output of the last layer. For instance, the network can
+        output the mean and covariance matrices which parametrizes a Gaussian probabilistic distribution"""
+        return False  # if False then it is deterministic
+
+    @staticmethod
+    def is_discriminative():
+        """A neural network is a discriminative model which given inputs predicts some outputs"""
+        return True
+
+    @staticmethod
+    def is_generative():  # unless VAE, GAN,...
+        """Standard neural networks are not generative, and thus we can not sample from it. This is different,
+        for instance, for generative adversarial networks (GANs) and variational auto-encoders (VAEs)."""
+        return False
+
+    ###########
+    # Methods #
+    ###########
+
+    def predict(self, x=None):
+        return self.model(x)
+
+    def get_input_dims(self):
+        return self.input_dims
+
+    def get_output_dims(self):
+        return self.output_dims
+
+    def parameters(self):
+        return self.model.parameters()
+
+    def get_params(self):
+        return list(self.parameters())
+
+    def get_hyperparams(self):
+        """
+        Return the number of units per layer, the number of layers, and the type of layers.
+        """
+        raise NotImplementedError
+
+    def hyperparameters(self):
+        raise NotImplementedError
+
+    #############
+    # Operators #
+    #############
+
+    def __str__(self):
+        """
+        Return string describing the NN model.
+        """
+        if self.framework == 'pytorch':
+            return str(self.model)
+        elif self.framework == 'keras':
+            summary = []
+            self.model.summary(print_fn=lambda s: summary.append(s))
+            return '\n'.join(summary)
+        else:
+            raise NotImplementedError
+
+    def __getitem__(self, key):
+        """
+        Return the corresponding layer(s).
+        """
+        return self.model[key]
+
+    def __rshift__(self, other):
+        """
+        Concatenate two NN models in sequence, and return the sequenced model.
+        Note that It doesn't modify the given models, but return a new one.
+
+        Args:
+            other (NN): other NN model
+
+        Returns:
+            NN: sequenced model
+        """
+        # copy current model
+        model = copy.deepcopy(self.model)
+
+        # concatenate the other model
+        for idx, item in enumerate(other.model, start=len(self.model)):
+            model.add_module(str(idx), item)
+
+        # return the concatenation
+        return NN(model)
+
+
+class NNTorch(NN):
+    r"""Neural Network written in PyTorch
+    """
+
+    def __init__(self, model, input_dims, output_dims):
+        super(NNTorch, self).__init__(model, input_dims, output_dims, framework='pytorch')
+
+    def save(self, filename):
+        """
+        Save the neural network to the specified file.
+
+        Args:
+            filename (str): filename to save the neural network
+        """
+        torch.save(self.model, filename)
+
+    def load(self, filename):
+        """
+        Load the neural network from the specified file.
+
+        Args:
+            filename (str): filename from which to load the neural network
+        """
+        self.model = torch.load(filename)
+        # check input and output dimensions
+
+    def __str__(self):
+        """
+        Return string describing the NN model.
+        """
+        return str(self.model)

pyrobolearn/models/nn/gan.py

@@ -0,0 +1,70 @@
+#!/usr/bin/env python
+"""Define the Generative Adversarial Network (GAN) learning model.
+
+This file provides the GAN model; a parametric, generally non-linear, non-recurrent, generative, and stochastic model.
+This is a generative model which works in a game theory setting by having a generator and discriminator compete
+between each other. The goal of the generator is to generate data samples that are similar to the provided dataset
+and fool the discriminator. The goal of the discriminator is to discriminate the given samples by identifying the fake
+ones from the true ones.
+
+We decided to use the `pytorch` framework because of its popularity in the research community field, flexibility,
+similarity with numpy (but with automatic differentiation: autograd), GPU capabilities, and more Pythonic approach.
+While we thought about using other frameworks (such as Keras, Tensorflow, and others) as well, it would have
+unnecessarily complexify the whole framework, as these frameworks would not only influence the learning models,
+but also the losses, optimizers, and other modules. While we could have written some interfaces that makes the bridge
+between these various frameworks and ours, we came to the conclusion that this would take a considerable amount of
+efforts and time that we do not have for the moment.
+
+References:
+    [1] "Deep Learning" (http://www.deeplearningbook.org/), Goodfellow et al., 2016
+    [2] PyTorch: https://pytorch.org/
+"""
+
+import copy
+import inspect
+import numpy as np
+import torch
+
+from dnn import NN
+
+__author__ = "Brian Delhaisse"
+__copyright__ = "Copyright 2018, PyRoboLearn"
+__credits__ = ["Brian Delhaisse"]
+__license__ = "(c) Brian Delhaisse"
+__version__ = "1.0.0"
+__maintainer__ = "Brian Delhaisse"
+__email__ = "briandelhaisse@gmail.com"
+__status__ = "Development"
+
+
+class GAN(NN):
+    r"""Generative Adversarial Network
+
+    Type: generative model
+
+    .. seealso:: Variational Auto-Encoders
+
+    References:
+        [1] "NIPS:
+        [2]
+    """
+
+    def __init__(self):
+        pass
+
+    @staticmethod
+    def isDiscriminative():
+        """A neural network is a discriminative model which given inputs predicts some outputs"""
+        return True
+
+    @staticmethod
+    def isGenerative():  # unless VAE, GAN,...
+        """Standard neural networks are not generative, and thus we can not sample from it. This is different,
+        for instance, for generative adversarial networks (GANs) and variational auto-encoders (VAEs)."""
+        return True
+
+
+class GANTorch(NNTorch):
+    r"""Generative Adversarial Network in PyTorch
+    """
+    pass

pyrobolearn/models/nn/mlp.py

@@ -0,0 +1,179 @@
+#!/usr/bin/env python
+"""Define the Multi-Layer Perceptron (MLP) learning model.
+
+This file provides the MLP model; a parametric, generally non-linear, non-recurrent, discriminative,
+and deterministic model.
+
+We decided to use the `pytorch` framework because of its popularity in the research community field, flexibility,
+similarity with numpy (but with automatic differentiation: autograd), GPU capabilities, and more Pythonic approach.
+While we thought about using other frameworks (such as Keras, Tensorflow, and others) as well, it would have
+unnecessarily complexify the whole framework, as these frameworks would not only influence the learning models,
+but also the losses, optimizers, and other modules. While we could have written some interfaces that makes the bridge
+between these various frameworks and ours, we came to the conclusion that this would take a considerable amount of
+efforts and time that we do not have for the moment.
+
+References:
+    [1] "Deep Learning" (http://www.deeplearningbook.org/), Goodfellow et al., 2016
+    [2] PyTorch: https://pytorch.org/
+"""
+
+import copy
+import inspect
+import numpy as np
+import torch
+
+from dnn import NN, NNTorch
+
+__author__ = "Brian Delhaisse"
+__copyright__ = "Copyright 2018, PyRoboLearn"
+__credits__ = ["Brian Delhaisse"]
+__license__ = "(c) Brian Delhaisse"
+__version__ = "1.0.0"
+__maintainer__ = "Brian Delhaisse"
+__email__ = "briandelhaisse@gmail.com"
+__status__ = "Development"
+
+
+class MLP(NN):
+    r"""Multi-Layer Perceptron
+
+    Feed-forward and fully-connected neural network, where linear layers are followed by non-linear activation
+    functions.
+
+    .. math::
+
+        h_{l} = f_{l}(W_{l} h_{l-1} + b_{l})
+
+    where :math:`l \in [1,...,L]` with :math:`L` is the total number of layers,
+    :math:`W_{l}` and :math:`b_{l}` are the weight matrix and bias vector at layer :math:`l`, :math:`f_{l}` is
+    the nonlinear activation function, and :math:`h_{0} = x` and :math:`y = h_{L}` are the input and output vectors.
+
+    The parameters of the neural networks are all the weight matrices and bias vectors.
+    """
+
+    def __init__(self, num_units=(), activation_fct='Linear', last_activation_fct=None, dropout_prob=None,
+                 framework='pytorch'):
+        """
+        Initialize a MLP network.
+
+        Args:
+            num_units (list/tuple of int): number of units in each layer (this includes the input and output layer)
+            activation_fct (None, str, or list/tuple of str/None): activation function to be applied after each layer.
+                                                                   If list/tuple, then it has to match the
+            last_activation_fct (None or str): last activation function to be applied. If not specified, it will check
+                                               if it is in the list/tuple of activation functions provided for the
+                                               previous argument.
+            dropout_prob (None, float, or list/tuple of float/None): dropout probability.
+            framework (str): specifies which framework we want to use between 'pytorch' and 'keras' (default: 'pytorch')
+        """
+
+        # check framework
+        framework = framework.lower()
+        if framework == 'pytorch':
+            model = MLPTorch(num_units, activation_fct, last_activation_fct, dropout_prob)
+        elif framework == 'keras':
+            model = MLPKeras(num_units, activation_fct, last_activation_fct, dropout_prob)
+        else:
+            raise ValueError("The current frameworks allowed are pytorch and keras")
+
+        # instantiate the super class
+        super(MLP, self).__init__(model.model, input_dims=num_units[0], output_dims=num_units[-1], framework=framework)
+
+        # rewrite methods
+        self.save = model.save
+        self.load = model.load
+        self.__str__ = model.__str__
+
+
+class MLPTorch(NNTorch):
+    r"""Multi-Layer Perceptron in PyTorch
+
+    Feed-forward and fully-connected neural network, where linear layers are followed by non-linear activation
+    functions.
+
+    .. math::
+
+        h_{l} = f_{l}(W_{l} h_{l-1} + b_{l})
+
+    where :math:`l \in [1,...,L]` with :math:`L` is the total number of layers,
+    :math:`W_{l}` and :math:`b_{l}` are the weight matrix and bias vector at layer :math:`l`, :math:`f_{l}` is
+    the nonlinear activation function, and :math:`h_{0} = x` and :math:`y = h_{L}` are the input and output vectors.
+
+    The parameters of the neural networks are all the weight matrices and bias vectors.
+    """
+
+    def __init__(self, num_units=(), activation_fct='Linear', last_activation_fct=None, dropout_prob=None):
+        """
+        Initialize a MLP network.
+
+        Args:
+            num_units (list/tuple of int): number of units in each layer (this includes the input and output layer)
+            activation_fct (None, str, or list/tuple of str/None): activation function to be applied after each layer.
+                                                                   If list/tuple, then it has to match the
+            last_activation_fct (None or str): last activation function to be applied. If not specified, it will check
+                                               if it is in the list/tuple of activation functions provided for the
+                                               previous argument.
+            dropout_prob (None, float, or list/tuple of float/None): dropout probability.
+        """
+        # check number of units
+        if len(num_units) < 2:
+            raise ValueError("The num_units list/tuple needs to have at least the input and output layers")
+
+        # set the dimensions of the input and output
+        self.input_dims = num_units[0]
+        self.output_dims = num_units[-1]
+
+        # check for activation fcts
+        activations = dir(torch.nn.modules.activation)
+        activations = {act: act for act in activations}
+        activations.update({act.lower(): act for act in activations})
+
+        def check_activation(activation):
+            if activation is None or activation.lower() == 'linear':
+                activation = None
+            else:
+                if activation not in activations:
+                    raise ValueError("The given activation function is not available")
+                activation = getattr(torch.nn, activations[activation])
+            return activation
+
+        activation_fct = check_activation(activation_fct)
+        last_activation_fct = check_activation(last_activation_fct)
+
+        # check dropout
+        dropout_layer = None
+        if dropout_prob is not None:
+            dropout_layer = torch.nn.Dropout(dropout_prob)
+
+        # build pytorch network
+        layers = []
+        for i in range(len(num_units[:-2])):
+            # add linear layer
+            layer = torch.nn.Linear(num_units[i], num_units[i + 1])
+            layers.append(layer)
+
+            # add activation layer
+            if activation_fct is not None:
+                layers.append(activation_fct())
+
+            # add dropout layer
+            if dropout_layer is not None:
+                layers.append(dropout_layer)
+
+        # last output layer
+        layers.append(torch.nn.Linear(num_units[-2], num_units[-1]))
+        if last_activation_fct is not None:
+            layers.append(last_activation_fct)
+
+        # create nn model
+        model = torch.nn.Sequential(*layers)
+
+        super(MLPTorch, self).__init__(model, input_dims=num_units[0], output_dims=num_units[-1])
+
+
+# Tests
+if __name__ == '__main__':
+
+    # create MLP network
+    mlp = MLPTorch(num_units=(2,10,3), activation_fct='relu')
+    print(mlp)

pyrobolearn/models/nn/neat_model.py

@@ -0,0 +1,329 @@
+#!/usr/bin/env python
+"""Define the NEAT model class.
+
+This uses the Neuro-Evolution through Augmenting topologies (NEAT) framework. It allows the evolution of not only the
+parameters/weights but also the topological structure of neural networks. Note that the model associated with
+this policy (i.e. the neural network) is tightly coupled with the algorithm that modifies it.
+"""
+
+import numpy as np
+import cPickle as pickle
+
+try:
+    import neat
+except ImportError as e:
+    raise ImportError(e.__str__() + "\n HINT: you can install NEAT directly via 'pip install neat-python'.")
+
+# from pyrobolearn.models import Model
+
+__author__ = "Brian Delhaisse"
+__copyright__ = "Copyright 2018, PyRoboLearn"
+__credits__ = ["Brian Delhaisse"]
+__license__ = "(c) Brian Delhaisse"
+__version__ = "1.0.0"
+__maintainer__ = "Brian Delhaisse"
+__email__ = "briandelhaisse@gmail.com"
+__status__ = "Development"
+
+
+class NEATModel(object):  # Model):
+    r"""NEAT Model
+
+    NEAT stands for "Neuro-Evolution through Augmenting Topologies" [1] and allows the evolution of not only the
+    parameters/weights but also the topological structure of neural networks. The model (i.e. the neural network)
+    is tightly coupled with the algorithm that modifies it. By the structure of the neural network, we mean
+    the type (i.e. forward or recurrent) and number of connection, as well as the type (i.e. using non-linearity
+    activation function) and number of nodes can change.
+
+    This model works in a Reinforcement Learning setting, where exploration is carried out in the parameter and
+    hyper-parameter spaces of the neural network.
+
+    Warnings: The associated algorithm is a little bit special and currently only works with the corresponding
+    learning model.
+
+    References:
+        [1] "Evolving Neural Networks through Augmenting Topologies", Stanley et al., 2002
+        [2] NEAT-Python
+            - documentation: https://neat-python.readthedocs.io/en/latest/index.html
+            - github repo: https://github.com/CodeReclaimers/neat-python
+        [3] PyTorch NEAT (built upon NEAT-Python): https://github.com/uber-research/PyTorch-NEAT
+    """
+
+    def __init__(self, num_inputs, num_outputs, num_hidden=0, activation_fct='relu', network_type='feedforward',
+                 aggregation='sum', weights_limits=(-20, 20), bias_limits=(-20, 20)):
+        # super(NEATModel, self).__init__()
+
+        if network_type != 'feedforward' and network_type != 'recurrent':
+            raise ValueError("Expecting the 'network_type' argument to be 'feedforward' or 'recurrent'. Received "
+                             "instead {}".format('network_type'))
+        self.network_type = network_type
+
+        # set config file
+        # more info about genome's config file: https://neat-python.readthedocs.io/en/latest/config_file.html
+        # more info about activation fct: https://neat-python.readthedocs.io/en/latest/activation.html
+        self.config_dict = {'[NEAT]': {'fitness_criterion': 'max',
+                                       'pop_size': 10,
+                                       'fitness_threshold': 100,
+                                       'no_fitness_termination': True,
+                                       'reset_on_extinction': True},
+                            '[DefaultSpeciesSet]': {'compatibility_threshold': 3.0},
+                            '[DefaultStagnation]': {'species_fitness_func': 'max',
+                                                    'max_stagnation': 15,
+                                                    'species_elitism': 2},
+                            '[DefaultReproduction]': {'elitism': 2,
+                                                      'survival_threshold': 0.2,
+                                                      'min_species_size': 2},
+                            '[DefaultGenome]': {'activation_default': activation_fct,
+                                                'activation_mutate_rate': 0.0,
+                                                'activation_options': activation_fct,
+                                                # node aggregation options
+                                                'aggregation_default': aggregation,
+                                                'aggregation_mutate_rate': 0.0,
+                                                'aggregation_options': aggregation,
+                                                # node bias options
+                                                'bias_init_mean': 0.0,
+                                                'bias_init_stdev': 1.0,
+                                                'bias_init_type': 'gaussian',
+                                                'bias_max_value': bias_limits[1],
+                                                'bias_min_value': bias_limits[0],
+                                                'bias_mutate_power': 0.5,
+                                                'bias_mutate_rate': 0.7,
+                                                'bias_replace_rate': 0.1,
+                                                # genome compatibility options
+                                                'compatibility_disjoint_coefficient': 1.0,
+                                                'compatibility_weight_coefficient': 0.5,
+                                                # connection add/remove rates
+                                                'conn_add_prob': 0.5,
+                                                'conn_delete_prob': 0.5,
+                                                # connection enable options
+                                                'enabled_default': True,
+                                                'enabled_mutate_rate': 0.01,
+                                                'feed_forward': self.network_type == 'feedforward',
+                                                'initial_connection': 'full_nodirect',
+                                                # node add/remove rates
+                                                'node_add_prob': 0.1,
+                                                'node_delete_prob': 0.1,
+                                                # network parameters
+                                                'num_hidden': num_hidden,
+                                                'num_inputs': num_inputs,
+                                                'num_outputs': num_outputs,
+                                                # node response options
+                                                'response_init_mean': 1.0,
+                                                'response_init_stdev': 0.0,
+                                                'response_max_value': 30.0,
+                                                'response_min_value': -30.0,
+                                                'response_mutate_power': 0.0,
+                                                'response_mutate_rate': 0.0,
+                                                'response_replace_rate': 0.0,
+                                                # connection weight options
+                                                'weight_init_mean': 0.0,
+                                                'weight_init_stdev': 1.0,
+                                                'weight_max_value': weights_limits[1],
+                                                'weight_min_value': weights_limits[0],
+                                                'weight_mutate_power': 0.5,
+                                                'weight_mutate_rate': 0.8,
+                                                'weight_replace_rate': 0.1}}
+
+        # create config
+        config_file = self._create_config_file()
+        self._config = neat.Config(neat.DefaultGenome, neat.DefaultReproduction, neat.DefaultSpeciesSet,
+                                   neat.DefaultStagnation, config_file)
+
+        # create initial population
+        self.population = neat.Population(self.config)
+
+        # set initial genome (first genome from the population)
+        self._genome = self.population.population[1]
+
+        # create network
+        self.model = self.set_network(self.genome, self.config)
+
+    ##############
+    # Properties #
+    ##############
+
+    @property
+    def input_dims(self):
+        """Return the input dimension of the model"""
+        return len(self.model.input_nodes)
+
+    @property
+    def output_dims(self):
+        """Return the output dimension of the model"""
+        return len(self.model.output_nodes)
+
+    @property
+    def input_shape(self):
+        """Return the input shape of the model"""
+        return tuple([self.input_dims])
+
+    @property
+    def output_shape(self):
+        """Return the output shape of the model"""
+        return tuple([self.output_dims])
+
+    @property
+    def config(self):
+        """Return the config object"""
+        return self._config
+
+    @config.setter
+    def config(self, config):
+        """Set the config file (str) or object."""
+        if not isinstance(config, neat.config.Config):
+            raise TypeError("Expecting genome to be an instance of neat.config.Config.")
+        self._config = config
+        # create population
+        self.population = neat.Population(self._config)
+
+    @property
+    def genome(self):
+        return self._genome
+
+    @genome.setter
+    def genome(self, genome):
+        if not isinstance(genome, neat.genome.DefaultGenome):
+            raise TypeError("Expecting genome to be an instance of neat.genome.DefaultGenome type")
+        self._genome = genome
+
+        # create network
+        self.model = self.set_network(self._genome, self.config)
+
+    @property
+    def network(self):
+        return self.model
+
+    ##################
+    # Static Methods #
+    ##################
+
+    @staticmethod
+    def is_parametric():
+        return True
+
+    @staticmethod
+    def is_linear():
+        return False
+
+    @staticmethod
+    def is_recurrent():
+        return True
+
+    @staticmethod
+    def is_probabilistic():
+        return False
+
+    @staticmethod
+    def is_discriminative():
+        return True
+
+    @staticmethod
+    def is_generative():
+        return False
+
+    @staticmethod
+    def load(filename):
+        """
+        Load a model from memory.
+
+        Args:
+            filename (str): file that contains the model.
+        """
+        return pickle.load(open(filename, 'rb'))
+
+    ###########
+    # Methods #
+    ###########
+
+    def _create_config_str(self, config_dict=None):
+        """Create string describing the config file from a config dictionary"""
+        if config_dict is None:
+            config_dict = self.config_dict
+
+        # create config file
+        config = []
+        for section, parameters in config_dict.items():
+            config.append(section)
+            for key, value in parameters.items():
+                config.append(key + ' = ' + str(value))
+            config.append('')
+
+        # return string describing the config file
+        return '\n'.join(config)
+
+    def _create_config_file(self, config_dict=None):
+        """Create config file from a config dictionary"""
+        config = self._create_config_str(config_dict)
+        filename = 'config.txt'
+
+        # create config file
+        with open(filename, 'w') as f:
+            f.write(config)
+
+        # return path to the config file
+        return filename
+
+    def set_network(self, genome=None, config=None):
+        # check arguments
+        if genome is None:
+            genome = self.genome
+        if config is None:
+            config = self.config
+
+        # Create the neural network
+        if self.network_type == 'feedforward':
+            self.model = neat.nn.FeedForwardNetwork.create(genome, config)
+        elif self.network_type == 'recurrent':
+            self.model = neat.nn.RecurrentNetwork.create(genome, config)
+        else:
+            raise TypeError('Choose between feedforward or recurrent or implement your own type of NN.')
+
+        return self.model
+
+    def update_config(self, config):
+        # update config (dict)
+        if isinstance(config, dict):
+            self.config_dict.update(config)
+            config_file = self._create_config_file()
+            self._config = neat.Config(neat.DefaultGenome, neat.DefaultReproduction, neat.DefaultSpeciesSet,
+                                       neat.DefaultStagnation, config_file)
+        elif isinstance(config, neat.Config):
+            self._config = config
+
+        # create new population
+        self.population = neat.Population(self.config)
+
+        # set new genome (first genome from the population)
+        self._genome = self.population.population[1]
+
+        # set new network
+        self.model = self.set_network(self.genome, self.config)
+
+    def parameters(self):
+        """Returns an iterator over the model parameters."""
+        return []
+
+    def named_parameters(self):
+        """Returns an iterator over the model parameters, yielding both the name and the parameter itself"""
+        return []
+
+    def hyperparameters(self):
+        pass
+
+    def get_params(self):
+        pass
+
+    def get_hyperparams(self):
+        pass
+
+    def predict(self, x=None):
+        return self.model.activate(x)
+
+    def save(self, filename):
+        """
+        Save the model in memory.
+
+        Args:
+            filename (str): file to save the model in.
+        """
+        pickle.dump(self, open(filename, 'wb'))

pyrobolearn/models/nn/rcnn.py

@@ -0,0 +1,47 @@
+#!/usr/bin/env python
+"""Define the Recurrent Convolutional Neural Network (RCNN) learning model.
+
+This file provides the RCNN model; a parametric, generally non-linear, recurrent, discriminative,
+and deterministic model. This model is convenient for sequential data arrays/tensors where at each instant, the cells
+in the data arrays/tensors have a spatial relationship between them. For instance, this model can be used with videos
+where there is a temporal and spatial relationship between the pixels.
+
+We decided to use the `pytorch` framework because of its popularity in the research community field, flexibility,
+similarity with numpy (but with automatic differentiation: autograd), GPU capabilities, and more Pythonic approach.
+While we thought about using other frameworks (such as Keras, Tensorflow, and others) as well, it would have
+unnecessarily complexify the whole framework, as these frameworks would not only influence the learning models,
+but also the losses, optimizers, and other modules. While we could have written some interfaces that makes the bridge
+between these various frameworks and ours, we came to the conclusion that this would take a considerable amount of
+efforts and time that we do not have for the moment.
+
+References:
+    [1] "Deep Learning" (http://www.deeplearningbook.org/), Goodfellow et al., 2016
+    [2] PyTorch: https://pytorch.org/
+"""
+
+import copy
+import inspect
+import numpy as np
+import torch
+
+from dnn import NN
+
+__author__ = "Brian Delhaisse"
+__copyright__ = "Copyright 2018, PyRoboLearn"
+__credits__ = ["Brian Delhaisse"]
+__license__ = "(c) Brian Delhaisse"
+__version__ = "1.0.0"
+__maintainer__ = "Brian Delhaisse"
+__email__ = "briandelhaisse@gmail.com"
+__status__ = "Development"
+
+
+class RCNN(NN):
+    r"""Recurrent CNN
+    """
+    pass
+
+class RCNN(NNTorch):
+    r"""Recurrent CNN in PyTorch
+    """
+    pass

pyrobolearn/models/nn/rnn.py

@@ -0,0 +1,46 @@
+#!/usr/bin/env python
+"""Define the Recurrent Neural Network (RNN) learning model.
+
+This file provides the RNN model; a parametric, generally non-linear, recurrent, discriminative,
+and deterministic model. This model is convenient for sequential data. For instance, they can be used with language,
+where each word in a sentence is conditioned on the previous words.
+
+We decided to use the `pytorch` framework because of its popularity in the research community field, flexibility,
+similarity with numpy (but with automatic differentiation: autograd), GPU capabilities, and more Pythonic approach.
+While we thought about using other frameworks (such as Keras, Tensorflow, and others) as well, it would have
+unnecessarily complexify the whole framework, as these frameworks would not only influence the learning models,
+but also the losses, optimizers, and other modules. While we could have written some interfaces that makes the bridge
+between these various frameworks and ours, we came to the conclusion that this would take a considerable amount of
+efforts and time that we do not have for the moment.
+
+References:
+    [1] "Deep Learning" (http://www.deeplearningbook.org/), Goodfellow et al., 2016
+    [2] PyTorch: https://pytorch.org/
+"""
+
+import copy
+import inspect
+import numpy as np
+import torch
+
+from dnn import NN
+
+__author__ = "Brian Delhaisse"
+__copyright__ = "Copyright 2018, PyRoboLearn"
+__credits__ = ["Brian Delhaisse"]
+__license__ = "(c) Brian Delhaisse"
+__version__ = "1.0.0"
+__maintainer__ = "Brian Delhaisse"
+__email__ = "briandelhaisse@gmail.com"
+__status__ = "Development"
+
+
+class RNN(NN):
+    r"""Recurrent Neural Network
+    """
+    pass
+
+class RNNTorch(NNTorch):
+    r"""Recurrent Neural Network in PyTorch
+    """
+    pass

pyrobolearn/models/nn/vae.py

@@ -0,0 +1,76 @@
+#!/usr/bin/env python
+"""Define the Variational Auto-Encoder (VAE) learning model.
+
+This file provides the VAE model; a parametric, generally non-linear, non-recurrent, generative, and stochastic model.
+This model is a generative latent variable model which projects the input data into a latent lower
+dimensional space through an encoder, and re-projects it to the original data space through the use of the decoder.
+
+We decided to use the `pytorch` framework because of its popularity in the research community field, flexibility,
+similarity with numpy (but with automatic differentiation: autograd), GPU capabilities, and more Pythonic approach.
+While we thought about using other frameworks (such as Keras, Tensorflow, and others) as well, it would have
+unnecessarily complexify the whole framework, as these frameworks would not only influence the learning models,
+but also the losses, optimizers, and other modules. While we could have written some interfaces that makes the bridge
+between these various frameworks and ours, we came to the conclusion that this would take a considerable amount of
+efforts and time that we do not have for the moment.
+
+References:
+    [1] "Deep Learning" (http://www.deeplearningbook.org/), Goodfellow et al., 2016
+    [2] PyTorch: https://pytorch.org/
+"""
+
+import copy
+import inspect
+import numpy as np
+import torch
+
+from dnn import NN
+
+__author__ = "Brian Delhaisse"
+__copyright__ = "Copyright 2018, PyRoboLearn"
+__credits__ = ["Brian Delhaisse"]
+__license__ = "(c) Brian Delhaisse"
+__version__ = "1.0.0"
+__maintainer__ = "Brian Delhaisse"
+__email__ = "briandelhaisse@gmail.com"
+__status__ = "Development"
+
+
+class VAE(NN):
+    r"""Variational AutoEncoder
+
+    Type: generative model
+
+    .. seealso:: Generative Adversarial Networks
+
+    References:
+        [1] "Deep Learning" (http://www.deeplearningbook.org/), Goodfellow et al., 2016
+        [2] "Tutorial on Variational Autoencoder"
+    """
+    def __init__(self, layer_sizes, activation_fct=None, dropout=None):
+        self.encoder = None
+        self.decoder = None
+
+    @staticmethod
+    def isDiscriminative():
+        """A neural network is a discriminative model which given inputs predicts some outputs"""
+        return True
+
+    @staticmethod
+    def isGenerative():  # unless VAE, GAN,...
+        """Standard neural networks are not generative, and thus we can not sample from it. This is different,
+        for instance, for generative adversarial networks (GANs) and variational auto-encoders (VAEs)."""
+        return True
+
+    ###########
+    # Methods #
+    ###########
+
+    def sample(self, size=None, seed=None):
+        """Sample from the VAE"""
+        pass
+
+
+class VAETorch(NNTorch):
+    r"""Variational AutoEncoder in PyTorch
+    """
+    pass

pyrobolearn/models/pca.py

@@ -0,0 +1,273 @@
+#!/usr/bin/env python
+"""Define the PCA model.
+"""
+
+import numpy as np
+
+__author__ = "Brian Delhaisse"
+__copyright__ = "Copyright 2018, PyRoboLearn"
+__credits__ = ["Brian Delhaisse"]
+__license__ = "MIT"
+__version__ = "1.0.0"
+__maintainer__ = "Brian Delhaisse"
+__email__ = "briandelhaisse@gmail.com"
+__status__ = "Development"
+
+
+class PCA(object):
+    r"""Principal Component Analysis (PCA)
+
+    This class describes PCA; a non-parametric, linear model which possesses the 3 following properties [1]:
+    * linearity
+    * orthogonality
+    * high signal-noise ratio
+
+    PCA can be achieved by performing Singular Value Decomposition (SVD) on the data, or performing an
+    eigendecomposition on the covariance data matrix.
+
+    Assuming the mean-centered data matrix is given by :math:`X \in \mathbb{R}^{N \times M}`, applying SVD on it
+    gives us:
+    .. math:: X = USV^T,
+    where :math:`U \in \mathbb{R}^{N \times N}` is an orthogonal matrix where its columns represent the eigenvectors
+    of :math:`XX^T` also known as the left-singular vectors of :math:`X`, :math:`S \in \mathbb{R}^{N \times M}`
+    contains the singular values ordered by descending order, and :math:`V \in \mathbb{R}^{M \times M}` is an
+    orthogonal matrix in which its columns represent the eigenvectors of :math:`X^TX` also known as the right-singular
+    vectors of :math:`X`. The columns of :math:`V` form a basis spanning :math:`\mathbb{R}^M`.
+
+    Applying eigendecomposition on the covariance matrix :math:`C_X \sim X^TX` gives us:
+    .. math:: X^TX = QLQ^T.
+
+    While applying SVD on :math:`X` and then computing the covariance gives us:
+    .. math:: X^TX = (USV^T)^T (USV^T) = VSSV^T = VS^2V^T
+
+    Thus, :math:`Q=V` (i.e. same orthogonal matrix containing the evecs) and :math:`L=S^2` (that is, the eigenvalues
+    are the square of the singular values). Note that the covariance :math:`C_X` is formally defined as
+    :math:`\frac{1}{(N-1)} X^TX` where :math:`N` is the number of data points, and not :math:`X^TX`, then
+    :math:`L=\frac{S^2}{(N-1)}`.
+
+    PCA can be formulated as an optimization process, which consists to find the dimensions that maximize
+    the projected variance. Specifically, this can be written as:
+
+    .. math::
+        max_{u_i} ||Xu_i||^2  & \mbox{ subj. to } u_i^Tu_i = 1; u_j^Tu_i = 0  \\
+        max_{u_i} (Xu_i)^T(Xu_i) & \mbox{ subj. to } u_i^Tu_i = 1; u_j^Tu_i = 0 \\
+        max_{u_i} u_i^TX^TXu_i & \mbox{ subj. to } u_i^Tu_i = 1; u_j^Tu_i = 0 \\
+
+    :math:`\forall i, \forall j < i`. The solution is given by the column vectors of the matrix :math:`V`, that is,
+    the eigenvectors, while the obtained values during the maximization process are given by the eigenvalues.
+    This is why PCA can be performed by applying SVD (on the data matrix) or Eigendecomposition (on the covariance
+    matrix).
+
+    Complexity:
+        * Spatial complexity: O(NM)
+        * Time complexity: O(min(NM^2, MN^2))
+
+    References:
+        [1] "A Tutorial on Principal Component Analysis", Shlens, 2014
+    """
+    def __init__(self, X=None, normalize_data=False):
+        if X is not None:
+            self.train(X, normalize=normalize_data)
+
+        self.evals, self.evecs = None, None
+
+    ##############
+    # Properties #
+    ##############
+
+    @property
+    def eigenvalues(self): # alias to evals
+        return self.evals
+
+    @property
+    def eigenvectors(self): # alias to evecs
+        return self.evecs
+
+    ##################
+    # Static Methods #
+    ##################
+
+    # TODO: think if PCA can be considered as a model
+
+    @staticmethod
+    def isParametric():
+        """PCA is a non-parametric approach"""
+        return False
+
+    @staticmethod
+    def isLinear():
+        """PCA does not have parameters, but it is a linear dimensionality reduction algo"""
+        return True
+
+    @staticmethod
+    def isRecurrent():
+        """PCA is not recurrent"""
+        return False
+
+    @staticmethod
+    def isLatent():
+        """PCA gives a latent model"""
+        return True
+
+    @staticmethod
+    def isProbabilistic():
+        """PCA is not a probabilistic approach but a deterministic one"""
+        return False
+
+    @staticmethod
+    def isDiscriminative():
+        """PCA is a discriminative model, which projects the given data into a lower space"""
+        return True
+
+    @staticmethod
+    def isGenerative():
+        """PCA is not a generative model from which you can sample from it"""
+        return False
+
+    ###########
+    # Methods #
+    ###########
+
+
+    def parameters(self):
+        raise RuntimeError("PCA doesn't have any parameters.")
+
+    def getParams(self):
+        raise RuntimeError("PCA doesn't have any parameters.")
+
+    def getHyperparams(self):
+        pass
+
+    def train(self, X, normalize=False, copy=True):
+        """
+        Compute PCA on the given data. This method will center the data.
+
+        Args:
+            X (float[N,D]): the matrix to apply PCA on (with shape(X)=NxD, where `N` is the nb of samples, and
+                `D` is the dimensionality of a sample).
+            normalize (bool): if True, it will normalize the data using the std dev. PCA will then be applied on
+                the correlation matrix instead of the covariance matrix.
+            copy (bool): if True, it will first copy the data before centering it, and possibly normalizing it.
+
+        Return:
+            float[D]: the sorted eigenvalues
+            float[D]: the sorted eigenvectors
+        """
+        if copy:
+            X = np.copy(X)
+
+        # 1. Center the data
+        mean = X.mean(axis=0)
+        X -= mean
+        N = X.shape[0]
+
+        # Normalize using the std dev
+        if normalize:
+            X /= X.std(axis=0)
+
+        # 2. Compute the covariance/correlation matrix
+        CovX = 1./(N-1) * X.T.dot(X) # TxT (same as np.cov(X, rowvar=False)))
+
+        # 3. Compute the eigenvectors of this covariance matrix
+        # np.linalg.eigh is more efficient than np.linalg.eig for symmetric matrix
+        evals, evecs = np.linalg.eigh(CovX)
+
+        # 4. Sort the eigenvalues (in decreasing order) and eigenvectors
+        idx = np.argsort(evals)[::-1]
+        evals, evecs = evals[idx], evecs[:,idx]
+
+        # save values
+        self.evals = evals
+        self.evecs = evecs
+
+        return evals, evecs
+
+    def predict(self, x):
+        # TODO
+        pass
+
+
+class RecursivePCA(PCA):
+    r"""Recursive PCA
+
+    Given a new data point, the
+
+    Assume the covariance matrix is given by :math:`C`, and the mean by :math:`m`, then the new covariance matrix
+    accounting for the new data point is given by:
+
+    .. math:: C' = C + \frac{N}{N+1} (m m^T - m x'^T - x' m^T + x'x'^T)
+
+    where :math:`C` can be reconstructed from the eigenvalues and eigenvectors using the eigendecomposition:
+
+    .. math:: C = V \Lambda V^T
+
+    and the new mean is given by:
+
+    .. math::
+
+    Finally, the total number of data points N is updated. This reduces the the spatial and time complexities to
+    O(T^2) and O(T^3), respectively.
+
+    Time complexity:
+    * O(min(NT^2, TN^2))
+    """
+
+    def __init__(self, X=None, normalize_data=False):
+        super(RecursivePCA, self).__init__(X, normalize_data)
+
+    def train_recursive(self, X):
+        if self.X is None:
+            # apply std PCA
+            self.X = self.train(X)
+            self.mean = np.mean(X, axis=0).reshape(-1, 1)
+            self.cov = np.cov(X, rowvar=False)
+            self.N = len(X)
+
+        else:  # recursive
+            for x in X:
+                x = x.reshape(-1, 1)
+                Y = self.mean.dot(x.T)
+                self.cov = self.cov + self.N / (self.N + 1.) * (self.mean.dot(self.mean.T) - Y - Y.T
+                                                            + self.mean.dot(self.mean.T))
+                self.mean = self.N / (self.N + 1.) * self.mean + 1 / (N+1) * x
+
+                self.N += 1
+
+
+class HierarchicalPCA(PCA):
+    r"""Hierarchical PCA
+
+    Assume :math:`X = [X_1, X_2] \in \mathbb{R}^{N \times 2D}`, where :math:`X_1, X_2 \in \mathbb{R \times D}`.
+
+    We first decompose the covariance matrix :math:`C` of :math:`X` in terms of :math:`X_1` and :math:`X_2`, and
+    apply SVD on each term as follows:
+
+    .. math::
+
+        C = X^T X
+        = \left[ \begin{array}{cc}
+            X_1^T X_1 & X_1^T X_2 \\
+            X_2^T X_1 & X_2^T X_2
+        \end{array} \right]
+        = \left[ \begin{array}{cc}
+            V_1 \Lambda_1 V_1^T & V_1 \Sigma_1 U_1^T U_2 \Sigma_2 V_2^T \\
+             V_2 \Sigma_2 U_2^T U_1 \Sigma_1 V_1^T & V_2 \Lambda_2 V_2^T
+        \end{array} \right]
+
+    Similarly, we apply the eigendecomposition on C which results in:
+
+    .. math::
+
+        C = X^T X = V \Lambda V^T
+        = \left[ \begin{array}{cc} V_{11} & V_{12} \\ V_{21} & V_{22} \end{array} \right]
+        \left[ \begin{array}{cc} \Lambda_{11} & 0 \\ 0 & \Lambda_{22} \end{array} \right]
+        \left[ \begin{array}{cc} V_{11}^T & V_{21}^T \\ V_{12}^T & V_{22}^T \end{array} \right]
+        = \left[ \begin{array}{cc}
+            V_{11} \Lambda_{11} V_{11}^T + V_{12} \Lambda_{22} V_{12}^T
+            & V_{11} \Lambda_{11} V_{21}^T + V_{12} \Lambda_{22} V_{22}^T \\
+            V_{21} \Lambda_{11} V_{11}^T + V_{22} \Lambda_{22} V_{12}^T
+            & V_{21} \Lambda_{11} V_{21}^T + V_{22} \Lambda_{22} V_{22}^T
+        \end{array} \right]
+
+    """
+    pass

pyrobolearn/models/polynomial.py

@@ -0,0 +1,292 @@
+#!/usr/bin/env python
+"""Define the polynomial learning model.
+
+The polynomial model is a discriminative deterministic model given by: :math:`y = f(x) = W \phi(x)`, where
+:math:`\phi` is a function that returns a transformed input vector (possibly of higher dimension).
+"""
+
+import copy
+# import inspect
+import collections
+import numpy as np
+import torch
+
+__author__ = "Brian Delhaisse"
+__copyright__ = "Copyright 2018, PyRoboLearn"
+__credits__ = ["Brian Delhaisse"]
+__license__ = "MIT"
+__version__ = "1.0.0"
+__maintainer__ = "Brian Delhaisse"
+__email__ = "briandelhaisse@gmail.com"
+__status__ = "Development"
+
+
+class Polynomial(object):
+    r"""Polynomial model
+
+    The polynomial model is a discriminative deterministic model expressed mathematically as
+    :math:`y = f(x) = W \phi(x)`, where :math:`x` is the input vector, :math:`y` is the output vector, :math:`W`
+    is the weight matrix, and :math:`\phi` is the polynomial function which returns the transformed input vector.
+    This transformed input vector is often of higher dimension, based on the idea that if it is not linear with
+    respect to the parameters in the current space, it might be in a higher dimensional space.
+    """
+
+    def __init__(self, num_inputs, num_outputs, polynomial_fct):
+        """
+        Initialize the polynomial model: :math:`y = W \phi(x)`
+
+        Args:
+            num_inputs (int): dimension of the input vector x
+            num_outputs (int): dimension of the output vector y
+            polynomial_fct (callable, PolynomialFunction): polynomial function to be applied on the input vector x.
+                It has to be callable, returns the output vector :math:`\phi(x)` when called, and has a `size`
+                attribute which returns the number of exponents in the polynomial.
+        """
+        self.phi = polynomial_fct
+        num_inputs = num_inputs * self.phi.size
+        self.model = torch.nn.Linear(num_inputs, num_outputs, bias=False)
+        self._num_parameters = len(self.get_vectorized_parameters())
+
+    ##############
+    # Properties #
+    ##############
+
+    @property
+    def phi(self):
+        """Return the polynomial function"""
+        return self._phi
+
+    @phi.setter
+    def phi(self, fct):
+        """Set the polynomial function"""
+        if not callable(fct):
+            raise TypeError("Expecting the polynomial function to be callable.")
+
+        if not hasattr(fct, 'size'):
+            raise ValueError("Expecting the polynomial function to have the 'size' attribute.")
+
+        # if len(inspect.getargspec(fct).args) < 1:
+        #     raise TypeError("Expecting the polynomial function to accept an argument (the state).")
+
+        self._phi = fct
+
+    @property
+    def polynomial_function(self):
+        """Return the polynomial function"""
+        return self._phi
+
+    @polynomial_function.setter
+    def polynomial_function(self, fct):
+        """Set the polynomial function"""
+        self.phi = fct
+
+    @property
+    def input_dims(self):
+        """Return the input dimension of the model"""
+        return self.model.weight.shape[1]
+
+    @property
+    def output_dims(self):
+        """Return the output dimension of the model"""
+        return self.model.weight.shape[0]
+
+    @property
+    def input_shape(self):
+        """Return the input shape of the model"""
+        return tuple([self.input_dims])
+
+    @property
+    def output_shape(self):
+        """Return the output shape of the model"""
+        return tuple([self.output_dims])
+
+    @property
+    def num_parameters(self):
+        """Return the total number of parameters"""
+        return self._num_parameters
+
+    ##################
+    # Static methods #
+    ##################
+
+    @staticmethod
+    def copy(other):
+        """Return another copy of the polynomial model"""
+        if not isinstance(other, Polynomial):
+            raise TypeError("Trying to copy an object which is not a Polynomial model")
+        return copy.copy(other)
+
+    @staticmethod
+    def is_parametric():
+        """The polynomial model is a parametric model"""
+        return True
+
+    @staticmethod
+    def is_linear():
+        """The polynomial model is linear with respect to its parameters"""
+        return True
+
+    @staticmethod
+    def is_recurrent():
+        """The polynomial model is not recurrent; current outputs do not depend on previous inputs"""
+        return False
+
+    @staticmethod
+    def is_deterministic():
+        """The polynomial model is a deterministic model"""
+        return True
+
+    @staticmethod
+    def is_probabilistic():     # same as is_stochastic()
+        """The polynomial model is not a probabilistic model; it is a deterministic one"""
+        return False
+
+    @staticmethod
+    def is_discriminative():
+        """The polynomial model is a discriminative model."""
+        return True
+
+    @staticmethod
+    def is_generative():
+        """The polynomial model is not a generative model, and thus we can not sample from it"""
+        return False
+
+    ###########
+    # Methods #
+    ###########
+
+    def parameters(self):
+        """Returns an iterator over the model parameters."""
+        return self.model.parameters()
+
+    def named_parameters(self):
+        """Returns an iterator over the model parameters, yielding both the name and the parameter itself"""
+        return self.model.named_parameters()
+
+    def list_parameters(self):
+        """Return a list of parameters"""
+        return list(self.parameters())
+
+    def get_vectorized_parameters(self, to_numpy=True):
+        """Return a vectorized form (1 dimensional array) of the parameters."""
+        parameters = self.parameters()
+        vector = torch.cat([parameter.reshape(-1) for parameter in parameters])  # np.concatenate = torch.cat
+        if to_numpy:
+            return vector.detach().numpy()
+        return vector
+
+    def set_vectorized_parameters(self, vector):
+        """Set the vector parameters."""
+        # convert the vector to torch array
+        if isinstance(vector, np.ndarray):
+            vector = torch.from_numpy(vector).float()
+
+        # set the parameters from the vectorized one
+        idx = 0
+        for parameter in self.parameters():
+            size = parameter.nelement()
+            parameter.data = vector[idx:idx+size].reshape(parameter.shape)
+            idx += size
+
+    def predict(self, x, to_numpy=True):
+        """
+        Predict output vector :math:`y` given input vector :math:`x`, using the formula: :math:`y = W \phi(x)`.
+
+        Args:
+            x (np.ndarray, torch.Tensor): input vector
+            to_numpy (bool): if True, return a np.array
+
+        Returns:
+            np.ndarray, torch.Tensor: output vector
+        """
+        # convert from numpy to pytorch if necessary
+        if isinstance(x, np.ndarray):
+            x = torch.from_numpy(x).float()
+
+        # predict the output
+        y = self.model(self.phi(x))
+
+        # return the output and convert it if necessary
+        if to_numpy:
+            if y.requires_grad:
+                return y.detach().numpy()
+            return y.numpy()
+        return y
+
+    def __call__(self, x, to_numpy=True):
+        return self.predict(x, to_numpy=to_numpy)
+
+
+class PolynomialFunction(object):
+    r"""Polynomial function
+
+    Polynomial function to be applied on the given inputs.
+    """
+
+    def __init__(self, degree=1):
+        """
+        Initialize the polynomial function.
+
+        Args:
+            degree (int, list of ints): degree(s) of the polynomial. Setting `degree=3`, will return `[1,x,x^2,x^3]`
+                as output, while setting `degree=[1,3]` will return `[x,x^3]` as output.
+        """
+        self.degree = degree
+
+    @property
+    def degree(self):
+        """Return the degree of the polynomial"""
+        return self._degree
+
+    @degree.setter
+    def degree(self, degree):
+        """Set the degree of the polynomial"""
+        # checks
+        if isinstance(degree, int):
+            degree = range(degree + 1)
+        elif isinstance(degree, collections.Iterable):
+            for d in degree:
+                if not isinstance(d, int):
+                    raise TypeError("Expecting the given degrees to be positive integers, but got {}".format(type(d)))
+                if d < 0:
+                    raise ValueError("Expecting the given degrees to be positive integers, but got {}".format(d))
+        else:
+            raise TypeError("Expecting the degree to be a positive integer or a list of positive integers.")
+
+        self._degree = degree
+
+    @property
+    def size(self):
+        """Return the number of exponents. The output vector is then of size = size * len(x)"""
+        return len(self.degree)
+
+    def predict(self, x):
+        """Return output polynomial vector"""
+        if isinstance(x, np.ndarray):
+            return np.concatenate([x**d for d in self.degree])
+        elif isinstance(x, torch.Tensor):
+            return torch.cat([x**d for d in self.degree])
+
+    def __call__(self, x):
+        return self.predict(x)
+
+
+# Tests
+if __name__ == '__main__':
+    # test with numpy
+    x = np.array(range(3))
+    fct = PolynomialFunction(degree=3)
+    model = Polynomial(num_inputs=len(x), num_outputs=2, polynomial_fct=fct)
+    print("Polynomial input size: {}".format(model.input_dims))
+    print("Polynomial output size: {}".format(model.output_dims))
+    y = model(x)
+    print("Polynomial input: {}".format(x))
+    print("Polynomial output: {}".format(y))
+
+    # test with pytorch
+    x = torch.from_numpy(x).float()
+    y = model(x, to_numpy=False)
+    print("Polynomial input: {}".format(x))
+    print("Polynomial torch output: {}".format(y))
+    y = model(x, to_numpy=True)
+    print("Polynomial numpy output: {}".format(y))

pyrobolearn/optim/__init__.py

@@ -0,0 +1,3 @@
+
+# import optimizers
+from optimizer import *

pyrobolearn/optim/cio.py

@@ -0,0 +1,65 @@
+# This file implements the 'Contact Invariant Optimization' framework developed by Igor Mordatch.
+# Ref: "Automated Discovery and Learning of Complex Movement Behaviors" (PhD thesis), Mordatch, 2015
+# See also: presentation given CS294
+
+import numpy as np
+from scipy.interpolate as interp1d
+
+
+class CIO(object):
+    r"""
+    The Contact Invariant Optimization (CIO) algorithm [1] consists to minimize the following cost:
+
+    .. math::
+
+        s* = argmin_s L(s)
+           = argmin_s L_{CI}(s) + L_{physics}(s) + L_{task}(s) + L_{hint}(s)
+
+    where :math:`s` is the state which contains :math:`x_k`, :math:`\dot{x}_k`, and :math:`c_k` for each phase/interval
+    :math:`k`. The vector :math:`x_k` contains the torso and end-effector position and orientation, while
+    the :math:`c_k` vector represents the auxiliary contact variables.
+
+    Here is what each term represents in the total cost:
+    - :math:`L_{CI}` is the contact invariant cost
+    - :math:`L_{physics}` penalizes physics violation
+    - :math:`L_{task}` describes the task objectives (i.e. high-level goals of the movement)
+    - :math:`L_{hint}` provides hints to accelerate the optimization. This term is optional.
+
+    The CIO consists of 3 phases:
+    1. only :math:`L_{task}` is enabled
+    2. All 4 terms (:math:`L_{task}`, :math:`L_{physics}`, :math:`L_{CI}`, :math:`L_{hint}`) are enabled but with
+       :math:`L_{physics}` down-weighted by 0.1
+    3. :math:`L_{task}`, :math:`L_{physics}`, and :math:`L_{CI}` are fully enabled
+
+    Note that the solution obtained at the end of each phase is perturbed with small zero-mean Gaussian noise to
+    break any symmetries, and used to initialize the next phase.
+
+    From the optimized state :math:`s^*`, the optimal joints :math:`q^*` at each time step can be computed (using IK).
+    Then, a PD controller can be used to move the joints to their desired configuration.
+
+    Note that the framework do not take into account any sensory feedbacks.
+
+    References:
+        [1] "Automated Discovery and Learning of Complex Movement Behaviors" (PhD thesis), Mordatch, 2015
+    """
+
+    def __init__(self, robot, T, num_interval=20):
+        # TODO: think about optimizing multiple actors
+        self.robot = robot
+        self.K = num_interval
+
+        x = np.linspace(0, T, self.K)
+        y = np.array(range(1, self.K+1))
+        self.phase = scipy.interpolate.interp1d(x, y, kind='zero')
+
+    def get_phase_index(self, t):
+        return self.phase(t)
+
+    def compute_state(self):
+        base_pos = self.robot.getBasePosition()
+        base_quat = self.robot.getBaseOrientation()
+        end_effector_pos = self.robot.getEndEffectorPositions()
+        end_effector_quat = self.robot.getEndEffectorOrientations()
+
+    def optimize(self):
+        pass

pyrobolearn/optim/optimizer.py

@@ -0,0 +1,786 @@
+#!/usr/bin/env python
+"""Provide various optimizers.
+
+Optimizers allows to optimize (i.e. minimize or maximize) a utility function (also known as objective function,
+fitness function, loss, etc.) with or without constraints and bounds. Notably, it assumes the functions have some
+parameters that the optimizer can update.
+
+Mathematically, this is described as:
+
+.. math:: \min_{x \in R^n} f(x)
+
+subject to
+
+.. math::
+
+    g_i(x) \geq 0, \quad i = 1,...,m
+    h_j(x)  =   0, \quad j = 1,...,p
+    x_l \leq x \leq x_u
+
+For trajectory optimization, check "An Introduction to Trajectory Optimization: How to do your own Direct Collocation".
+"""
+
+# TODO: trajectory optimization
+
+from abc import ABCMeta, abstractmethod
+
+# Numpy with autograd
+import autograd.numpy as np  # Thinly-wrapped numpy
+from autograd import grad    # The only autograd function you may ever need
+
+# Scipy optimizer
+import scipy
+
+# Pytorch optimizers
+import torch.nn as nn
+import torch.optim as optim
+
+# NLopt optimizers
+try:
+    import nlopt
+except ImportError as e:
+    raise ImportError(e.__str__() + "\n HINT: you can install nlopt via `pip install nlopt`.")
+
+# IPopt optimizer
+try:
+    import ipopt
+except ImportError as e:
+    raise ImportError(e.__str__() + "\n HINT: you can install ipopt using the `pyrobolearn/scripts/install_ipopt.sh`."
+                                    "If ipopt is already installed, you can install the python wrapper via "
+                                    "`pip install ipopt`.")
+
+# CVXOPT
+# import cvxopt
+# CVXPY: nice wrapper around cvxopt
+# import cvxpy
+# Quadprog
+# import quadprog
+
+# QPsolvers optimizers: unified Python interface for multiple QP solvers (cvxopt, cvxpy, quadprog,...)
+try:
+    import qpsolvers
+except ImportError as e:
+    raise ImportError(e.__str__() + "\n HINT: you can install qpsolvers directly via 'pip install qpsolvers'.")
+
+# Bayesian optimization
+try:
+    import GPy
+    import GPyOpt
+except ImportError as e:
+    raise ImportError(e.__str__() + "\n HINT: you can install GPy/GPyOpt directly via 'pip install GPy' and "
+                                    "'pip install GPyOpt'.")
+
+# CMA-ES
+try:
+    import cma
+except ImportError as e:
+    raise ImportError(e.__str__() + "\n HINT: you can install CMA-ES or `pycma` directly via 'pip install cma'.")
+
+
+__author__ = "Brian Delhaisse"
+__copyright__ = "Copyright 2018, PyRoboLearn"
+__credits__ = ["Brian Delhaisse"]
+__license__ = "MIT"
+__version__ = "1.0.0"
+__maintainer__ = "Brian Delhaisse"
+__email__ = "briandelhaisse@gmail.com"
+__status__ = "Development"
+
+
+##################################################################
+#                           OPTIMIZER                            #
+##################################################################
+
+
+class Optimizer(object):
+    r"""Optimizer abstract class
+
+    This is an abstract class from which all optimizers inherit from. Most of the child optimizer classes are wrappers
+    around the original optimizer. This is to provide the same common interface to all the optimizers, and convert
+    seamlessly to the correct data types.
+
+    Optimizers are often given as a parameter to the learning algorithms, but can also be used of out of the box
+    directly on models. Several original optimizers can be given to the learning algorithm which will automatically
+    wrap the optimizer with the corresponding wrapper to provide a common interface.
+
+    In their most natural form, optimizers are used to ... optimization process which can be described mathematically
+    by:
+
+    .. math::
+
+        \min_{\theta} J(\theta) \mbox{ subj. to constraints}
+
+        \max_{\theta} J(\theta)
+
+    Several optimizers provides
+
+    The list of optimizers available are from the following libraries:
+    * nlopt
+    * ipopt
+    * torch.optim
+    * GPy / GPyOpt
+    * scipy.optimize
+    * qpsolvers (which includes cvxpy)
+    * cmaes
+    * pso
+
+    Each one of them expect a certain kind of type of the parameters.
+
+    Optimizers can be divided into:
+    * global vs local
+    * derivative-free vs gradient-based
+    * without constraints vs with (equality and/or inequality) constraints
+
+    Many implemented optimizers that can be found online are specific to a certain type of learning model.
+    If necessary, a conversion or wrapping process is carried out to make the optimizer work with the given learning
+    model.
+    """
+    __metaclass__ = ABCMeta
+
+    def __init__(self, model, losses, hyperparameters):
+        """
+
+        :param model: a certain type of model. If the original model is given instead of an instance of `Model`, then
+                      the model will be wrapped appropriately.
+        :param losses:
+        :param hyperparameters:
+        """
+        pass
+
+
+##################################################################
+#                             SCIPY                              #
+##################################################################
+
+class Scipy(object):
+    r"""Scipy optimizer
+
+    This uses the `scipy.optimize.minimize` to optimize a given objective function under various bounds and
+    constraints. Specifically, it consists of the minimization of a scalar function of one or more variables.
+    In general, the optimization problems are of the form:
+
+    .. math::
+
+        \min_{x \in R^n} f(x)
+
+    subject to
+
+    .. math::
+
+        g_i(x) \geq 0, \quad i = 1,...,m
+        h_j(x)  =   0, \quad j = 1,...,p
+
+    where :math:`x` is a vector of one or more variables, :math:`g_i(x)` are the inequality constraints, and
+    :math:`h_j(x)` are the equality constrains.
+
+    Optionally, the lower and upper bounds for each element in :math:`x` can also be specified using the `bounds`
+    argument.
+
+    Several methods/optimizers are available:
+    -
+
+    Note that only 'COBYLA' and 'SLSQP' support constraints, where the former only supports inequality constraints.
+
+    References:
+        [1] scipy.optimize.minimize: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html
+    """
+
+    def __init__(self, method='SLSQP'):
+        """
+        Initialize the scipy method
+
+        Args:
+            method (str, callable):
+                - 'Nelder-Mead' :ref:`(see here) <scipy.optimize.minimize-neldermead>`
+                - 'Powell'      :ref:`(see here) <scipy.optimize.minimize-powell>`
+                - 'CG'          :ref:`(see here) <scipy.optimize.minimize-cg>`
+                - 'BFGS'        :ref:`(see here) <scipy.optimize.minimize-bfgs>`
+                - 'Newton-CG'   :ref:`(see here) <scipy.optimize.minimize-newtoncg>`
+                - 'L-BFGS-B'    :ref:`(see here) <scipy.optimize.minimize-lbfgsb>`
+                - 'TNC'         :ref:`(see here) <scipy.optimize.minimize-tnc>`
+                - 'COBYLA'      :ref:`(see here) <scipy.optimize.minimize-cobyla>`
+                - 'SLSQP'       :ref:`(see here) <scipy.optimize.minimize-slsqp>`
+                - 'dogleg'      :ref:`(see here) <scipy.optimize.minimize-dogleg>`
+                - 'trust-ncg'   :ref:`(see here) <scipy.optimize.minimize-trustncg>`
+                - custom - a callable object (added in version 0.14.0),
+
+        """
+        # define optimization method
+        # By default, it will be 'BFGS', 'L-BFGS-B', or 'SLSQP' depending on the constraints and bounds
+        # If constraints, it can only be 'COBYLA' or 'SLSQP'. COBYLA only supports inequality constraints.
+        self.method = method
+
+    def optimize(self, maxiter=1e6, verbose=True):
+        # define objective function to MINIMIZE
+        # f = lambda x: -(x.T.dot(C)).dot(x)
+        def f(x):
+            return -(x.T.dot(C)).dot(x)
+
+        # define initial guess
+        x0 = np.ones((M,))  # np.zeros((M,))
+
+        # define 1st constraints: norm of 1
+        constraints = [{'type': 'eq', 'fun': lambda x: x.T.dot(x) - 1, 'jac': None, 'args': ()}]
+
+        # define bounds: each vector u have a norm of 1 thus each parameter is between -1 and 1
+        bounds = [(-1., 1.)] * M
+
+        # optimize recursively
+        evals, evecs = [], []
+        messages = {}
+        options = {'maxiter': maxiter, 'disp': verbose}
+        for i in range(M):
+            if i != 0:
+                # add orthogonality constraint
+                constraints.append({'type': 'eq', 'fun': lambda u: u1.T.dot(u)})
+
+            # minimize --> it returns an instance of OptimizeResult
+            result = scipy.optimize.minimize(f, x0, args=(), method=self.method, jac=None, hess=None, bounds=bounds,
+                                         constraints=constraints, tol=None, callback=None, options=options)
+
+            print(result.success)
+            print(result.message)
+            print(result.fun)
+            print(result.x)
+
+
+##################################################################
+#                    Quadratic Programming                       #
+##################################################################
+
+# class CVXOPT(Optimizer):
+#     r"""Convex Optimizer
+#
+#     Note: cvxpy module is a nice wrapper around cvxopt that follows paradigm of a disciplined convex programming.
+#
+#     References:
+#         [1] Python Software for Convex Optimization: https://cvxopt.org/
+#         [2] Github repo: https://github.com/cvxopt/cvxopt
+#     """
+#     pass
+#
+#
+# class CVXPY(Optimizer):
+#     r"""Convex Optimizer
+#
+#     References:
+#         [1] CVXPY: http://www.cvxpy.org/
+#         [2] Github repo: https://github.com/cvxgrp/cvxpy
+#     """
+#     pass
+#
+#
+# class QuadProg(object):
+#     r"""Quadprog
+#
+#     References:
+#         [1] Github repo: https://github.com/rmcgibbo/quadprog
+#     """
+#     pass
+
+class QP(object):
+    r"""Quadratic Programming solvers
+
+    This class uses the `qpsolvers` which is a unified Python interface for multiple QP solvers [1,2].
+
+    .. math::
+
+        \min_{x \in R^n} \frac{1}{2} x^T P x + q^T x
+
+    subject to
+
+    .. math::
+
+        Gx \leq h
+        Ax = b
+
+    where :math:`x` is the vector of optimization variables, the matrix :math:`P` and vector :math:`q` are used to
+    define any quadratic objective function on these variables, while the matrix-vector couples :math:`(G,h)` and
+    :math:`(A,b)` respectively define inequality and equality constraints. Vector inequalities apply coordinate by
+    coordinate [1].
+
+    - Dense solvers:
+        - CVXOPT
+        - CVXPY
+        - qpOASES
+        - quadprog
+    - Sparse solvers:
+        - ECOS as wrapped by CVXPY
+        - Gurobi
+        - MOSEK
+        - OSQP
+
+    Check the available solvers by calling `print(qpsolvers.available_solvers)`.
+
+    Notes: Many solvers (including CVXOPT, OSQP and quadprog) assume that `P` is a symmetric matrix, and may return
+    erroneous results when that is not the case. You can set ``sym_proj=True`` to project `P` on its symmetric part,
+    at the cost of some computation time.
+
+    References:
+        [1] QP in Python: https://scaron.info/blog/quadratic-programming-in-python.html
+        [2] Github repo: https://github.com/stephane-caron/qpsolvers
+    """
+
+    def __init__(self, method='quadprog'):
+        """
+        Initialize the QP solver.
+
+        Args:
+            method (str): ['cvxopt', 'cvxpy', 'ecos', 'gurobi', 'mosek', 'osqp', 'qpoases', 'quadprog']
+        """
+        solvers = set(qpsolvers.available_solvers)
+        if len(solvers) == 0:
+            raise ValueError("No QP solvers have been found on this computer. Please install one of the QP modules")
+        if method not in solvers:
+            method = 'quadprog'
+        self.method = method
+
+        # check methods that require a symmetric matrix for P
+        methods = ['cvxopt', 'osqp', 'quadprog']
+        self.sym_proj = True if self.method in set(methods) else False
+
+    def is_symmetric(self, X, tol=1e-8):
+        return np.allclose(X, X.T, atol=tol)
+
+    def optimize(self, P, q, x0=None, G=None, h=None, A=None, b=None):
+        return qpsolvers.solve_qp(P, q, G, h, A, b, solver=self.method, initvals=x0, sym_proj=self.sym_proj)
+
+
+##################################################################
+#                            NLOPT                               #
+##################################################################
+
+class NLopt(object):
+    r"""Non-Linear Optimizer
+
+    Non-linear optimizers based on the `nlopt` libraries.
+
+    Here is a brief of lists of the current algorithms implemented:
+    *
+
+    Nonlinear optimization algos that can handle nonlinear inequality and EQUALITY constraints are:
+    - ISRES (Improved Stochastic Ranking Evolution Strategy) --> global derivative-free
+    - COBYLA (Constrained Optimization BY Linear Approximations) --> local derivative-free
+    - SLSQP (Sequential Least-SQuares Programming) --> local gradient-based
+    - AUGLAG (AUGmented LAGrangian) --> global/local derivative-free/gradient based (determined based on the
+        subsidiary algo)
+
+    More information about:
+    - algorithms: https://nlopt.readthedocs.io/en/latest/NLopt_Algorithms/
+
+    References:
+        [1] NLopt: https://nlopt.readthedocs.io/en/latest/
+        [2] NLopt with Python: with Python: https://nlopt.readthedocs.io/en/latest/NLopt_Python_Reference/
+        [3] Github repo: https://github.com/stevengj/nlopt
+    """
+
+    # def __init__(self, model, losses, hyperparameters, seed):
+    #     super(NLopt, self).__init__(model, losses, hyperparameters)
+
+    def __init__(self, method, submethod=None, seed=None):
+
+        # define useful variables
+        self.results = {1: 'success', 2: 'stop_val reached', 3: 'ftol reached', 4: 'xtol reached',
+                        5: 'maxeval reached', 6: 'maxtime reached', -1: 'failure', -2: 'invalid args',
+                        -3: 'out of memory', -4: 'roundoff limited', -5: 'forced stop'}
+
+        # define random seed
+        nlopt.srand(seed)
+
+        # define which solver to use
+        def get_opt(method):
+            if method == 'ISRES':
+                return nlopt.opt(nlopt.GN_ISRES, M)
+            elif method == 'COBYLA':
+                return nlopt.opt(nlopt.LN_COBYLA, M)
+            elif method == 'SLSQP':
+                return nlopt.opt(nlopt.LD_SLSQP, M)
+            elif method == 'AUGLAG':
+                return nlopt.opt(nlopt.AUGLAG, M)
+            else:
+                raise NotImplementedError("The given method has not been implemented")
+
+        if method is None:
+            method = 'SLSQP'
+        self.opt = get_opt(method)
+
+        # define subsolver to use (if we use the AUGLAG method)
+        if method == 'AUGLAG':
+            if submethod is None:
+                submethod = 'SLSQP'
+            elif submethod == 'AUGLAG':
+                raise ValueError("Submethod should be different from AUGLAG")
+            subopt = get_opt(submethod)
+            subopt.set_lower_bounds(-1)
+            subopt.set_upper_bounds(1)
+            # subopt.set_ftol_rel(1e-2)
+            # subopt.set_maxeval(100)
+            self.opt.set_local_optimizer(subopt)
+
+    def optimize(self):
+        # define objective function and its gradient
+        def f(x, grad):
+            if grad.size > 0:
+                grad[:] = 2 * x.T.dot(C)
+            return x.T.dot(C).dot(x)
+
+        # define objective function to maximize
+        self.opt.set_max_objective(f)
+
+        # if nlopt.GN_ISRES, we can define the population size
+        self.opt.set_population(0)  # by default for ISRES: pop=20*(M+1)
+
+        # define bound constraints (should be between -1 and 1 because the norm should be 1)
+        self.opt.set_lower_bounds(-1.)
+        self.opt.set_upper_bounds(1.)
+
+        # define norm constraint and its gradient
+        def c1(x, grad):
+            if grad.size > 0:
+                grad[:] = 2 * x
+            return (x.T.dot(x) - 1)
+
+        # define orthogonal constraint
+        class OrthogonalConstraint(object):
+
+            def __init__(self, v):
+                self.v = np.copy(v)
+
+            def constraint(self, x, grad):
+                if grad.size > 0:
+                    grad[:] = self.v
+                return (x.T.dot(self.v))
+
+        # define equality constraints
+        opt.add_equality_constraint(c1, 0)
+        # opt.add_equality_mconstraint(constraints, tol)
+
+        # define stopping criteria
+        # opt.set_stopval(stopval)
+        opt.set_ftol_rel(1e-8)
+        # opt.set_xtol_rel(1e-4)
+        opt.set_maxeval(100000)  # nb of iteration
+        opt.set_maxtime(2)  # time in secs
+
+        # define initial value
+        x0 = np.array([0.1] * M)  # important that the initial value != 0 for the computation of the grad!
+
+        evals, evecs, msgs = [], [], {}
+        for i in range(M):
+            # add constraint
+            if i > 0:
+                c = OrthogonalConstraint(x)
+                opt.add_equality_constraint(c.constraint, 0)
+
+            # optimize
+            try:
+                x = opt.optimize(x0)
+            except nlopt.RoundoffLimited as e:
+                pass
+
+            # save values
+            evecs.append(x)  # param vector
+            evals.append(opt.last_optimum_value())  # max value
+            msgs[i] = nlopt_results[opt.last_optimize_result()]
+
+
+##################################################################
+#                            IPOPT                               #
+##################################################################
+
+class NormConstraint(object):
+
+    def __init__(self):
+        pass
+
+    def constraint(self, x):
+        return x.T.dot(x)
+
+    def jacobian(self, x):
+        return 2 * x
+
+
+class OrthogonalConstraint(object):
+
+    def __init__(self, v):
+        self.v = np.copy(v)
+
+    def constraint(self, x):
+        return x.T.dot(self.v)
+
+    def jacobian(self, x):
+        return self.v
+
+
+class _IPopt(object):
+
+    def __init__(self, verbose=True):
+        self.verbose = verbose
+        self.iter_count = 0
+        self.constraints = []
+
+    def add_constraint(self, constraint):
+        self.constraints.append(constraint)
+
+    def objective(self, x):
+        # objective fct to minimize
+        return -x.T.dot(C).dot(x)
+
+    def gradient(self, x):
+        # grad of the objective fct
+        return -2 * x.T.dot(C)
+
+    def constraints(self, x):
+        return np.array([c.constraint(x) for c in self.constraints])
+
+    def jacobian(self, x):
+        return np.array([c.jacobian(x) for c in self.constraints])
+
+    # def hessian(self, x):
+    #    pass
+
+    def intermediate(self, alg_mod, iter_count, obj_value, inf_pr, inf_du, mu, d_norm,
+                     regularization_size, alpha_du, alpha_pr, ls_trials):
+        if self.verbose:
+            print("Objective value at iteration #%d: %g" % (iter_count, obj_value))
+        self.iter_count = iter_count
+
+
+class IPopt(Optimizer):
+    r"""Interior-Point optimizer
+
+    This is a wrapper around the `ipopt` library. It can be used to solve general nonlinear programming problems of
+    the form:
+
+    .. math::
+
+        \min_{x \in R^n} f(x)
+
+    subject to
+
+    .. math::
+
+        g_L \leq g(x) \leq g_U
+
+        x_L \leq  x  \leq x_U
+
+    where :math:`x` are the optimization variables (possibly with upper an lower bounds, :math:`x_U` and :math:`x_L`
+    respectively), :math:`f(x)` is the objective function and :math:`g(x)` are the general nonlinear constraints.
+    The constraints, :math:`g(x)`, have lower and upper bounds. Note that equality constraints can be specified
+    by setting :math:`g^i_L = g^i_U`.
+
+    More info:
+    - Check the documentation of the `ipopt.problem` method
+
+    References:
+        [1] "On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear
+            programming", Wachter and Biegler, 2004
+        [2] Ipopt: https://projects.coin-or.org/Ipopt
+        [3] Ipopt in Python: https://pythonhosted.org/ipopt/
+        [4] Repos: https://github.com/coin-or/Ipopt   and   https://pypi.org/project/ipopt/
+    """
+
+    def __init__(self, model, losses, hyperparameters):
+        super(IPopt, self).__init__(model, losses, hyperparameters)
+
+    def optimize(self):
+        # define initial value
+        x0 = np.array([0.1] * N)  # important that the initial value != 0 for the computation of the grad!
+
+        # define (lower and upper) bound constraints
+        lb = [-1] * N
+        ub = [1] * N
+
+        # define constraints; if upper and lower constraints (resp. cu and cl) are equal then equality constraint
+        cl = [1] + [0] * (N - 1)
+        cu = [1] + [0] * (N - 1)
+
+        # create ipopt (which contains the objective function, its gradients, and constraints)
+        opt = _IPopt(verbose=False)
+        opt.add_constraint(NormConstraint())
+        opt.add_constraint(OrthogonalConstraint(x))
+
+        # define the nonlinear optimization problem
+        nlp = ipopt.problem(n=N, m=len(cl[:i]), problem_obj=opt, lb=lb, ub=ub, cl=cl[:i], cu=cu[:i])
+
+        # solve problem
+        x, info = nlp.solve(x0)
+
+        return x
+
+
+##################################################################
+#                      PYTORCH OPTIMIZERS                        #
+##################################################################
+
+
+class PyTorchOpt(Optimizer):
+    r"""PyTorch Optimizers
+
+    This is a wrapper around the optimizers from pytorch.
+    """
+
+    def __init__(self, model, losses, hyperparameters):
+        super(PyTorchOpt, self).__init__(model, losses, hyperparameters)
+
+    def add_constraint(self):
+        # it will add a constraint as the augmented lagrangian
+        pass
+
+
+class Adam(object):
+    r"""Adam Optimizer
+
+    References:
+        [1] "Adam: A Method for Stochastic Optimization", Kingma et al., 2014
+    """
+
+    def __init__(self, learning_rate=1e-3, betas=(0.9, 0.999), eps=1e-08, weight_decay=0, amsgrad=False,
+                 max_grad_norm=None):  # 0.5
+        self.optimizer = None
+        self.learning_rate = learning_rate
+        self.betas = betas
+        self.eps = eps
+        self.weight_decay = weight_decay
+        self.amsgrad = amsgrad
+        self.max_grad_norm = max_grad_norm
+
+    def reset(self):
+        self.optimizer = None
+
+    def optimize(self, params, loss):
+        # create optimizer if necessary
+        if self.optimizer is None:
+            self.optimizer = optim.Adam(params, lr=self.learning_rate, betas=self.betas, eps=self.eps,
+                                        weight_decay=self.weight_decay, amsgrad=self.amsgrad)
+
+        # optimize
+        self.optimizer.zero_grad()
+        loss.backward()
+        if self.max_grad_norm is not None:
+            nn.utils.clip_grad_norm_(params, self.max_grad_norm)
+        self.optimizer.step()
+
+
+class Adadelta(object):
+    r"""Adadelta Optimizer
+
+    References:
+        [1] "ADADELTA: An Adaptive Learning Rate Method", Zeiler, 2012
+    """
+
+    def __init__(self, learning_rate=1., rho=0.9, eps=1e-6, weight_decay=0, max_grad_norm=None): #0.5
+        self.optimizer = None
+        self.learning_rate = learning_rate
+        self.rho = rho
+        self.eps = eps
+        self.weight_decay = weight_decay
+        self.max_grad_norm = max_grad_norm
+
+    def optimize(self, params, loss):
+        if self.optimizer is None:
+            self.optimizer = optim.Adadelta(params, lr=self.learning_rate, rho=self.rho, eps=self.eps,
+                                            weight_decay=self.weight_decay)
+
+        # optimize
+        self.optimizer.zero_grad()
+        loss.backward()
+        if self.max_grad_norm is not None:
+            nn.utils.clip_grad_norm_(params, self.max_grad_norm)
+        self.optimizer.step()
+
+
+class Adagrad(object):
+    r"""Adagrad Optimizer
+
+    References:
+        [1] "Adaptive Subgradient Methods for Online Learning and Stochastic Optimization", Duchi et al., 2011
+    """
+
+    def __init__(self, learning_rate=0.01, learning_rate_decay=0, weight_decay=0, initial_accumumaltor_value=0,
+                 max_grad_norm=None):  # 0.5
+        self.optimizer = None
+        self.learning_rate = learning_rate
+        self.learning_rate_decay = learning_rate_decay
+        self.weight_decay = weight_decay
+        self.initial_accumulator_value = initial_accumumaltor_value
+        self.max_grad_norm = max_grad_norm
+
+    def optimize(self, params, loss):
+        if self.optimizer is None:
+            self.optimizer = optim.Adagrad(params, lr=self.learning_rate, lr_decay=self.learning_rate_decay,
+                                           weight_decay=self.weight_decay,
+                                           initial_accumulator_value=self.initial_accumulator_value)
+
+        # optimize
+        self.optimizer.zero_grad()
+        loss.backward()
+        if self.max_grad_norm is not None:
+            nn.utils.clip_grad_norm_(params, self.max_grad_norm)
+        self.optimizer.step()
+
+
+class RMSprop(object):
+    r"""RMSprop
+
+    References:
+        [1] "RMSprop:  Divide the gradient by a running average of its recent magnitude" (lecture 6.5), Tieleman and
+            Hinton, 2012
+        [2] "Generating Sequences With Recurrent Neural Networks", Graves, 2014
+    """
+
+    def __init__(self, learning_rate=1e-2, alpha=0.99, eps=1e-8, weight_decay=0, momentum=0, centered=False,
+                 max_grad_norm=None):  # 0.5
+        self.optimizer = None
+        self.learning_rate = learning_rate
+        self.alpha = alpha
+        self.eps = eps
+        self.weight_decay = weight_decay
+        self.momentum = momentum
+        self.centered = centered
+        self.max_grad_norm = max_grad_norm
+
+    def optimize(self, params, loss):
+        if self.optimizer is None:
+            self.optimizer = optim.RMSprop(params, lr=self.learning_rate, alpha=self.alpha, eps=self.eps,
+                                           weight_decay=self.weight_decay, momentum=self.momentum,
+                                           centered=self.centered)
+
+        # optimize
+        self.optimizer.zero_grad()
+        loss.backward()
+        if self.max_grad_norm is not None:
+            nn.utils.clip_grad_norm_(params, self.max_grad_norm)
+        self.optimizer.step()
+
+
+class SGD(object):
+    r"""Stochastic Gradient Descent
+
+    References:
+        [1] "A Stochastic Approximation Method", Robbins and Monro, 1951
+        [2] "On the importance of initialization and momentum in deep learning", Sutskever et al., 2013
+    """
+
+    def __init__(self, learning_rate=1e-3, momentum=0, dampening=0, weight_decay=0, nesterov=False,
+                 max_grad_norm=None): #0.5
+        self.optimizer = None
+        self.learning_rate = learning_rate
+        self.momentum = momentum
+        self.dampening = dampening
+        self.weight_decay = weight_decay
+        self.nesterov = nesterov
+        self.max_grad_norm = max_grad_norm
+
+    def optimize(self, params, loss):
+        # create optimizer if necessary
+        if self.optimizer is None:
+            self.optimizer = optim.SGD(params, lr=self.learning_rate, momentum=self.momentum, dampening=self.dampening,
+                                       weight_decay=self.weight_decay, nesterov=self.nesterov)
+
+        # optimize
+        self.optimizer.zero_grad()
+        loss.backward()
+        if self.max_grad_norm is not None:
+            nn.utils.clip_grad_norm_(params, self.max_grad_norm)
+        self.optimizer.step()

pyrobolearn/rewards/reward.py

@@ -27,7 +27,7 @@ from pyrobolearn.actions import *
 __author__ = "Brian Delhaisse"
 __copyright__ = "Copyright 2018, PyRoboLearn"
 __credits__ = ["Brian Delhaisse"]
-__license__ = "(c) Brian Delhaisse"
+__license__ = "MIT"
 __version__ = "1.0.0"
 __maintainer__ = "Brian Delhaisse"
 __email__ = "briandelhaisse@gmail.com"

pyrobolearn/states/basic_states.py

@@ -12,7 +12,7 @@ from pyrobolearn.actions import Action
 __author__ = "Brian Delhaisse"
 __copyright__ = "Copyright 2018, PyRoboLearn"
 __credits__ = ["Brian Delhaisse"]
-__license__ = "(c) Brian Delhaisse"
+__license__ = "MIT"
 __version__ = "1.0.0"
 __maintainer__ = "Brian Delhaisse"
 __email__ = "briandelhaisse@gmail.com"

pyrobolearn/states/gym_states.py

@@ -13,7 +13,7 @@ from state import State
 __author__ = "Brian Delhaisse"
 __copyright__ = "Copyright 2018, PyRoboLearn"
 __credits__ = ["Brian Delhaisse"]
-__license__ = "(c) Brian Delhaisse"
+__license__ = "MIT"
 __version__ = "1.0.0"
 __maintainer__ = "Brian Delhaisse"
 __email__ = "briandelhaisse@gmail.com"

pyrobolearn/states/object_states.py

@@ -13,7 +13,7 @@ from pyrobolearn.robots import Object
 __author__ = "Brian Delhaisse"
 __copyright__ = "Copyright 2018, PyRoboLearn"
 __credits__ = ["Brian Delhaisse"]
-__license__ = "(c) Brian Delhaisse"
+__license__ = "MIT"
 __version__ = "1.0.0"
 __maintainer__ = "Brian Delhaisse"
 __email__ = "briandelhaisse@gmail.com"

pyrobolearn/states/robot_states/joint_states.py

@@ -11,7 +11,7 @@ from robot_states import RobotState
 __author__ = "Brian Delhaisse"
 __copyright__ = "Copyright 2018, PyRoboLearn"
 __credits__ = ["Brian Delhaisse"]
-__license__ = "(c) Brian Delhaisse"
+__license__ = "MIT"
 __version__ = "1.0.0"
 __maintainer__ = "Brian Delhaisse"
 __email__ = "briandelhaisse@gmail.com"

pyrobolearn/states/robot_states/link_states.py

@@ -11,7 +11,7 @@ from robot_states import RobotState
 __author__ = "Brian Delhaisse"
 __copyright__ = "Copyright 2018, PyRoboLearn"
 __credits__ = ["Brian Delhaisse"]
-__license__ = "(c) Brian Delhaisse"
+__license__ = "MIT"
 __version__ = "1.0.0"
 __maintainer__ = "Brian Delhaisse"
 __email__ = "briandelhaisse@gmail.com"

pyrobolearn/states/robot_states/robot_states.py

@@ -17,7 +17,7 @@ from pyrobolearn.robots import Robot
 __author__ = "Brian Delhaisse"
 __copyright__ = "Copyright 2018, PyRoboLearn"
 __credits__ = ["Brian Delhaisse"]
-__license__ = "(c) Brian Delhaisse"
+__license__ = "MIT"
 __version__ = "1.0.0"
 __maintainer__ = "Brian Delhaisse"
 __email__ = "briandelhaisse@gmail.com"

pyrobolearn/states/robot_states/sensor_states.py

@@ -11,7 +11,7 @@ from robot_states import RobotState
 __author__ = "Brian Delhaisse"
 __copyright__ = "Copyright 2018, PyRoboLearn"
 __credits__ = ["Brian Delhaisse"]
-__license__ = "(c) Brian Delhaisse"
+__license__ = "MIT"
 __version__ = "1.0.0"
 __maintainer__ = "Brian Delhaisse"
 __email__ = "briandelhaisse@gmail.com"

pyrobolearn/states/state.py

@@ -16,7 +16,7 @@ from pyrobolearn.utils.data_structures.orderedset import OrderedSet
 __author__ = "Brian Delhaisse"
 __copyright__ = "Copyright 2018, PyRoboLearn"
 __credits__ = ["Brian Delhaisse"]
-__license__ = "(c) Brian Delhaisse"
+__license__ = "MIT"
 __version__ = "1.0.0"
 __maintainer__ = "Brian Delhaisse"
 __email__ = "briandelhaisse@gmail.com"

pyrobolearn/states/time_states.py

@@ -13,7 +13,7 @@ from state import State
 __author__ = "Brian Delhaisse"
 __copyright__ = "Copyright 2018, PyRoboLearn"
 __credits__ = ["Brian Delhaisse"]
-__license__ = "(c) Brian Delhaisse"
+__license__ = "MIT"
 __version__ = "1.0.0"
 __maintainer__ = "Brian Delhaisse"
 __email__ = "briandelhaisse@gmail.com"