# -*- coding: utf-8 -*-
"""


@author: mlazaro
"""
import numpy as np
import math
#import time

#------------------------------------------------------------------------------
# definición de métodos útiles
#------------------------------------------------------------------------------
import numpy as np
from scipy import special

def find(datos,valor):
    indices = np.asarray([i for (i, val) in enumerate(datos) if val == valor])
    return indices

def findL(datos,valor):
    indices = np.asarray([i for (i, val) in enumerate(datos) if val < valor])
    return indices    

def findG(datos,valor):
    indices = np.asarray([i for (i, val) in enumerate(datos) if val > valor])
    return indices
    
def findLE(datos,valor):
    indices = np.asarray([i for (i, val) in enumerate(datos) if val <= valor])
    return indices    
    
def findGE(datos,valor):
    indices = np.asarray([i for (i, val) in enumerate(datos) if val >= valor])
    return indices    

def setdiff(a,b):
    bs=set(b)
    return [x for x in a if x not in bs]
    
def setintersect(a,b):
    bs=set(b)
    return [x for x in a if x in bs]

    
"""
function o = mlp0b(x,ws,tAct);
%
%  o = mlp0b(x,ws,tAct)
%
%  Función que evalúa un MLP de una capa oculta con términos de bias
%
%     x : patrones de entrada (Ne x Np)
%    ws : parámetros (pesos) de la capa de salida (Ns x Ne+1)
%  tAct : parámetro que define la función de activación de las neuronas
%         - param = 0: lineal
%         - param = 1: tanh
%         - param = 2: logistic
%         - param = 3: rectified linear (ReLU)
%
%     o : valor de salida de la red (Ns x Np)
%
%       Donde Ne, Ns y Np son respectivamente el número de entradas
% 		el número de salidas y el de patrones
%
%--------------------------------------------------------------------------
%         Autor: Marcelino Lázaro
%      Creación: enero 2014
% Actualización: diciembre 2015
%--------------------------------------------------------------------------
"""
# @profile
def mlp0b(x,ws,tAct=1):
    (Ne,Np)=x.shape
    (Ns,Nn)=ws.shape 
   
    o = np.matmul(ws,np.concatenate((x,np.ones((1,Np)))));
    
    if tAct == 1:
        o = np.tanh(o)
    elif tAct == 2:
        o = 1.0/(1+np.exp(-1*o))
    elif tAct[1] == 3:
        o[o<0]=0

    return (o)
    
"""
function [salS,salO] = mlp1b(x,wo,ws,tAct);
%
%  [y,salO] = mlp1b(x,wo,ws,tAct)
%
%  Función que evalúa un MLP de una capa oculta con términos de bias
%
%     x : patrones de entrada (Ne x Np)
%    wo : parámetros (pesos) de la capa oculta (Nn x Ne+1)
%    ws : parámetros (pesos) de la capa de salida (Ns x Nn+1)
%  tAct : parámetro que define el tipo de funciones de activación de cada 
%         - param = 0: lineal
%         - param = 1: tanh
%         - param = 2: logistic
%         - param = 3: rectified linear (ReLU)
%         Es un vector de dimensión 2, con los valores correspondientes a 
%         la capa oculta y a la capa de salida     
%
%     y : valor de salida de la red (Ns x Np)
%  salO : valores de salida de las neuronas (Nn x Np)
%
%
%       Donde Ne, Ns, Np, y Nn son respectivamente el número de entradas
% 		el número de salidas, el de patrones, y el de neuronas
%
%
%--------------------------------------------------------------------------
%         Autor: Marcelino Lázaro
%      Creación: enero 2014
% Actualización: diciembre 2015
%--------------------------------------------------------------------------
"""
# @profile
def mlp1b(x,wo,ws,tAct=['relu','relu']):
    (Ne,Np)=x.shape
    (Nn,Neb)=wo.shape
    (Ns,Nnb)=ws.shape
    
    salO=np.matmul(wo,np.concatenate((x,np.ones((1,Np)))))
    
    if tAct[0] == 'tanh':
        salO = np.tanh(salO)
    elif tAct[0] == 'sigmoid':
        salO = 1.0/(1+np.exp(-1*salO))
    elif tAct[0] == 'relu':
        salO[salO<0]=0
    
    salS = np.matmul(ws,np.concatenate((salO,np.ones((1,Np)))));
    
    if tAct[1] == 'tanh':
        salS = np.tanh(salS)
    elif tAct[1] == 'sigmoid':
        salS = 1.0/(1+np.exp(-1*salS))
    elif tAct[1] == 'relu':
        salS[salS<0]=0

    return (salS,salO)

    


"""
function [dwo,dws]=adapmlp1b_tied(x,y,wo,ws,tipoSalida);
%
% [dwo,dws]=adapmlp1b_tied(x,y,wo,ws,tSal); 
%
%     x : patrones de entrada (Ne x Np)
%     y : patrones de salida (Ns x Np)
%    wo : parámetros (pesos) de la capa oculta (Nn x Ne+1)
%    ws : parámetros (pesos) de la capa de salida (Ns x Nn+1)
%  tAct : parámetro que define el tipo de funciones de activación de cada 
%         - param = 0: lineal
%         - param = 1: tanh
%         - param = 2: logistic
%         - param = 3: rectified linear (ReLU)
%         Es un vector de dimensión 2, con los valores correspondientes a 
%         la capa oculta y a la capa de salida
%
%       Donde Ne, Ns, Np, y Nn son respectivamente el número de entradas
% 		el número de salidas, el de patrones, y el de neuronas.
%
%
%--------------------------------------------------------------------------
%         Autor: Marcelino Lázaro
%      Creación: noviembre 2014
% Actualización: noviembre 2014
%--------------------------------------------------------------------------
"""
# @profile
def adapmlp1b_tied(x,y,wo,ws,tAct):
    (Nn,Neb)=np.shape(wo)
    try: 
        (Ns,Np)=y.shape
    except:
        Ns=1
        Np=len(y)
    
    dw=np.zeros((Nn,Neb-1))
    
    (o2,o1)=mlp1b(x,wo,ws,tAct)
    d2=2*(o2-y)/Np

    if tAct[1]==1:
        d2 = d2*(1-np.square(o2))
    elif tAct[1]==2:
        d2 = d2*(o2-np.square(o2))
    elif tAct[1] == 3:
        d2[d2<0]=0
    
    dv=np.matmul(d2,np.ones((Np,1)))
    
    d1a=np.matmul(np.transpose(ws[:,0:Nn]), d2)
    if tAct[0]==1:
        d1a=np.multiply(d1a,1-np.square(o1))
    elif tAct[0]==2:
        d1a=np.multiply(d1a,o1-np.square(o1))
    elif tAct[0]==3:
        mascara=o1
        mascara[mascara<0]=0
        mascara[mascara>0]=1
        d1a=np.multiply(d1a,mascara)

    d1b=np.transpose(np.matmul(d2,o1.T))
    dw=np.matmul(d1a,x.T) + d1b
    db=np.matmul(d1a,np.ones((Np,1)))
    dwo=np.concatenate((dw,db),axis=1)
    dws=np.concatenate((dw.T,dv),axis=1) 

    return (dwo,dws)
        

"""
function dws=adapmlp0b(x,y,ws,tipoSalida);
%
% dws=adapmlp1b(x,y,wo,ws,tSal); 
%
%     x : patrones de entrada (Ne x Np)
%     y : patrones de salida (Ns x Np)
%    ws : parámetros (pesos) de la capa de salida (Ns x Ne+1)
%  tAct : parámetro que define el tipo de funciones de activación de cada 
%         - param = 0: lineal
%         - param = 1: tanh
%         - param = 2: logistic
%         - param = 3: rectified linear (ReLU)
%         Es un vector de dimensión 2, con los valores correspondientes a 
%         la capa oculta y a la capa de salida
%
%       Donde Ne, Ns, y Np son respectivamente el número de entradas
% 		el número de salidas, y el de patrones
%
%
%--------------------------------------------------------------------------
%         Autor: Marcelino Lázaro
%      Creación: noviembre 2014
% Actualización: noviembre 2014
%--------------------------------------------------------------------------
"""
# @profile
def adapmlp0b(x,y,ws,tAct):
    (Nn,Neb)=np.shape(ws)
    try: 
        (Ns,Np)=y.shape
    except:
        Ns=1
        Np=len(y)
    
    dws=np.zeros((Ns,Nn+1))
    
    o2=mlp0b(x,ws,tAct)
    
    d2=2*(o2-y)/Np
    if tAct==1:
        d2 = d2*(1-np.square(o2))
    elif tAct==2:
        d2 = d2*(o2-np.square(o2))
    elif tAct == 3:
        d2[d2<0]=0
    
    dws = np.matmul(d2,np.transpose(np.concatenate((x,np.ones((1,Np))))))
    
    return (dws)

