def dreal(file1,file2):

	halo=sc.genfromtxt(file1,dtype=float)
	sub=sc.genfromtxt(file2,dtype=float)
	d_real, Vcirc, Vreal=np.array([]),np.array([]), np.array([])
	Xoff_halo, Xoff_sub=np.array([]), np.array([])

	for i in np.arange(len(halo[:,0])):
	
		mask=sub[:,14]==halo[i,11]
		N=np.size(sub[mask,14])
	
		if N >= 1.:
		
			d_real=np.append(d_real,np.sqrt(((float(sub[mask,0][0])-float(halo[i,0]))**2)+((float(sub[mask,1][0])-float(halo[i,1]))**2)+((float(sub[mask,2][0])-float(halo[i,2]))**2))*1000.)
			Vcirc=np.append(Vcirc,float(sub[mask,10][0])/float(halo[i,10]))
			Vreal=np.append(Vreal,(np.sqrt((float(sub[mask,3][0])-float(halo[i,3]))**2+(float(sub[mask,4][0])-float(halo[i,4]))**2+(float(sub[mask,5][0])-float(halo[i,5]))**2)/float(halo[i,10])))
			Xoff_halo=np.append(Xoff_halo,float(halo[i,15]))
			Xoff_sub=np.append(Xoff_sub,float(sub[mask,15][0]))
			
			N=float(len(d_real))
			Y_=1.-(np.arange(N)/N)
			d_real=np.sort(d_real)
	
	return d_real, Y_, Vcirc, Vreal, Xoff_halo, Xoff_sub

def main():
    # Open the data files
    outward = datadir + "/85Outward.csv"
    inward = datadir + "/85Inward.csv"

    outdata = _sp.genfromtxt(outward, delimiter=",", skip_header=1)
    indata = _sp.genfromtxt(inward, delimiter=",", skip_header=1)

    # Create plots
    pcoeffs, ncoeffs = mkplot(outdata, "Current adjusted increasing", "b.")
    mkplot(indata, "Current adjusted decreasing", "g.")

    # Set display options
    _plt.legend()
    _plt.xlabel("Current (A)")
    _plt.ylabel("Resonance (MHz)")
    _plt.title(r"$^{85}Rb$")

    _plt.show()

    # Calculate nuclear spin
    fit_b, fit_m, freq_err, err_b, err_m = pcoeffs
    measured_spin, err_spin = nuc_spin(fit_m, err_m, 0.29, 0.005, 135)
    spin = get_spin(measured_spin)

    # Calculate the earth's field
    outfield, sigoutfield = earth_field(outdata, spin)
    infield, siginfield = earth_field(indata, spin)

    sigfield = 1.0 / sigoutfield ** 2 + 1.0 / siginfield ** 2
    sigfield = 1.0 / _sp.sqrt(sigfield)
    field = outfield / sigoutfield ** 2 + infield / siginfield ** 2
    field *= sigfield ** 2

    return (spin, measured_spin, field, sigfield, freq_err, ncoeffs[2])

def loadgs(num, r=False, mayavi=False, **kwargs):
    if isinstance(num,str):
        name = num
        if name.endswith('_r'):

            r = True
            num = gsdict[num[:-2]]
        else:
            num = gsdict[num]
 
    else:
        name = 'gs'
    
    if not mayavi:

        output = scipy.genfromtxt(cmd_folder+'/gs_colors.txt',
                                  skip_header=256*num,
                                  skip_footer=(14-num)*256)
        if r:
            output = output[::-1]
            
        return matplotlib.colors.LinearSegmentedColormap.from_list(name, output,**kwargs)
    else:
        output = scipy.ones((256,4),dtype=int)
        output[:,0:3] = scipy.genfromtxt(cmd_folder+'/gs_colors.txt',
                                         skip_header=256*num,
                                         skip_footer=(14-num)*256,dtype=int)
        if r:
            output = output[::-1]

        return output

