def surf(z,x=None,y=None,win=None,shade=0,edges=1,edge_color='fg',phi=-45.0,
         theta=30.0,zscale=1.0,palette=None,gnomon=0):
  '''Plot a three-dimensional wire-frame (surface): z=f(x,y)
  '''
  if win is None:
    pl3d.window3()
  else:
    pl3d.window3(win)
  pl3d.set_draw3_(0)
  phi0 = phi*numpy.pi/180.0
  theta0 = theta*numpy.pi/180.0
  pl3d.orient3(phi=phi0,theta=theta0)
  pl3d.light3()
  _change_palette(palette)
  sz = numpy.shape(z)
  if len(sz) != 2:
    raise ValueError('Input must be a 2-d array --- a surface.')
  N,M = sz
  if x is None:
    x = numpy.arange(0,N)
  if y is None:
    y = numpy.arange(0,M)
  x = numpy.squeeze(x)
  y = numpy.squeeze(y)
  if (len(numpy.shape(x)) == 1):
    x = x[:,newaxis]*numpy.ones((1,M))
  if (len(numpy.shape(y)) == 1):
    y = numpy.ones((N,1))*y[newaxis,:]
  plwf.plwf(z,y,x,shade=shade,edges=edges,ecolor=edge_color,scale=zscale)
  lims = pl3d.draw3(1)
  gist.limits(lims[0],lims[1],lims[2],lims[3])
  pl3d.gnomon(gnomon)

def plotForce():
    figure(size=3,aspect=0.5)
    subplot(1,2,1)
    from EvalTraj import plotFF
    plotFF(vp=351,t=28,f=900,cm=0.6,foffset=8)
    subplot_annotate()
    
    subplot(1,2,2)
    for i in [1,2,3,4]:
        R=np.squeeze(np.load('Rdpse%d.npy'%i))
        R=stats.nanmedian(R,axis=2)[:,1:,:]
        dps=np.linspace(-1,1,201)[1:]
        plt.plot(dps,R[:,:,2].mean(0));
    plt.legend([0,0.1,0.2,0.3],loc=3) 
    i=2
    R=np.squeeze(np.load('Rdpse%d.npy'%i))
    R=stats.nanmedian(R,axis=2)[:,1:,:]
    mn=np.argmin(R,axis=1)
    y=np.random.randn(mn.shape[0])*0.00002+0.0438
    plt.plot(np.sort(dps[mn[:,2]]),y,'+',mew=1,ms=6,mec=[ 0.39  ,  0.76,  0.64])
    plt.xlabel('Displacement of Force Origin')
    plt.ylabel('Average Net Force Magnitude')
    hh=dps[mn[:,2]]
    err=np.std(hh)/np.sqrt(hh.shape[0])*stats.t.ppf(0.975,hh.shape[0])
    err2=np.std(hh)/np.sqrt(hh.shape[0])*stats.t.ppf(0.75,hh.shape[0])
    m=np.mean(hh)
    print m, m-err,m+err
    np.save('force',[m, m-err,m+err,m-err2,m+err2])
    plt.xlim([-0.5,0.5])
    plt.ylim([0.0435,0.046])
    plt.grid(b=True,axis='x')
    subplot_annotate()

def log_diff_exp(x, axis=0):
    """ Calculates the logarithm of the diffs of e to the power of input 'x'. The method tries to avoid
        overflows by using the relationship: log(diff(exp(x))) = alpha + log(diff(exp(x-alpha))).
        
    :Parameter:
        x:    data.
             -type: float or numpy array 
          
        axis: Sums along the given axis.
             -type: int
        
    :Return:
        Logarithm of the sum of exp of x. 
       -type: float or numpy array.
        
    """
    alpha = x.max(axis) - numx.log(numx.finfo(numx.float64).max)/2.0
    if axis == 1:
        return numx.squeeze(alpha + numx.log(
                                             numx.diff(
                                                       numx.exp(x.T - alpha)
                                                       , n=1, axis=0)))
    else:
        return numx.squeeze(alpha + numx.log(
                                             numx.diff(
                                                       numx.exp(x - alpha)
                                                       , n=1, axis=0)))

