def test_1d_gaussian_slope_error(self):

        mirror_length=200.0
        step=1.0
        rms_slopes=1.3e-7
        rms_heights = 1e-7
        correlation_length = 10.0

        x, f = simulate_profile_1D_gaussian(step=step, \
                                              mirror_length=mirror_length, \
                                              rms_heights=rms_heights, \
                                              correlation_length=correlation_length,\
                                              renormalize_to_slopes_sd=None)
        slopes = numpy.gradient(f, x[1]-x[0])
        print("test_1d_gaussian: test function: %s, Stdev (not normalized): HEIGHTS=%g.SLOPES=%g"%("test_1d_gaussian_slope_error",f.std(),slopes.std()))

        x, f = simulate_profile_1D_gaussian(step=step, \
                                              mirror_length=mirror_length, \
                                              rms_heights=rms_heights, \
                                              correlation_length=correlation_length,\
                                              renormalize_to_slopes_sd=rms_slopes)
        slopes = numpy.gradient(f, x[1]-x[0])

        print("test_1d_gaussian: test function: %s, SLOPES Stdev (normalized to %g)=%g"%("test_1d_gaussian_slope_error",rms_slopes,slopes.std()))
        assert numpy.abs( rms_slopes - slopes.std() ) < 0.01 * numpy.abs(rms_slopes)

        if do_plot:
            from srxraylib.plot.gol import plot
            plot(x,slopes,title="test_1d_gaussian_slope_error",xtitle="Y",ytitle="slopes Z'")

def spline_do(x, y, ax, k=3, s=7, n=1, diff=True, plot=True, error=False):
    '''
    Finds  spline with degree k and smoothing factor s for x and y.
    Returns ys,d_ys: the y values of the spline and the derivative.

    '''
    # xs = np.array(x).astype(float)
    s = UnivariateSpline(x, y, k=k, s=s)
    xs = np.linspace(x.values.min(), x.values.max(), len(x))
    ys = s(xs)
    if plot:
        if not error:
            ax.plot(x, norm(y), '.', label='Data')
        ax.plot(xs, norm(ys), label='Spline fit')
        if error:
            ax.errorbar(x, norm(y), yerr=np.sqrt(norm(y)))

    d_ys = np.gradient(ys)
    dd_ys = np.gradient(d_ys)
    if diff:
        if plot:
            ax.plot(xs, norm(d_ys), '-.', label='1st derivative')
            #ax.plot(xs, norm(dd_ys),'--',label='2nd derivative')
            ax.plot(xs, norm(smooth(dd_ys, 81)),
                    '--', label='2nd derivative (smothed)')
    if plot:
        ax.legend(loc=0)
    return xs, ys, d_ys, dd_ys

def get_beta_deff(savepath, filename):
  filepath = os.path.join(savepath,filename)
  with open(filepath,'r') as f1:
    #col1 = [] #typically mean x-position
    #col2 = [] #typically mean (x-position)^2
    col3 = [] #typically MSD-x as calculated from the above quantities
    for line in f1:
      ls = line.split()
      #Analytical data output is currently set at 3 cols.
      #Make sure you are using the correct data. 
      #col1 += [float(ls[0])]
      #col2 += [float(ls[1])]
      col3 += [float(ls[2])]

  times = np.array(range(1, len(col3) + 1)) #create an array of times
  msd_x = np.array(col3) #create an array of msd_x

  #For effective diffusion:
  deff  = np.true_divide(msd_x, 2*times)

  #For beta(t) plot.
  log_times = np.log10(times)
  log_msd_x = np.log10(msd_x)
  dt = np.gradient(log_times)
  dd = np.gradient(log_msd_x, dt)

  return times, deff, log_times, dd, msd_x