def fun_(var1,var2):

	halo_=sc.genfromtxt(var1)
	sub_=sc.genfromtxt(var2)

	#Variables
	
	Xoff_new_, V_circ_,dist_3d_,dist_2d_XY_, dist_2d_YZ_, dist_2d_XZ_=np.array([]), np.array([]), np.array([]),np.array([]), np.array([]), np.array([])#, Xoff_old_, V_circ_, dist_3d_, dist_2d_XY_, dist_2d_YZ_, dist_2d_XZ_, Rvirial_host_=np.array([]), np.array([]), np.array([]), np.array([]), np.array([]), np.array([]), np.array([]), np.array([])
	angulo_=np.array([])	
	#h=1
	Parameter=0.
	
	for i in np.arange(len(halo_[:,0])):
	
		mask=sub_[:,14]==halo_[i,11]
		N_=np.size(sub_[mask,14])
	
	#	print '['+str(i)+'/'+str(len(halo_[:,0]))+']'
	
		if N_==1.:
	
	#Parameters
	
			V_circ=float(sub_[mask,10][0])/float(halo_[i,10])
			dist_3d=np.sqrt(((float(sub_[mask,0][0])-float(halo_[i,0]))**2)+((float(sub_[mask,1][0])-float(halo_[i,1]))**2)+((float(sub_[mask,2][0])-float(halo_[i,2]))**2))*1000.
			dist_2d_XY=np.sqrt(((float(sub_[mask,0][0])-float(halo_[i,0]))**2)+((float(sub_[mask,1][0])-float(halo_[i,1]))**2))*1000.
			dist_2d_YZ=np.sqrt(((float(sub_[mask,1][0])-float(halo_[i,1]))**2)+((float(sub_[mask,2][0])-float(halo_[i,2]))**2))*1000.
			dist_2d_XZ=np.sqrt(((float(sub_[mask,0][0])-float(halo_[i,0]))**2)+((float(sub_[mask,2][0])-float(halo_[i,2]))**2))*1000.
			Vx, Vy, Vz=float(sub_[mask,3][0])-float(halo_[i,3]), float(sub_[mask,4][0])-float(halo_[i,4]), float(sub_[mask,5][0])-float(halo_[i,5])
			X_sh, Y_sh, Z_sh=(float(sub_[mask,0][0])-float(halo_[i,0]))*1000., (float(sub_[mask,1][0])-float(halo_[i,1]))*1000., (float(sub_[mask,2][0])-float(halo_[i,2]))*1000.
			V_r=(Vx*X_sh)+(Vy*Y_sh)+(Vz*Z_sh) # Escalar Product
			d_norm=np.sqrt(X_sh**2+Y_sh**2+Z_sh**2)
			V_norm=np.sqrt(Vx**2+Vy**2+Vz**2)
			angulo=V_r/(d_norm*V_norm)
			Rvirial_host=float(halo_[i,8])
			Xoff_new=(dist_3d/Rvirial_host)
			Xoff_old=float(halo_[i,15])
	
	#Vector
	
			angulo_=np.append(angulo_,angulo)
			Xoff_new_=np.append(Xoff_new_,Xoff_new)
#			Xoff_old_=np.append(Xoff_old_,Xoff_old)
			V_circ_=np.append(V_circ_,V_circ)
			dist_3d_=np.append(dist_3d_,dist_3d)
			dist_2d_XY_=np.append(dist_2d_XY_,dist_2d_XY)
			dist_2d_YZ_=np.append(dist_2d_YZ_,dist_2d_YZ)
			dist_2d_XZ_=np.append(dist_2d_XZ_,dist_2d_XZ)
#			Rvirial_host_=np.append(Rvirial_host_,Rvirial_host)
	
	mask_v=(V_circ_>=Parameter)
#	N=float(len(V_circ_[mask_v]))
#	Y_=1.-(np.arange(N)/N)
#	XYZ=np.sort(dist_3d_[mask_v])
#	XY=np.sort(dist_2d_XY_[mask_v])	

	return angulo_[mask_v], Xoff_new_[mask_v], dist_3d_[mask_v], dist_2d_XY_[mask_v], dist_2d_YZ_[mask_v], dist_2d_XZ_[mask_v]

def import_data_rows(data_path, data_delimiter, rows_numb = 2, cols = None, data_type = None):
    with open (data_path) as f_in:
        if cols is None:
            return scipy.genfromtxt(itertools.islice(f_in, rows_numb),
                                 delimiter = data_delimiter, dtype = data_type, usemask = False, deletechars = '"')
        elif rows_numb is None:
           return scipy.genfromtxt(f_in,
                                 delimiter = data_delimiter, dtype = data_type, usecols = cols , deletechars = '"') 
        else:
            return np.genfromtxt(itertools.islice(f_in, rows_numb), delimiter = data_delimiter,
                                 usecols = cols, dtype = data_type , deletechars = '"')

