Add option to plot angle over mu_gamma

src/Plotq3-Angle.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 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
+# from subprocess import Popen, PIPE
+#import sys
+
+# unabhängig von alpha...
+def GetMuGamma(beta,theta,gamma,mu1,rho1, InputFilePath = os.path.dirname(os.getcwd()) +"/inputs/computeMuGamma.parset"):
+    # ------------------------------------ get mu_gamma ------------------------------
+    # ---Scenario 1.1: extreme regimes
+    if gamma == '0':
+        print('extreme regime: gamma = 0')
+        mu_gamma = (1.0/6.0)*arithmeticMean(mu1, beta, theta) # = q2
+        print("mu_gamma:", mu_gamma)
+    elif gamma == 'infinity':
+        print('extreme regime: gamma = infinity')
+        mu_gamma = (1.0/6.0)*harmonicMean(mu1, beta, theta)   # = q1
+        print("mu_gamma:", mu_gamma)
+    else:
+        # --- Scenario 1.2:  compute mu_gamma with 'Compute_MuGamma' (much faster than running full Cell-Problem)
+        # print("Run computeMuGamma for Gamma = ", gamma)
+        with open(InputFilePath, 'r') as file:
+            filedata = file.read()
+        filedata = re.sub('(?m)^gamma=.*','gamma='+str(gamma),filedata)
+        # filedata = re.sub('(?m)^alpha=.*','alpha='+str(alpha),filedata)
+        filedata = re.sub('(?m)^beta=.*','beta='+str(beta),filedata)
+        filedata = re.sub('(?m)^theta=.*','theta='+str(theta),filedata)
+        filedata = re.sub('(?m)^mu1=.*','mu1='+str(mu1),filedata)
+        filedata = re.sub('(?m)^rho1=.*','rho1='+str(rho1),filedata)
+        f = open(InputFilePath,'w')
+        f.write(filedata)
+        f.close()
+        # --- Run Cell-Problem
+
+        # Check Time
+        # t = time.time()
+        # subprocess.run(['./build-cmake/src/Cell-Problem', './inputs/cellsolver.parset'],
+        #                                      capture_output=True, text=True)
+        # --- Run Cell-Problem_muGama   -> faster
+        # subprocess.run(['./build-cmake/src/Cell-Problem_muGamma', './inputs/cellsolver.parset'],
+        #                                              capture_output=True, text=True)
+        # --- Run Compute_muGamma (2D Problem much much faster)
+
+        subprocess.run(['./build-cmake/src/Compute_MuGamma', './inputs/computeMuGamma.parset'],
+                                                             capture_output=True, text=True)
+        # print('elapsed time:', time.time() - t)
+
+        #Extract mu_gamma from Output-File                                           TODO: GENERALIZED THIS FOR QUANTITIES OF INTEREST
+        with open(OutputFilePath, 'r') as file:
+            output = file.read()
+        tmp = re.search(r'(?m)^mu_gamma=.*',output).group()                           # Not necessary for Intention of Program t output Minimizer etc.....
+        s = re.findall(r"[-+]?\d*\.\d+|\d+", tmp)
+        mu_gamma = float(s[0])
+        # print("mu_gamma:", mu_gammaValue)
+    # --------------------------------------------------------------------------------------
+    return mu_gamma
+
+
+
+
+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)
+
+
+#1. Define Inputs Gamma-Array..
+#2. for(i=0; i<length(array)) ..compute Q_hom, B_eff from Cell-Problem
+#3
+
+# matrix = np.loadtxt(path + 'Qmatrix.txt', usecols=range(3))
+# print(matrix)
+
+# Use Shell Commands:
+# subprocess.run('ls', shell=True)
+
+
+#---------------------------------------------------------------
+
+# -------------------------- Input Parameters --------------------
+mu1 = 10.0
+rho1 = 1.0
+alpha = 10   #1.05263158
+beta = 40.0  #5.0
+# theta = 1.0/4.0
+theta = 1.0/8.0  # 0.5
+
+# InterestingParameterSet :
+# mu1 = 10.0
+# rho1 = 1.0
+# alpha = 10
+# beta = 40.0
+# theta = 1.0/8.0
+
+
+#set gamma either to 1. '0' 2. 'infinity' or 3. a numerical positive value
+# gamma = '0'
+# gamma = 'infinity'
+# gamma = 0.5
+
+print('---- Input parameters: -----')
+print('mu1: ', mu1)
+print('rho1: ', rho1)
+print('alpha: ', alpha)
+print('beta: ', beta)
+print('theta: ', theta)
+# print('gamma:', gamma)
+print('----------------------------')
+# ----------------------------------------------------------------
+
+
+
+
+
+Gamma_Values = np.linspace(0.01, 5, num=15)    # TODO variable Input Parameters...alpha,beta...
+print('(Input) Gamma_Values:', Gamma_Values)
+# mu_gamma = []
+
+
+#
+# # --- Options
+# RUN = True
+# # RUN = False
+# # make_Plot = False
+make_Plot = True      # vll besser : Plot_muGamma
+#
+# if RUN:
+#     for gamma in Gamma_Values:
+#         print("Run Cell-Problem for Gamma = ", gamma)
+#         # print('gamma='+str(gamma))
+#         with open(InputFilePath, 'r') as file:
+#             filedata = file.read()
+#         filedata = re.sub('(?m)^gamma=.*','gamma='+str(gamma),filedata)
+#         f = open(InputFilePath,'w')
+#         f.write(filedata)
+#         f.close()
+#         # --- Run Cell-Problem
+#         t = time.time()
+#         # subprocess.run(['./build-cmake/src/Cell-Problem', './inputs/cellsolver.parset'],
+#         #                                      capture_output=True, text=True)
+#         # --- Run Cell-Problem_muGama   -> faster
+#         # subprocess.run(['./build-cmake/src/Cell-Problem_muGamma', './inputs/cellsolver.parset'],
+#         #                                              capture_output=True, text=True)
+#         # --- Run Compute_muGamma (2D Problem much much faster)
+#         subprocess.run(['./build-cmake/src/Compute_MuGamma', './inputs/computeMuGamma.parset'],
+#                                                              capture_output=True, text=True)
+#         print('elapsed time:', time.time() - t)
+#
+#         #Extract mu_gamma from Output-File                                           TODO: GENERALIZED THIS FOR QUANTITIES OF INTEREST
+#         with open(OutputFilePath, 'r') as file:
+#             output = file.read()
+#         tmp = re.search(r'(?m)^mu_gamma=.*',output).group()                           # Not necessary for Intention of Program t output Minimizer etc.....
+#         s = re.findall(r"[-+]?\d*\.\d+|\d+", tmp)
+#         mu_gammaValue = float(s[0])
+#         print("mu_gamma:", mu_gammaValue)
+#         mu_gamma.append(mu_gammaValue)
+#     # ------------end of for-loop -----------------
+#     print("(Output) Values of mu_gamma: ", mu_gamma)
+# # ----------------end of if-statement -------------
+#
+# # mu_gamma=[2.06099, 1.90567, 1.905]
+# # mu_gamma=[2.08306, 1.905, 1.90482, 1.90479, 1.90478, 1.90477]
+#
+# ##Gamma_Values = np.linspace(0.01, 20, num=12) :
+# #mu_gamma= [2.08306, 1.91108, 1.90648, 1.90554, 1.90521, 1.90505, 1.90496, 1.90491, 1.90487, 1.90485, 1.90483, 1.90482]
+#
+# ##Gamma_Values = np.linspace(0.01, 2.5, num=12)
+# # mu_gamma=[2.08306, 2.01137, 1.96113, 1.93772, 1.92592, 1.91937, 1.91541, 1.91286, 1.91112, 1.90988, 1.90897, 1.90828]
+#
+# # Gamma_Values = np.linspace(0.01, 2.5, num=6)
+# # mu_gamma=[2.08306, 1.95497, 1.92287, 1.91375, 1.9101, 1.90828]
+
+
+GetMuGammaVec = np.vectorize(GetMuGamma)
+muGammas = GetMuGammaVec(beta,theta,Gamma_Values,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)
+
+# Make Plot
+if make_Plot:
+    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()
+
+#
+# # ------------- RUN Matlab symbolic minimization program
+# eng = matlab.engine.start_matlab()
+# # s = eng.genpath(path + '/Matlab-Programs')
+# s = eng.genpath(path)
+# eng.addpath(s, nargout=0)
+# # print('current Matlab folder:', eng.pwd(nargout=1))
+# eng.cd('Matlab-Programs', nargout=0)  #switch to Matlab-Programs folder
+# # print('current Matlab folder:', eng.pwd(nargout=1))
+# Inp = False
+# print('Run symbolic Minimization...')
+# G, angle, type, kappa = eng.symMinimization(Inp,Inp,Inp,Inp, nargout=4)  #Name of program:symMinimization
+# # G, angle, type, kappa = eng.symMinimization(Inp,Inp,Inp,Inp,path + "/outputs", nargout=4)  #Optional: add Path
+# G = np.asarray(G) #cast Matlab Outout to numpy array
+#
+# # --- print Output ---
+# print('Minimizer G:')
+# print(G)
+# print('angle:', angle)
+# print('type:', type )
+# print('curvature:', kappa)