def vorticity(x, y, u, v, coord_type='geographic'):
    """
    USAGE
    -----
    zeta = vorticity(x, y, u, v, coord_type='geographic')

    Calculates the vertical component 'zeta' (dv/dx - du/dy, in 1/s) of the
    relative vorticity vector from the 'u' and 'v' velocity arrays (in m/s)
    specified in spherical coordinates by the 'lon' and 'lat' 2D meshgrid-type
    arrays (in degrees).
    """
    x, y, u, v = map(np.array, (x, y, u, v))

    if coord_type=='geographic':
        dx, dy = deg2m_dist(lon, lat)
    elif coord_type=='cartesian':
        dy, _ = np.gradient(y)
        _, dx = np.gradient(x)
    elif coord_type=='dxdy':
        dx, dy = x, y
        pass

    duy, _ = np.gradient(u)
    _, dvx = np.gradient(v)

    dvdx = dvx/dx
    dudy = duy/dy
    vrt = dvdx - dudy # [1/s]

    return vrt

 def airfoil(self):
     """Biconvex airfoil"""
     L = self.L
     N = self.N
     R = 210
     theta_arc = np.arcsin(L/2 / R) * 2
     pi = np.pi
     theta_up = np.linspace((pi+theta_arc) / 2,
                            +(pi-theta_arc) / 2, N)
     #theta_up = theta[::-1]
     theta_down =np.linspace((3*pi-theta_arc) / 2,
                            +(3*pi+theta_arc) / 2, N)
     
     X_up = R * np.cos(theta_up)
     Y_up = R * np.sin(theta_up)
     X_down = R * np.cos(theta_down)
     Y_down = R * np.sin(theta_down)
     
     shift_r = X_up[0]
     X_up -= shift_r
     X_down -= shift_r
     shift_up = Y_up[0]
     shift_down = Y_down[0]
     Y_up = Y_up - shift_up
     Y_down = Y_down - shift_down
     
     X = np.concatenate((X_up, X_down))
     Y = np.concatenate((Y_up, Y_down))
     slope_up = np.gradient(Y_up,1)/np.gradient(X_up,1)
     slope_down = np.gradient(Y_down,1)/np.gradient(X_down,1)
     angle = np.arctan(np.concatenate((slope_up, slope_down)))
     return X, Y, angle

def extractAllDescriptors(signal):
    """
    Extracts the descriptors expected for the analysis of a given audio file.
    """

    described = {}

    described['Silence'] = _silence = silence(signal)
    signal = signal[config.hopSize * _silence[0]:config.hopSize * _silence[1]] / np.max(
        signal)  # Tomo solo la parte con sonido y promedio para que todas las señales sean parejas.
    described['mfcc'] = mfccs(signal)
    described['Inharmonicity'] = inharmonicity_tesis(signal)
    described['Energy'] = energy(signal)
    described['LogAttackTime'] = log_attack_time(signal)
    described['Standard-Dev'] = standard_dev(signal)
    described['Variance'] = variance(signal)
    described['Skewness'] = skewness(signal)
    described['kurtosis'] = kurtosis(signal)

    # described['mfcc-1st'] = np.gradient(described['mfcc'])[1]
    # described['mfcc-2nd'] = np.gradient(described['mfcc-1st'])[1]
    described['Inharmonicity-1st'] = np.gradient(described['Inharmonicity'])
    described['Inharmonicity-2nd'] = np.gradient(described['Inharmonicity-1st'])

    described['mfcc-Std-f'], described['mfcc-Var-f'], described['mfcc-Skew-f'], described['mfcc-Kurt-f']\
        = mfcc_std_frequency(described)

    return described

def read_HRRS_data(ff):

	"""
	Read in a .dat file from SPARC high-res radiosonde data 
	Input ff is a string pointing to the full path of the desired file. 
	"""

	# here is a dict that gives bad values for different columns 
	# alert: this is still incomplete 
	badvals = {'Temp':['999.0'],'Alt':['99.0','99999.0'],'Lat':['999.000'],'Lon':['9999.000']}
	
	D= pd.read_csv(ff,skiprows=13,error_bad_lines=False,delim_whitespace=True,na_values=badvals)
	colnames=list(D.columns.values)

	# kick out the first two rows - they hold units and symbols 
	D.drop(D.index[[0,1]], inplace=True)

	# also make sure that lat, lon, pressure, altitude, and temp are numerics 
	vars_to_float = ['Press','Temp','Lat','Lon','Alt']
	D[vars_to_float] = D[vars_to_float].astype(float)

	# compute the vertical gradient of potential temp and, from that, buoyancy frequency 
	P0=1000.0
	Rd = 286.9968933                # Gas constant for dry air        J/degree/kg
	g = 9.80616                     # Acceleration due to gravity       m/s^2
	cp = 1005.0                     # heat capacity at constant pressure    m^2/s^2*K
	theta=(D['Temp']+273.15)*(P0/D['Press'])**(Rd/cp)		# note that this includes conversion of Celsius to Kelvin  
	dZ = np.gradient(D['Alt']) 
	dthetadZ = np.gradient(theta,dZ)
	D["N2"]=(g/theta)*dthetadZ
	

	return(D)