def fun_(var1,var2):

	halo_=sc.genfromtxt(var1)
	sub_=sc.genfromtxt(var2)

	#Variables
	
	Xoff_new_, Xoff_old_, V_circ_, dist_3d_, dist_2d_XY_, dist_2d_YZ_, dist_2d_XZ_, Rvirial_host_=np.array([]), np.array([]), np.array([]), np.array([]), np.array([]), np.array([]), np.array([]), np.array([])
	
	#h=1
	Parameter=0.5
	
	for i in np.arange(len(halo_[:,0])):
	
		mask=sub_[:,14]==halo_[i,11]
		N_=np.size(sub_[mask,14])
	
	#	print '['+str(i)+'/'+str(len(halo_[:,0]))+']'
	
		if N_==1.:
	
	#Parameters
	
			V_circ=float(sub_[mask,10][0])/float(halo_[i,10])
			dist_3d=np.sqrt(((float(sub_[mask,0][0])-float(halo_[i,0]))**2)+((float(sub_[mask,1][0])-float(halo_[i,1]))**2)+((float(sub_[mask,2][0])-float(halo_[i,2]))**2))*1000.
			dist_2d_XY=np.sqrt(((float(sub_[mask,0][0])-float(halo_[i,0]))**2)+((float(sub_[mask,1][0])-float(halo_[i,1]))**2))*1000.
			dist_2d_YZ=np.sqrt(((float(sub_[mask,1][0])-float(halo_[i,1]))**2)+((float(sub_[mask,2][0])-float(halo_[i,2]))**2))*1000.
			dist_2d_XZ=np.sqrt(((float(sub_[mask,0][0])-float(halo_[i,0]))**2)+((float(sub_[mask,2][0])-float(halo_[i,2]))**2))*1000.
			Rvirial_host=float(halo_[i,8])
			Xoff_new=(dist_3d/Rvirial_host)
			Xoff_old=float(halo_[i,15])
	
	#Vector
	
			Xoff_new_=np.append(Xoff_new_,Xoff_new)
			Xoff_old_=np.append(Xoff_old_,Xoff_old)
			V_circ_=np.append(V_circ_,V_circ)
			dist_3d_=np.append(dist_3d_,dist_3d)
			dist_2d_XY_=np.append(dist_2d_XY_,dist_2d_XY)
			dist_2d_YZ_=np.append(dist_2d_YZ_,dist_2d_YZ)
			dist_2d_XZ_=np.append(dist_2d_XZ_,dist_2d_XZ)
			Rvirial_host_=np.append(Rvirial_host_,Rvirial_host)
	
	mask_v=(V_circ_>=Parameter)
	N=float(len(V_circ_[mask_v]))
	Y_=1.-(np.arange(N)/N)
	XYZ=np.sort(dist_3d_[mask_v])
	XY=np.sort(dist_2d_XY_[mask_v])	

	return XYZ, XY, Y_, dist_2d_XY_, V_circ_

def load_data_set(file_location, delimiter, column_y, column_x1):
    """
    :param file_location: string for data file location
    :param delimiter: ',' for CSV, etc.
    :param column_y: the column containing the target values
    :param column_x1: input data column -- right now, I only take 1
    TODO: add a parameter that is a set so I can take multiple input columns
    :return: Numpy Array of target values, input values, and bias term (bias term always = 1.0)
    """
    data = sp.genfromtxt(file_location, delimiter=delimiter, dtype=None)

    # Need to get everything after the headers

    X = data[1:, column_x1]
    Y = data[1:, column_y]

    # we make the cases 1 and -1 to fit with the Likelihood Formula
    # P(y | x) -> h(x) for y = +1
    # P(y | x) -> 1 - h(x) for y = -1
    y_numeric = [1.0 if entry == 'Yes' else -1.0 for entry in Y]

    # Will use this for x0
    ones = [1.0 for x in X]

    return np.array(zip(y_numeric, X, ones), dtype='float_')

