# 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
        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)
            lmax = abs(c + model_value(g, H, xmax-xbase))  # evaluate Lagrange poly
            if overall_max is None or lmax > overall_max:
                overall_max = lmax
        return overall_max


def build_interpolation_matrix(Y, approx_delta=1.0):
    # Y has columns Y[:,t] = yt - xk
    n, p = Y.shape  # p = npt-1
    assert n + 1 <= p + 1 <= (n + 1) * (n + 2) // 2, "npt must be in range [n+1, (n+1)(n+2)/2]"

    # What scaling was applied to each part of the matrix?
    # A(scaled) = diag(left_scaling) * A(unscaled) * diag(right_scaling)

    if p == n:  # linear models
        A = Y.T / approx_delta
        left_scaling = np.ones((n,))  # no left scaling
        right_scaling = np.ones((n,)) / approx_delta
    elif p + 1 == (n+1)*(n+2)//2:  # fully quadratic models
        A = np.zeros((p, p))
        A[:,:n] = Y.T / approx_delta
        for i in range(p):

# 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

        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)
            lmax = abs(c + model_value(g, H, xmax-xbase))  # evaluate Lagrange poly
            if overall_max is None or lmax > overall_max:
                overall_max = lmax
        return overall_max


def build_interpolation_matrix(Y, approx_delta=1.0):
    # Y has columns Y[:,t] = yt - xk
    n, p = Y.shape  # p = npt-1
    assert n + 1 <= p + 1 <= (n + 1) * (n + 2) // 2, "npt must be in range [n+1, (n+1)(n+2)/2]"

    # What scaling was applied to each part of the matrix?
    # A(scaled) = diag(left_scaling) * A(unscaled) * diag(right_scaling)

    if p == n:  # linear models
        A = Y.T / approx_delta
        left_scaling = np.ones((n,))  # no left scaling
        right_scaling = np.ones((n,)) / approx_delta
    elif p + 1 == (n+1)*(n+2)//2:  # fully quadratic models
        A = np.zeros((p, p))
        A[:,:n] = Y.T / approx_delta
        for i in range(p):

            A[i, n:] = to_upper_triangular_vector(np.outer(Y[:,i], Y[:,i]) - 0.5*np.diag(np.square(Y[:,i]))) / (approx_delta**2)