def exportGradientImage(sigma=3.0):
	x, y = 5000, 3500
	size = 2000
	topleft = (x, y)
	bottomright = (x+size, y+size)
	image = getImageByName(imagefilename='../resources/images/registered-to-2008-07-24-09_55.tif',
						   topleft=topleft,
						   bottomright=bottomright)
	smooth = vigra.filters.gaussianSmoothing(image, sigma)
	smoothswap = smooth.swapaxes(0, 1)
	m, n = vigra.Image((2000,2000)), vigra.Image((2000,2000))
	
	for i in range(size):
		grad = np.gradient(smooth[i])
		for j in range(len(grad)):
		    m[i][j] = grad[j]
		    
	for i in range(size):
		grad = np.gradient(smoothswap[i])
		for j in range(len(grad)):
		    n[j][i] = grad[j]
		    
	out = m + n
	vigra.impex.writeImage(vigra.colors.linearRangeMapping(out), '/home/max/Desktop/out.png')
	vigra.impex.writeImage(vigra.colors.linearRangeMapping(m), '/home/max/Desktop/m.png')
	vigra.impex.writeImage(vigra.colors.linearRangeMapping(n), '/home/max/Desktop/n.png')
	
	return smooth

def _threshold_gradient(im):
    """Indicate pixel locations with gradient below the bottom 10th percentile

    Parameters
    ----------
    im : array
        The mean intensity images for each channel.
        Size: (num_channels, num_rows, num_columns).

    Returns
    -------
    array
        Binary values indicating whether the magnitude of the gradient is below
        the 10th percentile.  Same size as im.

    """

    if im.shape[0] > 1:
        # Calculate directional relative derivatives
        _, g_x, g_y = np.gradient(np.log(im))
    else:
        # Calculate directional relative derivatives
        g_x, g_y = np.gradient(np.log(im[0]))
        g_x = g_x.reshape([1, g_x.shape[0], g_x.shape[1]])
        g_y = g_y.reshape([1, g_y.shape[0], g_y.shape[1]])
    gradient_magnitudes = np.sqrt((g_x ** 2) + (g_y ** 2))
    below_threshold = []
    for chan in gradient_magnitudes:
        threshold = mquantiles(chan[np.isfinite(chan)].flatten(), [0.1])[0]
        below_threshold.append(chan < threshold)
    return np.array(below_threshold)

def rect_guess(self, data, x=None, **kwargs):
    if x is None:
        return
    ymin, ymax = min(data), max(data)
    xmin, xmax = min(x), max(x)

    ntest = min(2, len(data)/5)
    step_up =  (data[:ntest].mean() > data[-ntest:].mean())


    dydx = savitzky_golay(np.gradient(data)/np.gradient(x), 5, 2)
    cen1 = x[np.where(dydx==dydx.max())][0]
    cen2 = x[np.where(dydx==dydx.min())][0]
    if step_up:
        center1 = cen1 # + (xmax+xmin)/4.0)/2.
        center2 = cen2 # + 3*(xmax+xmin)/4.0)/2.
    else:
        center1 = cen2 # + (xmax+xmin)/4.0)/2.0
        center2 = cen1 # + 3*(xmax+xmin)/4.0)/2.0

    pars = self.make_params(amplitude=(ymax-ymin),
                            center1=center1, center2=center2)

    pars['%ssigma1' % self.prefix].set(value=(xmax-xmin)/5.0, min=0.0)
    pars['%ssigma2' % self.prefix].set(value=(xmax-xmin)/5.0, min=0.0)
    return update_param_vals(pars, self.prefix, **kwargs)

