From 26ec8c07edf89cf9cfea1798153a33d80b9fa978 Mon Sep 17 00:00:00 2001 From: Klaus <klaus.boehnlein@tu-dresden.de> Date: Fri, 17 Sep 2021 22:34:00 +0200 Subject: [PATCH] Update --- src/CellScript.py | 24 +++-- src/ClassifyMin.py | 28 +++++ src/PhaseDiagram.py | 163 ++++++++++++++++------------- src/Plotq3-Angle.py | 248 +++++++++++++++----------------------------- 4 files changed, 218 insertions(+), 245 deletions(-) diff --git a/src/CellScript.py b/src/CellScript.py index b7afd08b..73f1b046 100644 --- a/src/CellScript.py +++ b/src/CellScript.py @@ -17,6 +17,11 @@ import sys # ---------------------------------------------------------------------------------------------------------------------------- + + + + + # ----- Setup Paths ----- InputFile = "/inputs/cellsolver.parset" OutputFile = "/outputs/output.txt" @@ -50,14 +55,19 @@ print('----------------------------') +print('RunCellProblem...') +RunCellProblem(alpha,beta,theta,gamma,mu1,rho1,InputFilePath) + +print('Read effective quantities...') +Q, B = ReadEffectiveQuantities() +print('Q:', Q) +print('B:', B) + -# -# print('RunCellProblem...') -# RunCellProblem(alpha,beta,theta,gamma,mu1,rho1,InputFilePath) -# TEST: -print('Compare_Classification...') -Compare_Classification(alpha,beta,theta,gamma,mu1,rho1,InputFilePath) +# Compare symbolicMinimization with Classification 'ClassifyMin' : +# print('Compare_Classification...') +# Compare_Classification(alpha,beta,theta,gamma,mu1,rho1,InputFilePath) @@ -73,7 +83,7 @@ Inp = False Inp_T = True print('Run symbolic Minimization...') #Arguments: symMinization(print_Input,print_statPoint,print_Output,make_FunctionPlot, InputPath) -G, angle, type, kappa = eng.symMinimization(Inp_T,Inp,Inp,Inp, nargout=4) #Name of program:symMinimization +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 --- diff --git a/src/ClassifyMin.py b/src/ClassifyMin.py index 9f7e7d22..1f78c5ef 100644 --- a/src/ClassifyMin.py +++ b/src/ClassifyMin.py @@ -12,6 +12,24 @@ import sys # from subprocess import Popen, PIPE +# -------------------------------------------------- +# 'classifyMin' classifies Minimizers by utilizing the result of +# Lemma1.6 +# +# +# +# +# 'classifyMin_ana': (Special-Case : Lemma1.4) +# ..additionally assumes Poisson-ratio=0 => q12==0 +# +# +# +# Output : MinimizingMatrix, Angle, Type, Curvature + + + + + def harmonicMean(mu_1, beta, theta): return mu_1*(beta/(theta+(1-theta)*beta)) @@ -55,6 +73,9 @@ def classifyMin_ana(alpha,beta,theta,q3,mu_1,rho_1,print_Cases=False, print_Outp b2 = prestrain_b2(rho_1, beta, alpha,theta) return classifyMin(q1, q2, q3, q12, b1, b2, print_Cases, print_Output) + + + # -------------------------------------------------------------------- # Classify Type of minimizer 1 = R1 , 2 = R2 , 3 = R3 # before : destinction between which axis.. (4Types ) # where @@ -67,6 +88,13 @@ def classifyMin_ana(alpha,beta,theta,q3,mu_1,rho_1,print_Cases=False, print_Outp # R3 = E2 U E3 U P1.1 U P2 U H # ------------------------------------------------------------------- def classifyMin(q1, q2, q3, q12, b1, b2, print_Cases=False, print_Output=False): #ClassifyMin_hom? + # Assumption of Classification-Lemma1.6: + # 1. [b3 == 0] + # 2. Q is orthotropic i.e. q13 = q31 = q23 = q32 == 0 + + # TODO: check if Q is orthotropic here - assert() + + if print_Output: print("Run ClassifyMin...") CaseCount = 0 epsilon = sys.float_info.epsilon #Machine epsilon diff --git a/src/PhaseDiagram.py b/src/PhaseDiagram.py index 486380a8..35081eb7 100644 --- a/src/PhaseDiagram.py +++ b/src/PhaseDiagram.py @@ -9,11 +9,13 @@ 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 # print(sys.executable) @@ -34,67 +36,65 @@ import time # # OUTPUT: Minimizer G, angle , type, curvature # ----------------------------------------------------------------------- +# +# +# def GetMuGamma(beta,theta,gamma,mu1,rho1, InputFilePath = os.path.dirname(os.getcwd()) +"/inputs/computeMuGamma.parset", +# OutputFilePath = os.path.dirname(os.getcwd()) + "/outputs/outputMuGamma.txt" ): +# # ------------------------------------ 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 +# + -# 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 - - - - - -# --- SETUPS PATHS +# ----------- SETUP PATHS # InputFile = "/inputs/cellsolver.parset" # OutputFile = "/outputs/output.txt" - InputFile = "/inputs/computeMuGamma.parset" OutputFile = "/outputs/outputMuGamma.txt" # --------- Run from src folder: @@ -114,11 +114,14 @@ mu1 = 10.0 # TODO : here must be the same values as in the Parset rho1 = 1.0 alpha = 2.0 beta = 2.0 +# beta = 10.0 theta = 1.0/4.0 #set gamma either to 1. '0' 2. 'infinity' or 3. a numerical positive value gamma = '0' -# gamma = 'infinity' +gamma = 'infinity' # gamma = 0.5 +# gamma = 0.25 +# gamma = 1.0 print('---- Input parameters: -----') print('mu1: ', mu1) @@ -138,12 +141,14 @@ print('----------------------------') # print_Cases = True # print_Output = True + #TODO +# generalCase = False #Read Output from Cell-Problem instead of using Lemma1.4 (special case) -make_3D_plot = True +# make_3D_plot = True make_3D_PhaseDiagram = True make_2D_plot = False make_2D_PhaseDiagram = False -# make_3D_plot = False +make_3D_plot = False # make_3D_PhaseDiagram = False # make_2D_plot = True # make_2D_PhaseDiagram = True @@ -183,9 +188,9 @@ make_2D_PhaseDiagram = False # ---------------------- MAKE PLOT / Write to VTK------------------------------------------------------------------------------ -SamplePoints_3D = 10 # Number of sample points in each direction -SamplePoints_2D = 10 # Number of sample points in each direction -SamplePoints_3D = 20 # Number of sample points in each direction +# SamplePoints_3D = 10 # Number of sample points in each direction +# SamplePoints_2D = 10 # Number of sample points in each direction +SamplePoints_3D = 300 # Number of sample points in each direction SamplePoints_2D = 10 # Number of sample points in each direction @@ -216,23 +221,32 @@ if make_3D_PhaseDiagram: print('Written to VTK-File:', VTKOutputName ) if make_2D_PhaseDiagram: - alphas_ = np.linspace(-20, 20, SamplePoints_2D) + # alphas_ = np.linspace(-20, 20, SamplePoints_2D) + # thetas_ = np.linspace(0.01,0.99,SamplePoints_2D) + alphas_ = np.linspace(8, 12, SamplePoints_2D) + thetas_ = np.linspace(0.05,0.2,SamplePoints_2D) # betas_ = np.linspace(0.01,40.01,1) #fix to one value: betas_ = 2.0; - thetas_ = np.linspace(0.01,0.99,SamplePoints_2D) - # print('type of alphas', type(alphas_)) - # print('Test:', type(np.array([mu_gamma])) ) alphas, betas, thetas = np.meshgrid(alphas_, betas_, thetas_, indexing='ij') - classifyMin_anaVec = np.vectorize(classifyMin_ana) + # if generalCase: #TODO + # classifyMinVec = np.vectorize(classifyMin) + # GetCellOutputVec = np.vectorize(GetCellOutput) + # Q, B = GetCellOutputVec(alpha,betas,thetas,gamma,mu1,rho1,InputFilePath ,OutputFilePath ) + # + # print('type of Q:', type(Q)) + # print('Q:', Q) + # + # else: + classifyMin_anaVec = np.vectorize(classifyMin_ana) GetMuGammaVec = np.vectorize(GetMuGamma) - muGammas = GetMuGammaVec(betas,thetas,gamma,mu1,rho1) + muGammas = GetMuGammaVec(betas,thetas,gamma,mu1,rho1,InputFilePath ,OutputFilePath ) G, angles, Types, curvature = classifyMin_anaVec(alphas,betas,thetas, muGammas, mu1, rho1) # Sets q12 to zero!!! # print('size of G:', G.shape) # print('G:', G) - print('Types:', Types) + # print('Types:', Types) # Out = classifyMin_anaVec(alphas,betas,thetas) # T = Out[2] # --- Write to VTK @@ -265,8 +279,13 @@ if(make_3D_plot or make_2D_plot): # print('Type 2 occured here:', np.where(T == 2)) -print(alphas_) -print(betas_) +# print(alphas_) +# print(betas_) + + + + + # ALTERNATIVE # colors = ("red", "green", "blue") # groups = ("Type 1", "Type2", "Type3") diff --git a/src/Plotq3-Angle.py b/src/Plotq3-Angle.py index fe6b7df6..678b66bd 100644 --- a/src/Plotq3-Angle.py +++ b/src/Plotq3-Angle.py @@ -9,74 +9,29 @@ 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 -# 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 +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" -# 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) @@ -88,42 +43,29 @@ 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 : +print('---- Input parameters: -----') # mu1 = 10.0 # rho1 = 1.0 # alpha = 10 # beta = 40.0 -# theta = 1.0/8.0 +# theta = 0.125 +# +mu1 = 10.0 +rho1 = 1.0 +# alpha = 10.02021333 +alpha = 10.0 +beta = 2.0 +theta = 0.125 +# 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('---- Input parameters: -----') print('mu1: ', mu1) print('rho1: ', rho1) print('alpha: ', alpha) @@ -131,116 +73,90 @@ 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... +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 = [] -# # # --- 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] - +make_Plot = True +# Get values for mu_Gamma GetMuGammaVec = np.vectorize(GetMuGamma) -muGammas = GetMuGammaVec(beta,theta,Gamma_Values,mu1,rho1) +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) -classifyMin_anaVec = np.vectorize(classifyMin_ana) +# 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) + + + + + + # Make Plot if make_Plot: plt.figure() - plt.title(r'angle$-\mu_\gamma(\gamma)$-Plot') - plt.plot(muGammas, angles) - plt.scatter(muGammas, angles) + # 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 = 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.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^*$') + + + # Plot Gamma Value that is closest to q3_star + plt.axvline(x = gammaClose, color = 'b', linestyle = 'dashed', label='$\gamma^*$') + + plt.axvspan(gamma_min, gammaClose, color='red', alpha=0.5) + plt.axvspan(gammaClose, gamma_max, color='green', alpha=0.5) + + plt.xlabel("$\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) -- GitLab