def _field_gradient_jac(ref, target):
    """
    Given a reference field ref and a target field target
    compute the jacobian of the target with respect to ref

    Parameters
    ----------
    ref: Field instance that yields the topology of the space
    target array of shape(ref.V,dim)

    Results
    -------
    fgj: array of shape (ref.V) that gives the jacobian
         implied by the ref.field->target transformation.
    """
    import numpy.linalg as nl
    n = ref.V
    xyz = ref.field
    dim = xyz.shape[1]
    fgj = []
    ln = ref.list_of_neighbors()
    for i in range(n):
        j = ln[i]
        if np.size(j) > dim - 1:
            dx = np.squeeze(xyz[j] - xyz[i])
            df = np.squeeze(target[j] - target[i])
            FG = np.dot(nl.pinv(dx), df)
            fgj.append(nl.det(FG))
        else:
            fgj.append(1)

    fgj = np.array(fgj)
    return fgj

    def downsize(self, coefs, cut=None, verbose=True):
        """
        Given a set of coefs, sort the coefs and get rid of the bottom cut
        percent of variables with lowest cut coefs. Return the new coefs.
        """


        downsized_coefs = np.squeeze(np.array(coefs))

        if cut is None:
            cut = self.cut

        n_trash = int(floor(cut * self.n_features))

        if verbose:
            print("Downsampling...")
            print("Current shape:", self.Xview.shape)
            print("Removing {} columns... ".format(n_trash))


        self.tail_start -= n_trash

        if self.tail_start <= 0:
            raise ValueError("Trying to downsize more variables than present")

        # get sorted order of coefs
        csort = np.squeeze(np.argsort(np.argsort(np.absolute(coefs))))
        keep_feature = np.squeeze(csort >= n_trash)

        tail_start = self.tail_start

        # columns in the tail we want to keep
        keep_idx = np.squeeze(
            np.where(keep_feature[tail_start:tail_start+n_trash]))
        keep_idx += tail_start

        # columns we want to move to the tail
        trash_idx = np.squeeze(np.where(keep_feature[0:tail_start] == False))
        if len(trash_idx) != len(keep_idx):
            raise ValueError("trash_idx and keep_idx not the same length")

        # swap the columns
        for trash, keep in zip(trash_idx, keep_idx):
            #print(keep, trash)
            keep_col = self.X[:, keep].copy()
            self.X[:, keep] = self.X[:, trash]
            self.X[:, trash] = keep_col
            self.orig_feature_index[trash], self.orig_feature_index[keep] = self.orig_feature_index[keep], self.orig_feature_index[trash]
            downsized_coefs[trash], downsized_coefs[keep] = downsized_coefs[keep], downsized_coefs[trash]
            if self.test_subj is not None:
                self.X_test[:, (trash, keep)] = self.X_test[:, (keep, trash)]

        self.n_features -= n_trash
        self.Xview = self.X.view()[:, :self.n_features]
        if self.test_subj is not None:
            self.X_testview = self.X_test.view()[:, :self.n_features]

        print("New Xview shape:", self.Xview.shape)

        return downsized_coefs[:-n_trash]

 def visualize_ns_old(self, term, points=200):
     """
     Use randomly selected coordinates instead of most active
     """
     if term in self.no.term:
         term_index = self.no._ns['features_df'].columns.get_loc(term)
         rand_point_inds = np.random.random_integers(0, len(np.squeeze(zip(self.no._ns['mni_coords'].data))), points)
         rand_points = np.squeeze(zip(self.no._ns['mni_coords'].data))[rand_point_inds]
         weights = []
         inds_of_real_points_with_no_fucking_missing_study_ids = []
         for rand_point in range(len(rand_points)):
             if len(self.no.coord_to_ns_act(rand_points[rand_point].astype(list))) > 0:
                 inds_of_real_points_with_no_fucking_missing_study_ids.append(rand_point_inds[rand_point])
                 weights.append(self.no.coord_to_ns_act(rand_points[rand_point].astype(list))[term_index])
         fig = plt.figure()
         ax = fig.add_subplot(111, projection='3d')
         colors = cm.jet(weights/max(weights))
         color_map = cm.ScalarMappable(cmap=cm.jet)
         color_map.set_array(weights)
         fig.colorbar(color_map)
         x = self.no._ns['mni_coords'].data[inds_of_real_points_with_no_fucking_missing_study_ids, 0]
         y = self.no._ns['mni_coords'].data[inds_of_real_points_with_no_fucking_missing_study_ids, 1]
         z = self.no._ns['mni_coords'].data[inds_of_real_points_with_no_fucking_missing_study_ids, 2]
     else:
         raise ValueError('Term '+term + ' has not been initialized. '
                                         'Use get_ns_act(' + term + ')')
     ax.scatter(x, y, z, c=colors, alpha=0.4)
     ax.set_title('Estimation of ' + term)

 def cos_distance(self, strike1, dip1, strike2, dip2):
     """Angular distance betwen the poles of two planes."""
     xyz1 = sph2cart(*mplstereonet.pole(strike1, dip1))
     xyz2 = sph2cart(*mplstereonet.pole(strike2, dip2))
     r1, r2 = np.linalg.norm(xyz1), np.linalg.norm(xyz2)
     dot = np.dot(np.squeeze(xyz1), np.squeeze(xyz2)) / r1 / r2
     return np.abs(np.degrees(np.arccos(dot)))

    def __init__(self, timber_variable_bbq, beam=0):

        if not (beam == 1 or beam == 2):
            raise ValueError('You need to specify which beam! (1 or 2)')
        

        if type(timber_variable_bbq) is dict:
            dict_timber = timber_variable_bbq
        
        self.beam = beam
        
        self.amp_1 = np.squeeze(np.array(
            dict_timber['LHC.BQBBQ.CONTINUOUS_HS.B{:d}:EIGEN_AMPL_1'.format(beam)][1]))
        self.amp_2  = np.squeeze(np.array(
            dict_timber['LHC.BQBBQ.CONTINUOUS_HS.B{:d}:EIGEN_AMPL_2'.format(beam)][1]))
        
        self.xamp_1 = np.squeeze(np.array(
            dict_timber['LHC.BQBBQ.CONTINUOUS_HS.B{:d}:EIGEN_X_AMPL_1'.format(beam)][1]))
        self.xamp_2 = np.squeeze(np.array(
            dict_timber['LHC.BQBBQ.CONTINUOUS_HS.B{:d}:EIGEN_X_AMPL_2'.format(beam)][1]))
        
        self.qh  = dict_timber['LHC.BQBBQ.CONTINUOUS_HS.B{:d}:TUNE_H'.format(beam)][1]
        self.qv  = dict_timber['LHC.BQBBQ.CONTINUOUS_HS.B{:d}:TUNE_V'.format(beam)][1]
        
        self.q1  = dict_timber['LHC.BQBBQ.CONTINUOUS_HS.B{:d}:EIGEN_FREQ_1'.format(beam)][1]
        self.q2  = dict_timber['LHC.BQBBQ.CONTINUOUS_HS.B{:d}:EIGEN_FREQ_2'.format(beam)][1]
        
        self.t_stamps = np.ravel(np.squeeze(np.array(
            dict_timber['LHC.BQBBQ.CONTINUOUS_HS.B{:d}:EIGEN_AMPL_1'.format(beam)][0])))
        
        self.t_str=[datetime.datetime.fromtimestamp(self.t_stamps[ii]) for ii in np.arange(len(self.t_stamps))]

 def action_value(self, obs):
     # executes call() under the hood
     logits, value = self.predict(obs)
     action = self.dist.predict(logits)
     # a simpler option, will become clear later why we don't use it
     # action = tf.random.categorical(logits, 1)
     return np.squeeze(action, axis=-1), np.squeeze(value, axis=-1)

    def __init__(self, complete_path):

        if complete_path.endswith('.mat.gz'):
            temp_filename = complete_path.split('.gz')[0]
            with open(temp_filename, "wb") as tmp:
                shutil.copyfileobj(gzip.open(complete_path), tmp)
            dict_mr = sio.loadmat(temp_filename)
            os.remove(temp_filename)
        elif complete_path.endswith('.mat'):
            dict_mr = sio.loadmat(complete_path)
        else:
            print('Unknown file extension for MountainRange file. Should be ' +
                  '.mat or .mat.gz')
        self.value = dict_mr['value']
        self.trigger_stamp = dict_mr['triggerStamp']
        self.SC_numb = np.int(np.squeeze(dict_mr['superCycleNb']))
        self.first_trigger_t_stamp_unix = dict_mr['first_trigger_t_stamp_unix']
        self.sample_interval = float(np.squeeze(dict_mr['sampleInterval']))
        self.first_sample_time = dict_mr['firstSampleTime']
        self.sensitivity = dict_mr['sensitivity']
        self.offset = dict_mr['offset']
        self.SPSuser = dict_mr['SPSuser']
        self.t_stamp_unix = dict_mr['t_stamp_unix']

        self.time_axis = np.float_(range(self.value.shape[1]))*self.sample_interval-self.value.shape[1]*self.sample_interval/2.