def _optimize_num_overlap_pixels(data):
    """
    """
    num_projection, num_slices, num_pixels = data.shape
    if num_projection % 2 != 0: # if odd
        img_first_half = np.squeeze(data[1:num_projection/2 + 1, num_slices/2, :])
        img_second_half = np.squeeze(data[num_projection/2:num_projection - 1, num_slices/2, :])
    else:
        img_first_half = np.squeeze(data[1:num_projection/2 + 1, num_slices/2, :])
        img_second_half = np.squeeze(data[num_projection/2:num_projection, num_slices/2, :])
    ind = range(0, num_pixels)[::-1]
    img_second_half = img_second_half[:, ind]
    
    img_first_half = ndimage.filters.gaussian_filter(img_first_half, sigma=2)
    img_second_half = ndimage.filters.gaussian_filter(img_second_half, sigma=2)
    
    gx1, gy1 = np.gradient(img_first_half)
    gx2, gy2 = np.gradient(img_second_half)
    img_first_half = np.power(gx1, 2) + np.power(gy1, 2)
    img_second_half = np.power(gx2, 2) + np.power(gy2, 2)
    
    img1 = np.fft.fft(img_first_half)
    img2 = np.fft.fft(img_second_half)
    tmp = np.real(np.fft.ifft(np.multiply(np.conj(img2), img1)))
    return np.argmax(np.sum(np.abs(tmp), axis=0))

def propagate_tie(mu, delta, pixel_size, dist):
    """
    Propagate emitting x-ray wave based on Transport of Intensity.

    Parameters
    ----------
    mu : ndarray, optional
        3D tomographic data for attenuation.
    delta : ndarray
        3D tomographic data for refractive index.
    pixel_size : float
        Detector pixel size in cm.
    dist : float
        Propagation distance of the wavefront in cm.

    Returns
    -------
    ndarray
        3D propagated tomographic intensity.
    """
    i1 = np.exp(-mu)
    i2 = np.zeros(delta.shape)
    for m in range(delta.shape[0]):
        dx, dy = np.gradient(delta[m], pixel_size)
        d2x, _ = np.gradient(i1[m] * dx, pixel_size)
        _, d2y = np.gradient(i1[m] * dy, pixel_size)
        i2[m] = i1[m] + dist * (d2x + d2y)
    return i2

def deblurImage(blurImg,iteration,mode=0):
    niter=0
    edge=getEdge(blurImg*255,mode)
    stDesv=np.std(edge)
    grady, gradx=np.gradient((blurImg*255))
    deblurImg=np.copy(blurImg)
    normalizar=False
    #print "Gradiente antes", np.sum(gradx+grady) 
    desv=np.std(deblurImg)
    extraGain=1.0
    while(niter<iteration):
        desv=np.std(deblurImg)
        #Dprint "Desviacion estandar borrosa", desv
        for j in range(deblurImg.shape[0]-1):
            for k in range(deblurImg.shape[1]-1):
                gain=gradx[j,k]*stDesv+grady[j,k]*stDesv
                if(gain<extraGain):
                    extraGain=gain
                deblurImg[j,k]=deblurImg[j,k]+gain
                if(deblurImg[j,k]<0.0 or deblurImg[j,k]>255.0):
                    normalizar=True
        deblurImg=extraGain/10.0+deblurImg
        if normalizar:
            normalize(deblurImg)
        
        edge=getEdge(deblurImg,mode)
        stDesv=np.std(edge)
        niter=niter+1
    gradx2, grady2=np.gradient(deblurImg)
    return deblurImg

    def get_quantum_driving_parameters(self):
        """Return the adapted parameters (eps_prime, delta, theta_prime) to
        obtain adiabatic dynamics for arbitrary length.
        """
        eps, delta = self.get_cycle_parameters()
        eps_dot, delta_dot = [np.gradient(x, self.dt) for x in eps, delta]

        mixing_angle_dot = 2.*np.abs(self.B0)*(delta*eps_dot-delta_dot*eps)
        mixing_angle_dot /= (delta**2 + 4.*np.abs(self.B0)**2*eps**2)
        self.mixing_angle_dot = mixing_angle_dot

        self.mixing_angle = np.arctan(2.*np.abs(self.B0)*eps/delta)
        self.mixing_angle_dot_alt = np.gradient(self.mixing_angle, self.dt)

        theta_prime = -2.*np.arctan2(mixing_angle_dot, (2*np.abs(self.B0)*eps))

        B_prime = (-1j * (np.exp(1j*theta_prime) + 1.) * np.pi**2 /
                   self.W**3 / np.sqrt(self.k0*self.k1))

        eps_prime = np.sqrt(4.*np.abs(self.B0)**2*eps**2 + mixing_angle_dot**2)
        eps_prime /= 2.*np.abs(B_prime)

        # avoid divergencies
        for n in (0, -1):
            eps_prime[n] = 0.0

        self.eps_prime = eps_prime
        self.delta_prime = delta
        self.theta_prime = theta_prime

        return eps_prime, delta, theta_prime

