Z_continuo.ipynb

{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"name":"Z_continuo.ipynb","provenance":[],"collapsed_sections":["u6x8F1nEQ2qg","4NdDDTiHUlNf","V80Eh7BjUlOA","TOonSiZhUlQH","hbQb89YcONWQ","Z_ZEO3LPUlTW","9hOvKVL82FFH","bPhJF38kqsHD","gBKXLaiYUlPQ","_ErWa2tV7mYg"]},"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.5.2"}},"cells":[{"cell_type":"markdown","metadata":{"id":"u6x8F1nEQ2qg"},"source":["# Importo y defino todo "]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"z9zR7bT-UtUi","executionInfo":{"status":"ok","timestamp":1623426378088,"user_tz":180,"elapsed":3154,"user":{"displayName":"Nicolás Chehebar","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GjWN7wrOddIeQwbW6mzML-_Z0jTJLZ1cnE5tDuipw=s64","userId":"03611693109350494416"}},"outputId":"4c77fd36-2c00-4022-c170-44bcaac5c082"},"source":["!pip install fermat"],"execution_count":18,"outputs":[{"output_type":"stream","text":["Requirement already satisfied: fermat in /usr/local/lib/python3.7/dist-packages (0.2.5)\n","Requirement already satisfied: numpy in /usr/local/lib/python3.7/dist-packages (from fermat) (1.19.5)\n","Requirement already satisfied: scikit-learn in /usr/local/lib/python3.7/dist-packages (from fermat) (0.22.2.post1)\n","Requirement already satisfied: scipy in /usr/local/lib/python3.7/dist-packages (from fermat) (1.4.1)\n","Requirement already satisfied: joblib>=0.11 in /usr/local/lib/python3.7/dist-packages (from scikit-learn->fermat) (1.0.1)\n"],"name":"stdout"}]},{"cell_type":"code","metadata":{"collapsed":true,"id":"NZkzOPezUlKj","executionInfo":{"status":"ok","timestamp":1623426378103,"user_tz":180,"elapsed":77,"user":{"displayName":"Nicolás Chehebar","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GjWN7wrOddIeQwbW6mzML-_Z0jTJLZ1cnE5tDuipw=s64","userId":"03611693109350494416"}}},"source":["import fermat as F\n","import numpy as np\n","import scipy.sparse as sp\n","import time\n","from scipy.spatial import  distance_matrix\n","from sklearn.manifold import TSNE\n","from scipy.optimize import minimize\n","from sklearn.neighbors import KDTree\n","from scipy.optimize import LinearConstraint\n","from scipy.sparse import csr_matrix\n","from scipy.spatial.distance import pdist,cdist, squareform\n","import scipy\n","import warnings\n","warnings.filterwarnings('ignore')\n","\n","\n","from mpl_toolkits.mplot3d import Axes3D\n","import matplotlib.pyplot as plt\n","\n","def cero(a,b):\n","  return 0, np.zeros(a.shape)"],"execution_count":19,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"7TfnCVmuUlLz"},"source":["importe lo necesario y armo las funciones que me dicen el costo de transportar una serie de puntos con fermat dada la matriz de distancias de fermat (lo que usaba en transporte optimo antes), tambien hago su gradiente y pruebo"]},{"cell_type":"code","metadata":{"collapsed":true,"id":"SyVlktUeUlL8","executionInfo":{"status":"ok","timestamp":1623426379022,"user_tz":180,"elapsed":973,"user":{"displayName":"Nicolás Chehebar","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GjWN7wrOddIeQwbW6mzML-_Z0jTJLZ1cnE5tDuipw=s64","userId":"03611693109350494416"}}},"source":["#Builds fermat graph given data, alpha, method (and if 'D' is chosen with k_fermat neighbours)\n","#RETURNS: fermat distance matrix between points of data and a KD-Tree \n","\n","\n","def fermat_graph(data,alpha,method='FW',k_fermat=30):\n","    model=F.Fermat(alpha, path_method=method,k=k_fermat)\n","    model.fit(np.matrix(distance_matrix(data,data)))\n","    fermat_distances = model.get_distances()\n","    tree = KDTree(data)\n","    return fermat_distances,tree\n","              \n","\n","#Computes fermat_cost of transporting initialData to x with the given fermat_distances\n","#initial data will be points\n","\n","#It uses the KD-Tree built with fermat_graph() as well as a parameter 'k' which is \n","#the amount of neighbours to look at when getting the fermat distance for a point that\n","#is not in the graph (including it may be computationally expensive)\n","#It also recieves the original data and the alpha used to build fermat_distances\n","              \n","#RETURNS: fermat cost and gradient of that cost\n","def fermat_cost_cuad(x,initialData,dimension,fermat_distances,tree,data,k,alpha):\n","    assert(len(x)==(initialData.shape[0])*dimension)\n","    real_x=x.reshape(-1,dimension)     \n","    res=0\n","    gradient=np.zeros(x.shape)\n","    for i in range(initialData.shape[0]):\n","        dist,grad_idx=fermat_dist(initialData[i,:],real_x[i,:],fermat_distances,tree,k,alpha) \n","        q=data[grad_idx,:] \n","        local_grad=alpha*(np.linalg.norm(q-real_x[i,:])**(alpha-2))*(real_x[i,:]-q)*2*dist\n","        res+=dist**2              \n","        gradient[dimension*i:dimension*(i+1)]=local_grad\n","    return res,gradient\n","#ojoo cambie para la DISTANCIA AL CUADRADO!! \n","\n","\n","\n","#RETURNS: fermat cost and gradient of that cost\n","def fermat_cost(x,initialData,dimension,fermat_distances,tree,data,k,alpha):\n","    assert(len(x)==(initialData.shape[0])*dimension)\n","    real_x=x.reshape(-1,dimension)     \n","    res=0\n","    gradient=np.zeros(x.shape)\n","    for i in range(initialData.shape[0]):\n","        dist,grad_idx=fermat_dist(initialData[i,:],real_x[i,:],fermat_distances,tree,k,alpha) \n","        q=data[grad_idx,:] \n","        local_grad=alpha*(np.linalg.norm(q-real_x[i,:])**(alpha-2))*(real_x[i,:]-q)\n","        res+=dist              \n","        gradient[dimension*i:dimension*(i+1)]=local_grad\n","    return res,gradient\n","\n","\n","\n","#Computes fermat distance between start  and end using fermat_distances\n","              \n","#it also recieves the KD-tree, numbre of neighbours to look at, and alpha used in fermat\n","              \n","#RETURNS: fermat distance between start and end, as well as the index of the k-nearest neighbour\n","#of end with smallest fermat distance to start\n","def fermat_dist(start,end,fermat_distances,tree,k,alpha):\n","    \n","    start_euclid_dist,start_ind=tree.query([start],k=1)\n","    idx=start_ind[0,0]\n","    start_euclid_dist=np.power(start_euclid_dist,alpha)\n","    \n","    end_euclid_dist,end_ind=tree.query([end],k=k)\n","    end_euclid_dist=np.power(end_euclid_dist,alpha)\n","        \n","    dist=min([ fermat_distances[idx,end_ind[0,i]]+end_euclid_dist[0,i] for i in range(end_euclid_dist.shape[1])])\n","    \n","    minim=np.argmin([ fermat_distances[idx,end_ind[0,i]] for i in range(end_ind.shape[1])]) \n","    grad_idx=end_ind[0,minim] #index of k-nearest neighbour with smallest fermat_distance to idx.\n","    \n","    return dist+start_euclid_dist[0,0],grad_idx\n","              \n","\n","#Here initialData are the points, and x the trasnportation of those\n","##UNUSED\n","def euclidean_cost(x,initialData):\n","    dim=initialData.shape[1]\n","    assert(initialData.shape[0]*dim==len(x))\n","    real_x=x.reshape(-1,dim)\n","    res=0\n","    gradient=np.zeros(x.shape)\n","    for i in range(initialData.shape[0]):\n","        res+=np.linalg.norm(initialData[i,:]-real_x[i,:])**2\n","        local_gradient=2*(real_x[i,:]-initialData[i,:])\n","        gradient[dim*i:(i+1)*dim]=local_gradient\n","    return res,gradient\n","\n","\n"],"execution_count":20,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"YZMRgVdEUlMM"},"source":["Un fermat que paralieliza, sino es una tortuga. Idem euclidea"]},{"cell_type":"code","metadata":{"collapsed":true,"id":"0lQ1vgttUlMS","executionInfo":{"status":"ok","timestamp":1623426379031,"user_tz":180,"elapsed":283,"user":{"displayName":"Nicolás Chehebar","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GjWN7wrOddIeQwbW6mzML-_Z0jTJLZ1cnE5tDuipw=s64","userId":"03611693109350494416"}}},"source":["           \n","#RETURNS: fermat cost and gradient of that cost\n","def fermat_cost_cuad_fast(x,initialData,dimension,fermat_distances,tree,data,k,alpha):\n","    assert(len(x)==(initialData.shape[0])*dimension)\n","    real_x=x.reshape(-1,dimension)     \n","    res=0\n","    gradient=np.zeros(x.shape)\n","    \n","    start_euclid_dist,start_ind=tree.query(initialData,k=1)\n","    idx=start_ind[:,0]\n","    start_euclid_dist=np.power(start_euclid_dist,alpha)\n","    \n","    end_euclid_dist,end_ind=tree.query(real_x,k=k)\n","    end_euclid_dist=np.power(end_euclid_dist,alpha)\n","    \n","    dist=([ [fermat_distances[idx[j],end_ind[j,i]]+end_euclid_dist[j,i] for i in range(end_euclid_dist.shape[1])] for j in range(end_euclid_dist.shape[0])])\n","    dist=np.array(dist)\n","    dist=np.min(dist,axis=1)\n","    \n","    minim=[[ fermat_distances[idx[j],end_ind[j,i]] for i in range(end_ind.shape[1])] for j in range(end_ind.shape[0])]  \n","    minim=np.argmin(minim,axis=1)\n","    \n","    grad_idx=[end_ind[j,minim[j]] for j in range(end_ind.shape[0])]\n","        \n","    distances=dist+start_euclid_dist[:,0]\n","    distances=distances.reshape(-1,1)\n","    \n","              \n","    q=data[grad_idx,:] \n","    aux=alpha*(np.linalg.norm(q-real_x,axis=1)**(alpha-2)).reshape(-1,1)\n","    local_grad=aux*(real_x-q)*2*distances\n","    gradient=local_grad.flatten()\n","    res=np.linalg.norm(distances)**2   \n","    return res,gradient\n","\n","def euclidean_cost_fast(x,initialData):\n","    assert(initialData.shape[0]*initialData.shape[1]==len(x))\n","    real_x=x.reshape(-1,initialData.shape[1])\n","    res=np.linalg.norm(np.linalg.norm(initialData-real_x,axis=1))**2\n","    gradient=2*(real_x-initialData).flatten()\n","    return res,gradient"],"execution_count":21,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"82EH0L-4UlMa"},"source":["Creo funciones del algoritmo propuesto por Tabak, Zhao, Triglia"]},{"cell_type":"code","metadata":{"collapsed":true,"id":"hGAd3m0uUlMh","executionInfo":{"status":"ok","timestamp":1623426379041,"user_tz":180,"elapsed":285,"user":{"displayName":"Nicolás Chehebar","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GjWN7wrOddIeQwbW6mzML-_Z0jTJLZ1cnE5tDuipw=s64","userId":"03611693109350494416"}}},"source":["def kernel(x,centro,b,dim):\n","  return np.exp(-(np.linalg.norm(x-centro)**2)/(2*b*b))/( (np.sqrt(2*np.pi)*b)*dim)\n","\n","#ARMO LA MATRIZ C    \n","\n","#vector con los valores z que hay, una matriz por filas paso\n","#viene en el mismo orden que los y\n","def armar_C(z,b=1):\n","    z_matrix=np.zeros((z.shape[0],z.shape[0]))\n","    dim=z.shape[1]\n","    for j in range(z.shape[0]):\n","        sumita=0\n","        for i in range (z.shape[0]):\n","            aux=kernel(z[j],z[i],b,dim)\n","            z_matrix[i,j]=aux\n","            sumita+=aux\n","        z_matrix[:,j]=z_matrix[:,j]/sumita\n","    C=z_matrix-np.mean(z_matrix,axis=1).reshape(-1,1) #CHECK, cambio axis 0 por 1\n","    return C\n","\n","\n","#VERSION DONDE Z ES CATEGORICO, SON ETIQUETAS, z es vector columna\n","def armar_C_categorico(z):\n","    z_matrix=np.zeros((z.shape[0],z.shape[0]))\n","    for j in range(z.shape[0]):\n","        sumita=0\n","        for i in range (z.shape[0]):\n","            aux=0\n","            if(z[i]==z[j]):\n","                aux=1\n","            z_matrix[i,j]=aux\n","            sumita+=aux\n","        z_matrix[:,j]=z_matrix[:,j]/sumita\n","    C=z_matrix-np.mean(z_matrix,axis=1).reshape(-1,1) #CHECK, cambio axias 0 por 1\n","    return C"],"execution_count":22,"outputs":[]},{"cell_type":"code","metadata":{"collapsed":true,"id":"8CEmegZfUlMo","executionInfo":{"status":"ok","timestamp":1623426379050,"user_tz":180,"elapsed":219,"user":{"displayName":"Nicolás Chehebar","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GjWN7wrOddIeQwbW6mzML-_Z0jTJLZ1cnE5tDuipw=s64","userId":"03611693109350494416"}}},"source":["#matriz que tieneen el lugar ij al kernel(y_i,y_j) \n","def kernel_matrix_y(y,h):\n","    pairwise_dists = squareform(pdist(y, 'euclidean'))\n","    return np.exp( -pairwise_dists**2/(2*(h**2))) #/ (np.sqrt( (2*np.pi)**dimension) )\n","    \n","    \n","#idem pero en el lugar ij esta el gradiente    \n","def grad_matrix_y(y,h):\n","    pairwise_dists=squareform(pdist(y,'euclidean'))\n","    #armo la matriz que tiene el escalar por el que multiplico a la direccion\n","    aux=np.exp( -pairwise_dists**2/(2*(h**2)))# / (np.sqrt( (2*np.pi)**dimension))\n","    aux=aux*(-1/(h**2))\n","    \n","    #simplemente hago en cada coordenada (De las d dimensiones) la resta de los y correspondeintes\n","    #en el lugar ij tendre y_i-y_j y eso luego lo multiplico por aux.\n","    unos=np.ones((y.shape[0],1))\n","    res=np.zeros((y.shape[0],y.shape[0],dimension))\n","    for d in range(dimension):\n","        res[:,:,d]=(y[:,d].reshape(-1,1)) @unos.T -  unos @ (y[:,d].reshape(-1,1)).T  \n","        res[:,:,d]=res[:,:,d]*aux\n","    return res\n","\n","\n","def kernel_matrix_y_centros(y,centros,h):\n","    pairwise_dists = cdist(y,centros, 'euclidean')\n","    return np.exp( -pairwise_dists**2/(2*(h**2)))# / (np.sqrt( (2*np.pi)**dimension) )"],"execution_count":23,"outputs":[]},{"cell_type":"code","metadata":{"collapsed":true,"id":"Ay-YCn-lUlMu","executionInfo":{"status":"ok","timestamp":1623426379064,"user_tz":180,"elapsed":222,"user":{"displayName":"Nicolás Chehebar","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GjWN7wrOddIeQwbW6mzML-_Z0jTJLZ1cnE5tDuipw=s64","userId":"03611693109350494416"}}},"source":["#AHORA DEFINO LAGRANGIANO Y SU GRADIENTE                \n","def lagrangiano(landa,kernel_matrix,f_matrix,C,cost,alpha): \n","    return  cost+ landa * ( np.sum(    f_matrix.shape[1]*np.diag(kernel_matrix @ C)/ (np.sum(f_matrix,axis=1)**((alpha-1)/alpha) )  ) )\n","    \n","#el gradiente de la parte de la F del lagrangiano\n","def fgrad(y,grad_matrix,C): \n","    dimension=y.shape[1]\n","    aux=np.zeros(y.shape).flatten()\n","    for i in range (y.shape[0]):\n","        local_grad=C[i,:]@grad_matrix[i,:,:]\n","        aux[i*dimension:(i+1)*dimension]= local_grad\n","    return aux\n","\n","def fgrad_fast(y,grad_matrix,kernel_matrix,fgrad_matrix,f_matrix,C,alpha): \n","    dimension=y.shape[1]\n","    #local_grad=np.zeros((y.shape[0],dimension))\n","    #for d in range (dimension):\n","    #   local_grad[:,d]=np.sum(C[:,:]*grad_matrix[:,:,d].T,axis=0) \n","    #NO ES UN AXIS1 ACA????? deberia sumar sobre las filas sin el transpuesto TODO lo arreglo en la version mas rapida\n","    #CHECK, TODO!!\n","\n","    arriba_deriv=np.sum((grad_matrix[:,:,:].T*C.T).T,axis=1)\n","    abajo_sin_deriv=np.sum(f_matrix,axis=1)**((alpha-1)/alpha)\n","\n","    arriba_sin_deriv=np.diag(kernel_matrix@C)\n","    auxi=((alpha-1)/alpha) * (np.sum(f_matrix,axis=1)**(-1/alpha))\n","    abajo_deriv= ( np.sum(fgrad_matrix,axis=1).T*auxi.T  ).T\n","    local_grad=(    ( (arriba_deriv.T*abajo_sin_deriv.T).T - (abajo_deriv.T * arriba_sin_deriv.T).T ).T  /  (abajo_sin_deriv**2).T ).T\n","    #print(np.sum(f_matrix,axis=1))\n","    return (f_matrix.shape[1])*local_grad.flatten()\n","\n","def fgrad_fast2(y,grad_matrix,C): \n","    dimension=y.shape[1]\n","    #local_grad=np.zeros((y.shape[0],dimension))\n","    #for d in range (dimension):\n","    #   local_grad[:,d]=np.sum(C[:,:]*grad_matrix[:,:,d].T,axis=0) \n","    #NO ES UN AXIS1 ACA????? deberia sumar sobre las filas sin el transpuesto TODO lo arreglo en la version mas rapida\n","    #CHECK, TODO!!\n","    local_grad=np.sum((grad_matrix.T*C.T).T,axis=1)\n","    return local_grad.flatten()\n","\n","\n","# derivo en y_i\n","def lagrangiano_grad(landa,cost_grad,f_grad):\n","    return cost_grad+landa*f_grad"],"execution_count":24,"outputs":[]},{"cell_type":"code","metadata":{"collapsed":true,"id":"eGNXwaJXUlM1","executionInfo":{"status":"ok","timestamp":1623426379074,"user_tz":180,"elapsed":225,"user":{"displayName":"Nicolás Chehebar","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GjWN7wrOddIeQwbW6mzML-_Z0jTJLZ1cnE5tDuipw=s64","userId":"03611693109350494416"}}},"source":["#La cuenta de como actualizar el landamin,hace la cuenta de que desciendo en direccion de L_F (que aca llamo f_grad)\n","def landamin_actualizacion(cost_grad,f_grad):\n","    num=cost_grad@f_grad\n","    denom=f_grad@f_grad\n","    return num/denom"],"execution_count":25,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"-ak0iJQ7UlM8"},"source":["Hago ahora una función que corra el algoritmo!"]},{"cell_type":"code","metadata":{"collapsed":true,"id":"WtajC6Y8UlNM","executionInfo":{"status":"ok","timestamp":1623426379083,"user_tz":180,"elapsed":226,"user":{"displayName":"Nicolás Chehebar","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GjWN7wrOddIeQwbW6mzML-_Z0jTJLZ1cnE5tDuipw=s64","userId":"03611693109350494416"}}},"source":["def barycenter(yinitial,n_iter,samples,landa0,landamax,learningrate0,z,h,beta_adapt,cost_function,landa_monotono=True,factor=1,categorical=True,b=1,callback=True):\n","    h_adapt=False\n","    if h==-1:\n","        h_adapt=True\n","    \n","    y=yinitial\n","    if categorical: \n","        C=armar_C_categorico(z)\n","    else: \n","        armar_C(z,b)\n","    landa=landa0\n","    learningrate=learningrate0\n","    n=0\n","\n","    while (n<n_iter):\n","        if h_adapt:\n","            h=np.std(y)*1.06*(y.shape[0]**(-1/5))*factor\n","            \n","        learningrate=min(learningrate*2.01,learningrate0)\n","\n","        cost,cost_grad=cost_function(y.flatten(),samples)\n","        grad_matrix=grad_matrix_y(y,h)\n","        f_grad=fgrad_fast2(y,grad_matrix,C)\n","        \n","        \n","        beta=landa*beta_adapt\n","        \n","        if landa_monotono:\n","            landamin=beta-landamin_actualizacion(cost_grad,f_grad)\n","            if(landamin>landa and landamin<landamax): \n","                landa=landamin\n","            elif(landamin>landamax): \n","                landa=landamax\n","        else:\n","            landa=beta-landamin_actualizacion(cost_grad,f_grad)\n","\n","        gradienteLagrangiano=lagrangiano_grad(landa,cost_grad,f_grad)\n","        yviejo=y\n","        costviejo=cost\n","\n","\n","        while(True):\n","            y=(yviejo.flatten()-(learningrate*gradienteLagrangiano)).reshape(yviejo.shape)  \n","            \n","            kernel_matrix_nuevo=kernel_matrix_y(y,h) \n","            cost,cost_grad=cost_function(y.flatten(),samples)\n","            kernel_matrix_viejo=kernel_matrix_y_centros(y,yviejo,h)\n","            \n","            lag_nuevo=lagrangiano(landa,kernel_matrix_nuevo,C,cost)\n","            lag_viejo=lagrangiano(landa,kernel_matrix_viejo,C,costviejo)\n","            \n","            if(lag_nuevo<=lag_viejo): \n","                break #es que el lagrangiano nuevo disminuyo efectivamente (uso el y nuevo en ambos -kernel_matrix- pero el costo que a cada uno le corresponde)\n","            learningrate=learningrate/2             \n","        \n","        paso=learningrate*np.linalg.norm(gradienteLagrangiano)\n","        if callback: print(\"it: \"+str(n)+\", paso: \"+str(paso)+\", landa: \"+str(landa)+\", h:\"+str(h))\n","        if paso<1e-7: landa=landa*1.1 \n","        n=n+1\n","    return y"],"execution_count":26,"outputs":[]},{"cell_type":"markdown","metadata":{"collapsed":true,"id":"tje_SA-7UlNS"},"source":["-El h adaptativo achica mucho mas rapido la varianza!! Pero favorece tener outliers, si se me juntan varios puntos y unos quedan muy lejos, el h será muy chicos y estos puntos alejados no van a tener en cuenta a los otros.\n","\n","-Los puntos se empiezana expandir para mejorar al principio, como despues pasa que los que se alejan hacia afuera se vienen al baricentro? Por qué pasaria? Ahi juega muy fuertemente la eleccion del h\n","\n","-Cuan adaptativo elijo el beta (alpha del paper)? Cuan dependiente del caso e esto? hacer un adaptativo on 0.01 o 0.05 parece andar bien. Todo indicaria que cuanto mas leve el aumento de landa mejor.\n","\n","-landa max??\n","\n","-landa0, es muy relevante? Arranca igual si tomo uno muy chico.\n"]},{"cell_type":"markdown","metadata":{"id":"XrqjhUG8UlNW"},"source":["### A implementar"]},{"cell_type":"markdown","metadata":{"id":"-AQFFiOqUlNY"},"source":["Ir con un bandwith muy grande e ir achicandolo. NO pasarse de la rule of thumb, no pasarse de la del principio o ir actualizando!!No tenes que irse con el bandwidth muy chico, sensible a outliers. Recocido simulado de agrandar el h de la nada por si se te escapo algo?\n","\n","Actualizacion del lambda? Tiene que ir creciendo necesariamente?\n","\n","Dejar kernels gaussianos pero evaluarlos en la distancia de fermat! EL gradietne hay que actualizarlo tambien con el gradiente de fermat. ESo haria que los puntos se muevan siguiendo la geodesica (pesada por la densidad)\n","\n","Rellenar el espacio con poca densidad de fermat y que gradualmente vaya aumentando hacia la misma sueprficie"]},{"cell_type":"markdown","metadata":{"id":"540iptbNUlNc"},"source":["Ver el tema del h\n","\n","Ejemplo no singular de densidad\n","\n","enchufar fermat en donde estan los kernels.\n","\n","Ejemplos distintos: una curva menos pronunciada, otro ejemplo tipo el circulo? Una curva mas separada para afuera"]},{"cell_type":"markdown","metadata":{"collapsed":true,"id":"4NdDDTiHUlNf"},"source":["# Defino Fermat Kernels y la funcion de baricentro con ellos"]},{"cell_type":"code","metadata":{"collapsed":true,"id":"qe-EglYoUlNi","executionInfo":{"status":"ok","timestamp":1623426379093,"user_tz":180,"elapsed":227,"user":{"displayName":"Nicolás Chehebar","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GjWN7wrOddIeQwbW6mzML-_Z0jTJLZ1cnE5tDuipw=s64","userId":"03611693109350494416"}}},"source":["           \n","#RETURNS: fermat cost and gradient of that cost\n","def pairwise_fermat_dist(real_x,initialData,fermat_distances,tree,data,k,alpha,C):\n","    iguales=False if real_x.shape!=initialData.shape  else np.alltrue(np.abs(real_x-initialData)<1e-10)\n","    dimension=real_x.shape[1]\n","    originalshape0=real_x.shape[0]\n","    originalshape1=initialData.shape[0]\n","    \n","    res=0   \n","    \n","    start_euclid_dist,start_ind=tree.query(initialData,k=1)\n","    start_euclid_dist=np.power(start_euclid_dist,alpha)\n","\n","    start_euclid_dist=np.tile(start_euclid_dist,(originalshape0,1))\n","    start_ind=np.tile(start_ind,(originalshape0,1))\n","    idx=start_ind[:,0]\n","    \n","    end_euclid_dist,end_ind=tree.query(real_x,k=k)\n","    end_euclid_dist=np.power(end_euclid_dist,alpha)\n","    end_euclid_dist=np.repeat(end_euclid_dist,originalshape1,axis=0)\n","    end_ind=np.repeat(end_ind,originalshape1,axis=0)\n","    \n","    real_x=np.repeat(real_x,originalshape1,axis=0) #me queda el punto 1 n veces, el 2 otras n veces y asi\n","    initialData=np.tile(initialData,(originalshape0,1)) #\n","    \n","    '''\n","    dist=([ [fermat_distances[idx[j],end_ind[j,i]]+end_euclid_dist[j,i] for i in range(end_euclid_dist.shape[1])] for j in range(end_euclid_dist.shape[0])])\n","    dist=np.array(dist)\n","    dist=np.min(dist,axis=1)\n","    \n","    minim=[[ fermat_distances[idx[j],end_ind[j,i]] for i in range(end_ind.shape[1])] for j in range(end_ind.shape[0])]  \n","    minim=np.argmin(minim,axis=1)\n","    '''\n","    \n","    \n","    dist=np.array([ fermat_distances[idx[:],end_ind[:,i]] for i in range(end_euclid_dist.shape[1])]).T\n","    dist=dist+end_euclid_dist\n","    dist=np.array(dist)\n","    dist=np.min(dist,axis=1)\n","    \n","    minim2=np.array([ fermat_distances[idx[:],end_ind[:,i]] for i in range(end_ind.shape[1])]).T \n","    minim=np.argmin(minim2,axis=1)\n","    #maxim=np.argmax(minim2,axis=1)\n","\n","    grad_idx=[end_ind[j,minim[j]] for j in range(end_ind.shape[0])]\n","    #grad_idx2=[end_ind[j,maxim[j]] for j in range(end_ind.shape[0])]\n","        \n","    distances=dist+start_euclid_dist[:,0]\n","    distances=distances.reshape(originalshape1,originalshape0).T #traspongo porque calcule al reves\n","\n","    if iguales:\n","      np.fill_diagonal(distances,0)#cada uno con sigo mismo esta a distancia 0\n","      #FUNDAMENTAL ESTO PORQUE SINO BUSCO ALEJARME MUCHO DE TODO Y EL COSTO VA A DAR 0\n","              \n","    q=data[grad_idx,:] \n","    #q2=data[grad_idx2,:]\n","    \n","    aux=alpha*(np.linalg.norm(q-real_x,axis=1)**(alpha-2)).reshape(-1,1)\n","    #aux2=alpha*(np.linalg.norm(q2-real_x,axis=1)**(alpha-2)).reshape(-1,1)\n","    local_grad=aux*(real_x-q) #apunto en la direccion opuesta al q\n","    #local_grad2=-aux2*(real_x-q2) #apunto en la direccion del q2\n","    grad=local_grad.reshape(originalshape0,originalshape1,dimension)\n","    if iguales:\n","      np.einsum('iij->ij',grad)[...] = 0   #ME REEMPLAZA EL FOR!\n","\n","    #grad2=local_grad2.reshape(originalshape0,originalshape1,dimension)\n","    #np.einsum('iij->ij',grad2)[...]\n","    #grad=np.zeros((originalshape0,originalshape0,dimension))\n","    #for d in range(dimension):\n","    #    grad[:,:,d]=local_grad[:,d].reshape(originalshape0,originalshape0)\n","    #    np.fill_diagonal(grad[:,:,d],0)\n","    #print(np.alltrue(np.abs(grad2-grad)<1e-12))\n","    #gradposta=np.where(C<0,grad,grad2)\n","    \n","    return distances,grad #PARA INCLUIR GRADIENTE PONER ACA GRADPOSTA"],"execution_count":27,"outputs":[]},{"cell_type":"code","metadata":{"collapsed":true,"id":"frOzDt6wUlNn","executionInfo":{"status":"ok","timestamp":1623426379101,"user_tz":180,"elapsed":222,"user":{"displayName":"Nicolás Chehebar","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GjWN7wrOddIeQwbW6mzML-_Z0jTJLZ1cnE5tDuipw=s64","userId":"03611693109350494416"}}},"source":["#matriz que tieneen el lugar ij al kernel(y_i,y_j) \n","def kernel_matrix_y_centros_fermat(y,centros,h,fermat_distances,tree,data,k,alpha):\n","    dimension=2\n","    assert(y.shape==centros.shape)\n","    Cfalsa=np.zeros((y.shape[0],y.shape[0],y.shape[1])) #basicamente no importa el gradiente aca. \n","    pairwise_dists,trash = pairwise_fermat_dist(y,centros,fermat_distances,tree,data,k,alpha,Cfalsa)\n","    return np.exp( -pairwise_dists**2/(2*(h**2)))# / ((np.sqrt( 2*np.pi)*h)**dimension)\n","\n","\n","def kernel_matrix_y_fermat(y,h,fermat_distances,tree,data,k,alpha,PD):\n","    dimension=2\n","    return np.exp( -PD**2/(2*(h**2)))# / ((np.sqrt( 2*np.pi)*h)**dimension)\n","    \n","    \n","#idem pero en el lugar ij esta el gradiente    \n","def grad_matrix_y_fermat(y,h,fermat_distances,tree,data,k,alpha,PD,G_PD):\n","    dimension=2\n","    pairwise_dists=PD\n","    grad = G_PD\n","    aux=np.exp( -pairwise_dists**2/(2*(h**2)))# / ((np.sqrt( 2*np.pi)*h)**dimension)\n","    aux=aux*(-pairwise_dists/(h**2))\n","    \n","    #simplemente hago en cada coordenada (De las d dimensiones) multiplicar por lo que seria el gradiente\n","    #de pairwise_dists que es lo que falta para completar regla de la cadena.\n","\n","    #res=np.zeros(grad.shape)\n","    res2=(grad.T*aux.T).T    #ESTO ME SACA EL FOR!\n","    #for d in range(dimension):\n","    #    res[:,:,d]=grad[:,:,d]*aux\n","    #print(np.alltrue(np.abs(res2-res)<1e-12))\n","    return res2\n"],"execution_count":28,"outputs":[]},{"cell_type":"code","metadata":{"collapsed":true,"id":"aujU-5IvUlNr","executionInfo":{"status":"ok","timestamp":1623426379107,"user_tz":180,"elapsed":217,"user":{"displayName":"Nicolás Chehebar","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GjWN7wrOddIeQwbW6mzML-_Z0jTJLZ1cnE5tDuipw=s64","userId":"03611693109350494416"}}},"source":["def barycenter_fermat(yinitial,n_iter,samples,landa0,landamax,learningrate0,z,h,beta_adapt,cost_function,f_dist,tree,data,k,alpha,quant=0.95,quant2=0,tasaenfrio=0.8,factor=1,categorical=True,b=1,callback=True,landa_monotono=True):\n","    \n","    \n","    if categorical: \n","        C=armar_C_categorico(z)\n","    else: \n","        C=armar_C(z,b)\n","    y=yinitial\n","    h_adapt=False\n","\n","    Crep=np.repeat(C[:, :, np.newaxis], yinitial.shape[1], axis=2)\n","\n","    PD,G_PD=pairwise_fermat_dist(y,y,f_dist,tree,data,k,alpha,Crep)\n","\n","\n","    if h==-1:\n","        h_adapt=True\n","        h=np.quantile(PD.flatten(),quant)*factor\n","    control=np.quantile(PD.flatten(),quant2)    \n","    \n","    kernel_matrix_nuevo=kernel_matrix_y_fermat(y,h,f_dist,tree,data,k,alpha,PD)\n","\n","    fPD,fG_PD=pairwise_fermat_dist(y,data,f_dist,tree,data,k,alpha,Crep)\n","    f_matrix_nuevo=kernel_matrix_y_fermat(y,h,f_dist,tree,data,k,alpha,fPD)\n","\n","\n","    landa=landa0\n","    learningrate=learningrate0\n","    n=0\n","\n","    while (n<n_iter):\n","            \n","        learningrate=min(learningrate*2.01,learningrate0)\n","\n","        cost,cost_grad=cost_function(y.flatten(),samples)\n","                \n","\n","        fgrad_matrix=grad_matrix_y_fermat(y,h,f_dist,tree,data,k,alpha,fPD,fG_PD)\n","        grad_matrix=grad_matrix_y_fermat(y,h,f_dist,tree,data,k,alpha,PD,G_PD) \n","        f_grad=fgrad_fast(y,grad_matrix,kernel_matrix_nuevo,fgrad_matrix,f_matrix_nuevo,C,alpha)  #ACA HAY QUE CAMBIAR y dividir por la de costos de fermat lugar a lugar\n","        \n","        '''\n","        beta=landa*beta_adapt\n","        \n","        if landa_monotono:\n","            landamin=beta-landamin_actualizacion(cost_grad,f_grad)\n","            if(landamin>landa and landamin<landamax): \n","                landa=landamin\n","            elif(landamin>landamax): \n","                landa=landamax\n","        else:\n","            landa=beta-landamin_actualizacion(cost_grad,f_grad)\n","        '''\n","        gradienteLagrangiano=lagrangiano_grad(landa,cost_grad,f_grad)\n","        yviejo=y\n","        costviejo=cost\n","\n","\n","        while(True):\n","            y=(yviejo.flatten()-(learningrate*gradienteLagrangiano)).reshape(yviejo.shape)  \n","            \n","            PD,G_PD=pairwise_fermat_dist(y,y,f_dist,tree,data,k,alpha,Crep)\n","            kernel_matrix_nuevo=kernel_matrix_y_fermat(y,h,f_dist,tree,data,k,alpha,PD)\n","            fPD,fG_PD=pairwise_fermat_dist(y,data,f_dist,tree,data,k,alpha,Crep) \n","            f_matrix_nuevo=kernel_matrix_y_fermat(y,h,f_dist,tree,data,k,alpha,fPD)\n","            \n","            cost,cost_grad=cost_function(y.flatten(),samples)\n","            kernel_matrix_viejo=kernel_matrix_y_centros_fermat(y,yviejo,h,f_dist,tree,data,k,alpha)\n","            fPD_viejo,fG_PD_viejo=pairwise_fermat_dist(y,data,f_dist,tree,data,k,alpha,Crep)\n","            f_matrix_viejo=kernel_matrix_y_fermat(y,h,f_dist,tree,data,k,alpha,fPD_viejo)\n","            \n","            lag_nuevo=lagrangiano(landa,kernel_matrix_nuevo,f_matrix_nuevo,C,cost,alpha) #ACA HAY QUE CAMBIAR y dividir por la de costos de fermat lugar a lugar\n","            lag_viejo=lagrangiano(landa,kernel_matrix_viejo,f_matrix_viejo,C,costviejo,alpha) #ACA HAY QUE CAMBIAR y dividir por la de costos de fermat lugar a lugar\n","            control2=np.quantile(PD.flatten(),quant2)\n","            \n","            if (lag_nuevo<=lag_viejo):      \n","                control=control2\n","                if h_adapt:\n","                  h=np.quantile(PD.flatten(),quant)*factor \n","                break #es que el lagrangiano nuevo disminuyo efectivamente (uso el y nuevo en ambos -kernel_matrix- pero el costo que a cada uno le corresponde)\n","\n","            learningrate=learningrate/2             \n","        \n","        paso=learningrate*np.linalg.norm(gradienteLagrangiano)\n","        if callback: print(\"it: \"+str(n)+\", paso: \"+str(paso)+\", landa: \"+str(landa)+\", h:\"+str(h))\n","            \n","             \n","        \n","        #PARA NO QUEDARME ATORADO!\n","        if paso<1e-6: \n","            if landa>=landamax:\n","              return y\n","            landa=landa*beta_adapt\n","            learningrate=learningrate0\n","            quant=quant2+((quant-quant2)*tasaenfrio)\n","            '''\n","            bari2=y\n","            fig = plt.figure()\n","\n","            ax = fig.add_subplot(111, projection='3d')\n","            ax.view_init(20, 75)\n","            ax.scatter(xs=data[:,0], ys=data[:,1], zs=data[:,2],c='b', s=0.01)\n","            ax.scatter(xs=samples[:,0], ys=samples[:,1], zs=samples[:,2],c='r', s=1)\n","            ax.scatter(xs=bari2[:n_samples//2,0], ys=bari2[:n_samples//2,1], zs=bari2[:n_samples//2,2], s=10,c='g')\n","            ax.scatter(xs=bari2[n_samples//2:,0], ys=bari2[n_samples//2:,1], zs=bari2[n_samples//2:,2], s=10,c='y')\n","\n","            plt.title('landa_actual:'+str(landa)+\", bbeta_adapt:\"+str(beta_adapt)+\", quant:\"+str(quant)+\", quant2:\"+str(quant2)+\", h_actual:\"+str(h)+\",control:\"+str(control))\n","\n","            plt.show()\n","            '''\n","            '''\n","            if landa>=landamax: #ya termine, atorado y llegue a landamax\n","                return y\n","            else:\n","                landa=landa*1.1 \n","                #learningrate=learningrate0*landa0/landa\n","                learningrate=learningrate*10\n","            '''    \n","              \n","        n=n+1   \n","        \n","        '''\n","        bari=y\n","        plt.plot(W[:,0],W[:,1],'.',c='y')\n","        plt.plot(x1,x2,'.',c='r')\n","        plt.plot(y1,y2,'.',c='b')\n","        plt.plot(bari[:puntos,0],bari[:puntos,1],'x',c='r')\n","        plt.plot(bari[puntos:,0],bari[puntos:,1],'x',c='b')\n","        plt.axis('equal')\n","        plt.show()\n","        \n","        bari3=y\n","        \n","        plt.plot(data3[:densidad,0],data3[:densidad,1],'.',c='y')\n","        plt.plot(data3[densidad:densidad+puntos1,0],data3[densidad:densidad+puntos1,1],'.',c='r')\n","        plt.plot(data3[densidad+puntos1:,0],data3[densidad+puntos1:,1],'.',c='b')\n","        plt.plot(bari3[:puntos1,0],bari3[:puntos1,1],'x',c='r')\n","        plt.plot(bari3[puntos1:,0],bari3[puntos1:,1],'x',c='b')\n","        plt.axis('equal')\n","        plt.title('landa0: '+str(g)+\", beta_adapt: \"+str(gg))\n","        plt.show()\n","\n","        bari2=y\n","        plt.plot(data2[:densidad,0],data2[:densidad,1],'.',c='y')\n","        plt.plot(data2[densidad:densidad+puntos1,0],data2[densidad:densidad+puntos1,1],'.',c='r')\n","        plt.plot(data2[densidad+puntos1:,0],data2[densidad+puntos1:,1],'.',c='b')\n","        plt.plot(bari2[:puntos1,0],bari2[:puntos1,1],'x',c='r')\n","        plt.plot(bari2[puntos1:,0],bari2[puntos1:,1],'x',c='b')\n","        plt.axis('equal')\n","        plt.title('landa0: '+str(g)+\", beta_adapt: \"+str(gg))\n","        plt.show()\n","        '''\n","    return y"],"execution_count":29,"outputs":[]},{"cell_type":"code","metadata":{"collapsed":true,"id":"2Y311JMVUlNy","executionInfo":{"status":"ok","timestamp":1623426379114,"user_tz":180,"elapsed":216,"user":{"displayName":"Nicolás Chehebar","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GjWN7wrOddIeQwbW6mzML-_Z0jTJLZ1cnE5tDuipw=s64","userId":"03611693109350494416"}}},"source":[""],"execution_count":29,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"hbQb89YcONWQ"},"source":["# Ejemplo z continuo"]},{"cell_type":"code","metadata":{"id":"lbxgfY7uOQHh","colab":{"base_uri":"https://localhost:8080/","height":281},"executionInfo":{"status":"ok","timestamp":1623426379853,"user_tz":180,"elapsed":942,"user":{"displayName":"Nicolás Chehebar","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GjWN7wrOddIeQwbW6mzML-_Z0jTJLZ1cnE5tDuipw=s64","userId":"03611693109350494416"}},"outputId":"93f5f531-51f7-45eb-e9a8-da7102e12ea3"},"source":["n=100\n","eps=0.3\n","#creo que es que los kernels se hacen tremendamente chicos, habria que normalizar por la dimension de la variedad en la que trabajo, lo que ni idea\n","#mas facil volar esa normalizacion como volumen del abola d dimensional y ya\n","\n","data1=np.array([np.linspace(-1,1,n),np.linspace(-1,1,n)**3]).T+np.random.random((n,2))*eps\n","data2=np.array([np.linspace(-1,1,n),np.linspace(-1,1,n)**3]).T+np.random.random((n,2))*eps\n","data=np.concatenate((data1,data2),axis=0)\n","\n","n_samples=50\n","\n","z2=np.random.normal(loc=-0.5,scale=0.25,size=n_samples).reshape(-1,1)\n","z2=np.concatenate((z2,np.random.normal(loc=0.7,scale=0.1,size=n_samples).reshape(-1,1)),axis=0)\n","zaux=np.concatenate((z2,z2),axis=1)\n","zaux[:,1]=zaux[:,1]**3\n","\n","\n","samples=zaux\n","aux=np.random.random(samples.shape[0])*0.4-0.2\n","samples[:,1]=samples[:,1]+aux\n","\n","data=np.concatenate((data,samples),axis=0)\n","\n","plt.plot(data[:,0],data[:,1],'.',c='y')\n","plt.plot(samples[n_samples:,0],samples[n_samples:,1],'.',c='r')\n","plt.plot(samples[:n_samples,0],samples[:n_samples,1],'.',c='b')\n","plt.axis('equal')\n","plt.title('Rojo es una distribucion, azul otra y amarillo es la superficie')\n","plt.show()"],"execution_count":30,"outputs":[{"output_type":"display_data","data":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAYsAAAEICAYAAACuxNj9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de3hU1bn48e87yQQ0QIKBElC5KMFCmhIspmJBaXVU2qJUW+uprbVSsedoj5zSQ6ttLdW2tvaHBau9ULWKYi1VyyFWqhEbLzUWUUIjUCEoyC2IgQxEhUwy6/fH3jPZM9lzI2Euyft5Hh4yM3v2XrNnZr+z1rsuYoxBKaWUiseT6QIopZTKfhoslFJKJaTBQimlVEIaLJRSSiWkwUIppVRCGiyUUkollPXBQkQ2iMj0TJcjE0Rkm4icZ/99k4jc04P7bhWRU+y/7xeRH/fUvh3HyNn37lidk95ORKaJyBuO287P8AIReSiDZbtKRF7M1PF7gogcJyLVIuIXkT+LyBUi8nQSz+v29SO/O09OlohsA4YBHUAr8DfgemNMa6LnGmPKj23pcoMx5qfJbCcitcBDxpi4HwxjzICeKFeCY/S5905ERgNvAV5jTHtmS5N+xpgXgNMyXY5e7PNY19ISx+drWaInJXv9iCedNYuZ9gWqEpgE3JjGYyubiKTlB4KKrbe+B731dWULEckDRgGbM/FDJO3NUMaYJuAprKABgIhcZDdZtIhIrYiMdzzmrMb2E5FFIrLb/rdIRPrFOpaIXC0im0TkgIg8JSKj7PtFRH4pIu+IyEERaRCRj8TYR/j49u1wVVpERouIEZGvisjbIvKuiHzPsW2ViNTZr2uPiNwlIgVxyvsVEdkuIs3O/bgct7+IPGRv1yIir4jIMBH5CTANuMtuZrrL3t6IyHUisgXY4rhvrOMQQ0SkRkQOichzjnMVeo35jrLUisjXHbevsc/zIRHZKCKnp/Leich0EdkpIvPs92SPiHwt1nlyOW+LRWSH/V6+KiLTHI+12OeiVUTes1/LaLcmCZdzEut4HhH5vv1evSMiS0WkyH74efv/0HGn2Mf6h/2ZawYWiMipIvKs/R6+KyLLRKQ4xvHuFpGFUfetFJH/OYrzsUCs5ouH7PerQUTGiciN9mvZISLnO7b/muO9fVNErnU8FnrfviMiTcAfQvclOof282N+7122/bD9+dwvIm+IyGWOxz5tf+4OicguEfl2ksePeZ5ctnU9RqLPkVjNmb91+24l8bruF5HfiMiTIvIe1mfrZuCL9mdrdvTxRaTcsb+9InKTfX9EE6CInCkiL9nnfr0k01xsjDnm/4BtwHn23ycBDcBi+/Y44D3AB3iB+UAjUODy3FuAl4EPAUOBl4BbYxzzYns/47Ga274PvGQ/dgHwKlAMiL3N8ERlt28vwGrmARgNGOD3wHHAROAIMN5+/GPAmfbxRwObgLkxjjMBq4nubKAfcAfQ7njtzuNeC1QDxwN59nEG2Y/VAl+P2rcBaoATgOMc9421/74fOOQ49mLgxajXmO/YX/gYwBeAXcAZ9rkcC4xK5b0Dptuv9Rb7M/Bp4H1gcJKfry8DJfZ5ngc0Af1dtvsp1hfOC1wVeo1R58l5Tn4c43hXY322TgEGAI8DD8Y5X1fZr++bdhmPs8+Tzz7fQ+1yLYpxvCpgN+Cxbw+xz8+wVM+H/Tk6jPUdyAeWYjWbfc8+L9cAbzn29RngVPu9Pcc+7ulR79vP7ddxnH3fzhjf/QV0fobjfu+jXk8hsAP4ml3mScC7wAT78T3ANPvvwaHyuewn4j1P9nMT7xhJfo5ifbcSva77AT/wCawf9v2d5zD6+MBAu5zz7G0HAh93OfcnAs1Y3zOP/R40A0Pjfs+S+TJ295/9gWm1T5oBVgPF9mM/AJY7tvVgXXymu3zYtgKfdmx7AbAtxjFXAbOj9vs+VjXuU8BmrAu5J4myJwoWJzkeXwNcHmNfc4G/xHjsZuCRqC9IG+5ftKuxLrYfddlPLe7B4lMJPtDOYw/Ayi+dTOJg8RRwQ6JzF++9w7rAfBB1jHeAM4/y83YAmBh13xft8gxN4UseK1isBv7Lcfs0IEDnjwK3YPF2gjLPAtbFeXwT4LP/vh548mjOh/05qnE8NhPru5ln3x5ol784xr5WhN5v+31rw3GBJflgEfd77/LevRB13++AH9p/v431A2pQgvPQ5T1P9LlxPOZ6jCQ/R7G+W4le1/3A0qjHw+cw+vjAf8T6DEWd++9g/7hxPP4U8NV45y+dzVCzjDEDsT5MH8b6dQQwAtge2sgYE8SKtie67CNiW/vvETGONwpYbFezWoD9WL+OTjTGPAvcBdwNvCMiS0Rk0NG+MKxfJCHvY30gsKv3T4hIk4gcxPplO8RtB/br2BG6YYx5Dyvau3kQ6819RKwmndtFxJugjDuSfdxYHQ/2E/vcOp2MFQgSSfTeNZvIdtjweUxERL5tN5X47fe6CMd5FpFJWO/354wx+5LZZwJuryUfK/EYS8T5F6vZ8BG7SeMg8BCxPxsAD2D9Esb+/8FYGyY6H8Bex98fAO8aYzoct6HzMzxDRF62mzVasH6NOve1zxhzOE65Y0nlez8K+Hjou2yX4wqg1H78Urtc2+1mninJFCCJ8+R0VMewxfpuJXpdEc9NQrLfxVHAF6KOOxUYHu9JmchZPIcVMf+ffddurMIDVj4B60Xvcnl6xLbASPs+NzuAa40xxY5/xxljXrLLcacx5mNYzT/jgP+NsZ/3sJp7QkpjbOfmN8C/gTJjzCDgJqyA5WYP1usGQESOx6oid2GMCRhjfmSMmQCcBXwWuDL0cIz9x7o/xHnsAVhNVruxXj/EPgc7sJopEknlvUua3c48H7gMq9mqGKvqLvbjH8L6NXydMWad46kR76uIpPK+ur2WdqyLcLLn/6f2fRX2Z+PLxP5sgBVMLhaRiVjNpivcNkp0PlIhVk7pMazv6jB7X09G7SvR5yqWVL73O4Dnor7LA4wx/wlgjHnFGHMxVhPnCmB5ooOnep7iHCOZz1Gs71bc1xU6dKLX4rADq2k0me0ejDpuoTHmZ/GelKlxFosAn/3BXw58RkTOtX8dz8Nq93/J5Xl/BL4vIkNFZAhW002sftu/BW4UkXIAESkSkS/Yf58hIh+3j/ceVhtuMMZ+6oHLRcQrIpOxuq4layBwEGgVkQ8D/xln20eBz4rIVLGS4LcQ4/0RkU+KSIVYvSMOYjWBhMq/l+Q+MNE+7Tj2rcDLxpgd9i/xXcCXRSRPRK4mMjjcA3xbRD4mlrHOBJ5DKu9d9Ou9Sqzu124GYl2o9wH5InIzMMh+Xj7WeX3IGBN9AVkPlItIpYiE2oKT9Ufgf0RkjP3l/ynwJ7tmtA/rvUj0HgzEav7xi8iJxP6xAoAxZifwClaN4jFjzAcxNo15Po5CAVY7+z6gXURmAOfHf0rSUvnePwGME6sDiNf+d4aIjBeRArHGGhQZYwJY34dY32WnpM9TgmMk8zly/W7Fe11JlN/NE8BwEZkrVoeSgSLycZftHgJmisgF9ne6v1gdE06Kt/OMBAv7ArQUuNkY8wbWr6pfYSV3ZmJ1s21zeeqPgbXAv7CS5K/Z97kd4y9YibdH7Gr+68AM++FBWEnpA1hV4WbgFzGK+wOsi+MB4EfAwym81G8DX8LK1fwe+FOsDY0xG4Dr7P3vsY8Xq1dJKdZF8CBWW/ZzdDZLLAY+L1YPsDtTKOvDwA+xqsgfo7PJA6yk5/9inadyHF9oY8yfgZ/Yzz+E9avrBJf9J/3euTgZ+EeMx57CGrezGeu9PExn1f0krN5hc6WzR1SriIw0xmzGCsjPYPUQS2Ww1n1Y5/t5rOTwYazkNcaY97HOxz/sKv6ZMfbxI+B0rF+zf8VKkifyAFBBnCYo4p+PlBhjDgH/jXVhP4D1WV55NPty2XfS33u7HOcDl2P9Im+iM6kO8BVgm/09/wZWU04iqZ4n12Mk+Tly/W4l8bpSYu/Ph3Uum+zyfNJlux1YHYBuwgqWO7C+33HjgdjJjawlIm8DXzbGPJ9wY9UriTVC9QZjzKZMlyWTRORsrF+Fo0y2f3EVYHV/xUr4fz/TZemurB5EIyJDsboVbstwUVQGGWN6qukjZ9lNNTcA92igUJmQtXNDicgZWNWoXxlj3s50eZTKFLsNuwWrt8qiDBdH9VFZ3wyllFIq87K2ZqGUUip7ZG3OYsiQIWb06NGZLoZSSuWUV1999V1jzNCe3m/WBovRo0ezdu3aTBdDKaVyiohsT7xV6rQZSimlVEIaLJRSSiWkwUIppVRCGiyUUkolpMFCKaVUQhoslFJKJaTBQimVUX5/Hdu334bfX5fpoqg4snachVKq9/P761i//lyCwTY8ngImTlxNUVEqi9CpdNGahVIqY1paagkG24AOgsE2WlpqM10kFYMGC6VUxhQXT8fjKQDy8HgKKC6enukiqRi0GUoplTFFRVOYOHE1LS21FBdP1yaoLKbBQimVUUVFUzRI5ABthlJKpZ32gMo9WrNQSqWV9oDKTVqzUEqllfaAyk0aLJRSaaU9oHKTNkMppdJKe0DlJg0WSqm00x5QuUeboZRSSiWkwUIplRG7dy9h/foL2L17SaaLopKgzVBKqbTbvXsJmzdfC8CBA08DMGLEnEwWSSWgNQulVNrt2/dY3NuJ6KC+9NNgoZRKuwEDKiNuezzHJ33hDw3qe+utH7B+/bkaMNJEg4VSKq38/jp27foVIPa/PJqbq5O+8OugvszokWAhIveJyDsi8nqMx0VE7hSRRhH5l4ic3hPHVUplp3jNRJ0Xe2PfEySVC78O6suMnkpw3w/cBSyN8fgMoMz+93HgN/b/SqleJtHcT6GLfTDYhkg+YDCmI+kLvw7qy4weCRbGmOdFZHScTS4GlhpjDPCyiBSLyHBjzJ6eOL5SKnu4NRM5L+jRF/vQc1K58OugvvRLV9fZE4Edjts77fsigoWIzAHmAIwcOTJNRVNK9SRnzSFWbcF5sdcEdW7IqnEWxpglwBKAyZMnmwSbK6WyUCrNRDpdee5IV7DYBZzsuH2SfZ9SqhdKtpkossnqME1NSzVYZKl0dZ1dCVxp94o6E/BrvkIpVVw8HZE8+5ahqekP2iyVpXqq6+wfgTrgNBHZKSKzReQbIvINe5MngTeBRuD3wH/1xHGVUrmtqGgKpaVXY423AGPaddxEluqp3lD/keBxA1zXE8dSSmWn1polBGoew+u7lAG+5Od5Ki29kr17H4ibEFeZl1UJbqVUbmqtWcJxM6+lMADBO5+mtZqkA4aOm8gNGiyUUt0WqHmMwgBIECRg3SaF2oWOm8h+OjeUUqrbvL5LCXoh6AHjtW4nQ2ePzR1as1BKddsA3xxaq0kpZ6FjLHKLBgulVI8Y4JuTUtNTomlBVHbRZiilVEbo7LG5RYOFUiplrTVLODD/AlprlsS9L55QL6gxY26NaILSPEZ20mYopVRK3LrJAkfVdTa6F5TmMbKX1iyUUikJ1DyGJ6qbrNt98cSqPegqeNlLg4VSKiVu3WRT6Tobbw1tzWNkL22GUkqlJFY32WS6zvr9dWzbtoBg8AgQ7NILSkdzZy+xpm3KPpMnTzZr167NdDGUUj2kMx9hBQrw4PH007xEDxORV40xk3t6v9oMpZRKi858hBUoBg8+TwNFDtFgoZRKi8h8RD9Gj16ggSKHaM5CKZWUVKYg9/vruuQdNB+R2zRYKKUSijUFuVtQiDdWQmeXzV0aLJRSCblNQe6vqnANCj0155NbIFKZozkLpVRCbuMoYg2gKy6eTtFGYeQyKNooeL0lKU/fEW8shsoMrVkopRJyG1vR4a/D4ynoshxq3poGPvqtdjwBCHrb+RfX4Z9gUpq+Q2ekzT4aLJRSSYmegjxWwjq6yapoXTv+CaR00Q/1nNJ1ubOHBgul1FFzS1i3T60keOfTSMBqsmqpzANI6aKvPaeyjwYLpZSrVLrKOrd9f2IxOxYKxfWGlkqh8LxrKDyK42vPqeyiwUIp1UW4q2wbmDue5sBPn2fw/Ifib2t3qz1++Xy2V/TnYLnVhFQ6cBKNjXMJBtvYu/cBJk5cDaC1hhyjwUIp1UWg5jEK20AM0AFFNy2jddLZrjWMwL2Lw9tKAPJfrGfi9zqbkKKT1U1NS9m79wFdsyLHaNdZpVQXXt+lGA8YQLACQeDexV1Wwtu38jsMenwjGGtb8qznFhVNYdSoGykqmtJl2nFA16zIQVqzUEq5OnROKQNrmwAweTDo8Y14OjaGR3B3VFVw6In/x5AOK6AEBfZcCAOqKiL2E52sBiJqFtrTKTdosFBKRQjlIDwBCOaB/5IJgKH4z5siRnAfHNdMy0RD0Eu459PeC4R2l+6x0cnqsWMXsW/fYwwdeqk2QeUIDRZKqQgR4yQARp5kjeBecW04KHh9l1JcXMH2iv6sX/gBxfXQUgmtFf05tXh63Kk6/P66cMLb73+BwsKKiHmlNPGdnTRYKKUieH2XRoyTCHWddVsJb+LE1bSMqsX7uRIKD62jEHjvvQY7GBxBJI+ysrsYMaIzMR5rdHa8CQhV5mmwUEpFiBUYQiO4W2uWcGD+BVYi2zcnfKEP1RaamgRjOgCDMUG2bLkOgECgmeLi6Y7R2UcQseaOAp3iI9tpsFBKdRE9tUdIrKnKnRd6YzxYDVjWks3GdLBly3UY0zk/1Nixi9iy5XqM6aCxcS6FhRU6xUeW066zSqmkBWoewxM1VTl0XQXv5JO/jUg+4EEkH2OCOGsMgUCzfV8wohYxceJqxoy5VZugspDWLJRSMYWm8WifWsn7E4s5PmreJ6/vUsB9LqchQ2bR0lKL11sSbqJy1hjcahE6xUf2EmNMpsvgavLkyWbt2rWZLoZSfYpzjiegswutF9YvFFor+jO+5Zvkv1if1JxRIbFW1NOeTz1PRF41xkzu6f32SM1CRC4EFgN5wD3GmJ9FPX4V8Atgl33XXcaYe3ri2EqpnhGdjzg4a3zEVOPF9YaD5UfYPaqe0ecsYEAKF3i3GoPWInJLt3MWIpIH3A3MACYA/yEiE1w2/ZMxptL+p4FCqSwTnY8AiVgdr6USIMiBA88kvXqd31+X8ip5Kjv1RM2iCmg0xrwJICKPABcDG3tg30qpNOkyvmL2DXwwm3DOwpxYC4deITopHYuOm+hdeiJYnAjscNzeCXzcZbtLReRsYDPwP8aYHdEbiMgcYA7AyJEje6BoSqlk+P11tIxr5vjl87vmI3xz8PvraK1fTKg7rEhewq6tOm6id0lXb6hq4I/GmCMici3wAPCp6I2MMUuAJWAluNNUNqX6tIgaQHEBE7+3OpyPCCW8D55+PKa03X6GUFp6dcILv46b6F16IljsAk523D6JzkQ2AMaYZsfNe4Dbe+C4SqluCPVGOnz4bdcagDPhPcgL++/Ixz/Bmma8tPTKLvuJ7tWkS6P2Lj0RLF4BykRkDFaQuBz4knMDERlujNlj37wI2NQDx1VKJRDrQu6sTYjkI5KHMZHrZAfuXRSxqNHI2pNo+3c++eddEtH9NV5eQns89R7dDhbGmHYRuR54Cqvr7H3GmA0icguw1hizEvhvEbkIaAf2A1d197hKqfjiXcgjp+eA4cOviXhua80Sih7f1LmokUDxim14OiC45HZaq0/tMs2H5iV6tx7JWRhjngSejLrvZsffNwI39sSxlFLJiXchD+UTBjQcZthTHRQU/B87P7UP/wTD3r0PMH5VBYX2okZG4P0PF1K46b3IaT58czQv0YfodB9K9UJ+fx2HD7/t2rwEVvPQ+JZvcsL/3I4nALCHkpVQ/0uAw/D2W5g8q1ZhvBC4chbB7y9DjLVqXrxpPlTvpMFCqV7G76/jzWVnU7SunYGVeRQVfYJBrx0m70gD+Dov5vkv1uNptxc4AmiHYU9B6dMGT2AfwTxonV5KcOgJgL0Wt8vxNC/RN2iwUKqXaVl1Ox/9Vru9LGoHwvNIBwR/vSY8pTjYg/AWPW3XLMDkQ0HBUDyBfVZzk4GBzzUhpgnz6MbwfaajsxlK9R0aLJTqZY5fs7tz2g5j1whMZ66hlc6FjT746+8I3LsIELyzb+A4IPiEtXwq0rkPPGA8YEzkbLOq79BgoVQvc9yM2QR/vcaatiPPus90WBd5hgyNmCzwg+rfMfiRyJl5QqvkMWQog25eFp7+4+AtV8C7+1KabVb1HhoslMoxiab2di6L2j61kp0776BoXTv+Sfmc1Lil6+JFURd+5yp5rZPODtdCBmuA6NM0WCiVQzrHThxBxENZ2d2MGNH1Ih664G/ffhv+QQb/BADDCSeMIOily+JFznUsQrUGt/vilUt7RPVuGiyUyiHW2IkjQBBjgmzZcj2FhRUxL9BdxkHMmM8H1TPCQaCjqoLtj3yOk65aEbGuNuC61rYbnV22b9BgoVSWifcrvbh4OiIee/1qMKYj7qhp13EQvim0Akd+fxuHf7WdfIzrutrOhY8+WHUvzeOaXcuko7j7Bg0WSmWRZOZaKiu7my1brseYDjyefhQXT+8SYPat/A7tzzxO/nmXMPSin0fso7VmCcd9xqo1gNVlNphnj9Z2NE0517bYOfY1/G+96lomK4DlYUwwqanLVW7SYKEU2dHm7vfXsW3bgnAzU6xf6SNGzKGwsCJcXiAiwJyy9xKGf2UZngCY397OvouWcty1Pwo3IwVqHqMwajDeoQ9De8Vo+l1zY2fOwk6SHzz9ePyl1cSvOUjU/6q36fayqkrlutCv+bfe+kHSy4UeqzIcOPAMEAQ8Sc+1FN0MFHx2VUSz0pDHmzhu5rW01iwBoH1qJcF8eyoPex+D/g2DV2yL2O8A3xwG3/4UxTPm4/EUAHmuZWppqcWYdsBgTDstLbXdOBMqW2mwUH1a5K/5zl/Oqe4jlXWm3bZ3Jq5BGDz4vJiJYr+/jvr66bz11veor5+O11sScTH3fGoGQa81ASBEDsgDeH9iMet/Cbtnwnuj7G2wt7l3cZfyhfIeY8bc6lqmUBI9VjBRvYM2Q6mcFK/ZKNkmpc78wGGs39jierEL7c/rLSEQiEzyptoTKNb2Xm8JVqAAMAwdemnM/TQ1LcWYNmtL08ahQ+u6JLF3AoU/XkZRg71HRy7C6y3hYDkcLIeyO6Bwe+e+29p28+/6T2JMe0T54s3/pJMJ9g0aLFTOiXeBTuXi3flrPtQYI4wduyjGIkGhX/0ePJ5+4f02NS0NB5tkegLF6jkUCDRjVfStY1i3k9PW1tTlYl74FhT/y3HOvlpFAXBg/gUETj8eSq1j7b0ASv8G0m4luned24KxT0cqPZt0MsHeT4OFyjnxumq6PRa6P1RjCP3t7MUTEn2R7txfaJsgweCR8H6bmv5AKNgk0xMo1voP1v39kloXYuDASezZ03l7//6/4vfXRVysvUtXWGWyS9f/6Xr6PdC5ROrh6zwM2GJt2/zD8+nY+yZ7TtvKwfLOwKlNSspJg4XKOV5viT3WwHS5oEVfjL3eEsfyoXmARDSxlJXdxZYt12FMMNwN1alzf6GmKoAgXm+JI7ELIJSWXp3w13VR0RTGjl3Evn2PRTQ1pdKUc+jQuojboaSy8zkdwwbChvc6txE6k95tULYoiNjxL/i3p9n/6Hxk5x2MXNZOS6VQcM7FnHzyfK0tqDBNcKuc4vfX0dg4F2M6EPF0aTaKTsYGAs2O5UMDdlt/Z62jsLCC0tKvM3z4HNcmq9D+Bg/20dkt1BPOXYjk2/d7GDhwUtLlP3BgNY2NcyOS3EVFUxg16sajuEB7ugS5thuu7Fy8KA/avvEFgl4IegCPHTTskks75D9czUe/1c6Y+2DiPEPg+SdSLIPq7bRmoXKKs1nIGOHQoXVs335bxK/x6PbzUE0jumbhrHV4PAWUll7pesyioimMHr2AlpbnMOYIIHZC2hpBbV2SO9iy5b8pLKwIl/NYjXYuLb2Spqb77MDnoaRkZsTjfn8dm4p/xYDFQvF6YeBnv83Qi34enhSQIUMp+v4ycAzKA4noblu0rp2Wz+hIbNVJg4XKKc5mJpE8mpr+0KXnjlN08w50XshTvXBbgQFCgaG09GtAh+PxNpqalrJ37wMxE+w9sWZ1UdEUKitraWpaynvP/J7C+hW8OekJTrnieYqKpoRf18Fyw8FyD2PGFAMus8neuxjPvv0Ehw6G0ycRXLExPGLbPymfUzRfoRw0WKisF90VNnTxP3z4bfbs+T2JLvbRNQ23WkeophFdSwmxEtqRgQFAxOv4uwDAMQL7SJcyJZubSNT9t6hoCi2rbmfsvA5rRTxvOztPuJ2iy/+SVEAa4JtDK9ZkgZ5AE8EVm2j+wfkE9mzk8JmncMrMn2mtQkXQYKGyWnRX2LFjF4XzBUDEr/hUf6U7L9xebwmNjXPj1giiA0Np6ZV2k9BSwGoeeu+9Bvbs6ew5FWqu6s5rnjhxNRDZtOX315H/4tqIpqPj1+yGyztfV6hcsQRqHouYLLB1+zO8/SXB42lmSMqlVr2dBguV1SKbio7YPZdM+CIa6lm0efNF/OQnR/D5GvD5KpLef6jWsX37bXGbpJxNP2AFBmeOxFneeOMlkhkHEt08Ft20NXbsIhob5zJg/GFKHWtTHDdjdsR+Qs/Zu/cB1+N4fZdGTBbYMtEQb04q1bdpsFBZp6amgZqaZny+EqqqnDkKj503CEZcRBsaJjFv3tcIBAq48842qqutgJHK5IDJNN0kM/AserxEdNNWMnmS6LIAEc/Zt++xcE5i/UJh+BunMnDm/0asNxE6zqANHRTXf0DLbquJyil6Rb3W4l9BN3IpqnfTYKGySk1NAzNnnkogMN6+8G+lqiq6qchaJS4QaCIYbKO+/mwCgQKCwXwCAUNNTTNVValNw5FMLiGZ4JOoaSvZoBSdlLdqCdbrHjCgEr//BYLBNlorCjj+y0sZUDSF3buXhMdvFBdPp2ij8NF51viK4IMraC1Z0mUBI2fSu8A/S6fsUDFpsFBZpaammUBgfMSF3+eLvHhZTVEdNDc/iUg+lZXP4/W2EQgYvN4AQ4YIP/jBJk47rZLy8jqCwcNWs82A0yOaj6LFqzmkMo2Ie9PWEbZtW8DQoZcybERsG1UAAByASURBVNhXAWskdmgkuFsvLud9Y8cuCr/uXbt+FZG7KbIDRdNfrqW4Hpoqn4YpZzNi8zg8gY2dixrdu5gD9gp5QJclU3XKDhWPBguVVXy+Eu68s/PC7/NFJogDgWaMsdrWX399Mps2XcGZZ77NwoUXUF8/jUGD9nPzzXcQCHwCr/dyFi48l/Lylzl0aA2HDq2hqek+Kivjt8e71SBS6WbrnHjQqkVYvaMOHKjhwIGnAQ8i+TQ1Sdxuv06HDq0Lj+kIBtsIBJoZNerG8OMHn1rExFAtwgvrFz7P7rI8hto5DfJg0OMb8XRsJLjoaWswXkfiJVOVCtFgodLGmYtwS0KHHr/llv28+65x3S7UjGPlKWoIBAq45542Fi70ccUVt7Fs2Y0RTVL19dMpL385/HxjAgkv9G41iGTHR7j13tq37zHHOhVgDSgMhEqUVPCJNQdVqOmpZO2BiJ5RxfXw9hVBdt4/i0GvvQ9v76D4z5usx409cts5dbkGC5WABguVFtG5iPvvv50ZM6ZRVDSFmpoG7r13P48//nE6Osbj9Vq5CreAEmrPf/TRTRFBYf366ZSX/5NJk16MaJKqrHwu4vki3ogLffT044cPv+06TiLZ8RHRNZBAoJnRoxfYOYbOmWtD04SEahaxpkUPDR50m4Nq9+4lbN58LQAdE2C4o2dUS6X9WmfMp+jyKbTWLCG44lo8oVHbHjDByKnLlYpHg4VKi+hcxOrVLYwYcS4tLQ/z+c9fSCBQQGi2os5cRez9nXXWuyxZ0hkUPvvZMxgz5lYmTZpOWdlWamqamTq1hYkTZ+L1XhWefM+Zs3CbfhzycNYAnOMkErXp+/11HD78NiL5GEM4CEQnvZ3jRNyCT3Tt5MQTv4mIYIw1PXpoWpJ9+x4LP+dgOWz5zXhKGgZy8PT+FJ4xgVMdr3WAbw4HbnmeopuWIUEIChz8wgS8s2/QJiiVFA0WKi2icxGVlX8nGGxj6dIOAoF+dE6mHXTNVYSELqSlpW3ccccTNDZ+ixkzTsXnmwXMso+FI9DM6vL8UFdWt+nHO2eWtTQ3r2LEiDkRz090gRfJY/jwa7qMxYg1Ejta9NiSnTt/iTFBRPIiJk4cOvRSOwdiGXTBXIZ+bQ5D7dutNUs4ULOgM4n97j6r6clY+QpGnqSBQiVNg4VKC5+vgurqBlat2srYsXcwYcIr1hKgnsKI7crLN/LLX5qYA+ucF1JjOmhv3w+cmlQZ3PIJzgS0VbMQQtN6bNhwJvX145k5cwUXXTQrbo+oyHJB//4jj7pnUeT8VxIeWzLwdYM89gtaZ1o1hVAQC3WXdQa11polHDfTWr8ilMSOHoSnzU8qFRosVNr4fBX2YLlh4V/n11wzgBUrjtDe7iU/PxA3UIB7gnvJks6BePG45ROim4e83hI2b/4vNmw4g3nzVhMIFPDgg9b+x42L3SMq1QkCE43ZcHavbWycy4CGw0ycZ/AEGgkuuTbcg2nEiDkRQSLkg1X3Rkzl8cGqexl6xz/Dg/CcXWaVSoYGC5V2ziYZnw/++tf4vaSin+uW4K6ufoWqqta4v+bdLujOmoHz9rJl21wG+sUOCKkM6os3D1V07WXgwEkMG/ZVjn/kcTyBdyIu/vEu9u9XjSDoSHi/XzUCiByEp1QqeiRYiMiFwGKs7OA9xpifRT3eD1gKfAxoBr5ojNnWE8dWuc+qcSTezvlrfObMMyIS3Keddi/r169LOFgu+oIenWsoLb2a0tIrmTnzQzz4YOR4j6KiCtdV7pz7T6ZLrtW0FMRtHqZYc2ENKg8ywnHxP3h6/3Buwk3xjPn8644nKFrXbk03PmN+4hOsVBzdDhZirShzN+ADdgKviMhKY8xGx2azgQPGmLEicjnwc+CL3T22yj3OsRbW7eRqFNG/uKuqVrN8+dOsWtXIRz7yuD1SO891vEJ0k0/0xH/OHMiePb9j794HOOec1VRXb40oX2iVu2CwDb//BQoLK5LOS0Qex2Ov/d11nevIfEXnXFgHy4X1C63xEy2VUHjGhLjHKyqawilXPE/LZ2o5RafvUD2gJ2oWVUCjMeZNABF5BLgYcAaLi4EF9t+PAneJiBhrKK7qI5xjLRYtssYNdHRYf8+a9TeuuebEpBLboUkEi4sf4Itf7ExOOy+8oaA0dWoLxcVfSmIxotAa252D5Hy+GyNqPN1Z5S66CSx6uo6QWHNLieRx6COGg+XtiBRwaoxV/Zx0+g7Vk3oiWJwI7HDc3gl8PNY2xph2EfEDJcC7zo1EZA4wB2DkyJE9UDSVLfz+OqqrN4XHWhhjrWdtTB7BYB6PPupj5cojMRPVsWditQLF4MHnMXr0gvAgP+cAwIULJ1Fe/lLMqcdDaz84V91zS1B3Z5W7ZAf1hbYNPV5YWOG6yp8GAZVuWZXgNsYsAZYATJ48WWsdvUSoCem00ybh9V5OIGDIy7O6p7a3gzEejMkjEPDGHIwXeyZW68IdChTQdQBgaHR3oqnHS0uvTLg6XbIXfDdH80s/3ip/SqVTTwSLXcDJjtsn2fe5bbNTrHkOirAS3aoPCDXflJe/xMKFPlav/ibt7S2Ulb3Kli2n87e/fY2Ojjy83gBTp7bE3E/0hTPWhTt6AGBodHcqv+i7s41SvVFPBItXgDIRGYMVFC4HvhS1zUrgq0Ad8HngWc1X9D6xxg54vSV2stYgkseTT15EIFCA13slDz74By67bAHPPZdHZWUtxcXr8Psj8wqpLGIEnQMAO5PTnaO7e4PWmiU6VkKlXbeDhZ2DuB54Cqvr7H3GmA0icguw1hizErgXeFBEGoH9WAFF9SKxRjeHehAZ04GIh82br4kYv/DCC8dx3XV+Tjrp91iJ48geTX5/HcuWfYd166YyadJ3uOKKn3fp8uqWuE62O26ucRuZrQFDpUOP5CyMMU8CT0bdd7Pj78PAF3riWCo7rVr1AqtXz6Wy8u+Ul78SvuA7518yRjjzzD389reR4yOaml5lw4Yp1NdPY9KkF5k0aXp4v9XV/8e3vvU3uybSxqBBt/DlL8dfnjTVmkguCdQ8FjEyW6cXV+mSVQlulZtqahq46qrrwxf0O+64kEmTpjtmYc3DGBDJZ9Kkt1i+/GmeeWY/Xu/L1NefzZtvlnP33YvCzy8r24rPZ130n3suL6Im8txzecycWRezZ1IqK9rlIp3fSWWKBgvVbc7eR21tQm3tTVxxBY4Ry/mUlMxk//4n2bPn9xQXP8C0ab/gK1+xAoRIkGDQgzHOpVStxHhlZW3U+hS1tLQMYNSoG10T3N0ZC5ELBvjm6PxOKiM0WKhu8/lKWLSonWAwD2M8PP74J7nkksWUlnbOwhoMvm+PRrYu4i+8cFy4xiDSjsdjrR7nnJ68uHg6FRW3snDhedTXn0NlZS0f+chavF5rkj23nkndGQuRK3R+J5UJGixUt/l8FUybVsezz54JCB0dBTz+eAXTp59lJ6ZfZNasS+0V36x1GaZN+yBibqfrrptLIHAmM2eeER6UFxrXMGpULRde2MLOna9iTJDGxrkxp9ro7lgIpZQ7DRaqRwwd2hZxu7l5cERi+qSTnqaoSOxHDR/60LdYuHCZXWN4jo98ZC1lZZMZMSJy9Hao9rB9+20xJ9+LpmMhlOp5nkwXQOW+mpoGAPLzA4h04PUeYejQI47EtJdnntlvryNtMKYdYwKUl9dxxRU/p7z8nxjTQWPjXPz+OtdjhJqXIA+PpwCvt4Tt22+Lub1SqmdpzUJ1i3Mepry8di67rI7Zs08ATmDFis5mpvPOO8Exm2o+VtDoiJqu+wjbti2ImLojJNYEe72xx5NS2UiDhUqac3rxUF7B2RMKDCNHGqqqWmlpqWX58vHU1kJV1WbOOWcaEDm3U1PTUtramti/fxXGWGMxDhyowe9/wTUAOJukeqLHk9vrUUq5k2yddWPy5Mlm7dq1mS6GsoVqEG1tBXg8QX760zrmzz/HUbPw4vUGWL786fCU4M4aRLwV4azWUKuJyuJhzJgfM2rUja5liR5L0dLyMC++WJzSRb+z3FZOpbp6qwYM1SuIyKvGmMk9vV/NWaik1NQ009ZWgDH5dHR4uekmaypwax6mrcyb9w+qq7cyceImxyI/bRgTwFkDCGlpqaWhYRLLlv0vGzZ8jM5AASJ5cbu8hpqkxoy5lZaWh7nssvNZuHAqM2eeGs6fJPN6nDmVmhqd11KpeDRYqKT4fCV4PEGsi7oQDHrCF9iqqlauu66OqqrWiES0SAEiXiCPjRvPYvHi8eGL+fr145k3r4b77ruVefNWs2HDVMCDSD5lZXclNfvrqFE38uKLxUd10ff5SvB62/B4Isd2KKXcac5CJcXnq+CnP32Om26aQjDooaCgDZ+vxHV6jeh1J1ateoFvfcuaDuTXv26juroh4iLf1gYPPfQLbrjheWbMmJZS/iF6OvJkL/pdZ6bVJiil4tGchUpJdFJ4+/bbeOutHwAdQB5jxtwazjWEJvRbvHg8ixd/lmAwH48nwLx5/8DnKwnnQIzJQ6SDgoKjyx1oolqpTscqZ6E1C5WS6Km/k5nQb+zYs/B6z4/49R/6Zf/DHx6hrq4SY6waRqyV8lIpk1Kq52mwUN0Sa3qNpqalBIOHAcOECS9x//138dprVRG//n2+Ctate466ujysXlN5DBkisQ+mlMoYDRbqqEQ2/UROr+H319HU9AdCPZxE8pgxYxqXX941F/HuuwaRDozJx+Np5913s7NZVKm+TntDqZSFxijE6q5qTRjYbt8SSkuvjpm09vlKKCgI9Upq015JSmUpDRYqZfHGKDgXPLLmcepPaemVMfcVPU5DE9RKZSdthlIpi9Vd1ZnUFsln+PBrKC29MmFXWE1QK5X9NFiolIV6Mq1atZWqqs1UVU0DIpPaxkD//iN1gj+legkNFuqolJfXUVBwHcYEWbbsE2zefA1lZa9TXt6Z1O6Nq9Qp1VdpsFAp8/vr2LLleoxpZ8OGM5k3L7TI0aUsXHgu5eX/jJvUVkrlHk1wq5StWvUCDz30bTZsOJP6+ukRye76+k8mTGorpXKP1iwUkPyUGTU1DVx11fXhqb2vu24uXm9nsvvcc4tpaXmYn/zkCD5fg/ZuUqqX0GChIla7u/NOa6K/WBd552JHgYDBmE9TXb01HGhgRnidiOh9heaKco70VkrlBm2GUimt7RA9tfeMGafi81Xwve/1Y9y4J1i1amuXffn9dbzxxn9SX/9J3nrrB6xff66una1UjtGahUp6mm+/v45x46zlUh9+eLB97wlxJw2cOrWF9eu/FO5SC3RrKVSlVGZosFBJre3gDAg7d57FihVWD6gVK9q4//67KC21VseLnjRw3LhNvPVWG50r4UnE7LRKqdygwUIBiUdRt7TUhpdLXbduqqOpybBmzThmzeqcpnzGjGnMmAEtLU/g9ZaEpzAXyaO09OqkRnUrpbKLBgsVFi8B7Vy3YtKkFyN6QM2YcWqX1fGcq+eNHbuIQKBZE9tK5TANFgrAdXlU54XduW7FpEnTKSvb2qXZKrT99u23hWshwWAbgUBzePU8pVRu0mChgMhmplgJ6KKiznUrfD5iNlvFWj1PKZW7NFgooGcv8LFWz1NK5S4NFgrovMCvWvUCa9aM48iRAd2aNtxZC1FK5b5uBQsROQH4EzAa2AZcZow54LJdBxBaTu1tY8xF3TmuOjbWrBkQnsrj17+OP5JbKdW3dHcE93eB1caYMmC1fdvNB8aYSvufBoos5RzJ3dbWj5/97D3mz6/tsmxqNL+/ju3bb9NR2Ur1YmKMSbxVrCeLvAFMN8bsEZHhQK0x5jSX7VqNMQNS2ffkyZPN2rVrj7psKnU1NQ185jNlBAL9wveJdFBQYA28mzFjWpempUS9qJRS6SUirxpjJvf0frtbsxhmjNlj/90EDIuxXX8RWSsiL4vIrFg7E5E59nZr9+3b182iKadYv/5rahrCtQefr4JLLvknIkFAADDGmuNp9eoW1zmd3HpRKaV6n4Q5CxF5Bih1eeh7zhvGGCMisaopo4wxu0TkFOBZEWkwxmyN3sgYswRYAlbNImHpVVJi/fp3m2129uwTWLHiiN0clYdIO15vgMrKv7t2qdVuskr1DQmDhTHmvFiPicheERnuaIZ6J8Y+dtn/vykitcAkoEuwUMdGrDEU0dON19Q0c/vt01m+fAXPPLOfYcMG0dycz9ixdzBhwiuuwUC7ySrVN3S36+xK4KvAz+z//y96AxEZDLxvjDkiIkOATwC3d/O4KgWxfv27zTbr99dRXPwlLrmksxYCP48bDLSbrFK9X3eDxc+A5SIyG9gOXAYgIpOBbxhjvg6MB34nVkO4B/iZMWZjN4+rUhDv1/9FFz1PMPgeV16Zh883q8tUHS0ttYwadaMGA6X6uG4FC2NMM3Cuy/1rga/bf78EaGf9DIv+9e/MV3i9bXzqUxfi9w/THIRSypWO4O6jovMV69ZN5YQTXuC116qYOvVhJk7cpDkIpVSYLqvaR0Uvjzpo0H6uuup6Fi6cymWXnc/mzZ91DRQ6AE+pvklrFn1UaHW8Vau2UlW1mTVrLoxY0MiafjzyOToAT6m+S4NFH2atjmelk0pKGvj1r+Ovw53MNOZKqd5Jg0UfFb0qXlVVK/fffxdr1oxjxoxTXScQ1OS3Un2XBos+yNmctHHjWWzefA3jxv2eCRNeYtas0NiKrnQAnlJ9lwaLPijUnLRhwxnMm/c3AoECvN5LWbjwXMrLX4nbvKQD8JTqm7Q3VB8Uak6qr/+kI6ntpb5+OiKC19s1X6GU6ts0WPRBoeaks87aG9F9trLyOYwJ0tg4V7vGKqUiaDNUH1VUNIVLLoEPPriQdeumUln5POXl/wSC2tNJKdWFBos+rKhoCldc8XM+85lavN4raWx8TXs6KaVcabDoZWpqGuwBdSVJrZ/tTFgXFlZoTyellKtuLat6LOmyqqlbuXIFl112vt27qY3q6q1JBQylVO+Rrcuqqizh99fxxBOvRPRuqqlpPqr96NxPSqlo2gzVS7S01DJxYi1eb/wpO+LRuZ+UUrFozaKXKC6eTkXFOhYu9DF79gKWL3/atQkqXs3Bbe4npZQCrVn0GqGxE6NG1fLlL3cmqJ0J76qq1oiaw9ixiwgEmsMJbZ37SSkViya4ewm3XlCdq+FZCe/f/OabjBlzP2CtcCuShzHBiCan6AkGlVK55VgluLVm0Qs4l0i98842qqsb8PkquqyG99JLwxgzJhQoPBjTQfQgPJ37SSnlRnMWvYAVFLr2gopeDa+y8u+Ah8GDz6Os7G48nn5AnjY5KaUS0ppFL+DzlXDnnZG9oPz+OsaNq2X58vE8++whysp+R3n5GjyefowevYCioik6CE8plTTNWfQSsRLZIvmAwZh2RDyUld3NiBFzMl1cpdQxojkLFZe1RKr19/btt4W7wBoTtLcwGCMEAqkP1FNKKc1Z9EKhLrCQh0gBIl40N6GU6g6tWfRC0cufApqbUEp1iwaLXsrZBdbKZ0zB5xsQbqpSSqlUaLDo5WKNwVBKqVRozqKXizUGQymlUqHBopeLHpiX6ky0SikF2gzV6/l8FVRXp7Z6nlJKRdNg0Qc4x2AopdTR0GYopZRSCWmwUEoplVC3goWIfEFENohIUERizkUiIheKyBsi0igi3+3OMZVSSqVfd2sWrwOXAM/H2kBE8oC7gRnABOA/RGRCN4+rlFIqjbqV4DbGbAIQkXibVQGNxpg37W0fAS4GNnbn2EoppdInHTmLE4Edjts77fu6EJE5IrJWRNbu27cvDUVTSimVjITBQkSeEZHXXf5d3NOFMcYsMcZMNsZMHjp0aE/vvtepqWlg/vxaamoaMl0UpVQvl7AZyhhzXjePsQs42XH7JPs+1Q0655NSKp3S0Qz1ClAmImNEpAC4HFiZhuP2ajrnk1IqnbrbdfZzIrITmAL8VUSesu8fISJPAhhj2oHrgaeATcByY8yG7hVb6ZxPSql00jW4c5hz3W1tglJKga7BrWzRAULnfFJKpYMGixyiSW2lVKbo3FA5RJPaSqlM0WCRQzSprZTKFG2GyiG6kJFSKlM0WOQYTWorpTJBm6GUUkolpMFCKaVUQhoslFJKJaTBQimlVEIaLJRSSiWkwUIppVRCWTuRoIjsA7Zn6PBDgHczdOxUaVl7Xq6UE7Ssx0KulBPcyzrKGNPjq8dlbbDIJBFZeyxmbTwWtKw9L1fKCVrWYyFXygnpLas2QymllEpIg4VSSqmENFi4W5LpAqRAy9rzcqWcoGU9FnKlnJDGsmrOQimlVEJas1BKKZWQBgullFIJabAAROQLIrJBRIIiErMbmohsE5EGEakXkbXpLKOjDMmW9UIReUNEGkXku+kso6MMJ4hIjYhssf8fHGO7Dvuc1ovIyjSWL+45EpF+IvIn+/F/isjodJXNpSyJynqViOxznMevZ6ic94nIOyLyeozHRUTutF/Hv0Tk9HSX0VGWRGWdLiJ+xzm9Od1ltMtxsoj8XUQ22t/9G1y2Ofbn1RjT5/8B44HTgFpgcpzttgFDsr2sQB6wFTgFKADWAxMyUNbbge/af38X+HmM7VozULaE5wj4L+C39t+XA3/K0HueTFmvAu7KRPmiynE2cDrweozHPw2sAgQ4E/hnFpd1OvBEFpzT4cDp9t8Dgc0u7/8xP69aswCMMZuMMW9kuhzJSLKsVUCjMeZNY0wb8Ahw8bEvXRcXAw/Yfz8AzMpAGWJJ5hw5y/8ocK6ISBrLGJIt72dCxpjngf1xNrkYWGosLwPFIjI8PaWLlERZs4IxZo8x5jX770PAJuDEqM2O+XnVYJEaAzwtIq+KyJxMFyaOE4Edjts76frhSodhxpg99t9NwLAY2/UXkbUi8rKIpCugJHOOwtsYY9oBP5CJhc+TfT8vtZsgHhWRk9NTtJRly2czWVNEZL2IrBKR8kwXxm4KnQT8M+qhY35e+8yyqiLyDFDq8tD3jDH/l+RuphpjdonIh4AaEfm3/eukR/VQWdMiXlmdN4wxRkRi9dMeZZ/XU4BnRaTBGLO1p8vay1UDfzTGHBGRa7FqRJ/KcJly3WtYn81WEfk0sAIoy1RhRGQA8Bgw1xhzMN3H7zPBwhhzXg/sY5f9/zsi8hes5oEeDxY9UNZdgPOX5Un2fT0uXllFZK+IDDfG7LGrxO/E2EfovL4pIrVYv5yOdbBI5hyFttkpIvlAEdB8jMvlJmFZjTHOct2DlS/KRmn7bHaX84JsjHlSRH4tIkOMMWmfZFBEvFiBYpkx5nGXTY75edVmqCSJSKGIDAz9DZwPuPaiyAKvAGUiMkZECrCSs2nrZeSwEviq/fdXgS61IhEZLCL97L+HAJ8ANqahbMmcI2f5Pw88a+xsYpolLGtU+/RFWO3a2WglcKXde+dMwO9oqswqIlIaylGJSBXW9TLtPxbsMtwLbDLG3BFjs2N/XjOd6c+Gf8DnsNr4jgB7gafs+0cAT9p/n4LVC2U9sAGrSSgry2o6e0dsxvqFnqmylgCrgS3AM8AJ9v2TgXvsv88CGuzz2gDMTmP5upwj4BbgIvvv/sCfgUZgDXBKBj+jicp6m/25XA/8Hfhwhsr5R2APELA/p7OBbwDfsB8X4G77dTQQp/dhFpT1esc5fRk4K0PlnIqVL/0XUG//+3S6z6tO96GUUiohbYZSSimVkAYLpZRSCWmwUEoplZAGC6WUUglpsFBKKZWQBgullFIJabBQSimV0P8HZR1AEmy1isQAAAAASUVORK5CYII=\n","text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"tags":[],"needs_background":"light"}}]},{"cell_type":"code","metadata":{"id":"ygvR39vtbGSK","executionInfo":{"status":"ok","timestamp":1623426379862,"user_tz":180,"elapsed":59,"user":{"displayName":"Nicolás Chehebar","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GjWN7wrOddIeQwbW6mzML-_Z0jTJLZ1cnE5tDuipw=s64","userId":"03611693109350494416"}}},"source":["#z2=np.matrix([0 for i in range(n_samples)]+[1 for i in range(n_samples)]).T\n","#le pongo label 0 a la primera y 1 a la segunda\n","\n","alpha2=2\n","k2=20\n","dimension2=2\n","f_dist2,tree2=fermat_graph(data,alpha=alpha2)\n","\n","def my_fermat_cost2(x,initialData):\n","    return fermat_cost_cuad_fast(x,initialData,dimension2,f_dist2,tree2,data,k2,alpha2)"],"execution_count":31,"outputs":[]},{"cell_type":"code","metadata":{"id":"i50CqYihiq75","executionInfo":{"status":"ok","timestamp":1623426379871,"user_tz":180,"elapsed":60,"user":{"displayName":"Nicolás Chehebar","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GjWN7wrOddIeQwbW6mzML-_Z0jTJLZ1cnE5tDuipw=s64","userId":"03611693109350494416"}}},"source":["def prueba2(g,gg,q,qq,h):\n","    n_iter=400\n","    learningrate0=1\n","    landamax=1e2\n","    \n","    \n","    landa0=g\n","    beta_adapt=gg\n","\n","    \n","    h=h\n","    quantil2=q\n","    q2=qq\n","    fctr=1\n","    #cost_function=euclidean_cost_fast\n","    def cero(a,b):\n","      return 0, np.zeros(a.shape)\n","    cost_function=my_fermat_cost2\n","    #cost_function=cero\n","\n","    yinit=samples\n","    #yinit=np.zeros(samples.shape)+np.random.random(samples.shape)*0.00001\n","    \n","    bb=1\n","\n","    monoton=True\n","    clb=True\n","    #clb=False\n","    cat=False\n","    \n","    return barycenter_fermat(yinit,n_iter,samples,landa0,landamax,learningrate0,z2,h,beta_adapt,cost_function,f_dist2,tree2,data,k2,alpha2,landa_monotono=monoton,factor=fctr,quant=quantil2,quant2=q2,b=bb,categorical=cat,callback=clb)"],"execution_count":32,"outputs":[]},{"cell_type":"code","metadata":{"id":"0LFbx9kfjxQ_","executionInfo":{"status":"ok","timestamp":1623426379879,"user_tz":180,"elapsed":62,"user":{"displayName":"Nicolás Chehebar","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GjWN7wrOddIeQwbW6mzML-_Z0jTJLZ1cnE5tDuipw=s64","userId":"03611693109350494416"}}},"source":["def plot2(bari2,g,gg,q,qq,h):\n","  plt.plot(data[:,0],data[:,1],'.',c='y')\n","  #plt.plot(samples[n_samples:,0],samples[n_samples:,1],'.',c='r')\n","  #plt.plot(samples[:n_samples,0],samples[:n_samples,1],'.',c='b')\n","  plt.plot(bari2[:n_samples,0],bari2[:n_samples,1],'x',c='b')\n","  plt.plot(bari2[n_samples:,0],bari2[n_samples:,1],'x',c='r')\n","  plt.axis('equal')\n","  plt.title('Rojo es una distribucion, azul otra y amarillo es la superficie')\n","  plt.title('landa0:'+str(g)+\", bbeta_adapt:\"+str(gg)+\", quant:\"+str(q)+\", quant2:\"+str(qq)+\", h:\"+str(h))\n","  plt.show()"],"execution_count":33,"outputs":[]},{"cell_type":"code","metadata":{"id":"x1ksvbKcUhYA","colab":{"base_uri":"https://localhost:8080/","height":1000},"executionInfo":{"status":"ok","timestamp":1623426592850,"user_tz":180,"elapsed":213027,"user":{"displayName":"Nicolás Chehebar","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GjWN7wrOddIeQwbW6mzML-_Z0jTJLZ1cnE5tDuipw=s64","userId":"03611693109350494416"}},"outputId":"ee5af20b-48c8-4716-e89a-be8eb679bef4"},"source":["gg=1.5\n","h=-1\n","g=0.001\n","q=1\n","qq=0.7\n","plot2(prueba2(g,gg,q,qq,h),g,gg,q,qq,h)"],"execution_count":34,"outputs":[{"output_type":"stream","text":["it: 0, paso: 2.1252553677778524e-08, landa: 0.001, h:0.21201147385115363\n","it: 1, paso: 7.969707549949542e-09, landa: 0.0015, h:0.12956104785983666\n","it: 2, paso: 2.727350930589419e-09, landa: 0.0022500000000000003, h:0.11416772364598397\n","it: 3, paso: 1.0283244654252433e-08, landa: 0.0033750000000000004, h:0.10652928520269088\n","it: 4, paso: 0.2848833063698161, landa: 0.005062500000000001, h:0.09923260761931382\n","it: 5, paso: 0.29120371278408214, landa: 0.005062500000000001, h:0.09392107836423393\n","it: 6, paso: 0.2879564069284317, landa: 0.005062500000000001, h:0.09033367966166132\n","it: 7, paso: 0.2829723905876679, landa: 0.005062500000000001, h:0.08696573135193358\n","it: 8, paso: 0.269102700988659, landa: 0.005062500000000001, h:0.08448632717375708\n","it: 9, paso: 0.2718341749155948, landa: 0.005062500000000001, h:0.08125727421725266\n","it: 10, paso: 0.2747368185466185, landa: 0.005062500000000001, h:0.07865018893286116\n","it: 11, paso: 0.283269163595701, landa: 0.005062500000000001, h:0.07538855496754114\n","it: 12, paso: 0.29090310922387186, landa: 0.005062500000000001, h:0.07295214112169539\n","it: 13, paso: 0.2948923905508752, landa: 0.005062500000000001, h:0.07087681338748525\n","it: 14, paso: 0.3028144704552121, landa: 0.005062500000000001, h:0.06783590871142564\n","it: 15, paso: 0.3040306405584044, landa: 0.005062500000000001, h:0.06506423974503368\n","it: 16, paso: 0.31001319512722075, landa: 0.005062500000000001, h:0.06253928708649348\n","it: 17, paso: 0.31685133205476124, landa: 0.005062500000000001, h:0.06065713192527797\n","it: 18, paso: 0.3232044466035445, landa: 0.005062500000000001, h:0.058880261854470525\n","it: 19, paso: 0.31142452275550137, landa: 0.005062500000000001, h:0.0570400606028613\n","it: 20, paso: 0.30588504585312687, landa: 0.005062500000000001, h:0.05442438597279309\n","it: 21, paso: 0.30866973588349367, landa: 0.005062500000000001, h:0.05261872401724224\n","it: 22, paso: 0.3124144400156078, landa: 0.005062500000000001, h:0.05208601030991469\n","it: 23, paso: 0.2707717412953498, landa: 0.005062500000000001, h:0.05070425309462786\n","it: 24, paso: 0.2503279905757858, landa: 0.005062500000000001, h:0.049048331283382714\n","it: 25, paso: 0.24492525411116323, landa: 0.005062500000000001, h:0.04847959504191438\n","it: 26, paso: 0.2079492747595411, landa: 0.005062500000000001, h:0.047593810750686846\n","it: 27, paso: 0.18452335029997494, landa: 0.005062500000000001, h:0.048851254779541074\n","it: 28, paso: 0.1461409898881339, landa: 0.005062500000000001, h:0.0488856825496338\n","it: 29, paso: 1.0802590710566284e-16, landa: 0.005062500000000001, h:0.0488856825496338\n","it: 30, paso: 0.251346146682963, landa: 0.0075937500000000015, h:0.04483071106488682\n","it: 31, paso: 0.21702993228100043, landa: 0.0075937500000000015, h:0.044625666498281036\n","it: 32, paso: 0.16983172492210033, landa: 0.0075937500000000015, h:0.04412186308769469\n","it: 33, paso: 9.625717136179004e-16, landa: 0.0075937500000000015, h:0.044121863087694696\n","it: 34, paso: 0.2807728948478446, landa: 0.011390625000000001, h:0.04108376982393898\n","it: 35, paso: 1.724260108025255e-16, landa: 0.011390625000000001, h:0.04108376982393898\n","it: 36, paso: 0.20185440000908353, landa: 0.017085937500000002, h:0.039110030949045554\n","it: 37, paso: 1.0708021889121807e-15, landa: 0.017085937500000002, h:0.039110030949045554\n","it: 38, paso: 0.14187536254582903, landa: 0.025628906250000003, h:0.037731698490666724\n","it: 39, paso: 1.1502243538384529e-10, landa: 0.025628906250000003, h:0.03773169849080882\n","it: 40, paso: 1.98714382765989e-10, landa: 0.03844335937500001, h:0.036333236385205525\n","it: 41, paso: 1.5304275022579396e-10, landa: 0.057665039062500006, h:0.035042503981530956\n","it: 42, paso: 1.243122672209514e-10, landa: 0.08649755859375001, h:0.03445894490652396\n","it: 43, paso: 1.0077911905877697e-10, landa: 0.129746337890625, h:0.033996697919779745\n","it: 44, paso: 1.5685797491246867e-10, landa: 0.1946195068359375, h:0.033559952382724675\n","it: 45, paso: 1.2086945060957322e-10, landa: 0.29192926025390625, h:0.03323497422878886\n","it: 46, paso: 9.298836154844136e-11, landa: 0.43789389038085935, h:0.03287605284519098\n","it: 47, paso: 1.4190051869060453e-10, landa: 0.656840835571289, h:0.03276990437642721\n","it: 48, paso: 1.0885675601846533e-10, landa: 0.9852612533569336, h:0.03263414590367431\n","it: 49, paso: 1.642482319810252e-10, landa: 1.4778918800354004, h:0.03247561192281506\n","it: 50, paso: 1.06568238362767, landa: 2.2168378200531005, h:0.034180384085399515\n","it: 51, paso: 0.3987242954902087, landa: 2.2168378200531005, h:0.03422523137580678\n","it: 52, paso: 0.30829321669657017, landa: 2.2168378200531005, h:0.03429372119749452\n","it: 53, paso: 0.41116026402661465, landa: 2.2168378200531005, h:0.03400216633434501\n","it: 54, paso: 9.060481697271131e-11, landa: 2.2168378200531005, h:0.03400216633142384\n","it: 55, paso: 1.341907301298962e-10, landa: 3.3252567300796505, h:0.033945503199088144\n","it: 56, paso: 1.00934694492359e-10, landa: 4.987885095119475, h:0.0339210748215018\n","it: 57, paso: 1.5180362473989254e-10, landa: 7.481827642679213, h:0.03390133716524837\n","it: 58, paso: 1.140364950203278e-10, landa: 11.222741464018819, h:0.033875385412974426\n","it: 59, paso: 1.7135308299061063e-10, landa: 16.83411219602823, h:0.03384773735927792\n","it: 60, paso: 1.2876412601069795e-10, landa: 25.251168294042344, h:0.03383990829432874\n","it: 61, paso: 9.667458156036468e-11, landa: 37.876752441063516, h:0.033831641917493924\n","it: 62, paso: 1.4509757343187977e-10, landa: 56.815128661595274, h:0.03382420040656914\n","it: 63, paso: 1.0888354082122967e-10, landa: 85.22269299239291, h:0.0338191154875876\n","it: 64, paso: 1.6339408766377357e-10, landa: 127.83403948858935, h:0.03381433792798939\n"],"name":"stdout"},{"output_type":"display_data","data":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAYEAAAEICAYAAAC55kg0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOydaXgUVdaA39NJJyzBBAKSILsEBEZJWFQcddIKo+DghkZQQUFFUD91dECdEQFxxMGNcRlUFARX4oIiroDJqCMqyKIiKCiLQMISSADZkvT9flR1Ut3pTjokZD3v8/TTXVW3qk5Vdd1z77nnniPGGBRFUZT6iau6BVAURVGqD1UCiqIo9RhVAoqiKPUYVQKKoij1GFUCiqIo9RhVAoqiKPWYGqUERGSjiPQ7BsdtLyJGRCIr+9i1ARFJFZEtpWx/UUQeqEqZKkptlFlRSqOs9/RYUaOUQE1BRM4VkbUickBEMkSkXSllo0VkpojsFZFsEbmjjGP/1S63194v2rGtvX2+A/b5+zm2/UFEPhaRXSJSYyZ31PbKOJwGgogkish8Edlml21fxjE3ishBEdlvfz6pbLkrGxGZKCIvl1EmTUS+tP+fmVUkWoUI9nxF5AIR+UJEcu138XkRaVLKMcJ+x0XkGcdz3y8ih0VkX2VfV4hzR4nIm/b/z4hIajj7qRIIQESaA28D44FmwDJgbim7TASSgHaABxgnIueHOPZ5wN3AuXb5jsAkR5HXgBVAPPAP4E0RaWFvywfSgeuO5rqUCuEFPgIGl2OfQcaYGPvz52MkV1WzG5gGPFTdglSQWOABoBXQFTgBeLiU8hMJ8x03xox2PPcYrHf6jUqUvSy+AK4GssPewxhTYz7ARqCf/ftUYAmQC2QBTwFRjrIGGA2ss8s8DYi9LQJ4BNgF/ArcbJePtLePANYA++ztNzqOOwr40rHcGDgInBRC5m3Anx3Lk4HXQ5R9FXjQsXwukG3/7gwcBpo4tn8OjA44RifrsZXrvqYCW4C/2/dkI3CVY/uLwDPAQvue/Bdo59h+kr1tN/ATkOa4V/nAEWA/8J69/m7gF/tYPwKXhCHjicCnQI4t4ytAnGN7CrDcPuZc4HXgAXtbU2ABsBPYY/9u7dg3E5gCfAPsBd4FmtnbNtv/jf32p28pMkbaZduH+z8+inegtOu8FvgioLwBOtm/L8BqROwFfgMmOsq1t8teY1/zLuAf9rbz7WeYb9+DVWXIeD2QWc7ramj/z/bY/4mxwJZg1+H4T5bn+U4G/mfft0+A5uE+X+BS4PtSZA/7HQ/Yr7Etz5/K+Z7eCezAqvdGHOX/aAuQGk7ZmtwTKAT+CjQH+mJVmDcFlPkL0Ac4BUgDzrPX32BvSwF6A5cF7LfD3n4clkJ4XER62tu6A6t8BY0xv2NVaN0BRORuEVlg/24KJDrL27+7h7im7kHKthSReHvbr8aYfQHbQx2rvCRg3csTsCqC50Ski2P7VVh/7ubASqxKGBFpjKUAXgWOB4YA/xGRbsaY5+xyU43V8hlkH+sX4CysFtck4GURSSxDPsGqqH2tszZYLTBEJAp4B3gJq3f2Bv6tchcwC6ul1hZLaT8VcPzhwEis51UAPGGvP9v+jrOvYYmItLVNBW3LkLk0XhGRnSLyiYj0CGeHMK6zLH7Hus44LIUwRkQuDihzJtAF6326T0S6GmM+Ah4E5tr3oIctT9F/vRKYgKXoT8R6T68px77hPN8rsd7l44Eo4G/2+hLPN8jxzwZW+xYq+I47GYyluD4Lo6yPBKz35gSsXv/TtgyIyJUi8l05jhUeR6NljtWHUlpQwO3AvICWw5mO5XTgbvv3pzha0MCfcfQEghz7HeA2+/cLwEMB2/8HXBtkvzb2cRs41vUHNoY4zy/A+Y5lt71/e2AY8FVA+X8CLwasO9qeQAHQOOB+jbd/v4ijZQPEYCnhNsAVwOcBx3sWmODY94Eyzr8SuKicMl8MrLB/n43VGhPH9i9DnRdIBvY4ljOdzxTohtXyjaC4hRz0vxFw3HB7An/Eavk2Au7B6prHhXH8Uq+TMnoCQY43DXjc/u27TmcL+htgiP17IvBymM/maHoCvwb890cRZk8gzOd7r2P5JuCjgOsO9e73x+pddA6xvVzveMC+i3H0xsIon4ql4CId63YAp5fnXtv71f6egIh0FpEFvkFUrJZK84BiTrvXAazKC6zW5G+ObZsCjj1ARL4Skd0ikgsMdBx7P1YPwclxWN26QPY7tpdVNtixfb/3lfO8R8MeY/VqfGzCuk8+iu6XMWY/lumnFVbr6zS7ZZxr36+rsFosQRGR4SKy0lH+D5R8doH7tBSR10Vkq/28X3bs0wrYaux/t0N+376NRORZEdlk7/sZECciEcGuz97XXZZMR4sx5n/GmIPGmAPGmClY5sqzwti11OssCxE5zXYs2CkieVjm0nDfmWNNqe9kaYT5fMt9XSJyOlYP9zJjzM8hipX3Hfcduy1WpT6nLDkCyDHGFDiWg16L3VstGoAu5zn8qLFKAJgOrAWSjDHHYdmzJcx9s7A0uI+ibr1Y3jhvYY0ZtDTGxAEfOI69GujhKN8Yqwtb1F30YYzZY5/L2d3vEaxssGPbv7cbY3LsbR0DvBRKO1Z5aWpfi4+2WK1OH0X3S0RisMwR27Be3P8aY+IcnxhjzBi7uLPCQixPqhnALUC8fX9/oOxn96B9rJPt5321Y58s4AQRcR7Daaq5E8vEcZq9r88E4Cwf+H/Ix7KL+8l/jDCE998t6zp/x+pdACAigYr4VWA+0MYYE4s1zhPuO3Os70PId9LmAI5rw7+REc7zDUXQ6xKRFKx7NdIYszjkzuV/x30MA/5njPk1DBnLjTFms/EfgD5qarISaII1wLVfRE4CxpRR3kk6cKuItLbtaXc7tkUB0Vi2ugIRGYBlLvIxD/iDiAwWkQbAfcB3xpi1Ic41B7hXRJract6A1ZUNVfY6EekmInHAvb6ydktkJTBBRBqIyCVYYx1vAYhFA1t+7DJO99IXRSTUeX1Mst3IzsIaE3F6LQwUkTNtu/RkLNPUb1iDcJ1FZJiIuO1PHxHpau+3HcvLyUdjrBdvpy3XCKyeQFk0wWp15YnICVgDhz6WYJmzbrXPfymW44Bz34NArog0w7I/B3K1fd8bAfcDbxpjCm05vQHXUAL73vvud7S9HKxcWxH5o32fG4jIWKzW+P/s7akS2sW3rOtcBXQXkWT7/BMD9m8C7DbGHBKRU7Hs5OGyHWgvIiHrBBGJsM8bCbjs63M7tm8UkWtD7J4O3GO/J62B/wvYvhK40j7H+cCfAq6rrOcbihLPV0T+gOXt9X/GmPfCOEZ53nEfw4OVCfM9PWrEcmf1/Td9/8HSlWV5bU3H8oO/d9DZWD2B/VheMvfjsIdSujdBJPA4lqfJBkp6B92M9afPxRqEK/LAsLf3s899EMve2N6x7e/Ah47laGAmlsLaDtzh2NbWlr+tY90ddrm9WINd0Y5t7e3zHcTywukXsM0EfDY6ti8GbijF1rgFy+10F5bHxLCAe+fzDtqP1d3u4NjeBXgf64XKwRpzSba3JWG9wLnAO/a6f2KZk3YBj2F5G11fxrPvDnxrn38lVuvPaTPujeX54vOamet43q3s+7Yf+Bm4MeB5Z+LvHfQetveIvf1++9pygdNDPLfAe28c254BnnFcx3dYrfYc+7n0dpT1tRBD3YeQ12lv9z3D37B6S0XvAZYDxCZ73wVYg6cvB/x/nPbmTN9zwXJL/gLLPr48xH/92iD34UV7W5R93lBedI2wKtNcgnsH9cZqXe/DeidfK+fzvT5ATmddEfh8Z2Ephv2Oz+pKfMf72s+/SZD7UOZ7WkqdeJVTzlLq0MBn1L60fXwulUotxm69rwJOMcbkV7c8NQ2xJja9bIx5vgbI8jzwhjHm4zDLv4hVMdx7TAWrICJyJnCzMWZomOVTsZ5J62MqWA2ipr6n9TKMQl3DGHMEy61SqeEYY66vbhmOBcaYL7B6EkoIaup7WpPHBJQ6hpScUu/7PFPdsilKfUXNQYqiKPUY7QkoiqLUY2rsmEDz5s1N+/btq1sMRVGUWsW33367yxjTouySFjVWCbRv355ly5ZVtxiKoii1ChEJezY2qDlIURSlXqNKQFEUpR6jSkBRFKUeo0pAURSlHqNKQFEUpR6jSkBRFKUeo0pAUZRqJS9vCZs2TSEvL1jmR+VYU2PnCSiKUvfJy1vCqlXn4vUeweWKokePxcTG9q1useoV2hNQFKXayM3NxOs9AhTi9R4hNzezukWqd6gSUBSl2oiLS8XligIicLmiiItLrW6R6h1qDlIUpdqIje1Ljx6Lyc3NJC4uVU1B1YAqAUVRqpXY2L5a+Vcjag5SFKVymDoVMjL812VkWOsDUI+gmoMqAUVRKoc+fSAtrVgRZGRYy336+BXzeQRt2DCeVavOVUVQzagSUBSlcvB4ID3dqvjvu8/6Tk+31jtQj6CahSoBRVEqD48HxoyByZOt7wAFAOoRVNPQgWFFUSqPjAyYPh3Gj7e+PZ4SikA9gmoWqgQURakcfGMAPhOQxxPSJKQeQTUHNQcpilI5LF3qX+H7xgiWLq1euZRSEWNMdcsQlN69exvNMawodZdt255j5863aNFiMK1ajapuceoMIvKtMaZ3uOXVHKQoSpWzbdtz/PzzjQDs2fMJgCqCakLNQYqiVDk7d75V6nJZ6GSzykOVgKIoVU5MTLLfssvVKOwKXSebVS6qBBRFqVLy8pawdeuTgNifCHJy3gu7QtfJZpVLpSgBEZkpIjtE5IcQ20VEnhCR9SLynYj0rIzzKopSMynNXFNcifucUryUp0LXyWaVS2UNDL8IPAXMCbF9AJBkf04DptvfiqLUMcrKFuarxL3eI4hEAgZjCsOu0HWyWeVSKUrAGPOZiLQvpchFwBxj+aN+JSJxIpJojMmqjPMrilJzCGaucVbUgZW4b5/yVOg62azyqCoX0ROA3xzLW+x1fkpAREYBowDatm1bRaIpilKZOFv6oVr3zkpcB3arlxo1T8AY8xzwHFiTxapZHEVRjoLymGs00Xz1U1VKYCvQxrHc2l6nKEodJFxzjb/p6BDZ2XNUCVQxVeUiOh8YbnsJnQ7k6XiAoihxcamIRNhLhuzsWWoeqmIqy0X0NWAJ0EVEtojIdSIyWkRG20U+AH4F1gMzgJsq47yKotRuYmP7kpAwEmu+ABhToH7/VUxleQcNLWO7AW6ujHMpilIzyctbclRumwkJw9m+fXapA8nKsaNGDQwrilI7qcgAr/r9Vy+qBBRFqTBlzQ0oC/X7rz40dpCiKBXmaEM5aDTQ6kd7AoqiVJijMenoHIGagSoBRVEqhfKadCpqQlIqBzUHKYpSLWg00JqB9gQURSk3wdxBy+siGsqEdLSupsrRoUpAUZRyEcyWDxyVfT/QhKTjBFWPmoMURSkXwWz55c32FcorSLOGVT3aE1AUpVyEChVdVvhoH6W19sMJQ61ULqoEFEUpF6Fs+eG4iOblLWHjxol4vYcBbwmvIJ09XPWoElAUpdwEcwcty0W0uAdgKQBwBW3t6+zhqkXHBBRFqRKK7f2WAmjatJ8O/NYAVAkoilIl+M8LiKZ9+4mqAGoAag5SFCUsyuO/H6ys2vtrJqoEFEUpk1AePaEmjYXy/lF7f81DlYCiKGUSyn8/WGVfWTGBdOZw1aBjAoqilEmwOD+hFENx3mBBJAK3O77c4aJ9vYkNG8azatW5Gmr6GKI9AUVRyiSUPT/0xC6xvw3r1v0fxhSWKwyERhitOlQJKIoSFoH2/FCKITc3E2PyAYMxBXZpU67KXGcOVx2qBBRFOWqCDfS63fFYcwEADOAGvOWqzNWTqOpQJaAoSlCO1iU0Pz8Ha7jRmhSWmHjdUZ1fPYmqBlUCiqKUwBniQcRFUtLTtGo1qoyylummU6dpuFzRRctNmqSwfv3teL1H2L59dlHoaW3l1wxUCSiKUgJrYNaK8WOMl3XrbqFx45ODVtjZ2XPweg/hs/vn5+f4mXICB3mzs+ewfftszRlQQ1AloChKCSw3TxfGWLZ9YwrJzp5TovW+bdtzZGXNwLL9g4iraLuzYncO8gLq+VODUCWgKEpQmjX7Czk57wEGETfZ2TP9XD0B1q27GSgs2senNJwEDvICfj0B9fypXlQJKIrih9PGLxJBQsJIALvF7z8xrGSlb4K27AN7Bp06TWPnzrdo0WKw9gKqGVUCiqL44bThGwMNGrQlLi41aOvdGgC2xgNAcLmiiYtLLdWzKC9vSdFAcV7e535jDRoqoupRJaAoih/BJmqVlU3M7Y5n374VAPz++/d2JX8YkQiSkp7y8ywKNRtYk8xXD6oEFEXxI1SF7zPp+JLEOweAna377GzBmEKsGcNee9wA8vNziItLdSiZw4iIPblMQ0VUF6oEFEUpQaiJWqFa6/4mJBdW7CDLY8iYQtatuxljTNE+nTpNY926WzCmkPXrb6dx45M1VEQ1oVFEFUUJm9IihzqzhrVp8zdEIgEXIpH2AHLxPvn5OfY6/2TzPXospkOHyWoKqkK0J6AoSkh8A7Vudzz5+Tm43fFBW+vBTEjNm19ctK/PVOQ/qBz8OFr5Vy2qBBRFKcLpnQMUhY7wxQFyuaLp1GlakX0/MKpoqOXGjU8OOaisnkDVS6UoARE5H/g3EAE8b4x5KGD7tcDDwFZ71VPGmOcr49yKolQOgfb+li2vsU0/vrkAXrzew+zc+Va5k8QHa+Frq79mUOExAbFSCD0NDAC6AUNFpFuQonONMcn2RxWAotQwAu39gG3n91UTAnjZs2dR2Nm+fJ5Emhms5lIZPYFTgfXGmF8BROR14CLgx0o4tqIoVUSgd05CwnASEoYX2fWzsl5g376lBA7mhkL9/msHlaEETgB+cyxvAU4LUm6wiJwN/Az81RjzW2ABERkFjAJo27ZtJYimKEo4+MYCgtn7ffMA9u9fRXGguIgyXTjV7792UFUDw+8BrxljDovIjcBs4JzAQsaY54DnAHr37m2qSDZFqdeU1mL3KYdDhzY7UkUKCQkjy6zQ1e+/dlAZSmAr0Max3JriAWAAjDE5jsXngamVcF5FUSqAs4IvK4yDSCQiERhDkako8DjBvIXUA6jmUxlKYCmQJCIdsCr/IcCVzgIikmiMybIXLwTWVMJ5FUUpg1AVdGkVvK/F7kwWYww0bz6IwsIDfpE/y7L7qwdQzafCSsAYUyAitwAfY7mIzjTGrBaR+4Flxpj5wK0iciFQAOwGrq3oeRVFKZ3SKujASKGJiTeU2Dc7exa+MQBwkZPzAcYU+kX+VLt/7adSxgSMMR8AHwSsu8/x+x7gnso4l6Io4VFaBe0M4gZw4MCPxD77P/Z2MazqNZuWLa/BmALiVkCTtZA7KoV9+74tcSy1+9d+NHaQotRB8vKWcOjQZqxpPBHFFfTUqZCRQWxsXzp1mgYIcSsKOe6Zz9jbpZBuk7wc9+0h9u9fTtOVLrpNgv1do0hMvM4+lvh5Bmm8n9qPho1QlDpGXt4SVq5MxZh8IJL4+EFERSVYG/v0gbQ0SE8nv2MOcSsK6TYJfpwAuSnWd7dJhm0XfkOr+bDl0T/S4MyT7SNLwLeF2v1rN6oEFKWOkZ09B2OO2Ev5dp5gK69vj56LiU1Ph7Q0Wo4YROJzxQpAxE3h2Slsu/Ab2r8EG4fBb52+wWR9RXa2y5EjoEBt/3UIVQKKUufxAqY49HPPVMyVJxP38CwOjR1Bw4HRNAQSEoYT8dm3RM3/ho3DhFbzDbkpBeSmGIwxiLgwRtT2X8dQJaAodYyEhOFkZ8+yewORiIAxBYhE4nbHs3FWKl1nH2HT8EjazJxHlwFvg8cDGRkwchL7X3oSOWUfewftotsNj/HjBBd7e4WOHqrUblQJKEpNZOpUy37v8RSvy8iApUvJu/GsUidgxcb2JTk5oyjmz7p1/2dvMRQsXEDXiUdsE5Ah+vzBJNhjBCxdCunpxHg8xAC0g/0NOtDuy/dx9bhPK/46inoHKUpNxDeAm5FhLWdkQFoa+7s2YtWqc9mw4V5WrjybbdueC7p7bGxf2rW7x87g5bPlF9Lwhx2smRhFborlMdRw4HWQns7Bz99k0xWF5PVsUHSMvLwl5JyyD9fdZSsAjRZae9GegKLURDweq3WelgZjxsD06ZCeTk7Hr/BusJK8WEncbymauBWMQD/+qHsfpz0Q50gc81OrdLL/NAuzoaBoUhkQdgRQjRZau1EloCg1jKJQDz1TiR0zBiZPhvHjweMhLq+BPUBrJXoxprBUT53S4vds3jzV9hwqLFrnzBvsnGiWnT0npAlKZw3XblQJKEoNwtmqbroygpP/04ivzh1P7yemE+XxEOvxkJT0NOnpb7BmTU+uuupJ4uJSS8QI2rbtOXbufIsWLQbTqtWoEnGDrHkERwLOLkFzAItEkp09E2MKg7b04+JS7dhD3rBCTCs1C1UCikLoQGtVLcPGjRPxeg8Tt8JL10mF7Jg+jMMJ95N2sYf0S9KImpfOT4zigQdG8OST6UWmm5y7zmZvFy+bekVzwgn/x2+/TSVuBRxc+wkrRr9Cx44P+cUNsiaSORGaN7+INm3GlcgBfOjQZrKyZlB6Sz/4RDKl5qMDw0q9x9f63rBhfNhpE4+VDHv2LAK8NFkrrJloDdx6PHDbOx7STDqf/HMpaWkwa9bPnHHGZsCq1Pd28dohHw6zc+fbxK2AbpNg30mQl/cZK1d6iq7L7Y6nODCcD0NOzod+a3yDywkJw+00kxFB5whYSqUA50QypfagPQGlXuNsfYebNjHYMcrTiwhW3rKrWzKAsP+m/n7J3D0eyLjNw3mTPYwdu4XjjuvJhg35iLhJSnqSvb2i+XHCYbpN8vJZt96c8t0mfpycT26Kdc7ly8/go48OMmUK5OfnYLXY/RWBMYfJzp5TJI9PvrLyAmgQudqNKgGlVlJaxRtupVxsf7di5gfaxAOP53bHl5gs5bThv/76XVxwwXEMGtS9aF/btZ9x4wLP6e9JY7XOvfZexi9mv+8406db48NPPRVLYuIZpKRk8uqrt5Gb25Wrr15Ku0vn481aR7+HZ/FgxN9Z8mlb7kwZzYoVqUyaNJeXXtoBhOoJ+OT7rz3RrMBPvtLiA2nymNqNKgGl1lFWOsRw3RWLW9++ClHo1GlaiOQrvla6C5cruui4zsQrXbp8zbXX/oM33yyegOubh+V/zpKeNFbr3FV0DmvZwnkcjwdOPPF5brstnQkT0jjppKXcd999fPxxDOPP2sGYjMf4V9R4bjgynS8WzeTBw7P5/PNLmTHjFc4++xQ2bZrCoUObHecSx2/DgQPF+Z7K0yvSIHK1F1UCSq2jNJfEYNt8630tfN9vp1eLD2fl6yt7wiuH2HeSsU0rXrzewxz84AWi1r1M9p+KE6/07Pk5L764gbS07k7Xfr9Jv6FMJ9b66KAmFXsib9Fx+vdvzIQJaaxd24ehQx8mNfU1ChZ2YPiCoVxAOksbe2hyYQ4vzxvO4IVvc1b/t+nc+QlWrdpQ5O1jVfwGcNGmzZ3s37+yaDzCQmME1RdUCSi1Drc73vaVNyUqqsBK1u2Od6RRtOLhO00dSUlPsW7dzRjjxeWKLlHpxcWlsrGrm25FoRYgboWXFg/MZedTV5RIvv6nP1kKwOHa74cvjr/PfdOnvEozqfhMST727VtBSkomKSmZ9hovKQXLSCOdTDxEF8Bt70zjjcLBXHLCK9z1+TSWLfuUbt1+wsok5uv5GKCQLVumkZT0pMNrSGje/EI/TyGl7qJKQKlV5OUtYf362zGmEBFXCfNNYGXqn0ax2Obu6yXExaWSkHA9YAVeC6z0YmP70n5EJptjbqXb35ax7UJoNR92Tr+ChgOvQ1a+bB/XRZMmKX62++nTLSXgVAQ++b3eI35pGn3nOppKd9euVkwpHIXPPfPwYQA3PzRPYsX+KYwYMZHx46czfnw2ycmLEBE7lIR9N0w++/atwFIK1icn5wPatBlX8mRKnUNdRJVaRXGl7sUYw759K0rErPG5NjrTH0IEIm5Eil0dfb2ErKwZbN8+O+Q5Y2P7csK2U9l9mov2L8G2CwXvn04l4rNvaf1qPr4WdXr6Wyw+7yEev/BV8vM/5fbbN5CWBm/clEHmwKlkZMBDDx0Maq5yYif/8iMjw1oPlrKyrgNWrDiHlSv7AeB2+0pbFfmuXa256qqHuPfezrz5ZkP27n2WDh0mk5T0dNH+YOURABy9GksxqKtn/UB7AkqtwmnuEYkI6sniJLBnAATtJQQdBHVG8oyIoOVCL7t7Qes3DBub3MTxr7nZN754POGxx/7NZd0/5/xZt7JrkIeJ//0PVya+R+r0kWSOSeemNJg5s6Wfuer42VvgrIyi7sLUqdBhYwafTVoKC8bh8cBnf5nKw4v7MPYDT9E19dr7KAc/e5nX9p5Kv35ziIgwvPvujTgnbSUlfcsrr9zNmWd+wtCh4PG0w5fqu3Hjk8nOnsORI9lERSXQpEkK2dlRGGPlHBZx63hAPUGVgFLjCXT5LM9MVqse74vHU7x++fK+LF0KN95IifGDTZumFNvkf/mFQxP+yW8jJtBx7gtsG2SZgvY0b0jH/xxk512nktfzW7DDLwwaNJPH//MQ+d2zmfDev3E3vYe0te9yZ+t03prtYcEC8Hi6k5dXrJQaFh4qdv1ZupRLNkbS9JkpHB7xby677AC3nfIRLTN/4Z2o+4lacT947oCMDGJGTuL365IYx0PMPyeV8ePn4XJ58XojAHC7D7Ft24lcc82DrFt3c4l76rtPvvGS7dujaN36NnJzM4mKakXbtjoeUF9QJaDUaAJdPp2JTcBKmej0qAkMw9+nD1xyCVxxBTz7rL+7pVOhuN3xRbb6ol7FkCFEvDKX1tP/wbbU/iQsfo/CSKHZroP82DuRDo1SSE5+tGiC1YMPXsK+fT/y1Ev30iZ6D+N2P86jje7ipS0ehg0rOUgM+EcLHTCApJdf5ocRf+OWt/ozuNkMxmQ+wPTUZ5HoHzB33kn+NwuJWryMg7cPJX7K0xR64NNtQygsjCQ6Oh/L1CUY46ZPn5957bXxzPSGcvAAACAASURBVJ69Oei99e8JHWbLlsfswfbvadtWxwPqC6oElBpNYEVlefKYoora52kTE5NMbm4mXbseR1padz+XSmNg7lxo2bKk26ZvMHbTpiklTUOee3B370LE8hW0yXyP/Igofi9swP7uCZy06ldcUy+jYcBg7uTJU/h9wUpG7HmJByPu4oYDL5Dd6zxe/NBDRgb07BlkHoPHQ5FL0bBhdH3naZ5q9jPnbfgfd3e4l3mr+tNt6BQudUPU3I84fFonIqc8CcCOc4BPISKigMsv38/w4VZC+YsvNsTEfMP48Xfw/vuncfbZfyl1pq/lbVXI0c6aVmovqgSUGofT/FNaRZWdPcfuCRxmz55PABexsZOZOXMpaWndebHbP3huxZk81/lLfm12Dn+f7GHmsAw8S5eyYkUf8hYtJfUDq8UbMvSBy4WrwAq2Fll4hG2dz6frj+/Baaf5zeAaOBD69YOCRUk8s2cM1zR4ng8PXcRnUefw8vI0eoz4N2lpV/LEE7+SmBigbJYfKnYpeuwxNib05Opf3uXlDhfx1u4rmTrgEfbN6MqRiBVEeQuJ/no9Xjd89y/ITXFx4qYIRo3a6jdT+Zln5pKZuZ6UlMWkpHxKdvbuUsdLAntCOh5Qf1AloNQogs34LVlRHUbERX5+dpGnkIWlHE45ZT4jRsSy5eFdvOq+hA+j+3L9wmdo2+se/vLSJH7pehZt1j4MjxRP5Y2N7Utm5hqSkr5g4MCOVoV5440UfruyyIXOIJz083xrTu/ll1srx4yh//rpuN2Q/eFSAOaddCt/WLuWD7mQj4/0Z1L3W/Gsfp6ZM3uwYsWZnHBCsbKJ/64JjCxWJltW7aLj/On8cuKZDN2zkOTBD9Nm5rO4ow4zP/ksBq/IhEKQfIhZJ+ztFc2kSd2Jje3uFz564MCTadVqJNaUAEN29qyQLrC+dY0bn6yhH+ohUjxxpGbRu3dvs2zZsuoWQ6liNm2awoYN47ESnUTQocNk2rW7p2j7tm3PFU3usma+ij3ByQqBIBLFjh1vcOONZzG291+5aeHbAPzctTO91iyjwBVFodfF1xf+k9R37yg+cUYGv7y+lNPfHldkLtp5+iDiv36fF9yjGZn/LBF4McBa6cyJDbchEsGsNvdxxU+TMEa4hHkApJNmT9z6E23b/si2bUlMnTqQc86JpEWLwbZPPjRpkkLDJ+fhPuMCYgbdAhkZHBiUxu5zL6f1otkweDC8/DI7T72ApTEn0v+rfyPA1sGRtHnTBQWFHHhzGjGDbmHbtuf4+ecbiy7nuOPOJiqqGbt2vYvlMhpBYuINNGjQtoSXlFb4dQsR+dYY0zvc8toTUGoUZUWktHLmGqx5AgU0adKbqKhW5OQsAApZvvxMJk06nQkTLiUlJYN/7fsr9371DH3WLKNAXER5j/BZu8s4dfEUyEjxC/JzYno66UPg8svzueqqr2n2/Rl0Pd3DiJXjceUXu4K2jP+N/LwICvINp7TJJWrtEQpx4SGDMUwnjXTAcG/UQ8T/9TOMyWft2l6kpDxSZLayErUI5rwCXK7/0iOvF7FLl9LovXQaeTxwX/OiMYIWf/gDZ3w7CwF+mOxLED+MhDFvErPgexgEWVkv+N2nvXs/AyIRicKYAr/EMMFmTqsiqL/oZDGlyigrGblvLKBTp2l06DA5aOVUPPnLCnq2b98ydu9+H18AtLVrezFhwhWkpFizrdLSHsfttnzfI4yXLzqcQrfNmWwdfo9l07/vPj/bfs+eS7jggkd44okz+fUyGPSXjUQeOoAAWzrGkU8kzXYdZEnh6azpfxunL5xMXv/LcQH3MZkPGQBAOldgeu/i6hO6c8Wm3xg69FGKA9V5MSbfzuzlmDQ2blyxUvKNEXz4Ifu7NmJ303W2AgCRCCtB/Lx57G2xi1WrzrMD3AVSSELCCDp0mExCwgh7PKUw+LmVeouag5QqwWnrF4kkIWFEkY06L28J2dlzypz45TzWxo0THQHPXEWxhHyt3D+MO8yhBDh+MZh8FwVHoomOOIQYw4/XjiVh1iwK+g0gYeFL5N7q4cBdQ8jPzyEzswG33XY1F144nfnzx7CpUQ9is7Jg2DC8775N/oFCpKCACDHsN43ZM+w22s9/gv15+URSSDSH+Z3GDOJdABY0Gox5/QGWx44rikT62mvj6Np1BSkpnxVdb17eUtas6c64PlavZP/MCeScso/475rQ8Nq7+P7eQ+SmWCavxMQb6dJlegkTUHE0UAuRKJKTM4vucXE0VCnKU6w9gbpHec1B2hNQqgT/GD6Hycp6llWrzmXbtudYuTKVrKxn7Nmq4bVOGzToaIc7iMDliiYp6Wk6dJhMcnImyckZSP8BtJoPEQWRrOw+iM2eFFyFBomK4g+ReewZfQ/NMt9kT08h5oUMdsy9kVtuiefmm0cw48r+TGg4gRlX9iciK4/xjR/ljZgRFBZGcqggis8638ARE4m4hOffiiMvz1CAm3QuR4DG/M7dPEQ6Q7ij9SvEDLqFnnlT6fHxn+nc+VlSU0/hgQfeZ+/eb+nQYTJ5eUsZObI7ffoAS5eyf+YElseOY8OG8SyPHUf2vwdx3E9ghbFuQELCcAB27nzL7540adKbxMTRxMdfTGLi6CIFAMWB60R80UOFxMQbVAEoOiagVA3Ftn5fAhcriFtW1gsBCc9LD2Hs36OwBjuDeb3wwAew8i/I++9zasxuyFwCjz7K/qQoVo//iOQNb7B79pVkZ8+i+3jDH8YJZzRdTWHh51ww50d+ngxd0nfwesQQphy6gzumT2Vm9DyuGplF8xfnMJAPaOAWJjZ8mIsPvEMv17dM8k7gg16n0e/b5ZzHIuZwNX+OWw2P/UTMlClw6aU0/SmJoQlw6uDHOX3kOKYO2MHmt94nfUF3a+6CZxw5m6bg3VA8N2J9m7cxQywzkDNgXosWg+0xBovExOto1WqU371yzoD2H08ppEGDtqoAFO0JKFWDzyc9MfFGvyBu0dGt/Mo1adInjEQwvh5FAYcO/Rr6pAsWwFlnweefw5lnknddX5bHjuNIm5+4rPAFPo+8iL29onlneAqHChoybveTvM0lXMFcPvzwWlp9fYTYUct57bX3eCxiHB8dTuXa2ZczwPsB3W76iFZXncQZez7gz/3hLu+/mHHhQFpd9TXeBl7yI9wM42V67lkEf/sb3HMPDBkCF18Ml1zCiUP6MHVABn95KY22g/sEzTlgBb0rnhthTCE7d75VNKbSqtUoOnd+lqZNrR5GoAJYudLDhg3/KMov7DyuzgVQfKgSUKqM2Ni+dOkyneTkzKKB3zZtxiESjWWnji4RGjqQwIHhPXsWhU4O/9hj8MUXliL44gsKH56M13uEwjvXc87kpxk9+lwWLPiZ0XO/YPPgQUQUGqIkgus6PsH1Cxcwo98gWqSt5PjjB3PZZTsAobDQzcknf4ExhcyalcDDD8M9/Zay7oGpLF58JZ3Hu/lpSgSH5z2G9OxJh3ULrVlkU6ZYA74iYAw/P5vOoJfTePuKpxlnzyZ20rLlNSQm3kBS0tO4XNGO613olzS+VatR9OjxsZ8CAMjOnmOb10xR7mCfIg416K7UU4wxNfLTq1cvo9QPcnO/NBs3Pmhyc78suXHAAGMefdRv1YEH/s/k/THeZGS4TEYGJiPDZVau/LP//o8+aoxI8b6PPmq8ImbdTZEmIyPC/Pe/Dc3Ysb8ZMOaZK+aagmaNzMGxI0xBdENjwMw7+U8mNnaHeeyxVPPYY+eYRo0OG/AWfUTyzciR2/zkWnTeXeb2Hi9Zcnz6qTHNmxszbJgxjRsb07OnMWD23OoxqwdfYQyY1b1aml9ujDTz5/9gmje3dsnN/dL8978Ni2TcuvVZs3btaPP1193sa7U+a9eOLvWerl07ulzllboDsMyUo66tlDEBETkf+DcQATxvjHkoYHs0MAfoBeQAVxhjNlbGuZXaT6nJVPr1s8wpQN51fSl8eDJNH/wIJt+Cy/V8kcfNnj2LyMv7nD6ZI2h41mWwaBE88gjcYUXdpKAAeeQR2n30Ju4Og/juuwuZNas1z1yRzuD0m3h7dBfk9+/of9hN4yjDhRu/4flh53Pt+MUcORJNfr6bRx8VvvzyJ956qwvGRPD664lcfXWxV+eQbx+yAtMtD0gKnJyMufNOjBuaPJdB+8IofuzVkq7Lt/PLqcIpp8wnPb07S5dCx47BYyUFJoY/ciS71HuakDDcnheQj4i7aDBZUQKpsDlILJ+8p4EBQDdgqIh0Cyh2HbDHGNMJeBz4V0XPq9ROnHMFSswbCJZNJSUFLrgA87e/Yc7+I00f/JBfxkRw5JahdOo0jQYNOmLF0LdCRuQlHbYq4LFjoaAAHnsM7+WXkt3mZ/Ku64v7nn9S8GAEI0daQebSOswiY3Rb5swZT//nf+KTMR1Z/S8v+YP/zKVzf+Xu0+Zy/PEuHn1USEpazaJFzejf/yUiI4/QuvWBYFMN/JMCZ2TAlCns+ccAdvcC1yFowBE6r93OL6OFtq8a4r9rgsdjTRPwHw/w5T8uJFAJREUllHqfY2P72ma3f/p5CSlKIBWeJyAifYGJxpjz7OV7AIwxUxxlPrbLLBFrrn820MKUcnKdJ1D3cHr2FPu0W37+zZsPov2Gc4kZOQkuvdQaRIWi2vXImKFE/bSd3JNh5RMRdH+/N9tOWMGe5IKiYzRdGUnnvdfT8KzL8F5+KQfObk/jd1byy5hItlxuaLoygpMfaET64LdpOcSDx1MskzycyDf0odedcykKV/Hr6VaFPm4cGRlw2WUHuPfeC0lJWcyKFefywAPzueCCRrz0kjWv6/77g1y0HVku77q+5Nx1No1+LSBhIRS0b8HWzL8S/10TYtYc8Esk7Js05wzqZrW1jD37N4rk5Ayt2JWglHeeQGUogcuA840x19vLw4DTjDG3OMr8YJfZYi//YpfZFXCsUcAogLZt2/batGlThWRTag4lJ3iVRCSaXnsfIWbYvXDkCERFwbx5LP/7G6R8NZ3fO0LjDfDLmEh+aXIDZ/z7RX558CC5KS5ar+9Jx7vX43rjbfJ6NiD31rNpN6eA7D8Lzb42RbmBd0+/joShz5eQrazJalOnQteuq4mN7VMU0uLzz1czZUoH7rijZIjqIhwJDPbvX02jq8YiBQaJagDz5oVIMuAvW7CsaKoAlFDUaiXgRHsCdQf/2arWrNdA84aF0KHDP2n3wkErbk7DhnDmmZiFC5kZPYbj37iZP349lqYPfsQnkefTqdNS2m3eZQVUez8G19/Hw08/sW/fCqI/WFpc8Z8GCZ/ApuGRxD3xWakzkcuqZH1lvvvuwiKTkiP8UGhFcPHFlnkqMhLeecdaH3IHRTl6qmPG8FagjWO5tb0uaBnbHBSLNUCs1AOcyeHBRaNGJ2H5EATiovE3u+DRR2HYMMudcuFCpH9/ej54OcuHvM801weMb/gIp5+YQ8eth3DlR9BuTgGugYMse8zLLxPz/hrWTIxi48gINl8VScuFsP+SZNq8H2PF7g+BM0F9WWXWrPFPXONLELZ0aZCdPB4rxdmBA3DbbdZyqTsoStVRGUpgKZAkIh3EmgU0BJgfUGY+cI39+zLg09LGA5TaSagAcW53vB2uwIqeefDgr4BBxE2bNuOIj78YiCBuhZfjbniMg3+/Dubbf6FGjfD+73N6TBxIwqD2TJ4MrjvvIHbN1xwaNwLJL8RERsBLL0F+Plx9NfLOfNqPyOSkrGvpODcGeeRRYk4fiuuNt63Wd+Dg81Hgi/XmxDe4W4KMDFi1yho4mD69+PwhdyidsgLxKUp5qLCLqDGmQERuAT7Gat7NNMasFpH7sfxV5wMvAC+JyHpgN5aiUOoQwZLB+AKXrV9/ux3C2EV8/EB27XoPawYs7N+/0vbwgSZrDT9OcJGQcICGxoDbzZHz+xIx7yPyvREsWOBm7NgtTJ/emvOjZ9H74SfZ3g8SFhYWCzJkCHg8xAKxv30Ob7ztX1v7Wt9VZYIJtBN5PBUyA4W6z4pytFTKPAFjzAfABwHr7nP8PgRcXhnnUmom2dlziuICOXPUOk1Bxghud4IdQ6jYv98KBCf8NlQQiaTLf6Mtu3lGBlGTJ/NV/3b888vHGef5G3HDTmHAgHl8csFWjvTvxJmL11MYDSJuq1v7+uvWwHBuJnE3BrHt+yriqsLpLuo7fwUUkX/OZc0FrFQcDSCnVBjLu2YWvsFekQji4lLJy1vCoUObbX937Exg0KnTNKZOddGu3aukpGTYAeRcrFiRyk8/nU7yvwbRcPkhvP95ks3DhO5v7GD4iEkUpv1KTs4mevZcQvQ/TiD5/t8wblj9UBQnnvgoMcPuxfv6K2zs+iJ7kgtrRks5mLmnAoqorKQ7ilJeNHaQUmFyczMxpqBoOT5+IACrVp1LVtYMQIiPHwQYsrJmcP/9m2jUqCWTJs1lxYpUwJCefit///t8unT5hoMfvABpaex4cjAbRwq/PHiQi15dRdwKAENubiZnuHcSceX55My4jvYjMq30jPPm8fsF3YhZk09dTZii8X+UykZ7AkqFiYtLRSTSjnYJu3YtwO1OcET7BK/3AK1fzWffSYYuXb5m8uSxXHnl/Xw6/lZiWh3PM+sfYfTosfTq9SWxmV0gPZ2GPRvgWvUquSmH+HGCoclayOvpwu2Oh3GjiAb85s16PHh7/putq86FOtxSLjXMhqKUE1UCSoWJje1Ls2YDyMmx/d8p4MiR7KKwByIRvPHG/SRH5nDhpG+Rif+lxRNzefUGDy8dvpbL1r1Jv/4vc+ONP9C+/WIanm0nQgF69FhMbm4mBW1y2dLzcYwpZP3622nc+OSgFaGvpayTqhQlPFQJKJVCYCwbEbAmhVnfPXrs4PobPmTGlX/mookrOXjav5h9aB+X8SY5J0ew+ZsB/PxzU3r08K+0fa3eTZum2HF0vGUOiGpLWVHCR8cElArj81f3efmIRON2J9jjBFa8m86dn2LChKHc8OonfND2T5z+yQamMxp3/9/YvLkrV175EDfccDrvvbc66DkCE6K43fHqK68olYD2BJQKUTLd441FYYu3b59d5MXSosVgevW6nbtOvYs/LlzNA/J3/hr9KBcNPMynA1KZNOl1rrzyIRYv7srZZ+8t0ZJ3mnmcgdVqhAeQotRitCeghE2wmar+6R6tvLW+9Z06TSMx8QZatryGxo1PJvLzp7lh0Xv885SZPNxkEl/edi3d74/kHPmcCROuoLAwgosvHhUyU5gvZEN+fk4JX/nKuh5FqW9oT0AJC2cQOJEIkpKeolWrUSX81mfMOIfmzSeQnLzInhdgWL78LD79dDW9P91Jw0fSmXaHh3PfW80VwyYx48qv6LlmFSlDM0lJscIpeL2HS7X5B57TZxoKNRA8daoVusfpmv/ee6t5//0FDBnyL+1NKPUa7QkoYWG1+K1Zvsbks27dzeTlLSnht969+3omTnyJFSvOwpgjLF/+RyZNeg3w0vn5tqTcYdXEp5wynwkThrCosB+/DfXijCrqm2wWCuc5O3Waxvr1t7Nhw/iQPYg+ffxDBmVkwLXXdqBLl6+pq/MJFCVctCeghIU1F8CX6QqM8QZtrQ8c2JGJE4cxcWI6F100g3ffvYEJE4bSs+fnJCSMJC+vI7GxfXG740lJ+ZSUlE/tPSMBLyIukpKeKrNV7vQaKiuMgi9SQ1oajBljxXB78cUNxMZ+idcbUWfnEyhKOKgSUMIiNrYvSUlP2Tlvvbhc0UWhIQIDmo0YMYENG37iiSf+ztixWxg4MIns7C/IyprB9u2z6dFjMfn5ORRnFxOaNOlJTExPEhKGl8ssE24YBY/HUgCTJ1vBPAcN6k5ens4nUBRVAkrYtGo1isaNT/arOIO1xH/99R5efRXGjt3CCy80o3v3zrRrV+BXxqq8o4tMTPv2LeP3378vd0L0cCeHZWRYPQBfNGcrfI/OJ1AUVQJKuQiciBXYErcybsHMmVYqxsTEM7jttteYOPEjkpMXF7XWfZW3lXJyIdYksNIHhMOVKZBKjuasKHUKHRhWKkTgwLAv41bnzk/h9R4iJWUxEyYMZevWG0sEPYuN7UuLFoMpHhT2WnGBKpnSojkrSn2nwjmGjxWaY7hmU1o+3ry8Jaxc6cGYwwCIRJGcHLyFv2nTFDZsuBdf6skOHR6gXbt7quAKFKVuUt4cw2oOUspNWdmt/ENLCwkJI8vw+Y/W+PiKUk2oElDKTWnZrQITybhcUaUO9mrUT0WpXlQJKOUmlFvmwftvZuNxz7MnuRCRSBITb+CEn7sT8+znMC505a5RPxWl+tCBYaXc+FrvvrhAPna2+42uE48Qt6IQYwqIXZ5PzMhJ1pRdRVFqJNoTUI6K33//nuzs5zHGS3b2TOLjB7Kr3UfsngDdJkHWRYbj338L3nhb/TAVpQajPQGl3OTlLWHdulvswV8vxhxh1653gHxyU2DbhdBujhfXTf+nCkBRajiqBJRyk509xy+xvJO4FdBqPhwaO8KamuuL2qYoSo1ElYAChB9bPy9vCdnZsyie4OWyM4pF0HSlmz9MbsDeGXew/eYk9s+c4B++U1GUGoeOCShl+v07CZwDkJg4ioSE4eTmZnL8Z1s4NLsrP8aOw7vhCJtio+g5cyoxS5eCx1PqBDNFUaoHVQJKqX7/gQS6hzqjfu64JpNDh1bjzSo+Vs4p+yg8O5Xsn8aQnT0LYwo0iYui1CBUCShhh2P2teQ7dZrGvn0r/NY78wyLRBZNFHO74+1th/CZkMpSNIqiVB2qBJSwZu0GVvQgGFPA9u2zadnyGkeeYUhMvIEGDdoSF5fq6GX4xhBEw0MoSg1ClYAClD1r1z+hvNdea+x1lDAR+fZxu+OLtolEkJAwstyJYxRFOXaoElCKKG3g1mkycvYEfJW+b3DY18J3DjR36jSN/PwcHRBWlBqIKgEFKNtDKNBkBJRQGL7vwGxj+fk5Gh5aUWooqgQUIDwPoUCTUbgeRGr/V5SaiyoBBajcilvDQytK7UGVgAIUV9zZ2XMq7Xha+StKzadCYSNEpJmILBSRdfZ30xDlCkVkpf2ZX5FzKseW7dtnk5U1g1Wrzi0zhISiKLWfisYOuhtYbIxJAhbby8E4aIxJtj8XVvCcyjHCf1zgEJs3Tw07nlA45RRFqXlU1Bx0EZBq/54NZAJ3VfCYSjURF5dqp4UsBAw5Oe+QkzMfETcJCSOC+veXJ+6Qoig1j4r2BFoaY7Ls39lAyxDlGojIMhH5SkQuDnUwERlll1u2c+fOCoqmOAnVWneuj43tS0LCSEAcJbwYc5isrGeDmoiCeRUpilJ7KLMnICKLgIQgm/7hXDDGGBExQcoBtDPGbBWRjsCnIvK9MeaXwELGmOeA5wB69+4d6lhKOQnVWg+2PiFhONu3z8brPQx4HUcxQV1H1R1UUWo3ZSoBY0y/UNtEZLuIJBpjskQkEdgR4hhb7e9fRSQTSAFKKAHl2BBqDkCw9e3a3UOnTtPYufMtYmKSKSjYS3b2TIwpDFrJqzuootRuKjomMB+4BnjI/n43sIDtMXTAGHNYRJoDfwSmVvC8SjkI1VoPtj4vbwnr19+O13uEvLzPi3oHpVXy6g6qKLWXiiqBh4B0EbkO2ASkAYhIb2C0MeZ6oCvwrIh4scYgHjLG/FjB8yrloLTWerNm53H48DYSE68jNrZviZAPvt6BVvKKUjepkBIwxuQA5wZZvwy43v79JXByRc6jVJzA1npe3hJWrvRgzGEA9u9fSePGJ6uNX1HqGTpjuJ5ipYk8UrRsTD7Z2XNo0KCtRv1UlHqEKoF6ijUnIKqoJwARYaV/1DzBilK3UCVQT4mN7UtycoZfrKCsrBmUFkVUJ4YpSt1DlUA9xjlOkJe3xJ4fEHosoDwJ6RVFqR2oEqinBDPrtGx5DUDI9I86aKwodQ9VAvWQwKTx8fEDycn5oGhCmC9HcCA6MUxR6h6qBOoh/knjC9m1652ibWWZeXRimKLULSoaQE6phfjMOv6B4ixEBLc7vuqFUhSlWlAlUA/xmXWaNOkTsEUwxsv69bdrbgBFqSeoEqinxMb2pVOnaYj4egQR9rdXQ0IrSj1CxwTqMdZcgUxyczNxu+OLAsep54+i1B9UCdQxyjuj1znQ27jxyer5oyj1DFUCdYht255j3bpbbFfP6HLP6FXPH0Wpf+iYQB0hL28J69bdjDH5WHb9w0dl19ek8YpSv9CeQB3BigpanA5SJKLcdn2NDaQo9Q/tCdQRLN//aMCFSCRJSU+FjAIaqqWvSeMVpf6hPYE6QqiQDs6BYsCvpR+YN0BjAylK/UOVQB0hmFeQf4ygSGJieuD1HsY3ZmANInv9TD8aG0hR6heqBOoAoWz5gTGC9u1bChgsk5ELYwpxTg7zeQdp5a8o9QcdE6gDhLLll4wRZCmApk37kZT0tD2GEKGmH0Wpx2hPoA4QzJbvMw916jSN3bs/ZNeu9wCDyxVN+/YTiY3tq5PDFEVRJVAXCLTlA35jAVYPwCDiolOnaUUVvpp+FEVRJVBHcFbomzZNcYwF+OYOGIwR8vNzqk1GRVFqHjomUAcpHguIQCQKETdq+1cUJRjaE6iDBDMPqe1fUZRgqBKoozjNQxoHSFGUUKgSqONoPCBFUUpDxwTqOBoPSFGU0lAlUMdxDhLrwLCiKIGoOaiOo/GAFEUpDVUC9QCdFKYoSijUHKQoilKPUSWgKIpSj6mQEhCRy0VktYh4RaR3KeXOF5GfRGS9iNxdkXMqxUydChkZ/usyMqz1iqIo4VDRnsAPwKXAZ6EKiEgE8DQwAOgGDBWRbhU8rwL06QNpacWKICPDWu7Tp3rlUhSl9lChgWFjzBoAESmt2KnAemPMr3bZ14GLgB8rcm4FPB5IT7cq/jFjYPp0a9njqW7JFEWpLVTFmMAJwG+O5S32uhKIyCgRWSYiy3bu3FkFdI13KgAABjVJREFUotV+PB5LAUyebH2rAlAUpTyUqQREZJGI/BDkc1FlC2OMec4Y09sY07tFixaVffg6R17eEl577RX+8598xo+3egKBYwSKoiilUaY5yBjTr4Ln2Aq0cSy3ttcpFSAvbwmzZk1i4sSXmDhxECNGTMDj6UtampqEFEUJn6owBy0FkkSkg4hEAUOA+VVw3jpNbm4ma9YkM2FCGsnJi8jNzSwaI1i6tLqlUxSltlChgWERuQR4EmgBvC8iK40x54lIK+B5Y8xAY0yBiNwCfAxEADONMasrLHk9Jy4ulauuOtcvrzBYPQDtBSiKEi5ijKluGYLSu3dvs2zZsuoWo0bjSyavMYEURfEhIt8aY0LO2wpEYwfVMgIrfq38FUWpCKoEahGaIEZRlMpGYwfVIjRBjKIolY0qgVqEJohRFKWyUXNQLUITxCiKUtmoEqhl6GCwoiiViZqDFEVR6jGqBBRFUeoxqgQURVHqMaoEFEVR6jGqBBRFUeoxqgQURVHqMTU2gJyI7AQ2VdPpmwO7qunc5UVlrXxqi5ygsh4LaoucEFzWdsaYsLNy1VglUJ2IyLLyROGrTlTWyqe2yAkq67GgtsgJlSOrmoMURVHqMaoEFEVR6jGqBILzXHULUA5U1sqntsgJKuuxoLbICZUgq44JKIqi1GO0J6AoilKPUSWgKIpSj1ElAIjI5SKyWkS8IhLS3UpENorI9yKyUkSWVaWMDhnClfV8EflJRNaLyN1VKaNDhmYislBE1tnfTUOUK7Tv6UoRmV+F8pV6j0QkWkTm2tu/FpH2VSVbEFnKkvVaEdnpuI/XV5OcM0Vkh4j8EGK7iMgT9nV8JyI9q1pGhyxlyZoqInmOe3pfVctoy9FGRDJE5Ef73b8tSJmjv6/GmHr/AboCXYBMoHcp5TYCzWu6rEAE8AvQEYgCVgHdqkHWqcDd9u+7gX+FKLe/GmQr8x4BNwHP2L+HAHOr6ZmHI+u1wFPVIV+AHGcDPYEfQmwfCHwICHA68HUNljUVWFAD7mki0NP+3QT4OcjzP+r7qj0BwBizxhjzU3XLEQ5hynoqsN4Y86sx5gjwOnDRsZeuBBcBs+3fs4GLq0GGUIRzj5zyvwmcKyJShTL6qCnPs0yMMZ8Bu0spchEwx1h8BcSJSGLVSOdPGLLWCIwxWcaY5fbvfcAa4ISAYkd9X1UJlA8DfCIi34rIqOoWphROAH5zLG+h5J+mKmhpjMmyf2cDLUOUayAiy0TkKxGpKkURzj0qKmOMKQDygPgqkS6EHDahnudg2xTwpoi0qRrRyk1N+W+GS18RWSUiH4pI9+oWxjZJpgBfB2w66vtab9JLisgiICHIpn8YY94N8zBnGmO2isjxwEIRWWu3JiqVSpK1SihNVueCMcaISCh/5Hb2fe0IfCoi3xtjfqlsWes47wGvGWMOi8iNWD2Yc6pZptrOcqz/5n4RGQi8AyRVlzAiEgO8BdxujNlbWcetN0rAGNOvEo6x1f7eISLzsLrpla4EKkHWrYCzJdjaXlfplCariGwXkURjTJbdNd0R4hi++/qriGRitXSOtRII5x75ymwRkUggFsg5xnIFo0xZjTFOuZ7HGo+piVTZf7OiOCtaY8wHIvIfEWlujKny4HIi4sZSAK8YY94OUuSo76uag8JERBqLSBPfb+DPQFCvghrAUiBJRDqISBTWoGaVed04mA9cY/++BijRixGRpiISbf9uDvwR+LEKZAvnHjnlvwz41NijcFVMmbIG2H8vxLIb10TmA8Ntb5bTgTyHybBGISIJvjEgETkVq76s8kaALcMLwBpjzGMhih39fa3uke+a8AEuwbKhHQa2Ax/b61sBH9i/O2J5ZawCVmOZZmqkrKbYW+BnrBZ1dckaDywG1gGLgGb2+t7A8/bvM4Dv7fv6PXBdFcpX4h4B9wMX2r8bAG8A64FvgI7V+B8tS9Yp9v9yFZABnFRNcr4GZAH59v/0OmA0/H87d2wCIAxEAfS7ljiBk1i4jjtZii5kYQorK0Hh3oNUaS7XfLiEZGr7XZKlnePIw2u8H9Q633q6Juk/qnPIdR+5J9naGt/qq28jAAozDgIoTAgAFCYEAAoTAgCFCQGAwoQAQGFCAKCwEx+VF+7KvfJ2AAAAAElFTkSuQmCC\n","text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"tags":[],"needs_background":"light"}}]},{"cell_type":"code","metadata":{"id":"LdSNNeBSUhjY","executionInfo":{"status":"ok","timestamp":1623426592866,"user_tz":180,"elapsed":144,"user":{"displayName":"Nicolás Chehebar","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GjWN7wrOddIeQwbW6mzML-_Z0jTJLZ1cnE5tDuipw=s64","userId":"03611693109350494416"}}},"source":[""],"execution_count":34,"outputs":[]},{"cell_type":"code","metadata":{"id":"DYr4wB8v5v8-","executionInfo":{"status":"ok","timestamp":1623426592873,"user_tz":180,"elapsed":120,"user":{"displayName":"Nicolás Chehebar","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GjWN7wrOddIeQwbW6mzML-_Z0jTJLZ1cnE5tDuipw=s64","userId":"03611693109350494416"}}},"source":[""],"execution_count":34,"outputs":[]}]}