# -*- coding: utf-8 -*- ############################################################################################# ############################################################################################# ############################################################################################# ############################################################################################# ############ Routine for finding local minima in D-dimensional potentials ################ ############ This version uses several minimization techniques beginning with ############## ############ a GA to generate clusters. The package contains three programmes ############## ############ -- this file, the file GA_latex which contains the latex ###################### ############ and the file potential.py for the potential (and its negative) ################ ############################################################################################# ############################################################################################# ############################################################################################# ############################################################################################# from itertools import * from decimal import * from scipy import optimize from sklearn import cluster from collections import defaultdict from matplotlib import cm from sklearn.metrics.cluster import normalized_mutual_info_score from sklearn.metrics.cluster import adjusted_rand_score from sklearn.cluster import KMeans from neldermead import NelderMead from scipy.optimize import minimize import networkx as nx import seaborn as sns import pandas as pd import numpy as np import math as math import matplotlib.pyplot as plt import random import sys import subprocess import shlex ############################################################################################# ################ set paramaters for the GA (dimensions etc) here ########################### ############################################################################################# ############# include your potential and latex comments in two .py files ==> ################ """ ============================== your filenames (.py) for potential and latex ============================== """ from potential6 import * from GA_latex import * """ ========================================================================================================== """ global ddim,dim,ch_number,setvariables,setvalues,pop,alpha,range,origin,mutate_rate,creatures,amoeba_num,clusternum,meteor_after,domain,scale_factors,angles # creatures[i->ch_length][j->ch_number][k->pop] is a 3 dim array global maxima,minima """ =============================== parameters are initialized in this block ================================== """ verbose = True # does what it says on the tin pop =150 # population -- with this method a larger population than a usual GA is better because it generates more clusters. meteor_after = 10 # number of generations to run before running the clustering algorithm. Value of 0 turns off the GA. alpha=3.0 # parameter for the breeding likelihood - the fittest individual is chosen alpha * the least fit (normally 2 - 3) mutate_rate=0.1 # mutation rate as a proportion of the alleles to randomly flip (2% = 0.02) -- to get better clusters # it pays to make this number higher than you would in a normal GA -- say 10-15% ch_length=8 # significant places for the coordinates """ ======= adjust paremeters in here ======= """ setvariables = np.array([2,3,4,5,6,7,8,9]) # list of variable positions that are fixed -- the dimensionality should be dim-ch_number freevariables = np.array([0,1]) # list of variable positions that are free -- complement of setvariables angles = np.array([]) # list of positions that are angles -- these get adjusted to principle branch: normally leave empty #================= make sure the below have the right dimensions ==> 5 setvalues = np.array([0.0 ,0.0 ,0.0 ,0.0, 0.0, 0.0, 0.0 ,0.0, 0.0, 0.0 ]) # holds values of variables in the fixed positions specified in the previous line -- with the others being dummy (e.g. set to zero) domain = np.array([ 6 , 6 , 1 , 1 , 1, 1, 1. , 1 , 1, 1]) # size of starting domain of the search-space -- # define origin below -- so far origin is at zero. Note that the 2/3 entries are powers of 10 in the potential origin = np.array([-3, -3, 0, 0.0, 0.0, 0.0, 0, 0.0, 0.0, 0.0]) # the position of the bottom left corner (will be scaled along with the domain by factor below) factor = np.array([1.2, 1.2 ,1., 1.,1.,1.,1., 1.,1.,1.]) # a factor that inflates the size of the domain each time plotfactor = 0.02 # this factor multiplies the FINAL domain to give the size of contour plots around the minima """ ======== end of adjust paremeters ======= """ """ do not adjust below unless you know what you are doing """ ch_number=len(freevariables) # dimensionality of freevariables being adjusted ddim=len(freevariables) # the dimensionality of the gene space (this can be larger than the number of variables you want -- fix the superfluous ones below) dim=ddim+len(setvariables) # the dimensionality of the gene space (this can be larger than the number of variables you want) factorings = 10 # number of times you expand the domain clusternum = 5 # the number of clusters to use in the cluster algorithm -- must be less than pop/(ddim+1) amoeba_num = 100 # number of iterations in the nelder-mead algorithm tolerance_add = 0.005 # this factor multiplies the domain and is added to the minima to check for flat directions tolerance_mult = 0.01 # this factor multiplies the domain and multiplies the minima to check for flat directions """ ========================================================================================================== """ ############################################################################################# ############################## GET LATEX DEFN's FOR OUTPUT ################################ ############################################################################################# content1 = text1() content3 = text3() def writeout(stuff): output_file.write(stuff) def view(): # latex/make/update the pdf # return_value = subprocess.call(['pdflatex', 'output.tex'], shell=False) # shell should be set to False proc=subprocess.Popen(shlex.split('pdflatex output.tex')) proc.communicate() output_file = open('output.tex',"w") # open the output file - touch this if it isn't there already ################################## USAGE TO COMPILE LATEX ################################### writeout(content1) # latex preamble and title may as well leave this here #writeout(content3) # the stuff to end the document - also bibliography maybe #output_file.close() # put this right at the end of the file #view() # compile the latex and produce/update the pdf file after file closed ############################################################################################# ############################################################################################# ############################################################################################# ############################ SUBROUTINES START HERE #################################### ############################ edit in above for diffn't preamble ############################ ############################################################################################# ############################################################################################# ############################################################################################# """ stuff specific for the minimisation problem """ x = [] # initialize coordinates b = np.array([10**(-i) for i in range(0,ch_length) ]) scale_factors = np.array([p/10. for p in domain]) # converson factor to get to actual x,y,z from a 0-10 interval def phenotype(creaturess): global origin # in this case returns the position for each chromosomes as a 1-dim aray" # in long hand - return array([[ dot(creatures[k][j],b) for j in range(0,ch_number) ] for k in range(0,pop)]) ptype = np.dot(creaturess,b) tmp = scale_factors * ptype + origin output = tmp for i in range(0,pop): for j in angles: output[i][j] = np.array([ ( tmp[i][j] * math.pi % 1 ) /math.pi ]) for j in setvariables: output[i][j] = setvalues[j] return output def inrange(coords,basis): # returns true if a point is inside the domain otherwise false. basis holds the flat directions global domain,origin res = True for i in freevariables: iinbasis = False if i in basis: iinbasis=True if iinbasis == False: if ( coords[i] - origin[i] > domain[i] or coords[i] - origin[i] < 0.) : res = False return res ############################################################################################# ################ below is portable as long as phenotype written as above ################## ############################################################################################# likelihood = np.array([ 2.0/(1+alpha)/pop * ( (alpha-1)*i/(pop-1) + 1 ) for i in range(0,pop) ]) # in this programme the fittest is ordered, likelihood is a trivial function of position def build(): return np.array([[[ np.random.randint(0,10) for i in range(0,ch_length)] for j in range(0,dim) ] for k in range(0,pop)]) # return random np.array for population. Include random placeholder rows for fixed variables def mag(x): return math.sqrt(sum(p**2 for p in x)) # returns magnitude of a vector def reorder(): # reorders creatures with least fit first and fittest last global creatures,maxima #in long hand - fitness = [ f(c[k]) for k in range(0,pop)] c=phenotype(creatures) ct=np.transpose(c) fitness = f(ct) # this is provided by the file potential4.py # "crowding penalty" for this problem based on proximity to previous solutions along non-flat directions - need to extract the right coordinates for q in maxima: qcoords= [] qbasis = [] # basis of non-flat directions. This also includes fixed variables which are also "flat" for i in range(0,dim): if i not in q[1] and i not in setvariables: qbasis.append(i) # print qbasis if qbasis != []: # if all directions seem to be flat there is nothing to do qcoords = np.array([ q[3][j] for j in qbasis]) # print qcoords for i in range(0,pop): p = c[i] pcoords = np.array([ p[j] for j in qbasis]) r = mag(pcoords-qcoords)+1e-5 fitness[i] *= (1-0.0**dim/r**dim) newi = np.argsort(fitness) # returns indices for re_ordering creature with fittest goes last oldcreatures = creatures[:] # copying the np.array creatures = oldcreatures[newi] def phenoav(creatures): #returns the average list of physical data from the chromosomes as a 1-dim np.array tmp=phenotype(creatures) return tmp.sum(axis=0)/Decimal(pop) def breeding_pair(): global likelihood "Generates a set of breeding pairs randomly based on fitness" pair=[] for k in range(0,2): bonk = 0 while bonk == 0: candidate = np.random.randint(0,pop) ok = np.random.random()/pop if ok < likelihood[candidate]: bonk = 1 pair.append(candidate) return pair def dying_pair(): global likelihood "Generates a set of dying pairs randomly based on fitness" pair=[] for k in range(0,2): die = 0 while die == 0: # die as in dice that is candidate = np.random.randint(1,pop) ok = np.random.random()/pop if ok > likelihood[candidate]: die = 1 pair.append(candidate) return pair def shuffle(zwei): "construct a chromosome at random from two others with a spot of mutation" yy= np.array(zwei[0]) zz= np.array(zwei[1]) length = len(yy) splice = random.randint(0,length-1) vv= np.array([yy[0:splice],zz[splice:]]) xx = np.hstack(vv) return xx def breed(creatures2): p=breeding_pair() dp=dying_pair() creatures[dp[0]] = np.array([ shuffle( [creatures2[ p[0] ][j],creatures2[ p[1] ][j] ] ) for j in range(0,dim) ]) creatures[dp[1]] = np.array([ shuffle( [creatures2[ p[0] ][j],creatures2[ p[1] ][j] ] ) for j in range(0,dim) ]) def mutate(creatures2): "introduce mutations everywhere except the best creature" mutation_number=int(mutate_rate*ch_length*dim*pop) for l in range(0,mutation_number): creatures2[np.random.randint(1,pop)][np.random.randint(0,dim)][np.random.randint(0,ch_length)]=np.random.randint(0,10) return creatures2 ############################################################################################# #################### Some plotting routines with matplotlib ############################### ############################################################################################# def plotterxy(plotnum): # this is only in 2D -- not decided what to do about projecting higher dims yet print ("plot") x = np.linspace(origin[0],origin[0]+1*scale_factors[0],100) #tau4 y = np.linspace(origin[1],origin[1]+1*scale_factors[1],100) #tau5 (X,Y) = np.meshgrid(x,y) # a = f([X,Y]) # c = plt.pcolor(x,y,a, cmap=plt.cm.Reds) lx = plt.xlabel(names[0],fontsize=18) ly = plt.ylabel(names[1],fontsize=18) plt.title(r'Generation number = ' + str(breeds)) vals=phenotype(creatures) # CS = plt.contour(X, Y, -a, (0) ) plt.scatter(vals[:, 0], vals[:, 1]) # plt.show() filename="file"+str(plotnum)+".png" print("filename="+filename) plt.savefig(filename) plt.close() def plottermax(minima): # plots the positions of the final minima in coord1 coord2 with contour plot of V=0 print ("plot minima") x = np.linspace(origin[0],origin[0]+2*scale_factors[0],200) y = np.linspace(origin[1],origin[1]+2*scale_factors[1],200) (X,Y) = np.meshgrid(x,y) lx = plt.xlabel(names[0],fontsize=18) ly = plt.ylabel(names[1],fontsize=18) plt.title('Minima found') xvals=[] yvals=[] for p in minima: xvals.append(p[3][0]) yvals.append(p[3][1]) plt.scatter(xvals, yvals) #plt.show() filename="maxima.png" print("filename="+filename) plt.savefig(filename) plt.close() def plot_minima(minima): # produces a set of plots slicing through all 2D directions around each minimum # and adds them into the latex file global text_minima,angles,names text_minima=r''' ~ ''' figcount=0 print("plot around the minima") min_num = 0 for min in minima: print ("min",min) sub_plot_num = 0 q= min[3] # should be all dim coordinates d= min[2] # depth qbasis = [] # basis of the true coordinates for i in freevariables: if i not in min[1]: qbasis.append(i) # gets the actual indices of coords to be plotted that are not flat if len(qbasis)>1: figcount+=1 if figcount > 15: # stops the latex falling over by inserting \newpage if too many figures text_minima += '\n' + r''' \newpage ''' figcount=0 text_minima += '\n' + r''' \begin{figure} ''' min_num+=1 plt.title(r'Minimum number ' + str(min_num) + r' Depth =' + str(d)) # take pairs of true-coordinates for contour plots for i in range(0,len(qbasis)-1): ii = qbasis[i] plotrange_x = plotfactor * domain[ii] if ii in angles: plotrange_x = 3.142 for j in range(i+1,len(qbasis)): sub_plot_num += 1 jj = qbasis[j] plotrange_y = 4 * plotfactor * domain[ii] if jj in angles: plotrange_y = 3.142 arg=q # the dimensions = dim vector of points x = np.linspace(q[ii]-plotrange_x,q[ii]+plotrange_x,200) y = np.linspace(q[jj]-plotrange_y,q[jj]+plotrange_y,200) # construct the argument in V (X,Y) = np.meshgrid(x,y) print (ii,jj) # must be a more pythonic way to do the below but can't be bothered to figure it out yet if ii==0 and jj==1: a = f([X, Y , q[2] , q[3], q[4] , q[5] ])/d elif ii==0 and jj==2: a = f([X, q[1], Y , q[3], q[4], q[5] ])/d elif ii==0 and jj==3: a = f([X, q[1], q[2], Y, q[4], q[5] ])/d elif ii==0 and jj==4: a = f([X, q[1], q[2], q[3], Y, q[5] ])/d elif ii==0 and jj==4: a = f([X, q[1], q[2], q[3], q[4], Y ])/d elif ii==1 and jj==2: a = f([q[0] , X , Y , q[3], q[4], q[5] ])/d elif ii==1 and jj==3: a = f([q[0], X , q[2] , Y, q[4], q[5] ])/d elif ii==1 and jj==4: a = f([q[0], X , q[2] , q[3], Y, q[5] ])/d elif ii==1 and jj==5: a = f([q[0], X , q[2] , q[3], q[4], Y ])/d elif ii==2 and jj==3: a = f([q[0] , q[1] , X , Y, q[4], q[5] ])/d elif ii==2 and jj==4: a = f([q[0] , q[1] , X , q[3], Y, q[5] ])/d elif ii==2 and jj==5: a = f([q[0] , q[1] , X , q[3], q[4], Y ])/d elif ii==3 and jj==4: a = f([q[0] , q[1] , q[2] , X, Y, q[5] ])/d elif ii==3 and jj==5: a = f([q[0] , q[1] , q[2] , X, q[4], Y ])/d elif ii==4 and jj==5: a = f([q[0] , q[1] , q[2] , q[3], X, Y ])/d conts= np.array( [ (1.0 + (jjj-100)*0.02) for jjj in range(0,200)] ) c = plt.pcolor(x,y,a,cmap=plt.cm.Greys) CS = plt.contour(X, Y, -a, conts ) lx = plt.xlabel(names[ii],fontsize=18) ly = plt.ylabel(names[jj],fontsize=18) filename="minimum_"+str(min_num)+"_"+str(ii)+str(jj)+".png" text_minima+= '\n'+ r'''\includegraphics[scale=0.3]{''' + filename + r'''}''' print("filename="+filename) plt.savefig(filename) plt.close() text_minima += r'''\caption{\emph{Local contour plots around ''' for p in qbasis: text_minima += names[p]+r',' if min[1]==[] : text_minima += r' = ' +str(min[0])+r'.\label{fig:minimum'+str(min_num)+r' }. Note that values of fixed parameters also appear here.}}\end{figure}''' else: text_minima += r' = ' +str(min[0])+r'. The flat directions have been identified as ' for p in min[1]: text_minima += names[p]+r', ' text_minima+= r'while the other directions are genuinely fixed.\label{fig:minimum'+str(min_num)+r' }. Note that values of fixed parameters also appear here.} }\end{figure}' if len(qbasis)==1: figcount+=1 if figcount > 15: # stops the latex falling over by inserting \newpage if too many figures text_minima += '\n' + r''' \newline ''' figcount=0 text_minima += '\n' + r''' \begin{figure} ''' min_num+=1 plt.title(r'Minimum number ' + str(min_num) + r' Depth =' + str(d)) ii = qbasis[0] plotrange_x = 5*q[ii] X = np.linspace(q[ii]-plotrange_x,q[ii]+plotrange_x,200) #(X) = meshgrid(x) if ii==0: a = potential([X, q[1] , q[2] , q[3], q[4], q[5] ]) elif ii==1 : a = potential([q[0] , X , q[2] , q[3], q[4], q[5] ]) elif ii==2 : a = potential([q[0] , q[1] , X , q[3], q[4], q[5] ]) elif ii==3 : a = potential([q[0] , q[1] , a[2] , X, q[4], q[5] ]) elif ii==4 : a = potential([q[0] , q[1] , a[2] , q[3], X, q[5] ]) c = plt.plot(X,a,'ro') lx = plt.xlabel(names[ii],fontsize=18) ly = plt.ylabel(r'$V$',fontsize=18) filename="minimum_"+str(min_num)+"_"+str(ii)+".png" text_minima+= '\n'+ r'''\includegraphics[scale=0.3]{''' + filename + r'''}''' print("filename="+filename) plt.savefig(filename) plt.close() text_minima += r'''\caption{\emph{Local plot around ''' for p in qbasis: text_minima += names[p] text_minima += r' = ' +str(min[0])+r'.\label{fig:minimum'+str(min_num)+r' }} }\end{figure}''' ############################################################################################# ######################### Main GA algorithm -- used to be stand-alone ####################### ############################################################################################# creatures = build() def GArun(): global creatures,mutate_rate,breeds,alpha,mutate_rate,meteor_after creatures = build() if meteor_after==0: return for mutate_from in range(0,1): plotornot=0 plotnum=1 for breeds in range(0,meteor_after): reorder() # order in terms of fitness ranking with fittest last creaturestmp=creatures[pop-1][:] for sub in range(0,pop): breed(creatures) creatures=mutate(creatures) creatures[0]=creaturestmp[:] # copies the previous fittest individual over the least fit plotornot += 1 if (plotornot > 10): plotterxy(plotnum) plotornot=0 plotnum +=1 plotterxy(99) ############################################################################################# ############# clustering algorithm -- returns a set of simplices around minima ############## ############################################################################################# def simplices(clusternum): # makes clusters then returns an np.array of np.array of [x,y,..z] first dim+1 entries for each cluster # change now it forms as many simplices as possible from each cluster global simp_count X=phenotype(creatures) kmeans = KMeans(n_clusters=clusternum) kmeans.fit(X) labels=kmeans.predict(X) centres=kmeans.cluster_centers_ counters= [ np.where(labels == p) for p in range(0,clusternum) ] # make the big clusters counters_points = (np.array([ [ X[i] for i in p] for p in counters])) # choose the best dim+1 for each cluster counters_refined=[] for i in range(0,len(counters_points)): p = counters_points[i][0] tmp=[] for q in p: if len(tmp) < dim+1 : tmp.append(q) # print q,tmp,len(tmp) if len(tmp)==dim+1 : counters_refined.append(tmp) return counters_refined counters_refined=[] for i in range(0,len(counters_points)): p = counters_points[i][0] # p are the coodinates in each cluster tmp=[] for q in p: # collect into x d+1 dimensional simplices if len(tmp) < dim+1: tmp.append(q) if len(tmp)==dim+1 : counters_refined.append(tmp) return counters_refined # returns a set of simplices ############################################################################################# ##################### amoeba algorithm -- optimizes around minima ########################### ############################################################################################# def nelder(simps): # do nelder-mead for a set of simplices global dim,amoeba_num,domain,maxima,maxima_coords for i in range(0,len(counters_refined)): simplex = np.transpose(np.array(counters_refined)[i]) simplex_test=simplex simplex_test2=simplex simplex_add=0*simplex for j in range(0,len(domain)): row=tolerance_add*domain[j]*np.random.rand(len(simplex[0])) simplex_add[j] += row simplex_test = (1+tolerance_mult) * simplex # test to determine valleys simplex_test2 = simplex+simplex_add # test to determine valleys OK_min = True # local switch that gets set to False if the point fails various tests # The routine below is a Numpy routine that seems to work less well on this problem # nm = NelderMead(dim, negf, simplex) # best_point = nm.optimize(amoeba_num) # coords=np.transpose(best_point[0])[0] params = {'maxiter': amoeba_num, 'initial_simplex': simplex.T, 'xatol': 1e-40, 'fatol': 1e-40 } params2 = {'maxiter': amoeba_num, 'initial_simplex': simplex_test.T, 'xatol': 1e-40, 'fatol': 1e-40 } params3 = {'maxiter': amoeba_num, 'initial_simplex': simplex_test2.T, 'xatol': 1e-40, 'fatol': 1e-40 } scipy_nm = minimize(negf, x0=np.zeros(dim), method='Nelder-Mead', options=params) scipy2_nm = minimize(negf, x0=np.zeros(dim), method='Nelder-Mead', options=params2) scipy3_nm = minimize(negf, x0=np.zeros(dim), method='Nelder-Mead', options=params3) coords= scipy_nm.x best_point=scipy_nm.fun coords2=scipy2_nm.x best_point2=scipy2_nm.fun coords3=scipy3_nm.x best_point3=scipy3_nm.fun print ("Scipy Nelder coords",coords) print ("Adjusted Scipy Nelder coords",coords2) print ("Coordinate difference",coords - coords2) pole_test_depth = -5.0 try: if (best_point < pole_test_depth): OK_min = False # reject if pole print ("fail at pole_test",best_point) if (isNaN(best_point) == True): OK_min = False # reject if overflow print ("fail at nan_test") print (best_point) except RuntimeWarning: print ("fail at runtimewarning_test",best_point) OK_min = False basis = [] # np.array to hold flat directions flatnum=0 if OK_min == True: for i in freevariables: if abs( (coords[i] - coords2[i])/(coords[i]+coords2[i]+10.0/domain[i] )) > 0.03 or abs( (coords[i] - coords3[i])/(coords[i] + coords3[i]+10.0/domain[i] )) > 0.03 : if abs( (best_point - best_point2 )/(best_point+ best_point2) ) > 0.01 : print ("fail at 1st coord_test",best_point,best_point2) OK_min = False print ('Nelder Fails') if abs( (best_point - best_point3 )/(best_point + best_point3) ) > 0.01 : print ("fail at 2nd coord_test",best_point,best_point3) OK_min = False print ('Nelder Fails') else: basis.append(i) flatnum+=1 print ('Nelder finds flat directions along the directions',basis) if flatnum == dim : OK_min = False # reject if ONLY flat directions. print ("fail because all flat") if (inrange(coords,basis) == False): # OK_min = False # reject if it falls outside domain. This also puts angles on principle branch print ('wtf', coords,domain) print ("We have found a minimum = ",OK_min) if (OK_min==True): true_coords=[] # fixed directions for i in range(0,dim): if i not in basis: true_coords.append(coords[i]) maxima.append([true_coords,basis,best_point,coords]) # each entry has an np.array for the true coordinates and a basis for the flat directions, and the value at the minimum. Also need to keep a set of the entire coordinates for plotting ############################################################################################# ############################################################################################# ############################################################################################# ######################### Main programme follows ... three elements ######################### ############################################################################################# ############################################################################################# ############################################################################################# maxima=[] # initialize the set of local maxima for runs in range(0,factorings): domain = factor*domain origin = factor*origin scale_factors = np.array([p/10. for p in domain]) # converson factor to get to actual x,y,z from a 0-10 interval print ("origin=",origin,"Domain size=",domain) print (" ########################################### Doing the GA ######################################### ") GArun() # do the genetic algorithm print (" ########################################### Clustering ########################################### ") counters_refined=simplices(clusternum) #kmeans clustering print (" ####################################### Nelder Mead Algo ######################################### ") nelder(counters_refined) #Nelder-mead optimization -- appends any new minima found to maxima and maxima_coords print (" ############################################################################################# ") print (" ########################################### Results ######################################### ") print (" ############################################################################################# ") print (" ") no_minima=False if maxima==[]: no_minima=True else: # print "This is the final set of minima:", maxima print ("The unique minima are ") #for p in maxima: # print (p) # print (p[0]) minima = [maxima[0]] # take the first value for p in maxima: # normally would set add_entry = False straight away but possibly the coordinate length is completely new # so set it true first and set to False only if the length is found in the old minima add_entry = False len_found = False for q in minima: # test against previous list of minima for different length if (len(q[0]) == len(p[0])) : len_found=True if len_found==False: add_entry=True if add_entry == False: found=False for q in minima: # test against previous list of minima for any difference same = True if (len(q[0]) == len(p[0])): for i in range(0,len(p[0])): if (abs( (p[0][i] - q[0][i])/(q[0][i]+1.e-6) ) > 0.01 ) : same = False if same == True : found = True if found == False : add_entry = True if add_entry == True: minima.append(p) for p in minima: print ("Coordinates = ", p[0]," Flat directions = ", p[1], "Depth = ", p[2]) plottermax(minima) plot_minima(minima) #plot_minima_y_theta5(minima_coords,minima) print (" ") writeout(content2) if no_minima==False: writeout(text_minima) writeout(content3) output_file.close() view() if no_minima==True: print("No mimima found -- Sorry")