def march(x,u_e,nu):
    dx = numpy.diff(x)
    du_e = numpy.gradient(u_e,numpy.gradient(x))
    delta = numpy.full_like(x,0.)
    lam = numpy.full_like(x,lam0)

    # Initial conditions must be a stagnation point. If u_e[0]>0
    # assume stagnation is at x=0 and integrate from x=0..x[0].
    if u_e[0]<0.01:                     # stagnation point
        delta[0] = numpy.sqrt(lam0*nu/du_e[0])
    elif x[0]>0:                        # just downstream
        delta[0] = numpy.sqrt(lam0*nu*x[0]/u_e[0])
        delta[0] += 0.5*x[0]*g_pohl(delta[0],0,u_e,du_e,nu)
        lam[0] = delta[0]**2*du_e[0]/nu
    else:
        raise ValueError('x=0 must be stagnation point')

    # march!
    for i in range(len(x)-1):
        delta[i+1] = heun(g_pohl,delta[i],i,dx[i],
                          u_e,du_e,nu)  # ...additional arguments
        lam[i+1] = delta[i+1]**2*du_e[i+1]/nu

        if lam[i+1] < -12: i-=1; break  # separation condition

    return delta,lam,i+1                # return with separation index

def PlotTSUMU():
  import matplotlib.pyplot as plt
  domain = history[:, 0]

  plot_deriv = True

  if plot_deriv:
    # plot derivative of norm...
    deriv = np.gradient(history[:,3])
    deriv2 = np.gradient(deriv)

    # plt.plot(domain, history[:, 3], color='blue') # Norm of S0
    plt.plot(cjr_data, color='black')
    # plt.plot(domain, deriv2, color='red')
    plt.plot(cjr_detections, np.ones(len(cjr_detections)), marker='o', color='r', ls='')


    # from scipy.interpolate import Rbf, InterpolatedUnivariateSpline
    # rbf = Rbf(domain, deriv)


    # plt.show()
    # return

  plt.plot(domain, history[:, 1], color='red') # Minimum of first window
  plt.plot(domain, history[:, 2], color='green') # Maximum of second window
  plt.plot(domain, history[:, 3], color='blue') # Norm of S0

  # Plot the detections.
  plt.plot(detections, np.ones(len(detections)) * 500, marker='o', color='r', ls='')

  # Plot the thresholds.
  plt.plot(domain, np.ones(len(domain)) * T1, color='#aadddd') # Minimum must be lower
  plt.plot(domain, np.ones(len(domain)) * T2, color='#ddaadd') # Maximum must be higher
  plt.show()

def wrf_pv( U, V, F, THETA, PRES, MAPFAC_M, dx ):
    """Calculate the potential vorticity given the U and V vector components in m/s,
    the Coriolis sine latitude term (F) in s^-1, THETA potential temperature in degrees
    Kelvin, PRES pressure in hPa or mb, the map scale factor on mass grid and the gridspacing 
    dx in meters. U, V, F, THETA, and PRES must be 4D arrays.
    ---------------------
    U (numpy.ndarray): ndarray of U vector values in m/s
    V (numpy.ndarray): ndarray of V vector values in m/s
    F (numpy.ndarray): ndarray of Coriolis sine latitude values in s^-1
    THETA (numpy.ndarray): ndarray of potential temperature in degrees Kelvin
    PRES (numpy.ndarray): ndarray of pressure in hPa or mb same shape as THETA
    MAPFAC_M (numpy.ndarray): 2D of map scale factor on mass grid. 
    dx (float or int): float or integer of U and V grispacing in meters
    ---------------------
    returns:
        numpy.ndarray of potential vorticity values in ( K * m^2 * kg^-1 * s^-1 ) * 10^6 
        ( or 1 PVU * 10^6).
    """
    assert U.shape == V.shape == THETA.shape == PRES.shape, 'Arrays are different shapes. They must be the same shape.'
    ## pres in hPa needs to convert to Pa
    PRES = PRES * 100
    dx = dx * MAPFAC_M
    dy = dx
    grav = 9.8

    dVt,dVp,dVy,dVx = numpy.gradient( V )
    dUt,dUp,dUy,dUx = numpy.gradient( U )
    dTt,dTp,dTy,dTx = numpy.gradient( THETA )
    dPt,dp,dPy,dPx = numpy.gradient( PRES )
    return ( -grav * ( -dVp/dp * dTx/dx + dUp/dp * dTy/dy + ( dVx/dx - dUy/dy + F ) * dTp/dp ) ) * pow(10, 6)

