Rewrite Classification in Python & Add Mixed Term (Lemma1.6)

src/ClassifyMin.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 subprocess import Popen, PIPE
+
+
+
+def harmonicMean(mu_1,beta,theta):
+    return mu_1*(beta/(theta+(1-theta)*beta))
+
+
+def arithmeticMean(mu_1,beta,theta):
+    return mu_1*((1-theta)+theta*beta)
+
+
+def prestrain_b1(rho_1,beta,alpha):
+    return (3.0*rho_1/2.0)*beta*(1-(theta*(1+alpha)));
+
+def prestrain_b2(rho_1,beta,alpha):
+    return (3.0*rho_1/(4.0*((1.0-theta)+ theta*beta)))*(1-theta*(1+beta*alpha));
+
+
+# define function to be minimized
+def f(a1,a2,q1,q2,q3,q12,b1,b2):
+    A = np.array([[q1, q3 + q12/2.0], [q3 + q12/2.0, q2]])
+    print(A)
+    B = np.array([-2.0*q1*b1-q12*b2, -2.0*q2*b2-q12*b1])
+    print(B)
+    a = np.array([a1, a2])
+    print(a)
+    tmp = np.dot(A,a)
+    print(tmp)
+    tmp2 = np.dot(a,tmp)
+    print(tmp2)
+    tmpB = np.dot(B,a)
+    print(tmpB)
+    # print(q1*(b1**2))
+    # print(q2*(b2**2))
+    # print(q12*b1*b2)
+    return tmp2 + tmpB + q1*(b1**2) + q2*(b2**2) + q12*b1*b2
+
+
+# Classify Type of minimizer  1 = (I) , 2 = (II) , 3 = (III) , 4 = (IV)
+def classifyMin(q1, q2, q3, q12,  b1, b2,  print_output = "True"):
+    print("Run ClassifyMin...")
+    CaseCount = 0
+    epsilon = sys.float_info.epsilon
+
+    #TEST
+    # A = np.array([[1, 2], [3, 4]])
+    # kappa = 1.0
+    # angle = 1.57
+    # type = 2
+
+    determinant = q1*q2 - (q3**2 + 2*q3*q12 + q12**2)
+    print("determinant:", determinant)
+
+
+    if abs(determinant) < epsilon:
+        print('parabolic case (determinant equal zero)')
+
+
+    if determinant >= epsilon:
+        print('elliptic case (determinant greater zero)')
+        a1_star = -(2*(b1*(q12**2) + 2*b1*q3*q12 - 4*b1*q1*q2 + 4*b2*q2*q3))/(4*(q3**2) + 4*q3*q12 + (q12**2) - 4*q1*q2)
+        a2_star = -(2*(b2*(q12**2) + 2*b2*q3*q12 + 4*b1*q1*q3 - 4*b2*q1*q2))/(4*(q3**2) + 4*q3*q12 + (q12**2) - 4*q1*q2)
+        prod = a1_star*a2_star
+
+        b1_star = (2.0*q1*b1 + b2*q12)/(2*q1)
+        b2_star = (2.0*q2*b2 + b1*q12)/(2*q2)
+
+        if prod >= epsilon:
+            print('(E1) - inside Lambda ')
+            a1 = a1_star
+            a2 = a2_star
+            type = 3  # non-axial Minimizer
+            CaseCount += 1
+        if abs(prod) < epsilon:
+            print('(E2) - on the boundary of Lambda ')
+            # if q2*b2**2 < q1*b1**2: # Needs to be updated to Mixed Term!!! just define f as a function and check value?!
+            #check
+            if f(b1_star,0,q1,q2,mu_gamma,q12,b1,b2) < f(0, b2_star, q1,q2,mu_gamma,q12,b1,b2):
+                a1 = b1_star
+                a2 = 0.0
+                type = 1
+                CaseCount += 1
+            if f(b1_star,0,q1,q2,mu_gamma,q12,b1,b2) > f(0, b2_star, q1,q2,mu_gamma,q12,b1,b2):
+                a1 = 0
+                a2 = b2_star
+                type = 2
+                CaseCount += 1
+            if f(b1_star,0,q1,q2,mu_gamma,q12,b1,b2) = f(0, b2_star, q1,q2,mu_gamma,q12,b1,b2):
+                # Two Minimizers pick one
+                a1 = b1_star
+                a2 = 0.0
+                type = 4
+                CaseCount += 1
+
+        if prod < -1.0*epsilon:
+            print('(E3) - Outside Lambda ')
+            if f(b1_star,0,q1,q2,mu_gamma,q12,b1,b2) < f(0, b2_star, q1,q2,mu_gamma,q12,b1,b2):
+                a1 = b1_star
+                a2 = 0.0
+                type = 1
+                CaseCount += 1
+            if f(b1_star,0,q1,q2,mu_gamma,q12,b1,b2) > f(0, b2_star, q1,q2,mu_gamma,q12,b1,b2):
+                a1 = 0
+                a2 = b2_star
+                type = 2
+                CaseCount += 1
+            #TODO ...does type4 happen here? 
+
+
+
+
+    return A, angle ,type ,kappa
+
+
+
+
+
+# ---------------------------------------------- Main ---------------------
+print('Running Python Code')
+
+
+# --- Input Parameters ----
+
+mu_1 = 1.0
+rho_1 = 1.0
+alpha = 14.0
+beta = 20.0
+theta = 1.0/4.0
+
+
+# define q1, q2 , mu_gamma, q12
+# 1. read from Cell-Output
+# 2. define values from analytic formulas (expect for mu_gamma)
+q1 = harmonicMean(mu_1,beta,theta)
+q2 = arithmeticMean(mu_1,beta,theta)
+q12 = 0.0  #  (analytical example)
+
+# set mu_gamma  to value or read from Cell-Output
+mu_gamma = q1   # TODO read from Cell-Output
+b1 =  prestrain_b1(rho_1,beta,alpha)
+b2 =  prestrain_b2(rho_1,beta,alpha)
+
+
+print("q1: ", q1)
+print("q2: ", q2)
+print("mu_gamma : ", mu_gamma)
+print("b1: ", b1)
+print("b2: ", b2)
+
+print("machine epsilon", sys.float_info.epsilon)
+
+## ------- options --------
+# print_output = false;
+
+
+
+
+
+
+
+A, angle, type, kappa = classifyMin(q1, q2, mu_gamma, q12,  b1, b2)
+
+print("Matrix A")
+print(A)
+print("angle",angle)
+print("type",type)
+print("kappa", kappa)
+
+
+Test = f(1,2 ,q1,q2,mu_gamma,q12,b1,b2)
+print("Test", Test)
+
+
+
+# Run several times / make Array :
+# Gamma_Values = np.linspace(0.01, 2.5, num=12)
+# print(Gamma_Values)
+
+
+# Make Plot
+# plt.figure()
+# plt.title(r'$\mu_\gamma(\gamma)$-Plot')  # USE MATHEMATICAL EXPRESSIONS IN TITLE
+# plt.plot(Gamma_Values, mu_gamma)
+# plt.scatter(Gamma_Values, mu_gamma)
+# plt.xlabel("$\gamma$")
+# plt.ylabel("$\mu_\gamma$")
+# plt.legend()
+# plt.show()