"""
function [dwo,dws]=adapmlp1b(x,y,wo,ws,tipoSalida);
%
% [dwo,dws]=adapmlp1b(x,y,wo,ws,tSal); 
%
%     x : patrones de entrada (Ne x Np)
%     y : patrones de salida (Ns x Np)
%    wo : parámetros (pesos) de la capa oculta (Nn x Ne+1)
%    ws : parámetros (pesos) de la capa de salida (Ns x Nn+1)
%  tAct : parámetro que define el tipo de funciones de activación de cada 
%         - param = 0: lineal
%         - param = 1: tanh
%         - param = 2: logistic
%         - param = 3: rectified linear (ReLU)
%         Es un vector de dimensión 2, con los valores correspondientes a 
%         la capa oculta y a la capa de salida
%
%       Donde Ne, Ns, Np, y Nn son respectivamente el número de entradas
% 		el número de salidas, el de patrones, y el de neuronas.
%
%
%--------------------------------------------------------------------------
%         Autor: Marcelino Lázaro
%      Creación: noviembre 2014
% Actualización: noviembre 2014
%--------------------------------------------------------------------------
"""
# @profile
def adapmlp1b(x,y,wo,ws,tAct):
    (Nn,Neb)=np.shape(wo)
    try: 
        (Ns,Np)=y.shape
    except:
        Ns=1
        Np=len(y)
    
    dwo=np.zeros((Nn,Neb))
    dws=np.zeros((Ns,Nn+1))
    
    (o2,o1)=mlp1b(x,wo,ws,tAct)
    
    d2=2*(o2-y)/Np
    if tAct[1]==1:
        d2 = d2*(1-np.square(o2))
    elif tAct[1]==2:
        d2 = d2*(o2-np.square(o2))
    elif tAct[1] == 3:
        d2[d2<0]=0
    
    dws = np.matmul(d2,np.transpose(np.concatenate((o1,np.ones((1,Np))))))
    
    d1=np.matmul(np.transpose(ws[:,0:Nn]), d2)
    if tAct[0]==1:
        d1=np.multiply(d1,1-np.square(o1))
    elif tAct[0]==2:
        d1=np.multiply(d1,o1-np.square(o1))
    elif tAct[0]==3:
        mascara=o1
        mascara[mascara<0]=0
        mascara[mascara>0]=1
        d1=np.multiply(d1,mascara)
    
    dwo=np.matmul(d1,np.transpose(np.concatenate((x,np.ones((1,Np))))))
    return (dwo,dws)
  


"""
function [ws,Coste,mu]=entrena_mlp0b(x,y,ws,Niter,tAct,mu_ini);
%
%  [ws,coste,mu]=entrena_mlp1b(x,y,ws,Niter,tAct,mu_ini);
%
%     x : patrones de entrada (Ne x Np)
%     y : patrones de salida (Ns x Np)
%    ws : parámetros (pesos) de la capa de salida (Ns x Ne+1)
% Niter : número de iteraciones de gradiente
%    Nn : número de neuronas (si no se especifican wo y ws)
%  tAct : parámetro que define el tipo de funciones de activación de cada 
%         - param = 0: lineal
%         - param = 1: tanh
%         - param = 2: logistic
%         Es un vector de dimensión 2, con los valores correspondientes a 
%         la capa oculta y a la capa de salida   
%
%       Donde Ne, Ns, Np, y Nn son respectivamente el número de entradas
% 		el número de salidas, el de patrones, y el de neuronas.
%
%
%--------------------------------------------------------------------------
%         Autor: Marcelino Lázaro
%      Creación: noviembre 2014
% Actualización: abril 2015
%--------------------------------------------------------------------------
%    
"""
# @profile
def entrena_mlp0b(x,y,ws,Niter,tAct=1,mu_ini=1):
    MuCrec=1.05
    MuDec=2.0
    mu=mu_ini
    (Ne,Np)=x.shape
    try:
        (Ns,Nada)=y.shape
    except:
        Ns=1
  
   
    Coste=np.zeros(Niter+1)
    ye=mlp0b(x,ws,tAct)
    coste=np.mean(np.square(y-ye))
    Coste[0]=coste
    for kiter in range(Niter):
        dws=adapmlp0b(x,y,ws,tAct)
        wsn = ws-mu*dws
        
        ye=mlp0b(x,wsn,tAct)
        costen = np.mean(np.square(y-ye))
        if (costen>=coste):
            aumenta=True
            while aumenta:
                mu=mu/MuDec
                wsn = ws - mu*dws
                
                ye=mlp0b(x,wsn,tAct)
                costen = np.mean(np.square(y-ye))
                if ((costen<coste) |  (mu<1e-20)):
                    aumenta=False
        mu = mu*MuCrec
        ws=wsn
        coste=costen
        Coste[kiter+1]=coste
    return (ws,Coste,mu)


"""
function [wo,ws,Coste,mu]=entrena_mlp1b(x,y,wo,ws,Niter,tAct,mu_ini);
%
%  [wo,ws,coste,mu]=entrena_mlp1b(x,y,wo,ws,Niter,tAct,mu_ini);
%  [wo,ws,coste,mu]=entrena_mlp1b(x,y,Nn,[],Niter,tAct,mu_ini);
%
%     x : patrones de entrada (Ne x Np)
%     y : patrones de salida (Ns x Np)
%    wo : parámetros (pesos) de la capa oculta (Nn x Ne+1)
%    ws : parámetros (pesos) de la capa de salida (Ns x Nn+1)
% Niter : número de iteraciones de gradiente
%    Nn : número de neuronas (si no se especifican wo y ws)
%  tAct : parámetro que define el tipo de funciones de activación de cada 
%         - param = 0: lineal
%         - param = 1: tanh
%         - param = 2: logistic
%         Es un vector de dimensión 2, con los valores correspondientes a 
%         la capa oculta y a la capa de salida   
%
%       Donde Ne, Ns, Np, y Nn son respectivamente el número de entradas
% 		el número de salidas, el de patrones, y el de neuronas.
%
%
%--------------------------------------------------------------------------
%         Autor: Marcelino Lázaro
%      Creación: noviembre 2014
% Actualización: abril 2015
%--------------------------------------------------------------------------
%    
"""
# @profile
def entrena_mlp1b(x,y,wo,ws,Niter,tAct=['relu','tanh'],paso=1e-3):
    try: 
        if len(paso)==1:
            mu=paso[0]
            MuCrec=1.05
            MuDec=2.0
        else:
            mu=paso[0]
            MuCrec=paso[1]
            MuDec=paso[2]
    except:
        mu=paso
        MuCrec=1.05
        MuDec=2.0    

    (Ne,Np)=x.shape
    try:
        (Ns,Nada)=y.shape
    except:
        Ns=1
   
    if len(np.shape(wo))==1:
        Nn=wo[0]
        wo = 0.1*np.random.randn(Nn,Ne+1)
    else:
        (Nn,Ne_w)=np.shape(wo)
    if len(np.shape(ws))==1:
        ws = 0.2*np.random.rand(Ns,Nn+1)-0.1
    
    Coste=np.zeros(Niter+1)
    ye=mlp1b(x,wo,ws,tAct)[0]
    Coste[0]=np.mean(np.square(y-ye))
    coste=Coste[0]
    for kiter in range(Niter):
        (dwo,dws)=adapmlp1b(x,y,wo,ws,tAct)
        won = wo-mu*dwo
        wsn = ws-mu*dws
        
        ye=mlp1b(x,won,wsn,tAct)[0]
        costen = np.mean(np.square(y-ye))
        if (costen>=coste):
            aumenta=True
            while aumenta:
                mu=mu/MuDec
                won = wo - mu*dwo
                wsn = ws - mu*dws
                
                ye=mlp1b(x,won,wsn,tAct)[0]
                costen = np.mean(np.square(y-ye))
                if ((costen<coste) |  (mu<1e-20)):
                    aumenta=False
        mu = mu*MuCrec
        wo=won
        ws=wsn
        coste=costen
        Coste[kiter+1]=coste
    
    paso=[mu,MuCrec,MuDec]    
    return (wo,ws,Coste,paso)


    


    