def main():
    data = sp.genfromtxt('./data/web_traffic.tsv', delimiter='\t')
    x = data[:, 0]
    y = data[:, 1]
    x = x[~sp.isnan(y)]
    y = y[~sp.isnan(y)]
    fp1 = sp.polyfit(x, y, 1)
    print('Model parameters for fp1 %s' % fp1)
    f1 = sp.poly1d(fp1)
    print('This is the error rate for fp1 %f' % error(f1, x, y))

    fp2 = sp.polyfit(x, y, 2)
    print('Model parameters for fp2 %s' % fp2)
    f2 = sp.poly1d(fp2)
    print('This is the error rate for fp2 %f' % error(f2, x, y))

    plt.scatter(x, y,color= 'pink')
    plt.title('My first impression')
    plt.xlabel('Time')
    plt.ylabel('#Hits')
    plt.xticks([w * 7 * 24 for w in range(10)], ['week %i' % w for w in range(10)])
    fx = sp.linspace(0, x[-1], 1000)
    plt.plot(fx, f1(fx), linewidth=3,color='cyan')


    plt.plot(fx, f2(fx), linewidth=3, linestyle='--',color= 'red')
    plt.legend(['d = %i' %f1.order, 'd = %i' %f2.order], loc='upper left')
    plt.autoscale(tight=True)
    plt.grid()
    plt.show()

def plot_2():
    data=sp.genfromtxt('F:\EPAM\coursera\ML_COURSERA_GARVARD\machine-learning-ex1\machine-learning-ex1\ex1\ex1data1.txt',delimiter=',')
    x=data[:,0]
    y=data[:,1] 
    m=len(y)
    y=y.reshape(m,1)
    x1=np.array([])
    
    for xi in x:
        x1=np.append(x1,[1,xi])
    x=x1.reshape(m,2)
    
    theta=np.zeros((2,1))
    iterations = 1500;
    alpha = 0.01;
   
    cost=computerCost(x,y,theta)
    theta=Ggradient_Descent(x,y,theta,alpha,iterations)
    
    print cost
    print theta 
    pr1=np.array([1,3.5]).dot(theta)
    pr2=np.array([1,7]).dot(theta)
   
    print pr1
    print pr2
    y_1=x.dot(theta)
    y_1.shape=(m,1)
    plt.title('Linear regression')
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.plot(x[:,1],y,'b-')
    plt.plot(x[:,1],y_1,'r-')
    plt.show(block=True)

def loadct(num, **kwargs):
    file = os.path.join(ifigure.__path__[0], 'utils', 'idl_colors.txt',)
    output = scipy.genfromtxt(file,
                              skip_header=256*num,
                              skip_footer=(39-num)*256)/255.
    return matplotlib.colors.LinearSegmentedColormap.from_list('idl'+str(num), 
                                                               output, **kwargs)

def init_and_cleanup_data(path, delimiter):
    data = sp.genfromtxt(path, delimiter=delimiter)
    hours = data[:, 0] # contains the hours
    webhits = data[:, 1] # contains the number of web hits at a particular hour
    hours = hours[~sp.isnan(webhits)]
    webhits = webhits[~sp.isnan(webhits)]
    return (hours, webhits)

	def AtmosphereCondition(self,Height):
		data = sp.genfromtxt('AC.d',delimiter=',')		
		PresGrand = 0.1013   #[MPa]
		for i in range(0,520):
			if Height < data[i,0]:
				break;
		return data[i-1,1],data[i-1,2]*PresGrand,data[i-1,3]

def read_inputcat_for_mgc3(filename,pardic=None):

  #Open file
  if '.gz' in filename:
     inputfile=gzip.open(filename,'r')
     filename=filename.replace('.gz','')
  else: inputfile=open(filename,'r')

  obsdata = scipy.genfromtxt(inputfile,comments='#')
  
  #Deal with one-line files
  if np.ndim(obsdata)==1: obsdata=np.reshape(obsdata,(1,obsdata.size))

  #Do cuts
  mask = obsdata[:,0]==obsdata[:,0]  #Initialize mask to all-True-vector
  if pardic:
    for NAUX in range(1,pardic['NAUX']+1,1):
      mykey_col='AUX%d_col' % (NAUX)
      mykey_valo='AUX%d_o' % (NAUX)
      mykey_valf='AUX%d_f' % (NAUX)
      #Skip if col=998   
      if pardic[mykey_col]!=998:
       print(' Cutting input catalogue with %.1f<%s[%d]<%.1f' % (pardic[mykey_valo],mykey_col,pardic[mykey_col]+1,pardic[mykey_valf]))
       #Create mask 
       mask_i = (obsdata[:,pardic[mykey_col]]>pardic[mykey_valo]) & (obsdata[:,pardic[mykey_col]]<pardic[mykey_valf])
       #Combine masks
       mask = mask & mask_i

  #Apply mask
  obsdata=obsdata[mask,:]
    
  #Return data
  return (obsdata,filename)

