"""Tests for ot.smooth model """

# Author: Remi Flamary <remi.flamary@unice.fr>
#
# License: MIT License

import numpy as np
import ot
import pytest
from scipy.optimize import check_grad


def test_smooth_ot_dual():

    # get data
    n = 100
    rng = np.random.RandomState(0)

    x = rng.randn(n, 2)
    u = ot.utils.unif(n)

    M = ot.dist(x, x)

    with pytest.raises(NotImplementedError):
        Gl2, log = ot.smooth.smooth_ot_dual(u, u, M, 1, reg_type='none')

    # squared l2 regularisation
    Gl2, log = ot.smooth.smooth_ot_dual(u, u, M, 1, reg_type='l2', log=True, stopThr=1e-10)

    # check constraints
    np.testing.assert_allclose(
        u, Gl2.sum(1), atol=1e-05)  # cf convergence sinkhorn
    np.testing.assert_allclose(
        u, Gl2.sum(0), atol=1e-05)  # cf convergence sinkhorn

    # kl regularisation
    G = ot.smooth.smooth_ot_dual(u, u, M, 1, reg_type='kl', stopThr=1e-10)

    # check constraints
    np.testing.assert_allclose(
        u, G.sum(1), atol=1e-05)  # cf convergence sinkhorn
    np.testing.assert_allclose(
        u, G.sum(0), atol=1e-05)  # cf convergence sinkhorn

    G2 = ot.sinkhorn(u, u, M, 1, stopThr=1e-10)
    np.testing.assert_allclose(G, G2, atol=1e-05)

    # sparsity-constrained regularisation
    max_nz = 2
    Gsc, log = ot.smooth.smooth_ot_dual(
        u, u, M, 1,
        max_nz=max_nz,
        log=True,
        reg_type='sparsity_constrained',
        stopThr=1e-10)

    # check marginal constraints
    np.testing.assert_allclose(u, Gsc.sum(1), atol=1e-03)
    np.testing.assert_allclose(u, Gsc.sum(0), atol=1e-03)

    # check sparsity constraints
    np.testing.assert_array_less(
        np.sum(Gsc > 0, axis=0),
        np.ones(n) * max_nz + 1)


def test_smooth_ot_semi_dual():

    # get data
    n = 100
    rng = np.random.RandomState(0)

    x = rng.randn(n, 2)
    u = ot.utils.unif(n)

    M = ot.dist(x, x)

    with pytest.raises(NotImplementedError):
        Gl2, log = ot.smooth.smooth_ot_semi_dual(u, u, M, 1, reg_type='none')

    # squared l2 regularisation
    Gl2, log = ot.smooth.smooth_ot_semi_dual(u, u, M, 1, reg_type='l2', log=True, stopThr=1e-10)

    # check constraints
    np.testing.assert_allclose(
        u, Gl2.sum(1), atol=1e-05)  # cf convergence sinkhorn
    np.testing.assert_allclose(
        u, Gl2.sum(0), atol=1e-05)  # cf convergence sinkhorn

    # kl regularisation
    G = ot.smooth.smooth_ot_semi_dual(u, u, M, 1, reg_type='kl', stopThr=1e-10)

    # check constraints
    np.testing.assert_allclose(
        u, G.sum(1), atol=1e-05)  # cf convergence sinkhorn
    np.testing.assert_allclose(
        u, G.sum(0), atol=1e-05)  # cf convergence sinkhorn

    G2 = ot.sinkhorn(u, u, M, 1, stopThr=1e-10)
    np.testing.assert_allclose(G, G2, atol=1e-05)

    # sparsity-constrained regularisation
    max_nz = 2
    Gsc = ot.smooth.smooth_ot_semi_dual(
        u, u, M, 1, reg_type='sparsity_constrained',
        max_nz=max_nz, stopThr=1e-10)

    # check marginal constraints
    np.testing.assert_allclose(u, Gsc.sum(1), atol=1e-03)
    np.testing.assert_allclose(u, Gsc.sum(0), atol=1e-03)

    # check sparsity constraints
    np.testing.assert_array_less(np.sum(Gsc > 0, axis=0),
                                 np.ones(n) * max_nz + 1)


def test_sparsity_constrained_gradient():
    max_nz = 5
    regularizer = ot.smooth.SparsityConstrained(max_nz=max_nz)
    rng = np.random.RandomState(0)
    X = rng.randn(10,)
    b = 0.5

    def delta_omega_func(X):
        return regularizer.delta_Omega(X)[0]

    def delta_omega_grad(X):
        return regularizer.delta_Omega(X)[1]

    dual_grad_err = check_grad(delta_omega_func, delta_omega_grad, X)
    np.testing.assert_allclose(dual_grad_err, 0.0, atol=1e-07)

    def max_omega_func(X, b):
        return regularizer.max_Omega(X, b)[0]

    def max_omega_grad(X, b):
        return regularizer.max_Omega(X, b)[1]

    semi_dual_grad_err = check_grad(max_omega_func, max_omega_grad, X, b)
    np.testing.assert_allclose(semi_dual_grad_err, 0.0, atol=1e-07)