"""
function [we,wr,coste,paso]=entrena_AE(x,we,wr,Niter,tAct,paso);
%
%  [we,wr]=entrena_AE(x,we,wr,Niter,tAct,paso);
%  [we,wr]=entrena_AE(x,Nn,[],Niter,tAct,paso);
%
%     x : patrones de entrada (Ne x Np)
%    we : parámetros (pesos) del autocodificador (Nn x Ne+1)
%    wr : parámetros (pesos) de la capa de reconstrucción (Ne x Nn+1)
% Niter : número de iteraciones de gradiente
%   tAct: tipo de función de activación de cada capa, incluyendo la de salida (vector 1 x 4)
%         - 0: activación lineal
%         - 1: activacion tanh
%         - 2: activacion logística
%         - 3: activación lineal rectificada (ReLU)
%   paso: parámetro de paso del algoritmo de gradiente
%
%
%       Donde Ne, Ns, Np, Nn son respectivamente el número de entradas
%       el número de salidas, el de patrones, y los de neuronas del codificador
%
%
%--------------------------------------------------------------------------
%         Autor: Marcelino Lázaro
%      Creación: diciembre 2015
% Actualización: diciembre 2015
%--------------------------------------------------------------------------
%
"""

def entrena_AE(x,we,wr,Niter,tAct=['relu','relu'],paso=[1e-3]):
    # try: 
    #     if len(paso)==1:
    #         mu=paso[0]
    #         MuCrec=1.05
    #         MuDec=2.0
    #     else:
    #         mu=paso[0]
    #         MuCrec=paso[1]
    #         MuDec=paso[2]
    # except:
    #     mu=paso
    #     MuCrec=1.05
    #     MuDec=2.0
    
    (we, wr, coste, paso)=entrena_mlp1b(x.T,x.T,we.T,wr.T,Niter,tAct,paso)
    return (we.T, wr.T, coste, paso)
    
    #(W,paso,coste,dWm) = entrena_mlp(x,x,[we,wr],Nepochs=[Niter,0],fCoste='mmse',optimizador='gradiente',pesos=[],tAct=[],paso=paso,dWm=[],pDO=[],flagEvo=True)
    #we = W[0].copy()
    #wr = W[1].copy()       
    #return (W[0].copy(), W[1].copy(), coste, paso)

    
"""
function [we,wr,coste,paso]=entrena_DAE(x,we,wr,,Niter,tAct,paso,pRuido);
%
%  [we,wr]=entrena_DAE(x,we,wr,Niter,tAct,paso);
%  [we,wr]=entrena_DAE(x,Nn,[],Niter,tAct,paso);
%
%     x : patrones de entrada (Ne x Np)
%    we : parámetros (pesos) del autocodificador (Nn x Ne+1)
%    wr : parámetros (pesos) de la capa de reconstrucción (Ne x Nn+1)
% Niter : número de iteraciones de gradiente
%   tAct: tipo de función de activación de cada capa, incluyendo la de salida (vector 1 x 4)
%         - 0: activación lineal
%         - 1: activacion tanh
%         - 2: activacion logística
%         - 3: activación lineal rectificada (ReLU)
%   paso: parámetro de paso del algoritmo de gradiente
%
%       Donde Ne, Ns, Np, Nn son respectivamente el número de entradas
%       el número de salidas, el de patrones, y los de neuronas del codificador
%
%
%--------------------------------------------------------------------------
%         Autor: Marcelino Lázaro
%      Creación: diciembre 2015
% Actualización: diciembre 2015
%--------------------------------------------------------------------------
%
"""

def entrena_DAE(x,wo,ws,Niter,tAct=['tanh','tanh'],paso=1,pRuido=0):
    
    x=x.T
    wo=wo.T
    ws=ws.T
 
    try: 
        if len(paso)==1:
            mu=paso[0]
            MuCrec=1.05
            MuDec=2.0
        else:
            mu=paso[0]
            MuCrec=paso[1]
            MuDec=paso[2]
    except:
        mu=paso
        MuCrec=1.05
        MuDec=2.0
                        
    (Ne,Np)=x.shape
    Ns=Ne
        
   
    if len(np.shape(wo))==1:
        Nn=wo[0]
        wo = 0.1*np.random.randn(Nn,Ne+1)
    else:
        (Nn,Ne_w)=np.shape(wo)
    if len(np.shape(ws))==1:
        ws = 0.2*np.random.rand(Ns,Nn+1)-0.1
    
    evoCoste=np.zeros(Niter+1)
    ye=mlp1b(x,wo,ws,tAct)[0]
    coste=np.mean(np.square(x-ye))
    evoCoste[0]=coste
    for kiter in range(Niter):
        mascara=np.random.rand(Ne,Np)
        mascara[mascara<pRuido]=0
        mascara[mascara>=pRuido]=1
        xn=np.multiply(x,mascara)
        ye=mlp1b(xn,wo,ws,tAct)[0]
        coste = np.mean(np.square(x-ye))

        (dwo,dws)=adapmlp1b(xn,x,wo,ws,tAct)
        won = wo-mu*dwo
        wsn = ws-mu*dws
        
        ye=mlp1b(xn,won,wsn,tAct)[0]
        costen = np.mean(np.square(x-ye))
        if (costen>=coste):
            aumenta=True
            while aumenta:
                mu=mu/MuDec
                won = wo - mu*dwo
                wsn = ws - mu*dws
                
                ye=mlp1b(xn,won,wsn,tAct)[0]
                costen = np.mean(np.square(x-ye))
                if ((costen<coste) |  (mu<1e-20)):
                    aumenta=False
        mu = mu*MuCrec
        wo=won
        ws=wsn
        #coste=costen
        ye=mlp1b(x,wo,ws,tAct)[0]
        coste = np.mean(np.square(x-ye))
        evoCoste[kiter+1]=coste

    return (wo.T, ws.T, evoCoste, mu)
    

def entrena_DAEb(x,wo,ws,Niters,tAct=np.array([1,1]),paso=1,pRuido=0):
 
    try: 
        if len(paso)==1:
            mu=paso[0]
            MuCrec=1.05
            MuDec=2.0
        else:
            mu=paso[0]
            MuCrec=paso[1]
            MuDec=paso[2]
    except:
        mu=paso
        MuCrec=1.05
        MuDec=2.0
        
    try:
        if len(Niters)==1:
            Niter=Niters[0]
            Nbatch=0
        else:
            Niter=Niters[0]
            Nbatch=Niters[1]
    except:
        Niter=Niters
        NBatch=0
                        
            
    (Ne,Np)=x.shape
    Ns=Ne
    
    if Nbatch==0 or Nbatch>Np:
        Nbatch=Np
        indBatch=np.arange(Np)
    else:
        ordenBatch = np.random.permutation(Np)
        Nini=0
        
   
    if len(np.shape(wo))==1:
        Nn=wo[0]
        wo = 0.1*np.random.randn(Nn,Ne+1)
    else:
        (Nn,Ne_w)=np.shape(wo)
    if len(np.shape(ws))==1:
        ws = 0.2*np.random.rand(Ns,Nn+1)-0.1
    
    evoCoste=np.zeros(Niter+1)
    ye=mlp1b(x,wo,ws,tAct)[0]
    coste=np.mean(np.square(x-ye))
    evoCoste[0]=coste
    for kiter in range(Niter):
        if Nbatch<Np:
            if (Nini+Nbatch) > Np:
                ordenBatch = np.random.permutation(Np)
                indBatch=ordenBatch[np.arange(Nbatch)]
                Nini=Nbatch
            else:        
                indBatch=ordenBatch[Nini+np.arange(Nbatch)]
                Nini=Nini+Nbatch
            
        mascara=np.random.rand(Ne,Nbatch)
        mascara[mascara<pRuido]=0
        mascara[mascara>=pRuido]=1
        xn=np.multiply(x[:,indBatch],mascara)
        ye=mlp1b(xn,wo,ws,tAct)[0]
        coste = np.mean(np.square(x[:,indBatch]-ye))

        (dwo,dws)=adapmlp1b(xn,x[:,indBatch],wo,ws,tAct)
        won = wo-mu*dwo
        wsn = ws-mu*dws
        
        ye=mlp1b(xn,won,wsn,tAct)[0]
        costen = np.mean(np.square(x[:,indBatch]-ye))
        if (costen>=coste):
            aumenta=True
            while aumenta:
                mu=mu/MuDec
                won = wo - mu*dwo
                wsn = ws - mu*dws
                
                ye=mlp1b(xn,won,wsn,tAct)[0]
                costen = np.mean(np.square(x[:,indBatch]-ye))
                if ((costen<coste) |  (mu<1e-20)):
                    aumenta=False
        mu = mu*MuCrec
        wo=won
        ws=wsn
        #coste=costen
        ye=mlp1b(x,wo,ws,tAct)[0]
        coste = np.mean(np.square(x-ye))
        evoCoste[kiter+1]=coste

    return (wo, ws, evoCoste, mu)    
    