def kalman(x, u, P, A, B, C, W, V, z=np.NaN):
    """
    This function returns an optimal expected value of the state and covariance
    error matrix given an update and system parameters.

    x:   Estimate of state at time t-1.
    u:   Input at time t-1.
    P:   Estimate of error covariance matrix at time t-1.
    A:   Discrete time state tranistion matrix at time t-1.
    B:   Input to state model matrix at time t-1.
    C:   Observation model matrix at time t.
    W:   Process noise covariance at time t-1.
    V:   Measurement noise covariance at time t.
    z:   Measurement at time t.

    returns: (x,P) tuple
    x: Updated estimate of state at time t.
    P: Updated estimate of error covariance matrix at time t.

    """

    x = np.atleast_2d(x)
    u = np.atleast_2d(u)
    P = np.atleast_2d(P)
    A = np.atleast_2d(A)
    B = np.atleast_2d(B)
    x_p = np.dot(A, x) + np.dot(B, u)  # Prediction of estimated state vector
    P_p = np.dot(A, np.dot(P, A.T)) + W  # Prediction of error covariance matrix

    if np.any(np.isnan(z)):
        return (x_p, P_p)
    else:
        C = np.atleast_2d(C)
        W = np.atleast_2d(W)
        V = np.atleast_2d(V)
        z = np.atleast_2d(z)

        [M, N] = np.shape(C)

        if W.shape[0] == 1 or W.shape[1] == 1:
            W = np.diag(np.squeeze(W))

        if (V.shape[0] == 1 or V.shape[1] == 1) and not (V.shape[0] == 1 and V.shape[1] == 1):
            V = np.diag(np.squeeze(V))

        I = np.eye(N)  # N x N identity matrix

        S = np.dot(C, np.dot(P_p, C.T)) + V  # Sum of error variances
        S_inv = np.linalg.inv(S)  # Inverse of sum of error variances
        K = np.dot(P_p, np.dot(C.T, S_inv))  # Kalman gain
        r = z - np.dot(C, x_p)  # Prediction residual
        w = np.dot(-K, r)  # Process error
        x = x_p - w  # Update estimated state vector
        # v = z - np.dot(C, x)  # Measurement error
        if np.any(np.isnan(np.dot(K, V))):
            P = P_p
        else:
            # Updated error covariance matrix
            P = np.dot((I - np.dot(K, C)), np.dot(P_p, (I - np.dot(K, C)).T)) + np.dot(K, np.dot(V, K.T))
        return (x, P)

  def testTrainNetwork(self, distribution, optimizer_fn,
                       use_callable_loss=True):
    with distribution.scope():
      model_fn, dataset_fn, layer = minimize_loss_example(
          optimizer_fn, use_bias=True, use_callable_loss=use_callable_loss)
      iterator = distribution.make_input_fn_iterator(lambda _: dataset_fn())

      def run_step():
        return control_flow_ops.group(
            distribution.experimental_local_results(
                distribution.extended.call_for_each_replica(
                    model_fn, args=(iterator.get_next(),))))

      if not context.executing_eagerly():
        with self.cached_session() as sess:
          sess.run(iterator.initialize())
          run_step = sess.make_callable(run_step())
        self.evaluate(variables.global_variables_initializer())

      weights, biases = [], []
      for _ in range(10):
        run_step()

        weights.append(self.evaluate(layer.kernel))
        biases.append(self.evaluate(layer.bias))

      error = abs(numpy.add(numpy.squeeze(weights), numpy.squeeze(biases)) - 1)
      is_not_increasing = all(y <= x for x, y in zip(error, error[1:]))
      self.assertTrue(is_not_increasing)

