# Path: pybobyqa/hessian.py
# def to_upper_triangular_vector(A):
#     n = A.shape[0]
#     hq = np.zeros((n*(n+1)//2,))
#     ih = -1
#     for j in range(n):  # j = 0, ..., n-1
#         for i in range(j + 1):  # i = 0, ..., j
#             ih += 1
#             hq[ih] = A[i,j]
#     return hq
# 
# Path: pybobyqa/trust_region.py
# def trsbox_geometry(xbase, c, g, H, lower, upper, Delta, use_fortran=USE_FORTRAN):
#     # Given a Lagrange polynomial defined by: L(x) = c + g' * (x - xbase) + 0.5*(x-xbase)*H*(x-xbase)
#     # Maximise |L(x)| in a box + trust region - that is, solve:
#     #   max_x  abs(c + g' * (x - xbase) + 0.5*(x-xbase)*H*(x-xbase))
#     #    s.t.  lower <= x <= upper
#     #          ||x-xbase|| <= Delta
#     # Setting s = x-xbase (or x = xbase + s), this is equivalent to:
#     #   max_s  abs(c + g' * s + 0.5*s*H*s)
#     #   s.t.   lower <= xbase + s <= upper
#     #          ||s|| <= Delta
#     smin, gmin, crvmin = trsbox(xbase, g, H, lower, upper, Delta, use_fortran=use_fortran)  # minimise L(x)
#     smax, gmax, crvmax = trsbox(xbase, -g, -H, lower, upper, Delta, use_fortran=use_fortran)  # maximise L(x)
#     if abs(c + model_value(g, H, smin)) >= abs(c + model_value(g, H, smax)):  # take largest abs value
#         return xbase + smin
#     else:
#         return xbase + smax
# 
# Path: pybobyqa/util.py
# def sumsq(x):
#     # There are several ways to calculate sum of squares of a vector:
#     #   np.dot(x,x)
#     #   np.sum(x**2)
#     #   np.sum(np.square(x))
#     #   etc.
#     # Using the timeit routine, it seems like dot(x,x) is ~3-4x faster than other methods
#     return np.dot(x, x)
# 
# def model_value(g, H, s):
#     # Calculate model value (s^T * g + 0.5* s^T * H * s) = s^T * (gopt + 0.5 * H*s)
#     assert g.shape == s.shape, "g and s have incompatible sizes"
#     Hs = H.dot(s)
#     return np.dot(s, g + 0.5*Hs)


import logging
import numpy as np
import scipy.linalg as LA
from math import sqrt
from .hessian import to_upper_triangular_vector
from .trust_region import trsbox_geometry
from .util import sumsq, model_value
            H = np.zeros((self.n(), self.n()))
        elif self.npt() == (self.n() + 1) * (self.n() + 2) // 2:
            g = coeffs[:self.n()]
            H = build_symmetric_matrix_from_vector(self.n(), coeffs[self.n():])  # rest of coeffs are upper triangular part of Hess
        else:
            g = coeffs[self.npt() - 1:]  # last n values
            fval_row_idx = np.delete(np.arange(self.npt()), self.kopt)  # indices of all rows except kopt
            H = np.zeros((self.n(), self.n()))
            for i in range(self.npt() - 1):
                dx = self.xpt(fval_row_idx[i]) - self.xopt()
                H += coeffs[i] * np.outer(dx, dx)

        # (c, g, hess) currently based around xopt
        return c, g, H

    def poisedness_constant(self, delta, xbase=None, xbase_in_abs_coords=True):
        # Calculate the poisedness constant of the current interpolation set in B(xbase, delta)
        # if xbase is None, use self.xopt()
        overall_max = None
        if xbase is None:
            xbase = self.xopt()
        elif xbase_in_abs_coords:
            xbase = xbase - self.xbase  # shift to correct position
        for k in range(self.npt()):
            c, g, H = self.lagrange_polynomial(k, factorise_first=True)  # based at self.xopt()
            # Switch base of poly from xopt to xbase, as required by trsbox_geometry
            base_chg = self.xopt() - xbase
            Hx = H.dot(base_chg)
            c += np.dot(-g + 0.5 * Hx, base_chg)
            g += -Hx