"""
function [wo,ws,Coste,mu]=entrena_AEtied(x,wo,ws,Niter,tAct,mu_ini);
%
%  [wo,ws,coste,mu]=entrena_AEtied(x,wo,ws,Niter,tAct,mu_ini);
%  [wo,ws,coste,mu]=entrena_AEtied(x,Nn,[],Niter,tAct,mu_ini);
%
%     x : patrones de entrada (Ne x Np)
%    wo : parámetros (pesos) de la capa oculta (Nn x Ne+1)
%    ws : parámetros (pesos) de la capa de salida (Ns x Nn+1)
% Niter : número de iteraciones de gradiente
%    Nn : número de neuronas (si no se especifican wo y ws)
%  tAct : parámetro que define el tipo de funciones de activación de cada 
%         - param = 0: lineal
%         - param = 1: tanh
%         - param = 2: logistic
%         Es un vector de dimensión 2, con los valores correspondientes a 
%         la capa oculta y a la capa de salida   
%
%       Donde Ne, Ns, Np, y Nn son respectivamente el número de entradas
% 		el número de salidas, el de patrones, y el de neuronas.
%
%
%--------------------------------------------------------------------------
%         Autor: Marcelino Lázaro
%      Creación: noviembre 2014
% Actualización: abril 2015
%--------------------------------------------------------------------------
%    
"""
# @profile
def entrena_AEtied(xe,woe,wse,Niter,tAct=np.array([1,1]),mu_ini=1):
    
    x=xe.T
    wo=woe.T
    ws=wse.T
    
    
    MuCrec=1.05
    MuDec=2.0
    mu=mu_ini
    (Ne,Np)=x.shape
    Ns=Ne
   
    if len(np.shape(wo))==1:
        Nn=wo[0]
        wo = 0.1*np.random.randn(Nn,Ne+1)
    else:
        (Nn,Ne_w)=np.shape(wo)
    if len(np.shape(ws))==1:
        ws = 0.2*np.random.rand(Ns,Nn+1)-0.1
    # Se asegura que de inicio los pesos están "atados"
    ws=np.concatenate((np.transpose(wo[:,0:Ne]),ws[:,Nn-1:Nn]),axis=1)    

    evoCoste=np.zeros(Niter+1)
    ye=mlp1b(x,wo,ws,tAct)[0]
    coste=np.mean(np.square(x-ye))
    evoCoste[0]=coste
    for kiter in range(Niter):
        (dwo,dws)=adapmlp1b_tied(x,x,wo,ws,tAct)
        won = wo-mu*dwo
        wsn = ws-mu*dws
        
        ye=mlp1b(x,won,wsn,tAct)[0]
        costen = np.mean(np.square(x-ye))
        if (costen>=coste):
            aumenta=True
            while aumenta:
                mu=mu/MuDec
                won = wo - mu*dwo
                wsn = ws - mu*dws
                
                ye=mlp1b(x,won,wsn,tAct)[0]
                costen = np.mean(np.square(x-ye))
                if ((costen<coste) |  (mu<1e-20)):
                    aumenta=False
        mu = mu*MuCrec
        wo=won
        ws=wsn
        coste=costen
        evoCoste[kiter+1]=coste
    return (wo.T,ws.T,evoCoste,mu)


    
"""
function [we,wr,coste,paso]=entrena_DAEtied(x,we,wr,,Niter,tAct,paso,pRuido);
%
%  [we,wr]=entrena_DAEtied(x,we,wr,Niter,tAct,paso,pRuido);
%  [we,wr]=entrena_DAEtied(x,Nn,[],Niter,tAct,paso,pRuido);
%
%     x : patrones de entrada (Ne x Np)
%    we : parámetros (pesos) del autocodificador (Nn x Ne+1)
%    wr : parámetros (pesos) de la capa de reconstrucción (Ne x Nn+1)
% Niter : número de iteraciones de gradiente
%   tAct: tipo de función de activación de cada capa, incluyendo la de salida (vector 1 x 4)
%         - 0: activación lineal
%         - 1: activacion tanh
%         - 2: activacion logística
%         - 3: activación lineal rectificada (ReLU)
%   paso: parámetro de paso del algoritmo de gradiente
%
%       Donde Ne, Ns, Np, Nn son respectivamente el número de entradas
%       el número de salidas, el de patrones, y los de neuronas del codificador
%
%
%--------------------------------------------------------------------------
%         Autor: Marcelino Lázaro
%      Creación: diciembre 2015
% Actualización: diciembre 2015
%--------------------------------------------------------------------------
%
"""

def entrena_DAEtied(x,wo,ws,Niter,tAct=np.array([1,1]),paso=1,pRuido=0):
    
    x=x.T
    wo=wo.T
    ws=ws.T
 
    MuCrec=1.05
    MuDec=2.0
    mu=paso
    (Ne,Np)=x.shape
    Ns=Ne
   
    if len(np.shape(wo))==1:
        Nn=wo[0]
        wo = 0.1*np.random.randn(Nn,Ne+1)
    else:
        (Nn,Ne_w)=np.shape(wo)
    if len(np.shape(ws))==1:
        ws = 0.2*np.random.rand(Ns,Nn+1)-0.1
    # Se asegura que de inicio los pesos están "atados"
    ws=np.concatenate((np.transpose(wo[:,0:Ne]),ws[:,Nn-1:Nn]),axis=1)  

    evoCoste=np.zeros(Niter+1)
    ye=mlp1b(x,wo,ws,tAct)[0]
    coste=np.mean(np.square(x-ye))
    evoCoste[0]=coste
    for kiter in range(Niter):
        mascara=np.random.rand(Ne,Np)
        mascara[mascara<pRuido]=0
        mascara[mascara>=pRuido]=1
        xn=np.multiply(x,mascara)
        ye=mlp1b(xn,wo,ws,tAct)[0]
        coste = np.mean(np.square(x-ye))

        (dwo,dws)=adapmlp1b_tied(xn,x,wo,ws,tAct)
        won = wo-mu*dwo
        wsn = ws-mu*dws
        
        ye=mlp1b(xn,won,wsn,tAct)[0]
        costen = np.mean(np.square(x-ye))
        if (costen>=coste):
            aumenta=True
            while aumenta:
                mu=mu/MuDec
                won = wo - mu*dwo
                wsn = ws - mu*dws
                
                ye=mlp1b(xn,won,wsn,tAct)[0]
                costen = np.mean(np.square(x-ye))
                if ((costen<coste) |  (mu<1e-20)):
                    aumenta=False
        mu = mu*MuCrec
        wo=won
        ws=wsn
        #coste=costen
        ye=mlp1b(x,wo,ws,tAct)[0]
        coste = np.mean(np.square(x-ye))
        evoCoste[kiter+1]=coste

    return (wo.T, ws.T, evoCoste, mu)