def gray2rgb(image):
    """Create an RGB representation of a gray-level image.

    Parameters
    ----------
    image : array_like
        Input image of shape ``(M, N [, P])``.

    Returns
    -------
    rgb : ndarray
        RGB image of shape ``(M, N, [, P], 3)``.

    Raises
    ------
    ValueError
        If the input is not a 2- or 3-dimensional image.

    """
    if np.squeeze(image).ndim == 3 and image.shape[2] in (3, 4):
        return image
    elif image.ndim != 1 and np.squeeze(image).ndim in (1, 2, 3):
        image = image[..., np.newaxis]
        return np.concatenate(3 * (image,), axis=-1)
    else:
        raise ValueError("Input image expected to be RGB, RGBA or gray.")

def plsurf(z,x=None,y=None,win=None,shade=0,edges=1,edge_color='fg',phi=-45.0,
         theta=30.0,zscale=1.0,palette=None,gnomon=0,animate=False,limits=True, ireg=None):
  '''Plot a 3-D wire-frame surface z=f(x,y)
  '''
  if win is None:
    pass
    #pl3d.window3()
  else:
    pl3d.window3(win)
  pl3d.set_draw3_(0)
  phi0 = phi*numpy.pi/180.0
  theta0 = theta*numpy.pi/180.0
  pl3d.orient3(phi=phi0,theta=theta0)
  pl3d.light3()
  _change_palette(palette)
  sz = numpy.shape(z)
  if len(sz) != 2:
    raise ValueError('Input must be a 2-d array --- a surface.')
  N,M = sz
  if x is None:
    x = numpy.arange(0,N)
  if y is None:
    y = numpy.arange(0,M)
  x = numpy.squeeze(x)
  y = numpy.squeeze(y)
  if (len(numpy.shape(x)) == 1):
    x = x[:,newaxis]*numpy.ones((1,M))
  if (len(numpy.shape(y)) == 1):
    y = numpy.ones((N,1))*y[newaxis,:]
  plwf.plwf(z,y,x,shade=shade,edges=edges,ecolor=edge_color,scale=zscale, ireg=ireg)
  # if animate, the application is responsible to fma
  lims = pl3d.draw3(not animate)
  if limits:
    gist.limits(lims[0],lims[1],lims[2],lims[3])
  pl3d.gnomon(gnomon)