def LDA_batch_normalization(dataset, sample_table, batch_col, output_folder, ncomps): # this is actually the batch normalization method
   
    tmp_output_folder = os.path.join(output_folder, 'tmp')

    if not os.path.isdir(tmp_output_folder):
        os.makedirs(tmp_output_folder)
    
    barcodes, filtered_conditions, filtered_matrix, conditions, matrix = dataset
    
    # Remove any remaining NaNs and Infs from the filtered matrix - they would screw
    # up the LDA. 
    filtered_matrix[scipy.isnan(filtered_matrix)] = 0
    filtered_matrix[scipy.isinf(filtered_matrix)] = 0

    # For full matrix, also eliminate NaNs and Infs, BUT preserve the indices and values
    # so they can be added back into the matrix later (not implemented yet, and may never
    # be - there should no longer be NaNs and Infs in the dataset)
    # The NaNs and Infs will mess up the final step of the MATLAB LDA script, which uses
    # matrix multiplication to remove the specified number of components!
    matrix_nan_inds = scipy.isnan(matrix)
    matrix_nan_vals = matrix[matrix_nan_inds]
    matrix_inf_inds = scipy.isinf(matrix)
    matrix_inf_vals = matrix[matrix_inf_inds]

    matrix[matrix_nan_inds] = 0
    matrix[matrix_inf_inds] = 0

    # Save both the small matrix (for determining the components to remove) and the 
    # full matrix for the matlab script
    filtered_matrix_tmp_filename = os.path.join(tmp_output_folder, 'nonreplicating_matrix.txt')
    full_matrix_tmp_filename = os.path.join(tmp_output_folder, 'full_matrix.txt')
    
    np.savetxt(filtered_matrix_tmp_filename, filtered_matrix)
    np.savetxt(full_matrix_tmp_filename, matrix)

    # Map the batch to integers for matlab, and write out to a file so matlab can read
    # Note that yes, the batch_classes should match up with the filtered matrix, not
    # the full matrix
    batch_classes = get_batch_classes(dataset = [barcodes, filtered_conditions, filtered_matrix], sample_table = sample_table, batch_col = batch_col)
    class_tmp_filename = os.path.join(tmp_output_folder, 'classes.txt')
    writeList(batch_classes, class_tmp_filename)
   
    output_tmp_filename = os.path.join(tmp_output_folder, 'full_matrix_lda_normalized.txt')
    runLDAMatlabFunc(filtered_matrix_filename = filtered_matrix_tmp_filename, \
            matrix_filename = full_matrix_tmp_filename, \
            class_filename = class_tmp_filename, \
            ncomps = ncomps, \
            output_filename = output_tmp_filename)
    # The X norm that is returned is the full matrix. In the future, we could add in
    # returning the components to remove so they can be visualized or applied to other
    # one-off datasets
    Xnorm =  scipy.genfromtxt(output_tmp_filename)

    ## Dump the dataset out!
    #output_filename = os.path.join(mtag_effect_folder, 'scaleddeviation_full_mtag_lda_{}.dump.gz'.format(ncomps))
    #of = gzip.open(output_filename, 'wb')
    #cPickle.dump([barcodes, conditions, Xnorm], of)
    #of.close()

    return [barcodes, conditions, Xnorm]

def read_itcsimlib_exp( file, exp_args={} ):
	from scipy import genfromtxt
	from .itc_experiment import ITCExperiment

	ignore = ("itcsim","Date","Ivol","units")
	data,h = genfromtxt(file,unpack=True),open(file)
	kwargs = {'Cell':{},'Syringe':{}}
	for a in [l.split()[1:] for l in h.readlines() if l[0]=='#']:
		if a == [] or a[0] in ignore:
			continue
		elif a[0] == 'Cell' or a[0] == 'Syringe':
			kwargs[a[0]][a[1]] = float(a[2])
		elif a[0].lower() == 'skip':
			kwargs['skip'] = map(int,a[1:])
		else:
			kwargs[a[0]] = float(a[1])
	h.close()
	if not 'title' in kwargs:
		kwargs['title'] = os.path.splitext(os.path.basename(file))[0]
	
	# overwrite any file-obtained info with explicit values
	kwargs.update(exp_args)
	if len(data) == 2:
		return ITCExperiment(injections=data[0],dQ=data[1],**kwargs)
	elif len(data) == 3:
		return ITCExperiment(injections=data[0],dQ=data[1],dQ_err=data[2],**kwargs)
	else:
		return None # TODO : parser errors