#------------------------------------------------------------------------------
"""
%
%  [y,O] = mlpM(x,W,tAct)
%
%  Función que evalúa un MLP de M capas ocultas
%
%     x : patrones de entrada (Ne x Np)
%     W : cell array con los pesos de las capas ocultas y de salida, 1 x (M+1)
%         capa k : Nk x (Nv +1), siendo Nv en número de neuronas de la
%         capa anterior, y Nk el de la capa k
%
%  tAct : parámetro que define el tipo de funciones de activación de cada 
%         - 0: lineal
%         - 1: tanh
%         - 2: logistic
%         - 3: lineal rectificada (ReLU)
%         - 4: softmax
%         Es un vector de dimensión 6, con los valores correspondientes a 
%         las cinco capas ocultas y a la capa de salida     
%
%
%     y : valor de salida de la red (Ns x Np)
%     O : valores de salida de las distintas capas ocultas: cell 1 x M
%
%
%       Donde Ne, Ns, Np, son respectivamente el número de entradas
% 		el número de salidas, el de patrones, y el de neuronas de cada capa
%
%
%--------------------------------------------------------------------------
%         Autor: Marcelino Lázaro
%      Creación: diciembre 2015
% Actualización: febrero 2018
%--------------------------------------------------------------------------

"""
def mlpM(x,W,tAct):
    Ncapas=len(W)
    No=Ncapas-1
    
    (Ne,Np)=np.shape(x)
    
    #Os=[[] for k in range(Ncapas)]
    Os=list()
        
    for ko in range(Ncapas):
        o=np.matmul(W[ko],np.concatenate((x,np.ones((1,Np)))));
        
        if tAct[ko] == 'tanh':
            o = np.tanh(o)
        elif tAct[ko] == 'sigmoid':
            o = 1.0/(1+np.exp(-1*o))
        elif tAct[ko] == 'relu':
            o[o<0]=0.0
        elif tAct[ko] == 'softmax':
            o = np.exp(o-np.max(o,axis=0))
            o = o/np.sum(o,axis=0)
            
            
        #Os[ko]=o
        Os.append(o)
        x=o
        
    return (Os[-1],Os[0:No])    

#------------------------------------------------------------------------------
"""
%
%  [y,O] = mlp(x,W,tAct)
%
%  Función que evalúa un MLP de M capas ocultas
%
%     x : patrones de entrada (Ne x Np)
%     W : cell array con los pesos de las capas ocultas y de salida, 1 x (M+1)
%         capa k : Nk x (Nv +1), siendo Nv en número de neuronas de la
%         capa anterior, y Nk el de la capa k
%
%  tAct : parámetro que define el tipo de funciones de activación de cada 
%         - 0: lineal
%         - 1: tanh
%         - 2: logistic
%         - 3: lineal rectificada (ReLU)
%         - 4: softmax
%         Es un vector de dimensión 6, con los valores correspondientes a 
%         las cinco capas ocultas y a la capa de salida     
%
%
%     y : valor de salida de la red (Ns x Np)
%     O : valores de salida de las distintas capas ocultas: cell 1 x M
%
%
%       Donde Ne, Ns, Np, son respectivamente el número de entradas
% 		el número de salidas, el de patrones, y el de neuronas de cada capa
%
%
%--------------------------------------------------------------------------
%         Autor: Marcelino Lázaro
%      Creación: diciembre 2015
% Actualización: febrero 2018
%--------------------------------------------------------------------------

"""
def mlpN(xe,We,tAct):
    Ncapas=len(We)
    No=Ncapas-1
    
    x=xe.T
    W=We.copy()
    for ko in range(Ncapas):
        W[ko]=np.transpose(W[ko])
    
    (Ne,Np)=np.shape(x)
    
    #Os=[[] for k in range(Ncapas)]
    Os=list()
        
    for ko in range(Ncapas):
        o=np.matmul(W[ko],np.concatenate((x,np.ones((1,Np)))))
        
        if tAct[ko] == 1:
            o = np.tanh(o)
        elif tAct[ko] == 2:
            o = 1.0/(1+np.exp(-1*o))
        elif tAct[ko] == 3:
            o[o<0]=0.0
        elif tAct[ko] == 4:
            o = np.exp(o-np.max(o,axis=0))
            o = o/np.sum(o,axis=0)
            
            
        #Os[ko]=o
        Os.append(o)
        x=o
    
    for ko in range(Ncapas):
        Os[ko]=np.transpose(Os[ko])     
     
    return (Os[-1],Os[0:No])    

#------------------------------------------------------------------------------
"""
%
%  [y,O] = mlp(x,W,tAct)
%
%  Función que evalúa un MLP de M capas ocultas
%
%     x : patrones de entrada (Np x Ne)
%     W : cell array con los pesos de las capas ocultas y de salida, 1 x (M+1)
%         capa k : (Nv +1) x Nk, siendo Nv en número de neuronas de la
%         capa anterior, y Nk el de la capa k
%
%  tAct : parámetro que define el tipo de funciones de activación de cada 
%         - linear
%         - tanh
%         - sigmoid
%         - relu
%         - softmax
%         Es un vector de dimensión 6, con los valores correspondientes a 
%         las cinco capas ocultas y a la capa de salida     
%
%
%     y : valor de salida de la red (Ns x Np)
%     O : valores de salida de las distintas capas ocultas: cell 1 x M
%
%
%       Donde Ne, Ns, Np, son respectivamente el número de entradas
% 		el número de salidas, el de patrones, y el de neuronas de cada capa
%
%
%--------------------------------------------------------------------------
%         Autor: Marcelino Lázaro
%      Creación: diciembre 2015
% Actualización: febrero 2018
%--------------------------------------------------------------------------

"""
def mlp(x,W,tAct):
    Ncapas=len(W)
    No=Ncapas-1
    
    (Np,Ne)=np.shape(x)
    
    #Os=[[] for k in range(Ncapas)]
    Os=list()
        
    for ko in range(Ncapas):
        o=np.matmul(np.concatenate((x,np.ones((Np,1))),axis=1),W[ko]);
        
        if tAct[ko] == 'tanh':
            o = np.tanh(o)
        elif tAct[ko] == 'sigmoid':
            o = 1.0/(1+np.exp(-1*o))
        elif tAct[ko] == 'relu':
            o[o<0]=0.0
        elif tAct[ko] == 'softmax':
            o = np.exp(np.transpose(o.T-np.max(o,axis=1)))
            o = np.transpose(o.T/np.sum(o,axis=1))
            
            
        #Os[ko]=o
        Os.append(o)
        x=o
        
    return (Os[-1],Os[0:No])    
    
    
#------------------------------------------------------------------------------
"""
  (W,paso,evoCoste,dWm)=entrena_mlpN(x,y,W,Nepochs)

  Función que entrena una red MLP de una o varias capas ocultas con varias 
  funciones de coste y posibles métodos de optimización

     x : patrones de entrada (Ne x Np)
     y : patrones de salida (Ns x Np)
     W : lista con los pesos de las capas ocultas y de salida, 1 x (H+1)
         capa k : Nk x (Nv +1), siendo Nv en número de neuronas de la
         capa anterior, y Nk el de la capa k
 Nepochs: lista con dos parámetros, el número de iteraciones de adaptación y 
          el tamaño del mini-batch Nep=[Niter,Nbatch]. 
          Nbatch=0 implica un entrenamiento puramente batch
 
  Existen una serie de parámetros adicionales con valores por defecto

fCoste : función de coste a minimizar
         - 'mmse' : error cuadrático medio (valor por defecto)
         - 'entropia': entropía cruzada (a utilizar combinada con salida softmax)
         - 'wmmse' : error cuadrático medio ponderado
           En este caso hay que especificar el parámetro 'pesos' (NsxNp)            

optimizador : tipo de optimizador empleado para la minimización de la función
              de coste
         - 'gradiente' : descenso de gradiente con paso adaptativo
            * paso=[mu,muCrec,muDec] : define el paso inicial, y los factores de
                                       incremento y decremento ([1e-3,1.05,2.0])
         - 'momento' : gradiente con término de momento
            * paso=[mu,momento] : define el paso y el término de momento
            
  tAct : parámetro que define el tipo de funciones de activación de cada 
         - 0: lineal
         - 1: tanh
         - 2: logistic
         - 3: lineal rectificada (ReLU)
         - 4: softmax
         Es un vector de dimensión Ncapas, con los valores correspondientes a 
         las capas ocultas y a la capa de salida
                   
  pDO : probabilidades de Drop-Out. Lista de H+1 elementos (incluye la entrada) 
  
  dWm : término inicial para el término de momento en el proceso de adaptación
        lista con elementos del mismo tamaño de los parámetros W
        
flagEvo : indica si se calcula el coste tras cada iteración (True) o no (False)
          En caso afirmativo, se almacena en la variable de salida 'evoCoste'
         
       Donde Ne, Ns y Np son respectivamente el número de entradas
 		el número de salidas y el de patrones.

%--------------------------------------------------------------------------
%         Autor: Marcelino Lázaro
%      Creación: febrero 2017
% Actualización: mayo 2018
%--------------------------------------------------------------------------
"""    