def generate_ic_grid(dR=0.1*u.kpc, dRdot=5.*u.km/u.s):
    # spacing between IC's in R and Rdot
    dR = dR.decompose(usys).value
    dRdot = dRdot.decompose(usys).value
    max_Rdot = (50*10*u.km/u.s).decompose(usys).value
    max_R = (15*u.kpc).decompose(usys).value

    # from the paper
    E = (600*100*(u.km/u.s)**2).decompose(usys).value
    Lz = (10.*10.*u.km*u.kpc/u.s).decompose(usys).value # typo in paper? km/kpc instead of km*kpc
    z = 0.
    params = oblate_params

    w0s = []
    for R in np.arange(0, max_R, dR):
        # zero velocity curve
        V = zotos_potential(R, z, *params)
        ZVC_Rdot = np.squeeze(np.sqrt(2*(E-V) - Lz**2/R**2))
        for Rdot in np.arange(0, max_Rdot, dRdot):
            if Rdot > ZVC_Rdot or R < 0.2 or R >= 13: continue
            zdot = np.squeeze(np.sqrt(2*(E - V) - Lz**2/R**2 - Rdot**2))
            w0 = [R,z,Rdot,zdot]
            w0s.append(w0)
    w0s = np.array(w0s)
    return w0s, Lz