def TSP(stops, Alg, steps, param, seed = None,
                      coordfile = 'xycoords.txt'):
    '''A wrapper function that attempts to optimize the traveling 
    salesperson problem using a specified algorithm. If coordfile
    exists, a preexisting set of coordinates will be used. Otherwise,
    a new set of "stops" coordinates will be generated for the person to 
    traverse, and will be written to the specified file.'''
    
    ## Create the distance matrix, which will be used to calculate
    ## the fitness of a given path
    if os.path.isfile(coordfile):
        coords = scipy.genfromtxt(coordfile)
        distMat = DistanceMatrix(coords)
    else:
        distMat = GenerateMap(stops, fname = coordfile, seed = seed)

    if Alg == 'HC':
        ## param is the number of solutions to try per step
        bestSol, fitHistory = HillClimber(steps, param, distMat, seed)
    elif Alg == 'SA':
        ## param is a placeholder
        bestSol, fitHistory = SimulatedAnnealing(steps, param, distMat, seed)
    elif Alg == 'MC3':
        ## param is the number of chains
        bestSol, fitHistory = MCMCMC(steps, param, distMat, seed)
    elif Alg == 'GA':
        ## param is the population size
        bestSol, fitHistory = GeneticAlgorithm(steps, param, distMat, seed)
    else:
        raise ValueError('Algorithm must be "HC", "SA", "MC3", or "GA".')

    outfname = coordfile + '-' + Alg + '-' + str(steps) + '-' + str(param) + '.txt'
    scipy.savetxt(outfname, scipy.array(bestSol), fmt = '%i')
    return bestSol, fitHistory

def MainFrac(dire,MdotO,Time):
	#{{{
	Frac = sp.genfromtxt(dire+'Fraction.d',delimiter=',')

	#CH4 = Frac[:,0]
	#CO2 = Frac[:,1]
	#CO  = Frac[:,2]
	#H   = Frac[:,3]
	#H2  = Frac[:,4]
	#H2O = Frac[:,5]
	#O   = Frac[:,6]
	#O2  = Frac[:,7]
	#OH  = Frac[:,8]

	plt.plot(Time,Frac[:,1])
	plt.plot(Time,Frac[:,2])
	plt.plot(Time,Frac[:,3])
	plt.plot(Time,Frac[:,4])
	plt.plot(Time,Frac[:,5],'--')
	plt.plot(Time,Frac[:,6],'--')
	plt.plot(Time,Frac[:,7],'--')
	plt.plot(Time,Frac[:,8],'--')
	plt.grid()
	##plt.title('MdotO:%s[kg/s]'%MdotO)
	plt.legend(( 'CO2 ',  'CO ',  'H ',  'H2 ',  'H2O ',  'O ',  'O2 ', 'OH' ))
	plt.xlabel('Time[s]')
	plt.ylabel('Mole Fraction')
	#plt.ylim([0,100])
	plt.savefig(dire+"Frac_Mdot%s.png"%float(MdotO))
	plt.close()

def extremalPlot(V,c):
    datZM = genfromtxt('zero-mode-V0='+str(V)+'.dat')
    rh = datZM[:,-3]
    a = datZM[:,-2]
    w = a/(a*a+rh*rh)
    Qe = V * (a*a+rh*rh)/rh
    M = (a*a+rh*rh+Qe*Qe)/(2*rh)
    plot(w,M,c+'-',ms=2.5)
    extremalLine = lambda w: sqrt((1.0 - 4.0*V*V + sqrt(1.0+8.0*V*V))/w)/(2.0*sqrt(2.0*w))
    ws = linspace(0.75,1,100)
    plot(ws,extremalLine(ws),c+'-',ms=2.5)

    datExtremal = genfromtxt('gnuplot-extremal-V0='+str(V)+'.dat')
    wExtremal = datExtremal[:,0]
    MExtremal = datExtremal[:,2]
    plot(wExtremal,MExtremal,c+'--',ms=2.5)

def ReadDatas(dataDirectory):
    datas = []
    for fileName in os.listdir(dataDirectory):
        fullPath = os.path.join(dataDirectory, fileName)
        print("LOADING : " + fullPath)
        if os.path.isfile(fullPath):
            datas.append(sp.genfromtxt(fullPath, delimiter=","))
    return datas[::-1]  # because the original data is from late to early

    def load_single_data(self):
        single_data = sp.genfromtxt(os.path.join(DATA_DIR, "ex1data1.txt"), delimiter=",")
        x = single_data[:, 0]
        y = single_data[:, 1]

        x = np.transpose(np.atleast_2d(x))

        return (x, y)