def entrena_mlpN(x,y,W,Nepochs,fCoste='mmse',optimizador='gradiente',pesos=[],tAct=[],paso=[1e-3,1.05,2.0],dWm=[],pDO=[],flagEvo=True):
    
    (Ne,Np)=np.shape(x)
    try:
        (Ns,nada)=np.shape(y)
    except:
        Ns=1
        
    Ncapas=len(W)
    Nini=0    
    Niter=Nepochs[0]
    Nbatch=Nepochs[1]
    if Nbatch==0:
        Nbatch=Np
        
    esMetodoBayes=False
        
    if optimizador == 'gradiente':    
        if len(paso) == 1:
            mu=paso[0]
            muCrec=1.05
            muDec=2.0
        else:
            mu=paso[0]
            muCrec=paso[1]
            muDec=paso[2]
            
    elif optimizador=='momento':
        if len(paso) == 1:
            mu=paso[0]
            momento=0.9
        else:
            mu=paso[0]
            momento=paso[1]
                    
        
    if len(tAct)==0:
        tAct=[3 for k in range(Ncapas-1)]
        tAct.append(1)      
    if len(pDO)==0:
        pDO=[0 for k in range(Ncapas)]
        
        
    ws=W[-1]
    if len(ws)==0:
        NnS=[Ne]
        for ko in range(Ncapas-1):
            NnS.append(W[ko])
            
        NnS.append(Ns)

        for ko in range(Ncapas):
            W[ko]=np.sqrt(1.0/NnS[ko])*(2.0*np.random.rand(NnS[ko+1],NnS[ko]+1)-1.0)

    if len(dWm)==0:
        for ko in range(Ncapas):
            dWm.append(np.zeros((np.shape(W[ko]))))
            
    evoCoste=np.zeros(Niter+1)
    (ye,Os)=mlp(x,W,tAct)
    # Se evalúa el coste inicial
    if fCoste=='mmse':
        coste=np.mean(np.square(y-ye))
    elif fCoste=='wmmse':
        coste=np.mean(np.multiply(pesos,np.square(y-ye)))
    elif fCoste=='entropia':
        coste=np.mean(-np.sum(np.multiply(y,np.log(np.maximum(ye,1e-20)))))
                
    v0=find(y[0,:],0)        
    if len(v0) == 0:
        v0=find(y[0,:],-1)
                        
    N0=len(v0)
    v1=find(y[0,:],1)
    N1=len(v1)
    P1=1.0*len(v1)/(len(v0)+len(v1))
        
        
    if Nbatch>=Np:
        Nbatch=Np
        if esMetodoBayes:
            indBatch=np.concatenate((v0,v1),axis=0)
            v0b=v0
            v1b=v1
            N0b=N0
            N1b=N1
        else:
            indBatch=range(Np)                                
    else:
        ordenBatch = np.random.permutation(Np)
        if esMetodoBayes:
            N1b=int(math.ceil(P1*Nbatch))
            N0b=Nbatch-N1b
                            
    evoCoste[0]=coste

    dW=[[] for k in range(Ncapas)]        
    Wdo=[[] for k in range(Ncapas)]    
    Wdon=[[] for k in range(Ncapas)]   
    mX=np.diag(np.ones(Ne)).astype(float)
    mW=[[] for k in range(Ncapas)]        
    for ko in range(Ncapas):
        W[ko]=W[ko]/(1-pDO[ko])
        mW[ko]=np.ones(np.shape(W[ko]))
    
    #mW[Ncapas-1]=np.ones(np.shape(W[Ncapas-1]))      
    #--------------------------------------------------------------------------
    # Se inicia el procedimiento de entrenamiento
    #--------------------------------------------------------------------------
    for kiter in range(Niter):
        if Nbatch < Np:
            if esMetodoBayes:
                v0b=v0[list(np.random.choice(N0,N0b,replace=False))]
                v1b=v1[list(np.random.choice(N1,N1b,replace=False))]
                indBatch=np.concatenate((v0b,v1b),axis=0)       
            else:    
                if (Nini+Nbatch) > Np:
                    ordenBatch = np.random.permutation(Np)
                    indBatch=ordenBatch[np.arange(Nbatch)]
                    Nini=Nbatch
                else:        
                    indBatch=ordenBatch[Nini+np.arange(Nbatch)]
                    Nini=Nini+Nbatch
                    
        if np.sum(pDO)>0:                                    
            mX=np.diag(np.random.uniform(0,1,Ne)>pDO[0]).astype(float)
            for ko in range(Ncapas-1):
                (Nb,Na)=np.shape(W[ko])
                mW[ko]=np.matmul(np.diag((np.random.uniform(0,1,Nb)>pDO[ko+1])).astype(float), np.ones((Nb,Na)))
                Wdo[ko]=np.multiply(W[ko],mW[ko])
                
            Wdo[Ncapas-1]=W[Ncapas-1]
        else:
            Wdo=list(W)
            
        (ye,Os)=mlp(np.matmul(mX,x[:,indBatch]),Wdo,tAct)
        #----------------------------------------------------------------------
        # Cálculo coste inicial y gradientes
        #----------------------------------------------------------------------
        if fCoste=='mmse':
            coste=np.mean(np.square(y[:,indBatch]-ye))
            d=2.0*(ye-y[:,indBatch])
        elif fCoste=='wmmse':
            coste=np.mean(np.multiply(pesos[:,indBatch],np.square(y[:,indBatch]-ye)))
            d=2.0*np.multiply((ye-y[:,indBatch]),pesos[:,indBatch])/float(Nbatch)
        elif fCoste=='entropia':
            coste=np.mean(-np.sum(np.multiply(y[:,indBatch],np.log(np.maximum(ye,1e-20)))))
            d=ye-y[:,indBatch]
                                        
        if tAct[-1]==3:
            mascara=ye[:]
            mascara[mascara<=0]=0
            mascara[mascara>0]=1
            d=np.multiply(d,mascara)
        elif tAct[-1]==1:
            d=np.multiply(d,1-np.square(ye))    
        elif tAct[-1]==0:
            d=np.multiply(d,np.ones(np.shape(ye)))    
        elif tAct[-1]==2:
            d=np.multiply(d,ye-np.square(ye))    
            
                
        dW[Ncapas-1]=np.matmul(d,np.transpose(np.concatenate((Os[-1],np.ones((1,Nbatch))))))
                
        o=Os[-1]
        
        for ko in range(Ncapas-2,-1,-1):
            wp=W[ko+1]
            d=np.matmul(np.transpose(wp[:,0:-1]),d)
            
            if tAct[ko]==3:
                mascara=o[:]
                mascara[mascara<=0]=0
                mascara[mascara>0]=1
                d=np.multiply(d,mascara)
            elif tAct[ko]==1:
                d=np.multiply(d,1-np.square(o))
            elif tAct[ko]==2:
                d=np.multiply(d,o-np.square(o))
            elif tAct[ko]==0:
                d=np.multiply(d,np.ones(np.shape(ye)))    
                
            if ko == 0:
                o=x[:,indBatch]
            else:
                o=Os[ko-1]
             
            dW[ko]=np.matmul(d,np.transpose(np.concatenate((o,np.ones((1,Nbatch))))))
            #dp=d
        #----------------------------------------------------------------------    
        # Nuevos pesos iteración
        #----------------------------------------------------------------------
        for ko in range(Ncapas):
            if optimizador=='gradiente':
                Wdon[ko]=Wdo[ko]-mu*np.multiply(dW[ko],mW[ko])
            elif optimizador=='momento':
                dWm[ko]=dWm[ko]*momento+mu*np.multiply(dW[ko],mW[ko])
                W[ko]=W[ko]-np.multiply(dWm[ko],mW[ko])
                
        #----------------------------------------------------------------------
        # Estima de los costes y ajuste del parámetro de paso (GRADIENTE)
        #----------------------------------------------------------------------            
        if optimizador == 'gradiente':            
            (ye,Os)=mlp(np.matmul(mX,x[:,indBatch]),Wdon,tAct)
            if fCoste=='mmse':
                costen=np.mean(np.square(y[:,indBatch]-ye))
            elif fCoste=='wmmse':
                costen=np.mean(np.multiply(pesos[:,indBatch],np.square(y[:,indBatch]-ye)))
            elif fCoste=='entropia':
                costen=np.mean(-np.sum(np.multiply(y[:,indBatch],np.log(np.maximum(ye,1e-20)))))
            #------------------------------------------------------------------    
            # Ajuste del parámetro de paso
            #------------------------------------------------------------------
            if costen >= coste:
                aumenta=True
    
                while aumenta:
                    mu=mu/muDec
    
                    for ko in range(Ncapas):
                        #Wn[kb]=W[kb]-mu*dW[kb]
                        Wdon[ko]=Wdo[ko]-mu*np.multiply(dW[ko],mW[ko])
                        
                    (ye,Os)=mlp(np.matmul(mX,x[:,indBatch]),Wdon,tAct)
                    if fCoste=='mmse':
                        costen=np.mean(np.square(y[:,indBatch]-ye))
                    elif fCoste=='wmmse':
                        costen=np.mean(np.multiply(pesos[:,indBatch],np.square(y[:,indBatch]-ye)))
                    elif fCoste=='entropia':
                        costen=np.mean(-np.sum(np.multiply(y[:,indBatch],np.log(np.maximum(ye,1e-20)))))
                        
                    if (costen < coste):
                        aumenta = False
                    elif  (mu<1e-20):
                        for ko in range(Ncapas):
                            W[ko]=W[ko]*(1-pDO[ko])
                        
                        if optimizador == 'gradiente':
                            paso=[mu,muCrec,muDec]
                        
                        evoCoste[kiter+1:Niter]=coste    
                        print("Parámetro de paso en el límite (cero) ...")                        
                        return (W,paso,evoCoste,dWm)
            
            for ko in range(Ncapas):
                W[ko]=W[ko]-mu*np.multiply(dW[ko],mW[ko])
                
            mu=mu*muCrec                    
        #----------------------------------------------------------------------
        # Actualización del coste (global)
        #----------------------------------------------------------------------
        if flagEvo:
            for ko in range(Ncapas):
                Wdon[ko]=W[ko]*(1-pDO[ko])
            
            (ye,Os)=mlp(x,Wdon,tAct)
            if fCoste=='mmse':
                coste=np.mean(np.square(y-ye))
            elif fCoste=='wmmse':
                coste=np.mean(np.multiply(pesos,np.square(y-ye)))
            elif fCoste=='entropia':
                coste=np.mean(-np.sum(np.multiply(y,np.log(np.maximum(ye,1e-20)))))
                
            evoCoste[kiter+1]=coste
    
    for ko in range(Ncapas):
        W[ko]=W[ko]*(1-pDO[ko])
    
    if optimizador == 'gradiente':
        paso=[mu,muCrec,muDec]
    return (W,paso,evoCoste,dWm)    
    
