Update Angle-plot

src/Plot-Angle-Gamma.py

+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
+import sys
+from ClassifyMin import *
+from HelperFunctions import *
+# from CellScript import *
+from mpl_toolkits.mplot3d import Axes3D
+import matplotlib.cm as cm
+from vtk.util import numpy_support
+from pyevtk.hl import gridToVTK
+import time
+import matplotlib.ticker as ticker
+
+import matplotlib as mpl
+from matplotlib.ticker import MultipleLocator,FormatStrFormatter,MaxNLocator
+import pandas as pd
+
+# from matplotlib import rc
+# rc('text', usetex=True) # Use LaTeX font
+#
+# import seaborn as sns
+# sns.set(color_codes=True)
+
+
+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)
+
+
+
+
+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
+
+
+InputFile  = "/inputs/computeMuGamma.parset"
+OutputFile = "/outputs/outputMuGamma.txt"
+# --------- 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 = 10.0
+# rho1 = 1.0
+# alpha = 2.56140350877193 #2.56140350877193, 4.0852130325814535
+# beta = 2.0  #5.0
+# theta = 1.0/4.0
+# theta = 1.0/8.0  # 0.5
+# theta = 0.075  # 0.5
+# mu1 = 10.0
+# rho1 = 1.0
+alpha = 10.0
+# beta = 40.0
+# theta = 0.125
+
+mu1 = 10.0
+rho1 = 1.0
+# alpha = 2.0
+beta = 2.0  #5.0
+# theta = 1.0/4.0
+theta = 1.0/8.0
+#
+
+
+
+# mu1 = 10.0
+# rho1 = 1.0
+# alpha = 10.0
+# beta = 2.0  #5.0
+# theta = 1.0/8.0
+# #
+
+
+# mu1 = 10.0
+# rho1 = 1.0
+# # alpha = 10.02021333
+# alpha = 10.0
+# beta = 2.0
+# theta = 0.125
+
+
+
+# mu1 = 10.0
+# rho1 = 1.0
+# # alpha = 10.02021333
+# alpha = 9.0
+# beta = 2.0
+# theta = 0.075
+
+#
+# mu1 = 10.0
+# rho1 = 1.0
+# alpha = 4.0
+# beta = 10.0
+# theta = 1.0/4.0
+
+# alpha = 10   #1.05263158
+# beta = 40.0  #5.0
+# # theta = 1.0/4.0
+# theta = 1.0/8.0  # 0.5
+#
+#
+# alpha = 2.0
+# beta = 2.0 #5.0
+# theta = 1/4.0
+# rho1 = 1.0
+# mu1 = 10.0
+
+# InterestingParameterSet :
+# mu1 = 10.0
+# rho1 = 1.0
+# alpha = 10
+# beta = 40.0
+
+
+# theta = 0.124242
+# gamma = 0.75
+#set gamma either to 1. '0' 2. 'infinity' or 3. a numerical positive value
+# gamma = '0'
+# gamma = 'infinity'
+# gamma = 0.5
+# gamma = 0.25
+
+print('mu1: ', mu1)
+print('rho1: ', rho1)
+print('alpha: ', alpha)
+print('beta: ', beta)
+print('theta: ', theta)
+# print('gamma:', gamma)
+print('----------------------------')
+
+# ----------------------------------------------------------------
+
+
+gamma_min = 0.01
+gamma_max = 1.0
+Gamma_Values = np.linspace(gamma_min, gamma_max, num=100)    # TODO variable Input Parameters...alpha,beta...
+print('(Input) Gamma_Values:', Gamma_Values)
+# mu_gamma = []
+
+# Gamma_Values = '0'
+
+# # --- Options
+# # make_Plot = False
+make_Plot = True
+
+# Get values for mu_Gamma
+GetMuGammaVec = np.vectorize(GetMuGamma)
+muGammas = GetMuGammaVec(beta,theta,Gamma_Values,mu1,rho1, InputFilePath ,OutputFilePath )
+print('muGammas:', muGammas)
+
+q12 = 0.0
+q1 = (1.0/6.0)*harmonicMean(mu1, beta, theta)
+q2 = (1.0/6.0)*arithmeticMean(mu1, beta, theta)
+print('q1: ', q1)
+print('q2: ', q2)
+b1 = prestrain_b1(rho1, beta, alpha,theta)
+b2 = prestrain_b2(rho1, beta, alpha,theta)
+q3_star = math.sqrt(q1*q2)
+print('q3_star:', q3_star)
+
+# TODO these have to be compatible with input parameters!!!
+# compute certain ParameterValues that this makes sense
+# b1 = q3_star
+# b2 = q1
+print('b1: ', b1)
+print('b2: ', b2)
+
+# return classifyMin(q1, q2, q3, q12,  b1, b2,  print_Cases, print_Output)
+
+
+
+# classifyMin_anaVec = np.vectorize(classifyMin_ana)
+# G, angles, Types, curvature = classifyMin_anaVec(alpha, beta, theta, muGammas,  mu1, rho1)
+classifyMin_anaVec = np.vectorize(classifyMin_ana)
+G, angles, Types, curvature = classifyMin_anaVec(alpha, beta, theta, muGammas,  mu1, rho1)
+
+# _,angles,_,_ = classifyMin_anaVec(alpha, beta, theta, muGammas,  mu1, rho1)
+
+print('angles:', angles)
+
+
+
+
+idx = find_nearestIdx(muGammas, q3_star)
+print('GammaValue Idx closest to q_3^*', idx)
+gammaClose = Gamma_Values[idx]
+print('GammaValue(Idx) with mu_gamma closest to q_3^*', gammaClose)
+
+
+
+determinantVec = np.vectorize(determinant)
+
+detValues = determinantVec(q1,q2,muGammas,q12)
+print('detValues:', detValues)
+
+
+detZeroidx = find_nearestIdx(detValues, 0)
+print('idx where det nearest to zero', idx)
+gammaClose = Gamma_Values[detZeroidx]
+print('gammaClose:', gammaClose)
+
+
+# Make Plot
+if make_Plot:
+    # plt.figure()
+    #
+    #
+    #
+    #
+    # # plt.title(r'angle$-\mu_\gamma(\gamma)$-Plot')
+    # # plt.title(r'angle$-\gamma$-Plot')
+    plt.plot(Gamma_Values, angles)
+    plt.scatter(Gamma_Values, angles)
+    # plt.plot(muGammas, angles)
+    # plt.scatter(muGammas, angles)
+    # # plt.axis([0, 6, 0, 20])
+    # # plt.axhline(y = 1.90476, color = 'b', linestyle = ':', label='$q_1$')
+    # # plt.axhline(y = 2.08333, color = 'r', linestyle = 'dashed', label='$q_2$')
+    # plt.axvline(x = q1, color = 'b', linestyle = ':', label='$q_1$')
+    # plt.axvline(x = q2, color = 'r', linestyle = 'dashed', label='$q_2$')
+    # # plt.axvline(x = q3_star, color = 'r', linestyle = 'dashed', label='$\gamma^*$')
+
+
+
+    f,ax=plt.subplots(1)
+    # ax.plot(muGammas, angles)
+    # ax.scatter(muGammas, angles)
+
+    ax.plot(Gamma_Values, angles, zorder=3)
+    ax.scatter(Gamma_Values, angles)
+
+
+    plt.xlabel("$q_3$")
+    plt.xlabel("$\gamma$")
+    plt.ylabel("angle")
+    ax.grid(True)
+
+
+    # ax.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 2))
+    # ax.yaxis.set_minor_locator(plt.MultipleLocator(np.pi / 12))
+
+    # 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))
+
+    # ax.yaxis.set_major_formatter(ticker.FormatStrFormatter('%g $\pi$'))
+    # ax.yaxis.set_major_locator(ticker.MultipleLocator(base=0.25))
+
+
+    # ax.yaxis.set_major_formatter(ticker.FuncFormatter(
+    # lambda val,pos: '{:.0g}$\pi$'.format(2*val/np.pi) if val !=0 else '0'))
+    # ax.yaxis.set_major_locator(ticker.MultipleLocator(base=0.5*np.pi))
+
+    # ---------------------------- show pi values ------------------------------------
+    # ax.axvline(x = q1, color = 'b', linestyle = ':', label='$q_1$')
+    # ax.axvline(x = q2, color = 'r', linestyle = 'dashed', label='$q_2$')
+    # ax.legend()
+    # # ax.set(xlim=(1.750, 1.880), ylim=(0, math.pi/2.0))
+    # ax.set(xlim=(1.760, 1.880), ylim=(-0.1, np.pi/4.0))
+    # # ax.set_yticks([0,  np.pi/4 ,np.pi/2])
+    # # labels = ['$0$', r'$\pi/4$', r'$\pi/2$']
+    # ax.set_yticks([0, np.pi/8, np.pi/4 ])
+    # labels = ['$0$',r'$\pi/8$', r'$\pi/4$']
+    # ax.set_yticklabels(labels)
+    # ---------------------------------------------------------------
+
+    ax.legend()
+    # ax.set(xlim=(1.750, 1.880), ylim=(0, math.pi/2.0))
+
+    # ax.set(xlim=(1.760, 1.880), ylim=(-0.1, np.pi/4.0))
+    # ax.set(xlim=(1.760, 1.880), ylim=(-0.1, np.pi/4.0))
+    # ax.set_yticks([0,  np.pi/4 ,np.pi/2])
+    # labels = ['$0$', r'$\pi/4$', r'$\pi/2$']
+    ax.set_yticks([0, np.pi/8, np.pi/4 ])
+    labels = ['$0$',r'$\pi/8$', r'$\pi/4$']
+    ax.set_yticklabels(labels)
+
+
+    # Plot Gamma Value that is closest to q3_star
+    ax.axvline(x = gammaClose, color = 'b', linestyle = 'dashed', label='$\gamma^*$')
+    ax.axvspan(gamma_min, gammaClose, color='red', alpha=0.5)
+    ax.axvspan(gammaClose, gamma_max, color='green', alpha=0.5)
+
+
+
+
+    # plt.xlabel("$q_3(\gamma)$")
+    plt.xlabel("$\gamma$")
+    plt.ylabel("angle")
+    plt.legend(loc='upper center')
+
+    plt.show()
+
+
+
+
+    # plt.figure()
+    # plt.title(r'angle$-\mu_\gamma(\gamma)$-Plot')
+    # plt.plot(muGammas, angles)
+    # plt.scatter(muGammas, angles)
+    # # plt.axis([0, 6, 0, 20])
+    # # plt.axhline(y = 1.90476, color = 'b', linestyle = ':', label='$q_1$')
+    # # plt.axhline(y = 2.08333, color = 'r', linestyle = 'dashed', label='$q_2$')
+    # plt.axvline(x = 1.90476, color = 'b', linestyle = ':', label='$q_1$')
+    # plt.axvline(x = 2.08333, color = 'r', linestyle = 'dashed', label='$q_2$')
+    # plt.xlabel("$\mu_\gamma$")
+    # plt.ylabel("angle")
+    # plt.legend()
+    # plt.show()
+    #

