Skip to content
Snippets Groups Projects
PhaseDiagram_ContourSubPlots.py 20.5 KiB
Newer Older
  • Learn to ignore specific revisions
  • 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
    from matplotlib.ticker import MultipleLocator,FormatStrFormatter,MaxNLocator
    
    from mpl_toolkits.axes_grid1.inset_locator import inset_axes
    
    import time
    # print(sys.executable)
    
    # --------------------------------------------------------------------
    # START :
    # INPUT (Parameters):   alpha, beta, theta, gamma, mu1, rho1
    #
    # -Option 1 : (Case lambda = 0 => q12 = 0)
    #   compute q1,q2,b1,b2 from Formula
    #       Option 1.1 :
    #           set mu_gamma = 'q1' or 'q2' (extreme regimes: gamma \in {0,\infty})
    #       Option 1.2 :
    #           compute mu_gamma with 'Compute_MuGamma' (2D problem much faster then Cell-Problem)
    # -Option 2 :
    #   compute Q_hom & B_eff by running 'Cell-Problem'
    #
    # -> CLASSIFY ...
    #
    # 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
    #
    
    
    
    # ----------- SETUP PATHS
    # InputFile  = "/inputs/cellsolver.parset"
    # OutputFile = "/outputs/output.txt"
    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)
    
    
    # -------------------------- Input Parameters --------------------
    # mu1 = 10.0               # TODO : here must be the same values as in the Parset for computeMuGamma
    mu1 = 1.0
    rho1 = 1.0
    alpha = 2.0
    beta = 2.0
    # beta = 5.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 = 0.5
    # gamma = 0.25
    # gamma = 1.0
    
    # gamma = 5.0
    
    #added
    # lambda1 = 10.0
    lambda1 = 0.0
    
    #Test:
    # rho1 = -1.0
    
    
    
    print('---- Input parameters: -----')
    print('mu1: ', mu1)
    print('rho1: ', rho1)
    print('alpha: ', alpha)
    print('beta: ', beta)
    print('theta: ', theta)
    print('gamma:', gamma)
    
    print('lambda1: ', lambda1)
    print('----------------------------')
    # ----------------------------------------------------------------
    
    #
    # gamma_min = 0.5
    # gamma_max = 1.0
    #
    # # gamma_min = 1
    # # gamma_max = 1
    # Gamma_Values = np.linspace(gamma_min, gamma_max, num=3)
    # # #
    # # # Gamma_Values = np.linspace(gamma_min, gamma_max, num=13)    # TODO variable Input Parameters...alpha,beta...
    # print('(Input) Gamma_Values:', Gamma_Values)
    
    print('type of gamma:', type(gamma))
    # # #
    Gamma_Values = ['0', 'infinity']
    # Gamma_Values = ['infinity']
    # Gamma_Values = ['0']
    print('(Input) Gamma_Values:', Gamma_Values)
    
    for gamma in Gamma_Values:
    
        print('Run for gamma = ', gamma)
        print('type of gamma:', type(gamma))
            # muGamma = GetMuGamma(beta,theta,gamma,mu1,rho1,InputFilePath)
            # # muGamma = GetMuGamma(beta,theta,gamma,mu1,rho1)
            # print('Test MuGamma:', muGamma)
    
            # ------- Options --------
            # print_Cases = True
            # print_Output = True
    
                                #TODO
        # generalCase = True #Read Output from Cell-Problem instead of using Lemma1.4 (special case)
        generalCase = False
    
        # make_3D_plot = True
        # make_3D_PhaseDiagram = True
        make_2D_plot = False
        make_2D_PhaseDiagram = False
        make_3D_plot = False
        make_3D_PhaseDiagram = False
        make_2D_plot = True
        make_2D_PhaseDiagram = True
        #
    
        # --- Define effective quantities: q1, q2 , q3 = mu_gamma, q12 ---
        # q1 = harmonicMean(mu1, beta, theta)
        # q2 = arithmeticMean(mu1, beta, theta)
        # --- Set q12
        # q12 = 0.0  # (analytical example)              # TEST / TODO read from Cell-Output
    
    
    
    
    
        # b1 = prestrain_b1(rho1, beta, alpha, theta)
        # b2 = prestrain_b2(rho1, beta, alpha, theta)
        #
        # print('---- Input parameters: -----')
        # print('mu1: ', mu1)
        # print('rho1: ', rho1)
        # print('alpha: ', alpha)
        # print('beta: ', beta)
        # print('theta: ', theta)
        # print("q1: ", q1)
        # print("q2: ", q2)
        # print("mu_gamma: ", mu_gamma)
        # print("q12: ", q12)
        # print("b1: ", b1)
        # print("b2: ", b2)
        # print('----------------------------')
        # print("machine epsilon", sys.float_info.epsilon)
    
        # G, angle, type, kappa = classifyMin(q1, q2, mu_gamma, q12,  b1, b2, print_Cases, print_Output)
        # Test = f(1,2 ,q1,q2,mu_gamma,q12,b1,b2)
        # print("Test", Test)
    
        # ---------------------- 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 = 300 # Number of sample points in each direction
        # SamplePoints_3D = 150 # Number of sample points in each direction
        # SamplePoints_3D = 100 # Number of sample points in each direction
        # SamplePoints_3D = 200 # Number of sample points in each direction
        # SamplePoints_3D = 400 # Number of sample points in each direction
        # SamplePoints_2D = 7500 # Number of sample points in each direction
        # SamplePoints_2D = 4000 # 4000 # Number of sample points in each direction
    
        SamplePoints_2D = 400 # 4000  # Number of sample points in each direction
    
        # SamplePoints_2D = 500 # 4000    # Number of sample points in each direction
        # SamplePoints_2D = 100 # 4000  # Number of sample points in each direction
    
    Klaus Böhnlein's avatar
    Klaus Böhnlein committed
        # SamplePoints_2D = 2000 # 4000 # Number of sample points in each direction
    
        # SamplePoints_2D = 1000   # 4000 # Number of sample points in each direction
    
    
        if make_3D_PhaseDiagram:
            alphas_ = np.linspace(-20, 20, SamplePoints_3D)
            # alphas_ = np.linspace(-10, 10, SamplePoints_3D)
    
            # betas_  = np.linspace(0.01,40.01,SamplePoints_3D) # Full Range
            # betas_  = np.linspace(0.01,20.01,SamplePoints_3D) # FULL Range
    
    
    
            # betas_  = np.linspace(0.01,0.99,SamplePoints_3D)  # weird part
            betas_  = np.linspace(1.01,40.01,SamplePoints_3D)     #TEST !!!!!  For Beta <1 weird tings happen...
            thetas_ = np.linspace(0.01,0.99,SamplePoints_3D)
    
    
            # TEST
            # alphas_ = np.linspace(-2, 2, SamplePoints_3D)
            # betas_  = np.linspace(1.01,10.01,SamplePoints_3D)
            # print('betas:', betas_)
    
            # TEST :
            # alphas_ = np.linspace(-40, 40, SamplePoints_3D)
            # betas_  = np.linspace(0.01,80.01,SamplePoints_3D) # Full Range
    
            # 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)
    
            # Get MuGamma values ...
            GetMuGammaVec = np.vectorize(GetMuGamma)
            muGammas = GetMuGammaVec(betas, thetas, gamma, mu1, rho1)
            # Classify Minimizers....
            G, angles, Types, curvature = classifyMin_anaVec(alphas, betas, thetas, muGammas,  mu1, rho1)   # Sets q12 to zero!!!
    
            # G, angles, Types, curvature = classifyMin_anaVec(alphas, betas, thetas, muGammas,  mu1, rho1, True, True)
            # print('size of G:', G.shape)
            # print('G:', G)
    
            # Option to print angles
            # print('angles:', angles)
    
    
            # Out = classifyMin_anaVec(alphas,betas,thetas)
            # T = Out[2]
            # --- Write to VTK
    
            GammaString = str(gamma)
            VTKOutputName = "outputs/PhaseDiagram3D" + "Gamma" + GammaString
            gridToVTK(VTKOutputName , alphas, betas, thetas, pointData = {'Type': Types, 'angles': angles, 'curvature': curvature} )
            print('Written to VTK-File:', VTKOutputName )
    
        if make_2D_PhaseDiagram:
            # alphas_ = np.linspace(-20, 20, SamplePoints_2D)
            # alphas_ = np.linspace(0, 1, SamplePoints_2D)
            thetas_ = np.linspace(0.01,0.99,SamplePoints_2D)
            alphas_ = np.linspace(-5, 5, SamplePoints_2D)
            # alphas_ = np.linspace(-5, 15, SamplePoints_2D)
            # thetas_ = np.linspace(0.05,0.25,SamplePoints_2D)
    
    
            # good range:
            # alphas_ = np.linspace(9, 10, SamplePoints_2D)
            # thetas_ = np.linspace(0.075,0.14,SamplePoints_2D)
    
            # range used:
            # alphas_ = np.linspace(8, 10, SamplePoints_2D)
            # thetas_ = np.linspace(0.05,0.16,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;
            # betas_ = 10.0;
            # betas_ = 5.0;
            # betas_ = 0.5;
    
    
            #intermediate Values
            alphas_ = np.linspace(-2, 1, SamplePoints_2D)
            # thetas_ = np.linspace(0.4,0.6,SamplePoints_2D)
            # betas_ = 10.0;
    
            # TEST
            # alphas_ = np.linspace(-8, 8, SamplePoints_2D)
            # thetas_ = np.linspace(0.01,0.99,SamplePoints_2D)
            # betas_ = 1.0; #TEST Problem: disvison by zero if alpha = 9, theta = 0.1 !
            # betas_ = 0.9;
            # betas_ = 0.5;  #TEST!!!
            # alphas, betas, thetas = np.meshgrid(alphas_, betas_, thetas_, indexing='ij')
            betas = betas_
            alphas, thetas = np.meshgrid(alphas_, thetas_, indexing='ij')
    
            if generalCase:
                classifyMin_matVec = np.vectorize(classifyMin_mat)
                GetCellOutputVec = np.vectorize(GetCellOutput, otypes=[np.ndarray, np.ndarray])
                Q, B = GetCellOutputVec(alphas,betas,thetas,gamma,mu1,rho1,lambda1, InputFilePath ,OutputFilePath )
    
    
                # print('type of Q:', type(Q))
                # print('Q:', Q)
                G, angles, Types, curvature = classifyMin_matVec(Q,B)
    
            else:
                classifyMin_anaVec = np.vectorize(classifyMin_ana)
                GetMuGammaVec = np.vectorize(GetMuGamma)
                # 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!!!
                muGammas = GetMuGammaVec(betas,thetas,gamma,mu1,rho1,InputFilePath ,OutputFilePath )
    
                if gamma == '0':
                    G, angles_0, Types, curvature = classifyMin_anaVec(alphas,betas,thetas, muGammas,  mu1, rho1)    # Sets q12 to zero!!!
                if gamma == 'infinity':
                    G, angles_inf, 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)
                # Out = classifyMin_anaVec(alphas,betas,thetas)
                # T = Out[2]
                # --- Write to VTK
                # VTKOutputName = + path + "./PhaseDiagram2DNEW"
    
            # print('angles:',angles)
            # GammaString = str(gamma)
            # VTKOutputName = "outputs/PhaseDiagram2D" + "Gamma_" + GammaString
            # gridToVTK(VTKOutputName , alphas, betas, thetas, pointData = {'Type': Types, 'angles': angles, 'curvature': curvature} )
            # print('Written to VTK-File:', VTKOutputName )
    
    
    # --- Make 3D Scatter plot
    if(make_3D_plot or make_2D_plot):
        # fig = plt.figure()
        # ax = fig.add_subplot(111, projection='3d')
        # colors = cm.plasma(Types)
        colors = cm.coolwarm(angles_inf)
    
    
        width = 6.28
    
    Klaus Böhnlein's avatar
    Klaus Böhnlein committed
        # height = width / 1.618
        height = width / 2.5
    
        # fig, ax = plt.subplots()
        fig,ax = plt.subplots(nrows=1,ncols=2,figsize=(width,height), sharey=True)
        # ax = plt.axes((0.15,0.21 ,0.8,0.75))
    
    
    
    
        # if make_2D_plot: pnt3d=ax.scatter(alphas,thetas,c=Types.flat)
        # if make_3D_plot: pnt3d=ax.scatter(alphas,betas,thetas,c=Types.flat)
        #
        # if make_2D_plot: pnt3d=ax.scatter(alphas,thetas,c=angles.flat)
        # if make_3D_plot: pnt3d=ax.scatter(alphas,betas,thetas,c=angles.flat)
    
    
        # pnt=ax.scatter(alphas,thetas,c=angles,cmap='coolwarm')
        # # ax.colorbar()
        # CS = ax.contourf(alphas, thetas, angles,6, cmap=plt.cm.coolwarm, linestyle=dashed)
        # # CS = ax.contour(alphas, thetas, angles,6, colors='k')
        # ax.clabel(CS, inline=True, fontsize=7.5)
        # # ax.set_title('Simplest default with labels')
    
    
    Klaus Böhnlein's avatar
    Klaus Böhnlein committed
        # matplotlib.rcParams['contour.linestyles'] = 'dashed'
    
    
    
        CS_0 = ax[0].contourf(alphas, thetas, angles_0, 10, cmap=plt.cm.coolwarm)
        # CS = ax.contourf(alphas, thetas, angles, 10, cmap='RdBu')
        CS_02 = ax[0].contour(CS_0, levels=CS_0.levels[::2], colors='black',inline=True, linewidths=(0.5,))
        # ax.clabel(CS2, inline=True, fontsize=9, colors='black')
        # ax.clabel(CS2, inline=True, inline_spacing=3, rightside_up=True, colors='k', fontsize=8)
        # manual_locations = [
        #     (-0.5, 0.3), (-0.7, 0.4), (-0.8, 0.5), (-0.9, 0.6), (-1,0.7)]
        manual_locations = [
            (-0.4, 0.2),(-0.6, 0.3), (-0.7, 0.4), (-0.8, 0.5), (-0.9, 0.6), (-1,0.7)]
        # ax.clabel(CS2, inline=True, fontsize=6, colors='black', manual=manual_locations)
        # ax.clabel(CS2, inline=True, fontsize=6, colors='black')
        # ax.clabel(CS2, CS2.levels, inline=True, fontsize=10)
        # ax.clabel(CS,  fontsize=5, colors='black')
        # cbar = fig.colorbar(CS,label=r'angle $\alpha$', ticks=[0, np.pi/8, np.pi/4, 3*np.pi/8 , np.pi/2 ])
        # cbar = fig.colorbar(CS_0, ticks=[0, np.pi/2 ])
        # cbar.ax.set_yticklabels(['$0$', r'$\pi/2$'])
        # cbar.ax.set_title(r'angle $\alpha$')
    
    
        CS_1 = ax[1].contourf(alphas, thetas, angles_inf, 10, cmap=plt.cm.coolwarm)
        # CS = ax.contourf(alphas, thetas, angles, 10, cmap='RdBu')
        CS_12 = ax[1].contour(CS_1, levels=CS_1.levels[::2], colors='black',inline=True, linewidths=(0.5,))
        # ax.clabel(CS2, inline=True, fontsize=9, colors='black')
        # ax.clabel(CS2, inline=True, inline_spacing=3, rightside_up=True, colors='k', fontsize=8)
    
    Klaus Böhnlein's avatar
    Klaus Böhnlein committed
        # manual_locations = [
        #     (-0.5, 0.3), (-0.7, 0.4), (-0.8, 0.5), (-0.9, 0.6), (-1,0.7)]
    
        manual_locations = [
    
    Klaus Böhnlein's avatar
    Klaus Böhnlein committed
            (-0.5, 0.15), (-0.7, 0.3), (-0.8, 0.45), (-0.9, 0.6), (-1,0.75)]
        # ax[1].clabel(CS_12, inline=True, fontsize=10, colors='black', manual=manual_locations)
        ax[1].clabel(CS_12, inline=True, fontsize=8, colors='black', manual=manual_locations)
    
        # ax[1].clabel(CS_12, inline=True, fontsize=8, colors='black')
    
    
        axins1 = inset_axes(ax[1],
                           width="5%",  # width = 5% of parent_bbox width
                           height="100%",  # height : 50%
                           loc='lower left',
                           bbox_to_anchor=(1.05, 0., 1, 1),
                           bbox_transform=ax[1].transAxes,
                           borderpad=0,
                           )
    
        # ax.clabel(CS2, CS2.levels, inline=True, fontsize=10)
        # ax.clabel(CS,  fontsize=5, colors='black')
        # cbar = fig.colorbar(CS,label=r'angle $\alpha$', ticks=[0, np.pi/8, np.pi/4, 3*np.pi/8 , np.pi/2 ])
        # cbar = fig.colorbar(CS_1, ticks=[0, np.pi/8, np.pi/4, 3*np.pi/8 , np.pi/2 ])
    
        # cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
        cbar = fig.colorbar(CS_1, cax=axins1, ticks=[0, np.pi/8, np.pi/4, 3*np.pi/8 , np.pi/2 ])
        # cbar = fig.colorbar(CS_1, cax=cbar_ax, shrink=0.2, location='right', ticks=[0, np.pi/8, np.pi/4, 3*np.pi/8 , np.pi/2 ])
        # cbar = fig.colorbar(CS_1,  ax=ax[:], shrink=0.8, location='right', ticks=[0, np.pi/8, np.pi/4, 3*np.pi/8 , np.pi/2 ])
    
        cbar.ax.set_yticklabels(['$0$',r'$\pi/8$', r'$\pi/4$' ,r'$3\pi/8$' , r'$\pi/2$'])
    
    Klaus Böhnlein's avatar
    Klaus Böhnlein committed
        # cbar.ax.set_title(r'angle $\alpha$')
        cbar.ax.set_title(r'$\alpha$')
    
        # cbar=plt.colorbar(pnt3d)
        # cbar.set_label("Values (units)")
        # plt.axvline(x = 8, color = 'b', linestyle = ':', label='$q_1$')
        # plt.axhline(y = 0.083333333, color = 'b', linestyle = ':', label='$q_1$')
    
        ax[0].set_xlabel(r'$\theta_\rho$',fontsize=10)
        # ax[0].yaxis.set_major_locator(MultipleLocator(0.1))
        # ax[0].xaxis.set_major_locator(MultipleLocator(1))
        ax[0].yaxis.set_major_locator(MultipleLocator(0.1))
        ax[0].xaxis.set_major_locator(MultipleLocator(0.5))
        ax[0].set_ylabel(r'$\theta$   ',fontsize=10, rotation=0)
        ax[0].tick_params(axis='x', labelsize=7 )
        ax[0].tick_params(axis='y', labelsize=7 )
    
        ax[1].set_xlabel(r'$\theta_\rho$',fontsize=10)
        # ax.xaxis.set_minor_locator(MultipleLocator(0.5))
        # ax[1].yaxis.set_major_locator(MultipleLocator(0.1))
        # ax[1].xaxis.set_major_locator(MultipleLocator(1))
        ax[1].yaxis.set_major_locator(MultipleLocator(0.1))
        ax[1].xaxis.set_major_locator(MultipleLocator(0.5))
        ax[1].tick_params(axis='x', labelsize=7 )
        ax[1].tick_params(axis='y', labelsize=7 )
        # ax.set_ylabel('beta')
        # ax[1].set_ylabel(r'$\theta$   ',fontsize=10, rotation=0)
        # if make_3D_plot: ax.set_zlabel('theta')
        # plt.subplots_adjust(bottom=0.2)
        # plt.subplots_adjust(wspace=0.22, hspace=0.1)
        plt.subplots_adjust(hspace=0.15, wspace=0.1)
    
    Klaus Böhnlein's avatar
    Klaus Böhnlein committed
        plt.subplots_adjust(bottom=0.2)
    
        # fig.subplots_adjust(right=0.75)
    
    
        ax[0].grid( linestyle = '--', linewidth = 0.25)
        ax[1].grid( linestyle = '--', linewidth = 0.25)
    
    
    
        fig.set_size_inches(width, height)
        outputName = 'Plot-Contour_Gamma' +str(gamma) + '.pdf'
        fig.savefig(outputName)
        # fig.savefig('Plot-Contour.pdf')
        plt.show()
        # plt.savefig('common_labels.png', dpi=300)
        # print('T:', T)
        # print('Type 1 occured here:', np.where(T == 1))
        # print('Type 2 occured here:', np.where(T == 2))
    
    
        # print(alphas_)
        # print(betas_)
    
    
    
    
    
    # ALTERNATIVE
    # colors = ("red", "green", "blue")
    # groups = ("Type 1", "Type2", "Type3")
    #
    # # Create plot
    # fig = plt.figure()
    # ax = fig.add_subplot(1, 1, 1)
    #
    # for data, color, group in zip(Types, colors, groups):
    #     # x, y = data
    #     ax.scatter(alphas, thetas, alpha=0.8, c=color, edgecolors='none', label=group)
    #
    # plt.title('Matplot scatter plot')
    # plt.legend(loc=2)
    # plt.show()