#------------------------------------------------------------------------------
    
#------------------------------------------------------------------------------
"""
  (W,paso,evoCoste,dWm)=entrena_mlp(x,y,W,Nepochs)

  Función que entrena una red MLP de una o varias capas ocultas con varias 
  funciones de coste y posibles métodos de optimización

     x : patrones de entrada (Ne x Np)
     y : patrones de salida (Ns x Np)
     W : lista con los pesos de las capas ocultas y de salida, 1 x (H+1)
         capa k : Nk x (Nv +1), siendo Nv en número de neuronas de la
         capa anterior, y Nk el de la capa k
 Nepochs: lista con dos parámetros, el número de iteraciones de adaptación y 
          el tamaño del mini-batch Nep=[Niter,Nbatch]. 
          Nbatch=0 implica un entrenamiento puramente batch
 
  Existen una serie de parámetros adicionales con valores por defecto

fCoste : función de coste a minimizar
         - 'mmse' : error cuadrático medio (valor por defecto)
         - 'entropia': entropía cruzada (a utilizar combinada con salida softmax)
         - 'wmmse' : error cuadrático medio ponderado
           En este caso hay que especificar el parámetro 'pesos' (NsxNp)            

optimizador : tipo de optimizador empleado para la minimización de la función
              de coste
         - 'gradiente' : descenso de gradiente con paso adaptativo
            * paso=[mu,muCrec,muDec] : define el paso inicial, y los factores de
                                       incremento y decremento ([1e-3,1.05,2.0])
         - 'momento' : gradiente con término de momento
            * paso=[mu,momento] : define el paso y el término de momento
            
  tAct : parámetro que define el tipo de funciones de activación de cada 
         - 0: lineal
         - 1: tanh
         - 2: logistic
         - 3: lineal rectificada (ReLU)
         - 4: softmax
         Es un vector de dimensión Ncapas, con los valores correspondientes a 
         las capas ocultas y a la capa de salida
                   
  pDO : probabilidades de Drop-Out. Lista de H+1 elementos (incluye la entrada) 
  
  dWm : término inicial para el término de momento en el proceso de adaptación
        lista con elementos del mismo tamaño de los parámetros W
        
flagEvo : indica si se calcula el coste tras cada iteración (True) o no (False)
          En caso afirmativo, se almacena en la variable de salida 'evoCoste'
         
       Donde Ne, Ns y Np son respectivamente el número de entradas
 		el número de salidas y el de patrones.

%--------------------------------------------------------------------------
%         Autor: Marcelino Lázaro
%      Creación: febrero 2017
% Actualización: mayo 2018
%--------------------------------------------------------------------------
"""    