# Path: pybobyqa/hessian.py
# def to_upper_triangular_vector(A):
#     n = A.shape[0]
#     hq = np.zeros((n*(n+1)//2,))
#     ih = -1
#     for j in range(n):  # j = 0, ..., n-1
#         for i in range(j + 1):  # i = 0, ..., j
#             ih += 1
#             hq[ih] = A[i,j]
#     return hq
# 
# Path: pybobyqa/trust_region.py
# def trsbox_geometry(xbase, c, g, H, lower, upper, Delta, use_fortran=USE_FORTRAN):
#     # Given a Lagrange polynomial defined by: L(x) = c + g' * (x - xbase) + 0.5*(x-xbase)*H*(x-xbase)
#     # Maximise |L(x)| in a box + trust region - that is, solve:
#     #   max_x  abs(c + g' * (x - xbase) + 0.5*(x-xbase)*H*(x-xbase))
#     #    s.t.  lower <= x <= upper
#     #          ||x-xbase|| <= Delta
#     # Setting s = x-xbase (or x = xbase + s), this is equivalent to:
#     #   max_s  abs(c + g' * s + 0.5*s*H*s)
#     #   s.t.   lower <= xbase + s <= upper
#     #          ||s|| <= Delta
#     smin, gmin, crvmin = trsbox(xbase, g, H, lower, upper, Delta, use_fortran=use_fortran)  # minimise L(x)
#     smax, gmax, crvmax = trsbox(xbase, -g, -H, lower, upper, Delta, use_fortran=use_fortran)  # maximise L(x)
#     if abs(c + model_value(g, H, smin)) >= abs(c + model_value(g, H, smax)):  # take largest abs value
#         return xbase + smin
#     else:
#         return xbase + smax
# 
# Path: pybobyqa/util.py
# def sumsq(x):
#     # There are several ways to calculate sum of squares of a vector:
#     #   np.dot(x,x)
#     #   np.sum(x**2)
#     #   np.sum(np.square(x))
#     #   etc.
#     # Using the timeit routine, it seems like dot(x,x) is ~3-4x faster than other methods
#     return np.dot(x, x)
# 
# def model_value(g, H, s):
#     # Calculate model value (s^T * g + 0.5* s^T * H * s) = s^T * (gopt + 0.5 * H*s)
#     assert g.shape == s.shape, "g and s have incompatible sizes"
#     Hs = H.dot(s)
#     return np.dot(s, g + 0.5*Hs)


# Path: pybobyqa/model.py
import logging
import numpy as np
import scipy.linalg as LA
from math import sqrt
from .hessian import to_upper_triangular_vector
from .trust_region import trsbox_geometry
from .util import sumsq, model_value

            H = np.zeros((self.n(), self.n()))
        elif self.npt() == (self.n() + 1) * (self.n() + 2) // 2:
            g = coeffs[:self.n()]
            H = build_symmetric_matrix_from_vector(self.n(), coeffs[self.n():])  # rest of coeffs are upper triangular part of Hess
        else:
            g = coeffs[self.npt() - 1:]  # last n values
            fval_row_idx = np.delete(np.arange(self.npt()), self.kopt)  # indices of all rows except kopt
            H = np.zeros((self.n(), self.n()))
            for i in range(self.npt() - 1):
                dx = self.xpt(fval_row_idx[i]) - self.xopt()
                H += coeffs[i] * np.outer(dx, dx)

        # (c, g, hess) currently based around xopt
        return c, g, H

    def poisedness_constant(self, delta, xbase=None, xbase_in_abs_coords=True):
        # Calculate the poisedness constant of the current interpolation set in B(xbase, delta)
        # if xbase is None, use self.xopt()
        overall_max = None
        if xbase is None:
            xbase = self.xopt()
        elif xbase_in_abs_coords:
            xbase = xbase - self.xbase  # shift to correct position
        for k in range(self.npt()):
            c, g, H = self.lagrange_polynomial(k, factorise_first=True)  # based at self.xopt()
            # Switch base of poly from xopt to xbase, as required by trsbox_geometry
            base_chg = self.xopt() - xbase
            Hx = H.dot(base_chg)
            c += np.dot(-g + 0.5 * Hx, base_chg)
            g += -Hx

            xmax = trsbox_geometry(xbase, c, g, H, self.sl, self.su, delta)