src/Plot_Angle_Lemma1.4.py

@@ -10,7 +10,11 @@ import matlab.engine
 from HelperFunctions import *
 from ClassifyMin import *
 
-import matplotlib.ticker as ticker
+import matplotlib.ticker as tickers
+import matplotlib as mpl
+from matplotlib.ticker import MultipleLocator,FormatStrFormatter,MaxNLocator
+import pandas as pd
+
 # from subprocess import Popen, PIPE
 #import sys
 
@@ -27,14 +31,26 @@ import matplotlib.ticker as ticker
 ###########################################################
 
 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(2 * value / np.pi))
+    N = int(np.round(4 * value / np.pi))
     if N == 0:
         return "0"
     elif N == 1:
-        return r"$\pi/2$"
+        return r"$\pi/4$"
     elif N == 2:
-        return r"$\pi$"
+        return r"$\pi/2$"
     elif N % 2 > 0:
         return r"${0}\pi/2$".format(N)
     else:
@@ -304,10 +320,86 @@ print("x_plotValues:", x_plotValues)
 # 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 -------------------
-plt.figure()
 
-f,ax=plt.subplots(1)
+#Adjust gobal matplotlib variables
+mpl.rcParams['pdf.fonttype'] = 42
+mpl.rcParams['ps.fonttype'] = 42
+# mpl.rcParams['font.family'] = 'Arial'
+
+
+#setup canvas first
+fig = plt.figure()
+mpl.rcParams['figure.figsize']
+fig = plt.figure(figsize=(10,6)) # default is [6.4,4.8] 6.4 is the width, 4.8 is the height
+
+ax = plt.axes((0.1,0.1,0.5,0.8))
+# 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.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.set_yticks([0, np.pi/8, np.pi/4 ])
+
+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')
@@ -340,14 +432,13 @@ f,ax=plt.subplots(1)
 
 # jump_xValues[0]
 