def sdss_source_props_ota(img,ota):
    """
    Use photutils to get the elongation of all of the sdss sources
    can maybe use for point source filter
    Also fit a gaussian along a row and col of pixels passing
    through the center of the star
    """

    image = odi.reprojpath+'reproj_'+ota+'.'+str(img[16:])
    hdulist = odi.fits.open(image)
    data = hdulist[0].data

    sdss_source_file = odi.coordspath+'reproj_'+ota+'.'+str(img[16:-5])+'.sdssxy'

    x,y,ra,dec,g,g_err,r,r_err = np.loadtxt(sdss_source_file,usecols=(0,1,2,3,
                                                                      6,7,8,9),unpack=True)

    box_centers = zip(y,x)
    box_centers = np.reshape(box_centers,(len(box_centers),2))
    source_dict = {}
    total_fwhm = []
    for i,center in enumerate(box_centers):
        x1 = center[0]-50
        x2 = center[0]+50
        y1 = center[1]-50
        y2 = center[1]+50

        #print x1,x2,y1,y2,center
        box = data[x1:x2,y1:y2]
        col = data[x1:x2,int(center[1]-0.5):int(center[1]+0.5)]
        row = data[int(center[0]-0.5):int(center[0]+0.5),y1:y2]
        row = np.squeeze(row) - np.median(row)
        col = np.squeeze(col) - np.median(col)
        g_init = models.Gaussian1D(amplitude=250., mean=50, stddev=2.)
        fit_g = fitting.LevMarLSQFitter()
        pix = np.linspace(0,100,num=100)
        g_row = fit_g(g_init, pix, row)
        g_col = fit_g(g_init, pix, col)
        mean_fwhm = 0.5*(g_row.stddev*2.355+g_col.stddev*2.355)
        total_fwhm.append(mean_fwhm)
        #odi.plt.imshow(box)
        #odi.plt.plot(row)
        #odi.plt.plot(pix,g(pix))
        #plt.imshow(row2)
        #plt.show()
        mean, median, std = odi.sigma_clipped_stats(box, sigma=3.0)
        threshold = median + (std * 2.)
        segm_img = odi.detect_sources(box, threshold, npixels=20)
        source_props = odi.source_properties(box,segm_img)
        if len(source_props) > 0:
            columns = ['xcentroid', 'ycentroid','elongation','semimajor_axis_sigma','semiminor_axis_sigma']
            if i == 0:
                source_tbl = odi.properties_table(source_props,columns=columns)
            else:
                source_tbl.add_row((source_props[0].xcentroid,source_props[0].ycentroid,
                                    source_props[0].elongation,source_props[0].semimajor_axis_sigma,
                                    source_props[0].semiminor_axis_sigma))
    elong_med,elong_std = np.median(source_tbl['elongation']),np.std(source_tbl['elongation'])
    hdulist.close()
    return elong_med,elong_std,np.mean(total_fwhm),np.std(total_fwhm)

    def reparametrization_LS_assembler(self):
        """In this function we compute the arclength reparametrization by mean of a Least Square
        problem."""

        s_array = np.linspace(0, self.points_s[-1,1], self.arcfactor * self.n_dofs)
        #self.point_ls = list()
        self.point_ls = np.asmatrix(np.zeros([s_array.shape[0],self.dim + 2]))
        tval = np.zeros(s_array.shape[0])
        sval = np.linspace(0,1,s_array.shape[0])
        for i in range(0, s_array.shape[0]):
            tval[i] = self.find_s(s_array[i])
            #rint tval
            # The curve should have a value( or __call__ if prefered) method that we can query to know its value in space
        self.point_ls[:,0:self.dim] = np.squeeze(self.curve(tval).transpose())
        self.rhsinit = np.squeeze(self.curve(tval).transpose())
        #self.point_ls[:,self.dim] = s_array[:,np.newaxis]
        self.point_ls[:,self.dim] = tval[:,np.newaxis]
            # In point_ls we store the value in space, the s_array and the tval obtained
            #self.point_ls.append([cval, s_array[i], tval])

        #self.point_ls = np.array(self.point_ls)
        # We compute the number of elements in the system rectangular matrix(Bmatrix), it will have dim*s_array.shape[0] rows and dim*nknot columns.
        # We want it to be rectangular because we are approximating its resilution so we search for something that solves the reparametrization in a
        # least square sense.
        #Bmatrixnumelem = self.dim * s_array.shape[0] * self.n_dofs * self.dim
        #self.matrixB = np.zeros(Bmatrixnumelem).reshape(self.dim * s_array.shape[0], self.n_dofs * self.dim)
        #self.rhsinit = np.zeros(self.dim * s_array.shape[0])
        self.matrixB = interpolation_matrix(self.vector_space, sval)

