# Path: dev/multipolyfit.py
# def multipolyfit(xs, y, deg, full=False, model_out=False, powers_out=False):
#     """
#     Least squares multivariate polynomial fit
#     Fit a polynomial like ``y = a**2 + 3a - 2ab + 4b**2 - 1``
#     with many covariates a, b, c, ...
#     Parameters
#     ----------
#     xs : array_like, shape (M, k)
#          x-coordinates of the k covariates over the M sample points
#     y :  array_like, shape(M,)
#          y-coordinates of the sample points.
#     deg : int
#          Degree o fthe fitting polynomial
#     model_out : bool (defaults to True)
#          If True return a callable function
#          If False return an array of coefficients
#     powers_out : bool (defaults to False)
#          Returns the meaning of each of the coefficients in the form of an
#          iterator that gives the powers over the inputs and 1
#          For example if xs corresponds to the covariates a,b,c then the array
#          [1, 2, 1, 0] corresponds to 1**1 * a**2 * b**1 * c**0
#     See Also
#     --------
#         numpy.polyfit
#     """
#     y = asarray(y).squeeze()
#     rows = y.shape[0]
#     xs = asarray(xs)
#     num_covariates = xs.shape[1]
#     xs = hstack((ones((xs.shape[0], 1), dtype=xs.dtype) , xs))
# 
#     generators = [basis_vector(num_covariates+1, i)
#                             for i in range(num_covariates+1)]
# 
#     # All combinations of degrees
#     powers = map(sum, itertools.combinations_with_replacement(generators, deg))
#     if powers_out:
#         exponents = [p for p in dc(powers)]
#     
#     # Raise data to specified degree pattern, stack in order
#     A = hstack(asarray([as_tall((xs**p).prod(1)) for p in powers]))
# 
#     beta = linalg.lstsq(A, y)[0]
# 
#     if model_out:
#         return mk_model(beta, powers)
# 
#     if powers_out:
#         return beta, exponents
#     return beta


import numpy as np
import matplotlib.pyplot as plt
from dev import multipolyfit as mpf

# %% hago el fiteo
ruidoY = np.random.randn(y.shape[0]) * 0.1

coefs, exps = mpf.multipolyfit(deltaX, y + ruidoY, 2, powers_out=True)

constRan, jacRan, hesRan = coefs2mats(coefs)


testRan = [constRan + jacRan.dot(x) + x.dot(hesRan).dot(x) / 2 for x in deltaX]

plt.figure()
plt.plot(y, testRan, 'x')
emin, emax = [np.min(y), np.max(y)]
plt.plot([emin, emax],[emin, emax])


plt.figure()
plt.plot(hes[np.tril_indices_from(hes)], hesRan[np.tril_indices_from(hesRan)], '+')

plt.figure()
plt.plot(y, y+ruidoY, 'x')

# %%

hUli = np.array([[2,-10],[-10,2]])
xx = np.random.rand(1000, 2) * 2 - 1

yy = [x.dot(hUli).dot(x) / 2 for x in xx]


# Path: dev/multipolyfit.py
# def multipolyfit(xs, y, deg, full=False, model_out=False, powers_out=False):
#     """
#     Least squares multivariate polynomial fit
#     Fit a polynomial like ``y = a**2 + 3a - 2ab + 4b**2 - 1``
#     with many covariates a, b, c, ...
#     Parameters
#     ----------
#     xs : array_like, shape (M, k)
#          x-coordinates of the k covariates over the M sample points
#     y :  array_like, shape(M,)
#          y-coordinates of the sample points.
#     deg : int
#          Degree o fthe fitting polynomial
#     model_out : bool (defaults to True)
#          If True return a callable function
#          If False return an array of coefficients
#     powers_out : bool (defaults to False)
#          Returns the meaning of each of the coefficients in the form of an
#          iterator that gives the powers over the inputs and 1
#          For example if xs corresponds to the covariates a,b,c then the array
#          [1, 2, 1, 0] corresponds to 1**1 * a**2 * b**1 * c**0
#     See Also
#     --------
#         numpy.polyfit
#     """
#     y = asarray(y).squeeze()
#     rows = y.shape[0]
#     xs = asarray(xs)
#     num_covariates = xs.shape[1]
#     xs = hstack((ones((xs.shape[0], 1), dtype=xs.dtype) , xs))
# 
#     generators = [basis_vector(num_covariates+1, i)
#                             for i in range(num_covariates+1)]
# 
#     # All combinations of degrees
#     powers = map(sum, itertools.combinations_with_replacement(generators, deg))
#     if powers_out:
#         exponents = [p for p in dc(powers)]
#     
#     # Raise data to specified degree pattern, stack in order
#     A = hstack(asarray([as_tall((xs**p).prod(1)) for p in powers]))
# 
#     beta = linalg.lstsq(A, y)[0]
# 
#     if model_out:
#         return mk_model(beta, powers)
# 
#     if powers_out:
#         return beta, exponents
#     return beta


# Path: dev/testHiperparabolaFitting.py
import numpy as np
import matplotlib.pyplot as plt
from dev import multipolyfit as mpf


# %% hago el fiteo
ruidoY = np.random.randn(y.shape[0]) * 0.1

coefs, exps = mpf.multipolyfit(deltaX, y + ruidoY, 2, powers_out=True)

constRan, jacRan, hesRan = coefs2mats(coefs)


testRan = [constRan + jacRan.dot(x) + x.dot(hesRan).dot(x) / 2 for x in deltaX]

plt.figure()
plt.plot(y, testRan, 'x')
emin, emax = [np.min(y), np.max(y)]
plt.plot([emin, emax],[emin, emax])


plt.figure()
plt.plot(hes[np.tril_indices_from(hes)], hesRan[np.tril_indices_from(hesRan)], '+')

plt.figure()
plt.plot(y, y+ruidoY, 'x')

# %%

hUli = np.array([[2,-10],[-10,2]])
xx = np.random.rand(1000, 2) * 2 - 1

yy = [x.dot(hUli).dot(x) / 2 for x in xx]


coefs, exps = mpf.multipolyfit(xx, yy, 2, powers_out=True)