-Y_Values = np.array(Y_Values)  #convert the np array
-X_Values = np.array(X_Values)  #convert the np array
 
-# leave out jumps:
+
+# --- leave out jumps:
 # ax.scatter(X_Values, Y_Values)
 # leave out jumps:
-# ax.plot(X_Values[X_Values>jump_xValues[0]], Y_Values[X_Values>jump_xValues[0]] , 'b')
-# ax.plot(X_Values[X_Values<jump_xValues[0]], Y_Values[X_Values<jump_xValues[0]], 'b')
+ax.plot(X_Values[X_Values>=jump_xValues[0]], Y_Values[X_Values>=jump_xValues[0]] , 'blue')
+ax.plot(X_Values[X_Values<jump_xValues[0]], Y_Values[X_Values<jump_xValues[0]], 'blue')
 
 
 
@@ -355,8 +446,8 @@ X_Values = np.array(X_Values)  #convert the np array
 
 # 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)
+# 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)
@@ -365,19 +456,16 @@ ax.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])
 
-plt.xlabel(xName)
-plt.ylabel(yName)
+ax.set_xlabel(r"$\theta$")
+ax.set_ylabel('angle')
 
 # plt.ylabel('$\kappa$')
 
 # ax.yaxis.set_major_formatter(ticker.FormatStrFormatter('%g $\pi$'))
 # ax.yaxis.set_major_locator(ticker.MultipleLocator(base=0.1))
 
-ax.grid(True)
-ax.yaxis.set_major_locator(plt.MultipleLocator(np.pi / 2))
-ax.yaxis.set_minor_locator(plt.MultipleLocator(np.pi / 12))
 
-ax.yaxis.set_major_formatter(plt.FuncFormatter(format_func))
+
 
 # Plot every other line.. not the jumps...
 tmp = 1
@@ -396,12 +484,15 @@ tmp = 1
 
 
 for x in jump_xValues:
-    plt.axvline(x,ymin=0, ymax= 1, color = 'g',alpha=0.5, linestyle = 'dashed')
+    plt.axvline(x,ymin=0, ymax= 1, color = 'orange',alpha=0.5, linestyle = 'dashed',  label='$?$')
 
 # 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()
+
+ax.legend()
+
 plt.show()
 # #---------------------------------------------------------------