"""

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.

The development of this software was sponsored by NAG Ltd. (http://www.nag.co.uk)
and the EPSRC Centre For Doctoral Training in Industrially Focused Mathematical
Modelling (EP/L015803/1) at the University of Oxford. Please contact NAG for
alternative licensing.

"""

# Ensure compatibility with Python 2
from __future__ import absolute_import, division, print_function, unicode_literals

from math import sqrt
import numpy as np
import unittest

from dfogn.linear_altmov import max_step_in_box_and_ball 
from dfogn.util import calculate_model_value


def cauchy_pt(g, hess, delta):
    # General expression for the Cauchy point
    crv = np.dot(g, hess.vec_mul(g))
    gnorm = np.linalg.norm(g)
    if crv <= 0.0:
        alpha = delta / gnorm
    else:
        alpha = min(delta / gnorm, gnorm**2 / crv)
    s = -alpha * g
    red = calculate_model_value(g, hess, s)
    crvmin = np.dot(s, hess.vec_mul(s)) / np.dot(s, s)
    if crvmin < 0.0:
        crvmin = -1.0
    return s, red, crvmin


def cauchy_pt_box(g, hess, delta, lower, upper):
    # General expression for the Cauchy point, lower <= s <= upper
    crv = np.dot(g, hess.vec_mul(g))
    gnorm = np.linalg.norm(g)
    if crv <= 0.0:
        alpha = delta / gnorm
    else:
        alpha = min(delta / gnorm, gnorm**2 / crv)
    # print("alpha = %g" % alpha)
    # Then cap with bounds:
    for i in range(len(g)):
        if g[i] > 0:  # s[i] negative, will hit lower
            alpha = min(alpha, -lower[i] / g[i])
        elif g[i] < 0:  # s[i] positive, will hit upper
            alpha = min(alpha, -upper[i] / g[i])
        # print("alpha = %g after i=%g" % (alpha, i))
    s = -alpha * g
    red = calculate_model_value(g, hess, s)
    crvmin = np.dot(s, hess.vec_mul(s)) / np.dot(s, s)
    if crvmin < 0.0:
        crvmin = -1.0
    return s, red, crvmin


class TestGeom(unittest.TestCase):
    def runTest(self):
        xbase = np.array([0.0, 0.0])
        g = np.array([1.0, -1.0])
        a = np.array([-2.0, -2.0])
        b = np.array([1.0, 2.0])
        delta = 2.0
        c = -1.0
        x = max_step_in_box_and_ball(xbase, c, g, a, b, delta)
        xtrue = np.array([-sqrt(2.0), sqrt(2.0)])
        self.assertTrue(np.max(np.abs(x - xtrue)) < 1e-10, 'Wrong step')


class TestGeom2(unittest.TestCase):
    def runTest(self):
        xbase = np.array([0.0, 0.0])
        g = np.array([1.0, -1.0])
        a = np.array([-2.0, -2.0])
        b = np.array([1.0, 2.0])
        delta = 5.0
        c = -1.0
        x = max_step_in_box_and_ball(xbase, c, g, a, b, delta)
        xtrue = np.array([-2.0, 2.0])
        self.assertTrue(np.max(np.abs(x - xtrue)) < 1e-10, 'Wrong step')


class TestGeom3(unittest.TestCase):
    def runTest(self):
        xbase = np.array([0.0, 0.0]) + 1
        g = np.array([1.0, -1.0])
        a = np.array([-2.0, -2.0]) + 1
        b = np.array([1.0, 2.0]) + 1
        delta = 5.0
        c = 3.0  # may want to max instead
        x = max_step_in_box_and_ball(xbase, c, g, a, b, delta)
        xtrue = np.array([1.0, -2.0]) + 1
        self.assertTrue(np.max(np.abs(x - xtrue)) < 1e-10, 'Wrong step')


class TestGeomOldBug(unittest.TestCase):
    def runTest(self):
        xbase = np.array([0.0, 0.0])
        g = np.array([-1.0, -1.0])
        a = np.array([-2.0, -2.0])
        b = np.array([0.1, 0.9])
        delta = sqrt(2.0)
        c = -1.0  # may want to max instead
        x = max_step_in_box_and_ball(xbase, c, g, a, b, delta)
        xtrue = b
        print(x)
        self.assertTrue(np.max(np.abs(x - xtrue)) < 1e-10, 'Wrong step')
        # self.assertFalse(True, "bad")


class TestGeomOldBug2(unittest.TestCase):
    def runTest(self):
        xbase = np.array([0.0, 0.0, 0.0])
        g = np.array([-1.0, -1.0, -1.0])
        a = np.array([-2.0, -2.0, -2.0])
        b = np.array([0.9, 0.1, 5.0])
        delta = sqrt(3.0)
        c = -1.0  # may want to max instead
        x = max_step_in_box_and_ball(xbase, c, g, a, b, delta)
        xtrue = np.array([0.9, 0.1, sqrt(3.0 - 0.81 - 0.01)])
        print(x)
        self.assertTrue(np.max(np.abs(x - xtrue)) < 1e-10, 'Wrong step')
        # self.assertFalse(True, "bad")


class TestGeom2WithZeros(unittest.TestCase):
    def runTest(self):
        xbase = np.array([0.0, 0.0])
        g = np.array([0.0, -1.0])
        a = np.array([-2.0, -2.0])
        b = np.array([1.0, 2.0])
        delta = 5.0
        c = 0.0
        x = max_step_in_box_and_ball(xbase, c, g, a, b, delta)
        xtrue = np.array([0.0, 2.0])
        self.assertTrue(np.max(np.abs(x - xtrue)) < 1e-10, 'Wrong step')


class TestGeom2WithAlmostZeros(unittest.TestCase):
    def runTest(self):
        xbase = np.array([0.0, 0.0])
        g = np.array([1e-15, -1.0])
        a = np.array([-2.0, -2.0])
        b = np.array([1.0, 2.0])
        delta = 5.0
        c = 0.0
        x = max_step_in_box_and_ball(xbase, c, g, a, b, delta)
        xtrue = np.array([0.0, 2.0])
        self.assertTrue(np.max(np.abs(x - xtrue)) < 1e-10, 'Wrong step')


class TestGeom2WithAlmostZeros2(unittest.TestCase):
    def runTest(self):
        xbase = np.array([0.0, 0.0])
        g = np.array([1e-15, 0.0])
        a = np.array([-2.0, -2.0])
        b = np.array([1.0, 2.0])
        delta = 5.0
        c = 0.0
        x = max_step_in_box_and_ball(xbase, c, g, a, b, delta)
        xtrue = np.array([0.0, 0.0])
        # print(x)
        # print(xtrue)
        self.assertTrue(np.max(np.abs(x - xtrue)) < 1e-10, 'Wrong step')
