## module error
''' 

'''
import sys
def err(string):
    print(string)
    input('Press return to exit')
    sys.exit()

## modulo swap
''' 

'''
def swapRows(v,i,j):
    if len(v.shape) == 1:
        v[i],v[j] = v[j],v[i]
    else:
        v[[i,j],:] = v[[j,i],:]
        
def swapCols(v,i,j):
    v[:,[i,j]] = v[:,[j,i]]

### module gaussPivot
''' x = gaussPivot(a,b,tol=1.0e-12).
    Resuelve [a]{x} = {b} mediante la eliminación de Gauss
    con renglón pivote escalado
'''
import numpy as np

def gaussPivot(a,b,tol=1.0e-12):
    n = len(b)
    
#   Definiendo los factores escalables
    s = np.zeros(n)
    for i in range(n):
        s[i] = max(np.abs(a[i,:]))
        
    for k in range(0, n-1):
        
    #   Intercambio de renglones, si se necesita
        p = np.argmax(np.abs(a[k:n,k])/s[k:n]) + k
        if abs(a[p,k]) < tol: error.err('Matrix is singular')
        if p != k:
            swapRows(b,k,p)
            swapRows(s,k,p)
            swapRows(a,k,p)
            
    #   Eliminacion
        for i in range(k+1,n):
            if a[i,k] != 0.0:
                lam = a[i,k]/a[k,k]
                a[i,k+1:n] = a[i,k+1:n] - lam*a[k,k+1:n]
                b[i] = b[i] - lam*b[k]
    if abs(a[n-1,n-1]) < tol: error.err('Matrix is singlar')
        
#   Substitution en reversa
    b[n-1] = b[n-1]/a[n-1,n-1]
    for k in range(n-2,-1,-1):
        b[k] = (b[k] - np.dot(a[k,k+1:n],b[k+1:n]))/a[k,k]
    return b

## module newtonRaphson2
''' soln = newtonRaphson2(f,x,tol=1.0e-9).
    Resuelve las excuaciones simultáneas f(x) = 0 by
    El metodo Newton-Raphson utilizando {x} como la variable
    inicial. Notar que {f} y {x} son vectores.
'''

import numpy as np
import math
def newtonRaphson2(f,x,tol=1.0e-9):
    
    def jacobian(f,x):
        h = 1.0e-4
        n = len(x)
        jac = np.zeros((n,n))
        f0 = f(x)
        for i in range(n):
            temp = x[i]
            x[i] = temp + h
            f1 = f(x)
            x[i] = temp
            jac[:,i] = (f1 - f0)/h
        return jac,f0
    
    for i in range(30):
        jac, f0 = jacobian(f,x)
        if math.sqrt(np.dot(f0,f0)/len(x)) < tol: return x
        dx = gaussPivot(jac, -f0)
        x = x + dx
        if math.sqrt(np.dot(dx,dx)) < tol*max(max(abs(x)),1.0):
            return x
    print('Demasiadas Iteraciones')

import numpy as np
import math

def f(x):
    f = np.zeros(len(x))
    f[0] = x[0]**3 + x[0]**2*x[1] - x[0]*x[2] + 6.0
    f[1] = math.exp(x[0]) + math.exp(x[1]) - x[2]
    f[2] = x[1]**2 - 2*x[0]*x[2] - 4.0
    return f

x = np.array([-1.0, -2.0, 1.0])
print('Los valores para [x, y, z], respectivamente:\n ')
print(newtonRaphson2(f,x))