import numpy as np import matplotlib.pyplot as plt import sympy as sym import math import os import subprocess import fileinput import re import matlab.engine from HelperFunctions import * from ClassifyMin import * import matplotlib.ticker as tickers import matplotlib as mpl from matplotlib.ticker import MultipleLocator,FormatStrFormatter,MaxNLocator import pandas as pd # import tikzplotlib # # from pylab import * # from tikzplotlib import save as tikz_save # Needed ? mpl.use('pdf') # from subprocess import Popen, PIPE #import sys ###################### makePlot.py ######################### # Generalized Plot-Script giving the option to define # quantity of interest and the parameter it depends on # to create a plot # # Input: Define y & x for "x-y plot" as Strings # - Run the 'Cell-Problem' for the different Parameter-Points # (alternatively run 'Compute_MuGamma' if quantity of interest # is q3=muGamma for a significant Speedup) ########################################################### # figsize argument takes inputs in inches # and we have the width of our document in pts. # To set the figure size we construct a function # to convert from pts to inches and to determine # an aesthetic figure height using the golden ratio: # def set_size(width, fraction=1): # """Set figure dimensions to avoid scaling in LaTeX. # # Parameters # ---------- # width: float # Document textwidth or columnwidth in pts # fraction: float, optional # Fraction of the width which you wish the figure to occupy # # Returns # ------- # fig_dim: tuple # Dimensions of figure in inches # """ # # Width of figure (in pts) # fig_width_pt = width * fraction # # # Convert from pt to inches # inches_per_pt = 1 / 72.27 # # # Golden ratio to set aesthetic figure height # # https://disq.us/p/2940ij3 # golden_ratio = (5**.5 - 1) / 2 # # # Figure width in inches # fig_width_in = fig_width_pt * inches_per_pt # # Figure height in inches # fig_height_in = fig_width_in * golden_ratio # # fig_dim = (fig_width_in, fig_height_in) # # return fig_dim # def format_func(value, tick_number): # # find number of multiples of pi/2 # N = int(np.round(2 * value / np.pi)) # if N == 0: # return "0" # elif N == 1: # return r"$\pi/2$" # elif N == 2: # return r"$\pi$" # elif N % 2 > 0: # return r"${0}\pi/2$".format(N) # else: # return r"${0}\pi$".format(N // 2) # find number of multiples of pi/2 N = int(np.round(4 * value / np.pi)) if N == 0: return "0" elif N == 1: return r"$\pi/4$" elif N == 2: return r"$\pi/2$" elif N % 2 > 0: return r"${0}\pi/2$".format(N) else: return r"${0}\pi$".format(N // 2) def find_nearest(array, value): array = np.asarray(array) idx = (np.abs(array - value)).argmin() return array[idx] def find_nearestIdx(array, value): array = np.asarray(array) idx = (np.abs(array - value)).argmin() return idx # TODO # - Fallunterscheidung (Speedup) falls gesuchter value mu_gamma = q3 # - Also Add option to plot Minimization Output # ----- Setup Paths ----- # InputFile = "/inputs/cellsolver.parset" # OutputFile = "/outputs/output.txt" InputFile = "/inputs/computeMuGamma.parset" OutputFile = "/outputs/outputMuGamma.txt" # path = os.getcwd() # InputFilePath = os.getcwd()+InputFile # OutputFilePath = os.getcwd()+OutputFile # --------- Run from src folder: path_parent = os.path.dirname(os.getcwd()) os.chdir(path_parent) path = os.getcwd() print(path) InputFilePath = os.getcwd()+InputFile OutputFilePath = os.getcwd()+OutputFile print("InputFilepath: ", InputFilePath) print("OutputFilepath: ", OutputFilePath) print("Path: ", path) #--------------------------------------------------------------- print('---- Input parameters: -----') mu1 = 1.0 #10.0 # lambda1 = 10.0 rho1 = 1.0 alpha = 5.0 beta = 10.0 # alpha = 2.0 # beta = 2.0 theta = 1.0/8.0 #1.0/4.0 lambda1 = 0.0 # gamma = 1.0/4.0 # TEST: alpha=3.0; # # INTERESTING!: alpha = 5.0 beta = 10.0 theta= 1/8 alpha_Values gamma = 'infinity' #Elliptic Setting # gamma = '0' #Hyperbolic Setting # gamma = 0.5 print('mu1: ', mu1) print('rho1: ', rho1) print('alpha: ', alpha) print('beta: ', beta) print('theta: ', theta) print('gamma:', gamma) print('----------------------------') # --- define Interval of x-va1ues: xmin = 0.01 xmax = 0.41 xmax = 0.99 Jumps = False numPoints = 2000 numPoints = 1000 X_Values = np.linspace(xmin, xmax, num=numPoints) print(X_Values) Y_Values = [] for theta in X_Values: print('Situation of Lemma1.4') q12 = 0.0 q1 = (1.0/6.0)*harmonicMean(mu1, beta, theta) q2 = (1.0/6.0)*arithmeticMean(mu1, beta, theta) b1 = prestrain_b1(rho1, beta, alpha,theta) b2 = prestrain_b2(rho1, beta, alpha,theta) b3 = 0.0 q3 = GetMuGamma(beta,theta,gamma,mu1,rho1,InputFilePath ,OutputFilePath) G, angle, Type, curvature = classifyMin_ana(alpha,beta,theta, q3, mu1, rho1) Y_Values.append(angle) print("(Output) Values of " + yName + ": ", Y_Values) idx = find_nearestIdx(Y_Values, 0) print(' Idx of value closest to 0', idx) ValueClose = Y_Values[idx] print('GammaValue(Idx) with mu_gamma closest to q_3^*', ValueClose) # Find Indices where the difference between the next one is larger than epsilon... jump_idx = [] jump_xValues = [] jump_yValues = [] tmp = X_Values[0] for idx, x in enumerate(X_Values): print(idx, x) if idx > 0: if abs(Y_Values[idx]-Y_Values[idx-1]) > 1: print('jump candidate') jump_idx.append(idx) jump_xValues.append(x) jump_yValues.append(Y_Values[idx]) print("Jump Indices", jump_idx) print("Jump X-values:", jump_xValues) print("Jump Y-values:", jump_yValues) y_plotValues = [Y_Values[0]] x_plotValues = [X_Values[0]] # y_plotValues.extend(jump_yValues) for i in jump_idx: y_plotValues.extend([Y_Values[i-1], Y_Values[i]]) x_plotValues.extend([X_Values[i-1], X_Values[i]]) y_plotValues.append(Y_Values[-1]) # x_plotValues = [X_Values[0]] # x_plotValues.extend(jump_xValues) x_plotValues.append(X_Values[-1]) print("y_plotValues:", y_plotValues) print("x_plotValues:", x_plotValues) # Y_Values[np.diff(y) >= 0.5] = np.nan #get values bigger than jump position # gamma = infty # x_rest = X_Values[X_Values>x_plotValues[1]] # Y_Values = np.array(Y_Values) #convert the np array # y_rest = Y_Values[X_Values>x_plotValues[1]] # # # # gamma = 0 # x_rest = X_Values[X_Values>x_plotValues[3]] # Y_Values = np.array(Y_Values) #convert the np array # y_rest = Y_Values[X_Values>x_plotValues[3]] # gamma between # Y_Values = np.array(Y_Values) #convert the np array # X_Values = np.array(X_Values) #convert the np array # # x_one = X_Values[X_Values>x_plotValues[3]] # # ax.scatter(X_Values, Y_Values) # y_rest = Y_Values[X_Values>x_plotValues[3]] # ax.plot(X_Values[X_Values>0.135], Y_Values[X_Values<0.135]) # # # # y_rest = Y_Values[np.nonzero(X_Values>x_plotValues[1]] # print('X_Values:', X_Values) # print('Y_Values:', Y_Values) # print('x_rest:', x_rest) # print('y_rest:', y_rest) # print('np.nonzero(X_Values>x_plotValues[1]', np.nonzero(X_Values>x_plotValues[1]) ) # --- Convert to numpy array Y_Values = np.array(Y_Values) X_Values = np.array(X_Values) # ---------------- Create Plot ------------------- #--- change plot style: SEABORN # plt.style.use("seaborn-paper") #--- Adjust gobal matplotlib variables # mpl.rcParams['pdf.fonttype'] = 42 # mpl.rcParams['ps.fonttype'] = 42 mpl.rcParams['text.usetex'] = True mpl.rcParams["font.family"] = "serif" mpl.rcParams["font.size"] = "9" # plt.rc('font', family='serif', serif='Times') # plt.rc('font', family='serif') # # plt.rc('text', usetex=True) #also works... # plt.rc('xtick', labelsize=8) # plt.rc('ytick', labelsize=8) # plt.rc('axes', labelsize=8) #---- Scale Figure apropriately to fit tex-File Width # width = 452.9679 # width as measured in inkscape width = 6.28 *0.5 height = width / 1.618 #setup canvas first fig = plt.figure() #main # fig, ax = plt.subplots() # fig, (ax, ax2) = plt.subplots(ncols=2) # fig,axes = plt.subplots(nrows=1,ncols=2,figsize=(width,height)) # more than one plot # fig.subplots_adjust(left=.15, bottom=.16, right=.99, top=.97) #TEST # TEST # mpl.rcParams['figure.figsize'] = (width+0.1,height+0.1) # fig = plt.figure(figsize=(width+0.1,height+0.1)) # mpl.rcParams['figure.figsize'] = (width,height) # fig = plt.figure(figsize=(10,6)) # default is [6.4,4.8] 6.4 is the width, 4.8 is the height # fig = plt.figure(figsize=(width,height)) # default is [6.4,4.8] 6.4 is the width, 4.8 is the height # fig = plt.figure(figsize=set_size(width)) # fig = plt.subplots(1, 1, figsize=set_size(width)) # --- To create a figure half the width of your document:# # fig = plt.figure(figsize=set_size(width, fraction=0.5)) #--- You must select the correct size of the plot in advance # fig.set_size_inches(3.54,3.54) ax = plt.axes((0.15,0.18,0.8,0.8)) # ax = plt.axes((0.1,0.1,0.5,0.8)) # ax = plt.axes((0.1,0.1,1,1)) # ax = plt.axes() # ax.spines['right'].set_visible(False) # ax.spines['left'].set_visible(False) # ax.spines['bottom'].set_visible(False) # ax.spines['top'].set_visible(False) # ax.tick_params(axis='x',which='major',direction='out',length=10,width=5,color='red',pad=15,labelsize=15,labelcolor='green', # labelrotation=15) # ax.tick_params(axis='x',which='major', direction='out',pad=5,labelsize=10) # ax.tick_params(axis='y',which='major', length=5, width=1, direction='out',pad=5,labelsize=10) ax.tick_params(axis='x',which='major', direction='out',pad=3) ax.tick_params(axis='y',which='major', length=3, width=1, direction='out',pad=3) ax.xaxis.set_major_locator(MultipleLocator(0.05)) ax.xaxis.set_minor_locator(MultipleLocator(0.025)) #---- print data-types print(ax.xaxis.get_major_locator()) print(ax.xaxis.get_minor_locator()) print(ax.xaxis.get_major_formatter()) print(ax.xaxis.get_minor_formatter()) #---- Hide Ticks or Labels # ax.yaxis.set_major_locator(plt.NullLocator()) # ax.xaxis.set_major_formatter(plt.NullFormatter()) #---- Reducing or Increasing the Number of Ticks # ax.xaxis.set_major_locator(plt.MaxNLocator(3)) # ax.yaxis.set_major_locator(plt.MaxNLocator(3)) #----- Fancy Tick Formats ax.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 4)) ax.yaxis.set_minor_locator(plt.MultipleLocator(np.pi / 12)) ax.yaxis.set_major_formatter(plt.FuncFormatter(format_func)) # --- manually change ticks&labels: # ax.set_xticks([0.2,1]) # ax.set_xticklabels(['pos1','pos2']) # ax.set_yticks([0, np.pi/8, np.pi/4 ]) # labels = ['$0$',r'$\pi/8$', r'$\pi/4$'] # ax.set_yticklabels(labels) a=ax.yaxis.get_major_locator() b=ax.yaxis.get_major_formatter() c = ax.get_xticks() d = ax.get_xticklabels() print('xticks:',c) print('xticklabels:',d) ax.grid(True,which='major',axis='both',alpha=0.3) # plt.figure() # f,ax=plt.subplots(1) # plt.title(r''+ yName + '-Plot') # plt.plot(X_Values, Y_Values,linewidth=2, '.k') # plt.plot(X_Values, Y_Values,'.k',markersize=1) # plt.plot(X_Values, Y_Values,'.',markersize=0.8) # plt.plot(X_Values, Y_Values) # ax.plot([[0],X_Values[-1]], [Y_Values[0],Y_Values[-1]]) # Gamma = '0' # ax.plot([x_plotValues[0],x_plotValues[1]], [y_plotValues[0],y_plotValues[1]] , 'b') # # ax.plot([x_plotValues[1],x_plotValues[3]], [y_plotValues[2],y_plotValues[3]] , 'b') # # ax.plot(x_rest, y_rest, 'b') # Gamma between # x jump values (gamma 0): [0.13606060606060608, 0.21090909090909093] # ax.plot([[0,jump_xValues[0]], [0, 0]] , 'b') # ax.plot([jump_xValues[0],xmin], [y_plotValues[2],y_plotValues[2]] , 'b') # ax.plot([[0,0.13606060606060608], [0, 0]] , 'b') # ax.plot([[0.13606060606060608,xmin], [(math.pi/2),(math.pi/2)]], 'b') # jump_xValues[0] # --- leave out jumps: # ax.scatter(X_Values, Y_Values) ax.set_xlabel(r"volume fraction $\theta$") ax.set_ylabel(r"angle $\alpha$") if Jumps: # --- leave out jumps: if gamma == 'infinity': ax.plot(X_Values[X_Values>=jump_xValues[0]], Y_Values[X_Values>=jump_xValues[0]] , 'royalblue') ax.plot(X_Values[X_Values<jump_xValues[0]], Y_Values[X_Values<jump_xValues[0]], 'royalblue') # ax.plot(X_Values[X_Values>=jump_xValues[0]], Y_Values[X_Values>=jump_xValues[0]]) # ax.plot(X_Values[X_Values<jump_xValues[0]], Y_Values[X_Values<jump_xValues[0]]) # ax.plot(X_Values[X_Values>0.136], Y_Values[X_Values>0.136]) # ax.plot(X_Values[X_Values<0.135], Y_Values[X_Values<0.135]) # ax.scatter(X_Values, Y_Values) # ax.plot(X_Values, Y_Values) # plt.plot(x_plotValues, y_plotValues,'.') # plt.scatter(X_Values, Y_Values, alpha=0.3) # plt.scatter(X_Values, Y_Values) # plt.plot(X_Values, Y_Values,'.') # plt.plot([X_Values[0],X_Values[-1]], [Y_Values[0],Y_Values[-1]]) # plt.axis([0, 6, 0, 20]) # ax.set_xlabel(r"volume fraction $\theta$", size=11) # ax.set_ylabel(r"angle $\angle$", size=11) # ax.set_xlabel(r"volume fraction $\theta$") # # ax.set_ylabel(r"angle $\angle$") # ax.set_ylabel(r"angle $\alpha$") # plt.ylabel('$\kappa$') # ax.yaxis.set_major_formatter(ticker.FormatStrFormatter('%g $\pi$')) # ax.yaxis.set_major_locator(ticker.MultipleLocator(base=0.1)) # Plot every other line.. not the jumps... if gamma == '0': tmp = 1 for idx, x in enumerate(x_plotValues): if idx > 0 and tmp == 1: # plt.plot([x_plotValues[idx-1],x_plotValues[idx]] ,[y_plotValues[idx-1],y_plotValues[idx]] ) ax.plot([x_plotValues[idx-1],x_plotValues[idx]] ,[y_plotValues[idx-1],y_plotValues[idx]], 'royalblue', zorder=2) tmp = 0 else: tmp = 1 # plt.plot([x_plotValues[0],x_plotValues[1]] ,[y_plotValues[0],y_plotValues[1]] ) # plt.plot([x_plotValues[2],x_plotValues[3]] ,[y_plotValues[2],y_plotValues[3]] ) # plt.plot([x_plotValues[4],x_plotValues[5]] ,[y_plotValues[4],y_plotValues[5]] ) # plt.plot([x_plotValues[6],x_plotValues[7]] ,[y_plotValues[6],y_plotValues[7]] ) for x in jump_xValues: plt.axvline(x,ymin=0, ymax= 1, color = 'orange',alpha=0.5, linestyle = 'dashed', linewidth=1, zorder=1) # plt.axvline(x,ymin=0, ymax= 1, color = 'orange',alpha=0.5, linestyle = 'dashed', label=r'$\theta_*$') # plt.axvline(x_plotValues[1],ymin=0, ymax= 1, color = 'g',alpha=0.5, linestyle = 'dashed') # plt.axhline(y = 1.90476, color = 'b', linestyle = ':', label='$q_1$') # plt.axhline(y = 2.08333, color = 'r', linestyle = 'dashed', label='$q_2$') # plt.legend() # -- SETUP LEGEND # ax.legend(prop={'size': 11}) # ax.legend() # ------------------ SAVE FIGURE # tikzplotlib.save("TesTout.tex") # plt.close() # mpl.rcParams.update(mpl.rcParamsDefault) # plt.savefig("graph.pdf", # #This is simple recomendation for publication plots # dpi=1000, # # Plot will be occupy a maximum of available space # bbox_inches='tight', # ) # plt.savefig("graph.pdf") # ---- ADD additional scatter: # ax.scatter(X_Values,Y_Values,s=1,c='black',zorder=4) # Find transition point lastIdx = len(Y_Values)-1 for idx, y in enumerate(Y_Values): if idx != lastIdx: if abs(y-0) < 0.01 and abs(Y_Values[idx+1] - 0) > 0.05: transition_point1 = X_Values[idx+1] print('transition point1:', transition_point1 ) if abs(y-0.5*np.pi) < 0.01 and abs(Y_Values[idx+1] -0.5*np.pi)>0.01: transition_point2 = X_Values[idx] print('transition point2:', transition_point2 ) if abs(y-0) > 0.01 and abs(Y_Values[idx+1] - 0) < 0.01: transition_point3 = X_Values[idx+1] print('transition point3:', transition_point3 ) # Add transition Points: if gamma == '0': ax.scatter([transition_point1, transition_point2],[np.pi/2,np.pi/2],s=6, marker='o', cmap=None, norm=None, facecolor = 'black', edgecolor = 'black', vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3) ax.text(transition_point1-0.02 , np.pi/2-0.02, r"$1$", size=6, bbox=dict(boxstyle="circle",facecolor='white', alpha=1.0, pad=0.1, linewidth=0.5) ) ax.text(transition_point2+0.012 , np.pi/2-0.02, r"$2$", size=6, bbox=dict(boxstyle="circle",facecolor='white', alpha=1.0, pad=0.1, linewidth=0.5) ) else: ax.scatter([transition_point1, transition_point2, transition_point3 ],[np.pi/2,np.pi/2,0 ],s=6, marker='o', cmap=None, norm=None, facecolor = 'black', edgecolor = 'black', vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3) ax.text(transition_point1-0.02 , np.pi/2-0.02, r"$1$", size=6, bbox=dict(boxstyle="circle",facecolor='white', alpha=1.0, pad=0.1, linewidth=0.5) ) ax.text(transition_point2 +0.011 , np.pi/2-0.02, r"$2$", size=6, bbox=dict(boxstyle="circle",facecolor='white', alpha=1.0, pad=0.1, linewidth=0.5) ) ax.text(transition_point3 +0.009 , 0+0.08, r"$3$", size=6, bbox=dict(boxstyle="circle",facecolor='white', alpha=1.0, pad=0.1, linewidth=0.5) ) else: ax.scatter(X_Values,Y_Values,s=1, marker='o', cmap=None, norm=None, facecolor = 'blue', edgecolor = 'none', vmin=None, vmax=None, alpha=None, linewidths=None, zorder=3) ax.set_yticks([0, np.pi/8, np.pi/4, 3*np.pi/8 , np.pi/2, 5*np.pi/8 ]) labels = ['$0$',r'$\pi/8$', r'$\pi/4$' ,r'$3\pi/8$' , r'$\pi/2$',r'$5\pi/8$'] ax.set_yticklabels(labels) # ax.set_yticks([1.570786327, np.pi/2 ]) # labels = [r'$\pi/2-0.0005 $' , r'$\pi/2$'] # ax.set_yticklabels(labels) fig.set_size_inches(width, height) fig.savefig('Plot-Angle-Theta.pdf') # tikz_save('someplot.tex', figureheight='5cm', figurewidth='9cm') # tikz_save('fig.tikz', # figureheight = '\\figureheight', # figurewidth = '\\figurewidth') # ---------------------------------------- plt.show() # #---------------------------------------------------------------