def plot_animation_each_neuron(name_s, save_name, print_loss=False):
    """Plot the movie for all the networks in the information plane"""
    # If we want to print the loss function also
    #The bins that we extened the x axis of the accuracy each time
    epochs_bins = [0, 500, 1500, 3000, 6000, 10000, 20000]
    data_array = utils.get_data(name_s[0][0])
    data = np.squeeze(data_array['information'])

    f, (axes) = plt.subplots(1, 1)
    axes = [axes]
    f.subplots_adjust(left=0.14, bottom=0.1, right=.928, top=0.94, wspace=0.13, hspace=0.55)
    colors = LAYERS_COLORS
    #new/old version
    Ix = np.squeeze(data[0,:, :, :])
    Iy = np.squeeze(data[1,:, :, :])
    #Interploation of the samplings (because we don't cauclaute the infomration in each epoch)
    #interp_data_x = interp1d(epochsInds,  Ix, axis=1)
    #interp_data_y = interp1d(epochsInds,  Iy, axis=1)
    #new_x = np.arange(0,epochsInds[-1])
    #new_data  = np.array([interp_data_x(new_x), interp_data_y(new_x)])
    """"
    train_data = interp1d(epochsInds,  np.squeeze(train_data), axis=1)(new_x)
    test_data = interp1d(epochsInds,  np.squeeze(test_data), axis=1)(new_x)

    if print_loss:
        loss_train_data =  interp1d(epochsInds,  np.squeeze(loss_train_data), axis=1)(new_x)
        loss_test_data=interp1d(epochsInds,  np.squeeze(loss_test_data), axis=1)(new_x)
    """
    line_ani = animation.FuncAnimation(f, update_line_each_neuron, Ix.shape[1], repeat=False,
                                       interval=1, blit=False, fargs=(print_loss, Ix, axes,Iy,train_data,test_data,epochs_bins, loss_train_data,loss_test_data, colors,epochsInds))
    Writer = animation.writers['ffmpeg']
    writer = Writer(fps=100)
    #Save the movie
    line_ani.save(save_name+'_movie.mp4',writer=writer,dpi=250)
    plt.show()

def main():
    x = np.loadtxt(sys.argv[1], skiprows=1, delimiter=",")
    x = np.array(x)
    # separating the class 1 rows from class 2 rows
    indic = np.where(x[:,1] == 1)
    x_one = np.squeeze(x[indic,:])
    indic = np.where(x[:,1] == 2)
    x_two = np.squeeze(x[indic,:])
    
    m = np.loadtxt(sys.argv[2])
    m = np.array(m)
    var = np.loadtxt(sys.argv[3])
    var = np.array(var)
    w = np.loadtxt(sys.argv[4])
    w = np.array(w)
    its = sys.argv[5]
    its = np.int32(its)
    
    ll = gmmest(x_two[:,0], m, var, w, its)
    
    plt.plot(ll[3])
    plt.ylabel('log likelihood')
    plt.xlabel('iterations')
    plt.show()
    print "mu: ", ll[0], " sigmasq: ", ll[1], " wt: ", ll[2], " ll: ", ll[3]

 def ll2ij(self, lon, lat, mask=None, cutoff=None, nei=1, all_nei=False,
           return_dist=False):
     """Reproject a lat-lon vector to i-j grid coordinates"""
     self.add_ij()
     if mask is not None:
         self.add_kd(mask)
     elif not hasattr(self, 'kd'):
         self.add_kd()
     dist,ij = self.kd.query(list(np.vstack((lon,lat)).T), nei)
     #if cutoff is not None:
     #    ij[dist>cutoff] = 0
     if nei == 1 :
         ivec = self.kdijvec[ij[:] - 1][:, 0]
         jvec = self.kdijvec[ij[:] - 1][:, 1]
         if cutoff is not None:
             ivec[dist>cutoff] = -999
             jvec[dist>cutoff] = -999
     elif all_nei == False:
         ivec = np.squeeze(self.kdijvec[ij[:,:]-1])[:, nei-1, 0]
         jvec = np.squeeze(self.kdijvec[ij[:,:]-1])[:, nei-1, 1]
         dist = np.squeeze(dist[:, nei-1])
     else:
         ivec = np.squeeze(self.kdijvec[ij[:,:]-1])[:, :, 0]
         jvec = np.squeeze(self.kdijvec[ij[:,:]-1])[:, :, 1]
     if return_dist == True:
         return ivec,jvec,dist
     else:
         return ivec,jvec