def V(U):
    """ Spatial Viscous Fluxes at cell centres """
    rho, momx, momy, E = rollaxis(U, U.ndim-1) # shape to (neq,nx,ny)
    rho = rho.copy() # GASMODEL needs contiguous arrays
    u = momx/rho 
    v = momy/rho
    e = E/rho - 0.5*(u**2+v**2) # Gas internal energy (J/kg)
    T,Xi = GASMODEL.decode_conserved(rho, e)

    mu = 8.6412e-6*(T/288.16)**1.5*(398.16)/(T+110) # Sutherland viscosity law
    k = mu*14320.0/0.70 # This thing is important and hard to calculate

    dudx,dudy = gradient(u,dx,dy)
    dvdx,dvdy = gradient(v,dx,dy)
    dTdx,dTdy = gradient(T,dx,dy)
    tauxx = 2.0/3.0*mu*(2*dudx-dvdy)
    tauxy = mu*(dudy + dvdx)
    qx = -k*dTdx

    nx,ny,neq = U.shape
    C = zeros((nx,ny,neq)) # Input is cells, output is faces
    Q = rollaxis(C,C.ndim-1) # View of C with nicer indexing
    Q[1] = -tauxx
    Q[2] = -tauxy
    Q[3] = -u*tauxx -v*tauxy + qx
    return C

def kalman_smooth_dataframe(df, arena=None, smooth=True):
    if arena:
        fsx = arena.scale_x
        fsy = arena.scale_y
    else:
        fsx = fsy = lambda _v: _v

    #we need dt in seconds to calculate velocity. numpy returns nanoseconds here
    #because this is an array of type datetime64[ns] and I guess it retains the
    #nano part when casting
    dt = np.gradient(df.index.values.astype('float64')/SECOND_TO_NANOSEC)

    if smooth:
        print "\tsmoothing (%r)" % arena
        #smooth the positions, and recalculate the velocitys based on this.
        kf = Kalman()
        smoothed = kf.smooth(df['x'].values, df['y'].values)
        _x = fsx(smoothed[:,0])
        _y = fsy(smoothed[:,1])
    else:
        _x = fsx(df['x'].values)
        _y = fsy(df['y'].values)

    _vx = np.gradient(_x) / dt
    _vy = np.gradient(_y) / dt
    _v = np.sqrt( (_vx**2) + (_vy**2) )

    df['x'] = _x
    df['y'] = _y
    df['vx'] = _vx
    df['vy'] = _vy
    df['v'] = _v

    return dt

def draw(x, y, x_, x_fit, y_interp, group, samplenum, result):
    #
    #start to draw the figs
    y_df = np.gradient(y)
    y_df2 = np.gradient(y_df)
    fig = plt.figure(figsize=(20, 10), dpi=50, facecolor='white')
    fig.suptitle(samplenum)
    ax = fig.add_subplot(121)
    ax.plot(np.log10(x), y)
    ax.plot(np.log10(x), y_df)
    ax.plot(np.log10(x), y_df2)
    #ax.set_xscale('log')
    #ax.set_xlim(0.1, 1000)
    ax1 = fig.add_subplot(122)
    ax1.scatter(x, y, c='red')
    #ax1.plot(10**x_fit, y_interp)
    ax1.plot(x_, func(x_fit, result.params))
    ax1.set_xscale('log')
    ax1.set_xlim(0.1, 10**3)
    ax1.set_ylim(0, y.max()*1.2)

    for i in np.arange(group):
        p = 'p'+str(i)
        a = 'a'+str(i)
        b = 'b'+str(i)
        para = [result.params[p].value, result.params[a].value, result.params[b].value]
        ax1.plot(x_,single(x_fit, para))
    #ax1.plot(10**x_fit, y_-sum(y_test))
    return fig