def entrena_mlp(xe,ye,We,Nepochs,fCoste='mmse',optimizador='gradiente',pesos=[],tAct=[],paso=[1e-3,1.05,2.0],dWm=[],pDO=[],flagEvo=True):
    
    x=xe.T
    y=ye.T
    W=We.copy()
    
    (Ne,Np)=np.shape(x)
    try:
        (Ns,nada)=np.shape(y)
    except:
        Ns=1
        
    Ncapas=len(W)
    
    for ko in range(Ncapas):
        W[ko]=np.transpose(W[ko])
    
    
    Nini=0    
    Niter=Nepochs[0]
    Nbatch=Nepochs[1]
    if Nbatch==0:
        Nbatch=Np
        
    esMetodoBayes=False
        
    if optimizador == 'gradiente':    
        if len(paso) == 1:
            mu=paso[0]
            muCrec=1.05
            muDec=2.0
        else:
            mu=paso[0]
            muCrec=paso[1]
            muDec=paso[2]
            
    elif optimizador=='momento':
        if len(paso) == 1:
            mu=paso[0]
            momento=0.9
        else:
            mu=paso[0]
            momento=paso[1]
                    
        
    if len(tAct)==0:
        tAct=['relu' for k in range(Ncapas-1)]
        tAct.append('tanh')      
    if len(pDO)==0:
        pDO=[0 for k in range(Ncapas)]
        
        
    ws=W[-1]
    if len(ws)==0:
        NnS=[Ne]
        for ko in range(Ncapas-1):
            NnS.append(W[ko])
            
        NnS.append(Ns)

        for ko in range(Ncapas):
            W[ko]=np.sqrt(1.0/NnS[ko])*(2.0*np.random.rand(NnS[ko+1],NnS[ko]+1)-1.0)


    if len(dWm)<len(W):
        dWm=[]
        for ko in range(Ncapas):
            dWm.append(np.zeros((np.shape(W[ko]))))
            
            
    evoCoste=np.zeros(Niter+1)
    (ye,Os)=mlpM(x,W,tAct)
    # Se evalúa el coste inicial
    if fCoste=='mmse':
        coste=np.mean(np.square(y-ye))
    elif fCoste=='wmmse':
        coste=np.mean(np.multiply(pesos,np.square(y-ye)))
    elif fCoste=='entropia':
        coste=np.mean(-np.sum(np.multiply(y,np.log(np.maximum(ye,1e-20)))))
                
    v0=find(y[0,:],0)        
    if len(v0) == 0:
        v0=find(y[0,:],-1)
                        
    N0=len(v0)
    v1=find(y[0,:],1)
    N1=len(v1)
    P1=1.0*len(v1)/(len(v0)+len(v1))
        
        
    if Nbatch>=Np:
        Nbatch=Np
        if esMetodoBayes:
            indBatch=np.concatenate((v0,v1),axis=0)
            v0b=v0
            v1b=v1
            N0b=N0
            N1b=N1
        else:
            indBatch=range(Np)                                
    else:
        ordenBatch = np.random.permutation(Np)
        if esMetodoBayes:
            N1b=int(math.ceil(P1*Nbatch))
            N0b=Nbatch-N1b
                            
    evoCoste[0]=coste

    dW=[[] for k in range(Ncapas)]        
    Wdo=[[] for k in range(Ncapas)]    
    Wdon=[[] for k in range(Ncapas)]   
    mX=np.diag(np.ones(Ne)).astype(float)
    mW=[[] for k in range(Ncapas)]        
    for ko in range(Ncapas):
        W[ko]=W[ko]/(1-pDO[ko])
        mW[ko]=np.ones(np.shape(W[ko]))
    
    #mW[Ncapas-1]=np.ones(np.shape(W[Ncapas-1]))      
    #--------------------------------------------------------------------------
    # Se inicia el procedimiento de entrenamiento
    #--------------------------------------------------------------------------
    for kiter in range(Niter):
        if Nbatch < Np:
            if esMetodoBayes:
                v0b=v0[list(np.random.choice(N0,N0b,replace=False))]
                v1b=v1[list(np.random.choice(N1,N1b,replace=False))]
                indBatch=np.concatenate((v0b,v1b),axis=0)       
            else:    
                if (Nini+Nbatch) > Np:
                    ordenBatch = np.random.permutation(Np)
                    indBatch=ordenBatch[np.arange(Nbatch)]
                    Nini=Nbatch
                else:        
                    indBatch=ordenBatch[Nini+np.arange(Nbatch)]
                    Nini=Nini+Nbatch
                    
        if np.sum(pDO)>0:                                    
            mX=np.diag(np.random.uniform(0,1,Ne)>pDO[0]).astype(float)
            for ko in range(Ncapas-1):
                (Nb,Na)=np.shape(W[ko])
                mW[ko]=np.matmul(np.diag((np.random.uniform(0,1,Nb)>pDO[ko+1])).astype(float), np.ones((Nb,Na)))
                Wdo[ko]=np.multiply(W[ko],mW[ko])
                
            Wdo[Ncapas-1]=W[Ncapas-1]
        else:
            Wdo=list(W)
            
        (ye,Os)=mlpM(np.matmul(mX,x[:,indBatch]),Wdo,tAct)
        #----------------------------------------------------------------------
        # Cálculo coste inicial y gradientes
        #----------------------------------------------------------------------
        if fCoste=='mmse':
            coste=np.mean(np.square(y[:,indBatch]-ye))
            d=2.0*(ye-y[:,indBatch])#/float(Nbatch)
        elif fCoste=='wmmse':
            coste=np.mean(np.multiply(pesos[:,indBatch],np.square(y[:,indBatch]-ye)))
            d=2.0*np.multiply((ye-y[:,indBatch]),pesos[:,indBatch])/float(Nbatch)
        elif fCoste=='entropia':
            coste=np.mean(-np.sum(np.multiply(y[:,indBatch],np.log(np.maximum(ye,1e-20)))))
            d=ye-y[:,indBatch]
                                        
        if tAct[-1]=='relu':
            mascara=ye[:]
            mascara[mascara<=0]=0
            mascara[mascara>0]=1
            d=np.multiply(d,mascara)
        elif tAct[-1]=='tanh':
            d=np.multiply(d,1-np.square(ye))    
        elif tAct[-1]=='linear':
            d=np.multiply(d,np.ones(np.shape(ye)))    
        elif tAct[-1]=='sigmoid':
            d=np.multiply(d,ye-np.square(ye))    
            
                
        dW[Ncapas-1]=np.matmul(d,np.transpose(np.concatenate((Os[-1],np.ones((1,Nbatch))))))
                
        o=Os[-1]
        
        for ko in range(Ncapas-2,-1,-1):
            wp=W[ko+1]
            d=np.matmul(np.transpose(wp[:,0:-1]),d)
            
            if tAct[ko]=='relu':
                mascara=o[:]
                mascara[mascara<=0]=0
                mascara[mascara>0]=1
                d=np.multiply(d,mascara)
            elif tAct[ko]=='tanh':
                d=np.multiply(d,1-np.square(o))
            elif tAct[ko]=='sigmoid':
                d=np.multiply(d,o-np.square(o))
            elif tAct[ko]=='linear':
                d=np.multiply(d,np.ones(np.shape(ye)))    
                
            if ko == 0:
                o=x[:,indBatch]
            else:
                o=Os[ko-1]
             
            dW[ko]=np.matmul(d,np.transpose(np.concatenate((o,np.ones((1,Nbatch))))))
            #dp=d
        #----------------------------------------------------------------------    
        # Nuevos pesos iteración
        #----------------------------------------------------------------------
        for ko in range(Ncapas):
            if optimizador=='gradiente':
                Wdon[ko]=Wdo[ko]-mu*np.multiply(dW[ko],mW[ko])
            elif optimizador=='momento':
                dWm[ko]=dWm[ko]*momento+mu*np.multiply(dW[ko],mW[ko])
                W[ko]=W[ko]-np.multiply(dWm[ko],mW[ko])
                
        #----------------------------------------------------------------------
        # Estima de los costes y ajuste del parámetro de paso (GRADIENTE)
        #----------------------------------------------------------------------            
        if optimizador == 'gradiente':            
            (ye,Os)=mlpM(np.matmul(mX,x[:,indBatch]),Wdon,tAct)
            if fCoste=='mmse':
                costen=np.mean(np.square(y[:,indBatch]-ye))
            elif fCoste=='wmmse':
                costen=np.mean(np.multiply(pesos[:,indBatch],np.square(y[:,indBatch]-ye)))
            elif fCoste=='entropia':
                costen=np.mean(-np.sum(np.multiply(y[:,indBatch],np.log(np.maximum(ye,1e-20)))))
            #------------------------------------------------------------------    
            # Ajuste del parámetro de paso
            #------------------------------------------------------------------
            if costen >= coste:
                aumenta=True
    
                while aumenta:
                    mu=mu/muDec
    
                    for ko in range(Ncapas):
                        #Wn[kb]=W[kb]-mu*dW[kb]
                        Wdon[ko]=Wdo[ko]-mu*np.multiply(dW[ko],mW[ko])
                        
                    (ye,Os)=mlpM(np.matmul(mX,x[:,indBatch]),Wdon,tAct)
                    if fCoste=='mmse':
                        costen=np.mean(np.square(y[:,indBatch]-ye))
                    elif fCoste=='wmmse':
                        costen=np.mean(np.multiply(pesos[:,indBatch],np.square(y[:,indBatch]-ye)))
                    elif fCoste=='entropia':
                        costen=np.mean(-np.sum(np.multiply(y[:,indBatch],np.log(np.maximum(ye,1e-20)))))
                        
                    if (costen < coste):
                        aumenta = False
                    elif  (mu<1e-20):
                        for ko in range(Ncapas):
                            W[ko]=W[ko]*(1-pDO[ko])
                        
                        if optimizador == 'gradiente':
                            paso=[mu,muCrec,muDec]
                        
                        evoCoste[kiter+1:Niter]=coste    
                        print("Parámetro de paso en el límite (cero) ...")
                        
                        for ko in range(Ncapas):
                            W[ko]=W[ko]*(1-pDO[ko])
                        
                        if optimizador == 'gradiente':
                            paso=[mu,muCrec,muDec]
                            
                        for ko in range(Ncapas):
                            W[ko]=np.transpose(W[ko])
                            dWm[ko]=np.transpose(dWm[ko])

                        
                        return (W,paso,evoCoste,dWm)
            
            for ko in range(Ncapas):
                W[ko]=W[ko]-mu*np.multiply(dW[ko],mW[ko])
                
            mu=mu*muCrec                    
        #----------------------------------------------------------------------
        # Actualización del coste (global)
        #----------------------------------------------------------------------
        if flagEvo:
            for ko in range(Ncapas):
                Wdon[ko]=W[ko]*(1-pDO[ko])
            
            (ye,Os)=mlpM(x,Wdon,tAct)
            if fCoste=='mmse':
                coste=np.mean(np.square(y-ye))
            elif fCoste=='wmmse':
                coste=np.mean(np.multiply(pesos,np.square(y-ye)))
            elif fCoste=='entropia':
                coste=np.mean(-np.sum(np.multiply(y,np.log(np.maximum(ye,1e-20)))))
                
            evoCoste[kiter+1]=coste
    
    for ko in range(Ncapas):
        W[ko]=W[ko]*(1-pDO[ko])
    
    if optimizador == 'gradiente':
        paso=[mu,muCrec,muDec]
        
    for ko in range(Ncapas):
        W[ko]=np.transpose(W[ko])
        dWm[ko]=np.transpose(dWm[ko])
    
    return (W,paso,evoCoste,dWm)    
    
#------------------------------------------------------------------------------