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 ClassifyMin_New 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

import seaborn as sns
import matplotlib.colors as mcolors

from chart_studio import plotly
import plotly.graph_objs as go
import plotly.express as px
import plotly.colors
# from matplotlib import rc
# rc('text', usetex=True) # Use LaTeX font
#
# import seaborn as sns
# sns.set(color_codes=True)


def show(fig):
    import io
    import plotly.io as pio
    from PIL import Image
    buf = io.BytesIO()
    pio.write_image(fig, buf)
    img = Image.open(buf)
    img.show()




# set the colormap and centre the colorbar
class MidpointNormalize(mcolors.Normalize):
	"""
	Normalise the colorbar so that diverging bars work there way either side from a prescribed midpoint value)

	e.g. im=ax1.imshow(array, norm=MidpointNormalize(midpoint=0.,vmin=-100, vmax=100))
	"""
	def __init__(self, vmin=None, vmax=None, midpoint=None, clip=False):
		self.midpoint = midpoint
		mcolors.Normalize.__init__(self, vmin, vmax, clip)

	def __call__(self, value, clip=None):
		# I'm ignoring masked values and all kinds of edge cases to make a
		# simple example...
		x, y = [self.vmin, self.midpoint, self.vmax], [0, 0.5, 1]
		return np.ma.masked_array(np.interp(value, x, y), np.isnan(value))



def set_size(width, fraction=1):
    """Set figure dimensions to avoid scaling in LaTeX.

    Parameters
    ----------
    width: float
            Document textwidth or columnwidth in pts
    fraction: float, optional
            Fraction of the width which you wish the figure to occupy

    Returns
    -------
    fig_dim: tuple
            Dimensions of figure in inches
    """
    # Width of figure (in pts)
    fig_width_pt = width * fraction

    # Convert from pt to inches
    inches_per_pt = 1 / 72.27

    # Golden ratio to set aesthetic figure height
    # https://disq.us/p/2940ij3
    golden_ratio = (5**.5 - 1) / 2

    # Figure width in inches
    fig_width_in = fig_width_pt * inches_per_pt
    # Figure height in inches
    fig_height_in = fig_width_in * golden_ratio

    fig_dim = (fig_width_in, fig_height_in)

    return fig_dim



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 == -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(4 * value / np.pi))
    if N == 0:
        return "0"
    elif N == 1:
        return r"$\pi/4$"
    elif N == -1:
        return r"$-\pi/4$"
    elif N == 2:
        return r"$\pi/2$"
    elif N == -2:
        return r"$-\pi/2$"
    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




def check_case(q1,q2,q12,q3,b1,b2):

    # a = np.array([a1,a2])
    b = np.array([b1,b2])
    H = np.array([[2*q1, q12+2*q3], [q12+2*q3,2*q2] ])
    A = np.array([[q1,1/2*q12], [1/2*q12,q2] ])


    print('det(H)=', np.linalg.det(H))

    # check if g* is in G^*_R^2
    tmp = A.dot(b)

    ## compute inverse of H :
    inv_H = np.linalg.inv(H)

    g_star = 2*inv_H.dot(tmp)

    print('g_star=', g_star)


    return g_star


def get_gstar(q1,q2,q12,q3,b1,b2):

    # a = np.array([a1,a2])
    b = np.array([b1,b2])
    H = np.array([[2*q1, q12+2*q3], [q12+2*q3,2*q2] ])
    A = np.array([[q1,1/2*q12], [1/2*q12,q2] ])
    print('det(H)=', np.linalg.det(H))
    # check if g* is in G^*_R^2
    tmp = A.dot(b)
    ## compute inverse of H :
    inv_H = np.linalg.inv(H)
    g_star = 2*inv_H.dot(tmp)
    # print('g_star=', g_star)

    return g_star


def determine_b(q1,q2,q12,q3,g_star):
    ##
    # Input: g_star
    # Output : b  such that g_star is minimizer, i.e A*b = (1/2)*H*g_star
    # q1=1;
    # q2=2;
    # q12=1/2;
    # q3=((4*q1*q2)**0.5-q12)/2;
    # H=[2*q1,q12+2*q3;q12+2*q3,2*q2];



    H = np.array([[2*q1, q12+2*q3], [q12+2*q3,2*q2] ])
    A = np.array([[q1,1/2*q12], [1/2*q12,q2] ])
    abar = np.array([q12+2*q3, 2*q2])
    rhs = (1/2)*H.dot(g_star)
    print('rhs:', rhs)
    b = np.linalg.lstsq(A, rhs)[0]
    print('b',b)
    b1=b[0]
    b2=b[1]

    return b
    ## ---------------


def get_minimizer(q1,q2,q3,b1,b2):

    # In the case if q12 == 0:
    quotient = (q1*q2-q3**2)
    g_star = np.array([(q1*q2*b1-q3*q2*b2)/quotient, (q1*q2*b2-q3*q1*b1)/quotient])
    print('g_star=', g_star)
    return g_star


def energy(a1,a2,q1,q2,q12,q3,b1,b2):


    a = np.array([a1,a2])
    b = np.array([b1,b2])
    H = np.array([[2*q1, q12+2*q3], [q12+2*q3,2*q2] ])
    A = np.array([[q1,1/2*q12], [1/2*q12,q2] ])


    tmp = H.dot(a)

    # print('H',H)
    # print('A',A)
    # print('b',b)
    # print('a',a)
    # print('tmp',tmp)

    tmp = (1/2)*a.dot(tmp)
    # print('tmp',tmp)

    tmp2 = A.dot(b)
    # print('tmp2',tmp2)
    tmp2 = 2*a.dot(tmp2)

    # print('tmp2',tmp2)
    energy = tmp - tmp2
    # print('energy',energy)


    # energy_axial1.append(energy_1)

    return energy



# def energy(a1,a2,q1,q2,q12,q3,b1,b2):
#
#
#     b = np.array([b1,b2])
#     H = np.array([[2*q1, q12+2*q3], [q12+2*q3,2*q2] ])
#     A = np.array([[q1,1/2*q12], [1/2*q12,q2] ])
#
#
#     tmp = H.dot(a)
#
#     print('H',H)
#     print('A',A)
#     print('b',b)
#     print('a',a)
#     print('tmp',tmp)
#
#     tmp = (1/2)*a.dot(tmp)
#     print('tmp',tmp)
#
#     tmp2 = A.dot(b)
#     print('tmp2',tmp2)
#     tmp2 = 2*a.dot(tmp2)
#
#     print('tmp2',tmp2)
#     energy = tmp - tmp2
#     print('energy',energy)
#
#
#     # energy_axial1.append(energy_1)
#
#     return energy
#


def draw_1ParameterLine(q1,q2,q12,q3):

    # q1=1;
    # q2=2;
    # q12=1/2;
    # q3=((4*q1*q2)**0.5-q12)/2;
    # H=[2*q1,q12+2*q3;q12+2*q3,2*q2];

    H = np.array([[2*q1, q12+2*q3], [q12+2*q3,2*q2] ])
    A = np.array([[q1,1/2*q12], [1/2*q12,q2] ])
    abar = np.array([q12+2*q3, 2*q2])
    abar = (1.0/math.sqrt((q12+2*q3)**2+(2*q2)**2))*abar

    print('abar:',abar)

    b = np.linalg.lstsq(A, abar)[0]
    print('b',b)

    b1=b[0]
    b2=b[1]


    # print('abar:',np.shape(abar))
    # print('np.transpose(abar):',np.shape(np.transpose(abar)))
    sstar = (1/(q1+q2))*abar.dot(A.dot(b))
    # sstar = (1/(q1+q2))*abar.dot(tmp)
    print('sstar', sstar)
    abarperp= np.array([abar[1],-abar[0]])
    print('abarperp:',abarperp)

    N=500;
    scale_domain = 5
    translate_startpoint = -1.8
    scale_domain = 4.65
    translate_startpoint = -1.85
    # --- FUll LINE:
    scale_domain = 4.5
    translate_startpoint = -1.4
    T_line = np.linspace(-sstar*(q12+2*q3)/(2*q2)*scale_domain + translate_startpoint, sstar*(2*q2)/(q12+2*q3)*scale_domain , num=N)
    line_values = []
    for t in T_line :
        print('sstar*abar+t*abarperp', sstar*abar+t*abarperp)
        line_values.append(sstar*abar+t*abarperp)
    # ------


    T = np.linspace(-sstar*(q12+2*q3)/(2*q2), sstar*(2*q2)/(q12+2*q3), num=N)
    # T = np.linspace(-sstar*(q12+2*q3)/(2*q2)*scale_domain + translate_startpoint, sstar*(2*q2)/(q12+2*q3)*scale_domain , num=N)
    # T = np.linspace(-2,2, num=N)
    # print('T:', T)

    print('T.min():', T.min())
    print('T.max():', T.max())

    kappas = []
    alphas = []
    # G.append(float(s[0]))

    G_container = []
    abar_container = []


    abar_tmp = abar

    for t in T :
        abar_current = sstar*abar+t*abarperp;
        abar_current[abar_current < 1e-10] = 0 # Project to x-y-axis!!
        current_energy = energy(abar_current[0],abar_current[1],q1,q2,q12,q3,b1,b2)
        print('energy along line:', current_energy )
        G = [abar_current[0], abar_current[1] , (2*abar_current[0]*abar_current[1])**0.5 ]
        G_container.append(G)
        abar_container.append(abar_current)
        e = [(abar_current[0]/(abar_current[0]+abar_current[1]))**0.5, (abar_current[1]/(abar_current[0]+abar_current[1]))**0.5]
        kappa = abar_current[0]+abar_current[1]
        alpha = math.atan2(e[1], e[0])
        kappas.append(kappa)
        alphas.append(alpha)


    alphas = np.array(alphas)
    kappas = np.array(kappas)

    # print('G_container', G_container)
    G = np.array(G_container)
    abar = np.array(abar_container)

    # print('G', G)
    # print('abar', abar)
    # print('abar.shape',abar.shape)
    line = ax.plot(abar[:,0],abar[:,1], linewidth=2, color='royalblue', linestyle='-', zorder=4)

    # Plot 1-Parameter Line
    line_values= np.array(line_values)
    ax.plot(line_values[:,0],line_values[:,1],'k--', linewidth=1,color='royalblue',alpha=0.8,zorder=4)

    return line




################################################################################################################
################################################################################################################
################################################################################################################

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: -----')

# q1=1;
# q2=2;
# q12=1/2;
# q3=((4*q1*q2)**0.5-q12)/2;
# # H=[2*q1,q12+2*q3;q12+2*q3,2*q2];
#
# H = np.array([[2*q1, q12+2*q3], [q12+2*q3,2*q2] ])
# A = np.array([[q1,1/2*q12], [1/2*q12,q2] ])
# abar = np.array([q12+2*q3, 2*q2])
# abar = (1.0/math.sqrt((q12+2*q3)**2+(2*q2)**2))*abar
#
# print('abar:',abar)
#
# b = np.linalg.lstsq(A, abar)[0]
# print('b',b)
#
#
# # print('abar:',np.shape(abar))
# # print('np.transpose(abar):',np.shape(np.transpose(abar)))
# sstar = (1/(q1+q2))*abar.dot(A.dot(b))
# # sstar = (1/(q1+q2))*abar.dot(tmp)
# print('sstar', sstar)
# abarperp= np.array([abar[1],-abar[0]])
# print('abarperp:',abarperp)


# -------------------------- Input Parameters --------------------
#
# mu1 = 1.0
# rho1 = 1.0
# alpha = 5.0
# theta = 1.0/2
# # theta= 0.1
# beta = 5.0
#



# mu1 = 1.0
# rho1 = 1.0
# alpha = -0.75
# theta = 1.0/2
# # theta= 0.1
# beta = 5.0


# mu1 = 1.0
# rho1 = 1.0
# alpha = 2.0
# theta = 1.0/2
# # theta= 0.1
# beta = 5.0


#Figure3:
# mu1 = 1.0
# rho1 = 1.0
# alpha = 2.0
# theta = 1.0/8
# # theta= 0.1
# beta = 2.0


# alpha= -5


#set gamma either to 1. '0' 2. 'infinity' or 3. a numerical positive value
# gamma = '0'
# # gamma = 'infinity'
#
#
# lambda1 = 0.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('----------------------------')
# ----------------------------------------------------------------
print('----------------------------')

# ----------------------------------------------------------------






# q1 = (1.0/6.0)*harmonicMean(mu1, beta, theta)
# q2 = (1.0/6.0)*arithmeticMean(mu1, beta, theta)
# q12 = 0.0
# q3 = GetMuGamma(beta, theta,gamma,mu1,rho1,InputFilePath ,OutputFilePath )
# b1 = prestrain_b1(rho1,beta, alpha, theta )
# b2 = prestrain_b2(rho1,beta, alpha, theta )
#


case = 'a'
# case = 'b'
# case = 'c'



if case == 'a':
    print('---- CASE (a) ----')
    q1 = 1
    q2 = 2
    q12 = 0
    q3 = 1
    g_star = np.array([0,1])
    g_star = np.array([-1,1])
    # ToDO
    b = determine_b(q1,q2,q12,q3, g_star)
    b1 = b[0]
    b2 = b[1]


if case == 'b':
    print('---- CASE (b) ----')
    q1 = 1
    q2 = 2
    q12 = 0
    q3 = 1
    g_star = np.array([1,1])
    # ToDO
    b = determine_b(q1,q2,q12,q3, g_star)
    b1 = b[0]
    b2 = b[1]


if case == 'c':
    print('---- CASE (c) ----')
    ## ---- 1-ParameterFamilyCase:
    q1=1;
    q2=2;
    q12=1/2;
    q3=((4*q1*q2)**0.5-q12)/2;
    # H=[2*q1,q12+2*q3;q12+2*q3,2*q2];
    H = np.array([[2*q1, q12+2*q3], [q12+2*q3,2*q2] ])
    A = np.array([[q1,1/2*q12], [1/2*q12,q2] ])
    abar = np.array([q12+2*q3, 2*q2])
    abar = (1.0/math.sqrt((q12+2*q3)**2+(2*q2)**2))*abar
    print('abar:',abar)
    b = np.linalg.lstsq(A, abar)[0]
    print('b',b)
    b1=b[0]
    b2=b[1]
    ## ---------------




print('q1 = ', q1)
print('q2 = ', q2)
print('q3 = ', q3)
print('q12 = ', q12)
print('b1 = ', b1)
print('b2 = ', b2)

## --- CHECK CASE ---


if case == 'a':
    # Check Minimizers:
    vec1 = np.array([(2*q1*b1+q12*b2)/(2*q1), 0])
    vec2 = np.array([0, (2*q12*b1+q2*b2)/(2*q2)])
    print('vec1:',vec1)
    print('vec2:',vec2)
    energy_1 = energy(vec1[0],vec1[1],q1,q2,q12,q3,b1,b2)
    energy_2 = energy(vec2[0],vec2[1],q1,q2,q12,q3,b1,b2)
    print('energy_1:',energy_1)
    print('energy_2:',energy_2)
    if energy_1 >= energy_2:
        minvec = vec2
    else:
        minvec = vec1
    print('minvec:', minvec)




check_case(q1,q2,q12,q3,b1,b2)

g_star = get_gstar(q1,q2,q12,q3,b1,b2)
print('g_star:',g_star)

# TEST
print('special case g_star:')
get_minimizer(q1,q2,q3,b1,b2)




## -------------------------------------

num_Points = 400
# num_Points = 600
# num_Points = 200
# num_Points = 100
# num_Points = 20


# Creating dataset
x = np.linspace(-5,5,num_Points)
y = np.linspace(-5,5,num_Points)

x = np.linspace(-20,20,num_Points)
y = np.linspace(-20,20,num_Points)


x = np.linspace(-2,2,num_Points)
y = np.linspace(-2,2,num_Points)
# x = np.linspace(-10,10,num_Points)
# y = np.linspace(-10,10,num_Points)

# x = np.linspace(-60,60,num_Points)
# y = np.linspace(-60,60,num_Points)
#
#
# x = np.linspace(-40,40,num_Points)
# y = np.linspace(-40,40,num_Points)


a1, a2 = np.meshgrid(x,y)

# geyser = sns.load_dataset("geyser")
# print('type of geyser:', type(geyser))
# print('geyser:',geyser)

# ContourRange=20
ContourRange=2

# x_in = np.linspace(-ContourRange,ContourRange,num_Points)
# y_in = np.linspace(-ContourRange,ContourRange,num_Points)
# a1_in, a2_in = np.meshgrid(x_in,y_in)

print('a1:', a1)
print('a2:',a2 )
print('a1.shape', a1.shape)

#-- FILTER OUT VALUES for G+ :

# tmp1 = a1[np.where(a1*a2 >= 0)]
# tmp2 = a2[np.where(a1*a2 >= 0)]
#
# np.take(a, np.where(a>100)[0], axis=0)
# tmp1 = np.take(a1, np.where(a1*a2 >= 0)[0], axis=0)
# tmp2 = np.take(a1, np.where(a1*a2 >= 0)[0], axis=0)
# tmp2 = a2[np.where(a1*a2 >= 0)]

tmp1 = a1[a1*a2 >= 0]
tmp2 = a2[a1*a2 >= 0]
tmp1 = tmp1.reshape(-1,5)
tmp2 = tmp2.reshape(-1,5)


# tmp1_pos = a1[np.where(a1*a2 >= 0) ]
# tmp2_pos = a2[np.where(a1*a2 >= 0) ]
# tmp1_pos = tmp1_pos[np.where(tmp1_pos  >= 0)]
# tmp2_pos = tmp2_pos[np.where(tmp2_pos  >= 0)]
#
# tmp1_neg = a1[a1*a2 >= 0 ]
# tmp2_neg = a2[a1*a2 >= 0 ]
# tmp1_neg = tmp1_neg[tmp1_neg  < 0]
# tmp2_neg = tmp2_neg[tmp2_neg  < 0]
# a1 = tmp1
# a2 = tmp2
#
# a1 = a1.reshape(-1,5)
# a2 = a2.reshape(-1,5)
#
# tmp1_pos = tmp1_pos.reshape(-1,5)
# tmp2_pos = tmp2_pos.reshape(-1,5)
# tmp1_neg = tmp1_neg.reshape(-1,5)
# tmp2_neg = tmp2_neg.reshape(-1,5)


# print('a1:', a1)
# print('a2:',a2 )
# print('a1.shape', a1.shape)





energyVec = np.vectorize(energy)

# Z = energyVec(np.array([a1,a2]),q1,q2,q12,q3,b1,b2)
Z = energyVec(a1,a2,q1,q2,q12,q3,b1,b2)

# Z_in = energyVec(a1_in,a2_in,q1,q2,q12,q3,b1,b2)




print('Z:', Z)
print('any', np.any(Z<0))







# negZ_a1 = a1[np.where(Z<0)]
# negZ_a2 = a2[np.where(Z<0)]
# negativeValues = Z[np.where(Z<0)]
# print('negativeValues:',negativeValues)
#
# print('negZ_a1',negZ_a1)
# print('negZ_a2',negZ_a2)
#
#
# negZ_a1 = negZ_a1.reshape(-1,5)
# negZ_a2 = negZ_a2.reshape(-1,5)
# negativeValues = negativeValues.reshape(-1,5)
#
# Z_pos =  energyVec(tmp1_pos,tmp2_pos,q1,q2,q12,q3,b1,b2)
# Z_neg = energyVec(tmp1_neg,tmp2_neg,q1,q2,q12,q3,b1,b2)





# print('Test energy:' , energy(np.array([1,1]),q1,q2,q12,q3,b1,b2))




# print('Z_pos.shape', Z_pos.shape)







## -- PLOT :
# mpl.rcParams['text.usetex'] = True
# mpl.rcParams["font.family"] = "serif"
# mpl.rcParams["font.size"] = "9"
#
# label_size = 8
# mpl.rcParams['xtick.labelsize'] = label_size
# mpl.rcParams['ytick.labelsize'] = label_size
#
# label_size = 9
# mpl.rcParams['xtick.labelsize'] = label_size
# mpl.rcParams['ytick.labelsize'] = label_size

label_size = 10
#--- change plot style:  SEABORN
# plt.style.use("seaborn-paper")

# plt.style.use("seaborn-darkgrid")
# plt.style.use("seaborn-whitegrid")
# plt.style.use("seaborn")
# plt.style.use("seaborn-paper")
# plt.style.use('ggplot')
# plt.rcParams["font.family"] = "Avenir"
# plt.rcParams["font.size"] = 16

# plt.style.use("seaborn-darkgrid")
mpl.rcParams['text.usetex'] = True
mpl.rcParams["font.family"] = "serif"
mpl.rcParams["font.size"] = "10"
# mpl.rcParams['xtick.labelsize'] = 16mpl.rcParams['xtick.major.size'] = 2.5
# mpl.rcParams['xtick.bottom'] = True
# mpl.rcParams['ticks'] = True
mpl.rcParams['xtick.bottom'] = True
mpl.rcParams['xtick.major.size'] = 3
mpl.rcParams['xtick.minor.size'] = 1.5
mpl.rcParams['xtick.major.width'] = 0.75
mpl.rcParams['ytick.left'] = True
mpl.rcParams['ytick.major.size'] = 3
mpl.rcParams['ytick.minor.size'] = 1.5
mpl.rcParams['ytick.major.width'] = 0.75

mpl.rcParams.update({'font.size': 10})
mpl.rcParams['axes.labelpad'] = 3.0

#--- Adjust gobal matplotlib variables
# mpl.rcParams['pdf.fonttype'] = 42
# mpl.rcParams['ps.fonttype'] = 42
mpl.rcParams['text.usetex'] = True
mpl.rcParams["font.family"] = "serif"
mpl.rcParams["font.size"] = "10"




width = 6.28 *0.5
# width = 6.28 *0.33
# width = 6.28
height = width #/ 1.618

# width = 452.9579/2
# size= set_size(width, fraction=0.5)
# print('set_size(width, fraction=0.5)', set_size(width, fraction=1))
# print('size[0]',size[0])

f_size = 10
# fig= plt.figure()
fig, ax = plt.subplots()
# fig.set_size_inches(width, height)
# fig.set_size_inches(set_size(width, fraction=0.5))
# ax = plt.axes(projection ='3d', adjustable='box')
# ax = plt.axes((0.17,0.21 ,0.75,0.75))

# ax = plt.axes((0.17,0.23 ,0.7,0.7))
# ax = plt.axes((0.17,0.23 ,1.0,1.0))
# ax=plt.axes()
# ax = plt.axes((0.15,0.18,0.8,0.8))
# ax.tick_params(axis='x',which='major', direction='out',pad=5)
# ax.tick_params(axis='y',which='major', length=3, width=1, direction='out',pad=3)
# ax.xaxis.set_major_locator(MultipleLocator(0.1))
# ax.xaxis.set_minor_locator(MultipleLocator(0.05))
# ax.xaxis.set_major_locator(plt.MultipleLocator(np.pi / 8))
# ax.xaxis.set_minor_locator(plt.MultipleLocator(np.pi / 16))
# ax.xaxis.set_major_locator(plt.MultipleLocator(np.pi / 2))
# ax.xaxis.set_minor_locator(plt.MultipleLocator(np.pi / 4))
# ax.xaxis.set_major_formatter(plt.FuncFormatter(format_func))
# ax.grid(True,which='major',axis='both',alpha=0.3)

# colorfunction=(B*kappa)
# print('colofunction',colorfunction)

#translate Data
# Z = Z - (Z.max()-Z.min())/2
# Z = Z - 50
# Z = Z - 500
#
# Z = Z.T


# Substract constant:
# c = (b1**2)*q1+b1*b2*q12+(b2**2)*q2
# Z = Z-c
#
# print('Value of c:', c)



print('Z.min()', Z.min())
print('Z.max()', Z.max())
norm=mcolors.Normalize(Z.min(),Z.max())
# facecolors=cm.brg(norm)



print('norm:', norm)
print('type of norm', type(norm))
print('norm(0):', norm(0))
print('norm(Z):', norm(Z))

# ax.plot(theta_rho, theta_values, 'royalblue', zorder=3, )

# ax.scatter(a1,a2, s=0.5)

# ax.scatter(tmp1_pos,tmp2_pos, s=0.5)
# ax.scatter(tmp1_neg,tmp2_neg, s=0.5)

# CS = ax.contour(a1, a2, Z,10, cmap=plt.cm.gnuplot, levels=100 )
# CS = ax.contour(a1, a2, Z,10, cmap=plt.cm.gnuplot,  levels=20 )


# sns.kdeplot(np.array([a1, a2, Z]))
# sns.kdeplot(tmp1_pos,tmp2_pos,Z_pos)

# levels = [-5.0, -4, -3, 0.0, 1.5, 2.5, 3.5]
# CS = ax.contour(a1, a2, Z,10, cmap=plt.cm.gnuplot, corner_mask=True,levels=levels)
# CS = ax.contour(a1, a2, Z, cmap=plt.cm.gnuplot(norm(Z)), corner_mask=True)
# CS = ax.contour(a1, a2, Z, cm.brg(norm(Z)), levels=20)
# CS = ax.contour(a1, a2, Z, cmap=plt.cm.gnuplot, levels=20)
# CS = ax.contour(a1, a2, Z, colors='k', levels=14,  linewidths=(0.5,))
# CS = ax.contour(a1, a2, Z, colors='k', levels=18,  linewidths=(0.5,))

# ax.contour(negZ_a1, negZ_a2, negativeValues, colors='k',  linewidths=(0.5,))
# CS =  ax.contour(a1_in, a2_in, Z_in, colors='k',  linewidths=(0.75) ,zorder=5)


step = (-1)*Z.min()
# step = 3.0

levels = np.arange(Z.min(),Z.max()/2,step)
# evels = np.arange(-4, 25, 4)
# levels = np.arange(-3, 23, 4)

# levels = np.arange(Z.min(),Z.max()/12,step)
# levels = np.arange(Z.min(),Z.max(),step)
# levels = np.arange(0,Z.max()/2,200)
# CS =  ax.contour(a1_in, a2_in, Z_in, colors='k',  linewidths=(0.75), levels=levels ,zorder=5)

# CS =  ax.contour(a1, a2, Z, colors='k',  linewidths=(0.75), levels=levels ,zorder=5)
# CS =  ax.contour(a1, a2, Z,5,  colors='k',  linewidths=(0.75) ,zorder=5, vmin=Z.min(), vmax=1)

# CS =  ax.contour(a1, a2, Z,5,  colors='k',  linewidths=(0.75) ,zorder=5, vmin=Z.min(), vmax=1)

if case == 'b':
    levels = np.arange(-4, 25, 4)


if case == 'a':
    levels = np.arange(-1, 25, 5)  # if g_star = [0 1]
    levels = np.arange(0, 24, 4)  # if g_star = [-1 1]
if case == 'c':
    levels = np.arange(0, 25, 4)

CS =  ax.contour(a1, a2, Z,levels,  colors='k',linewidths=(0.75), extent=(-2, 2, -2, 2), zorder=5)


# CS =  ax.contour(a1 , a2 , Z, colors='k',  linewidths=(0.75), zorder=5)
# CS =  ax.contour(a1_in, a2_in, Z_in, colors='k',  linewidths=(1) , vmin= Z.min()+0.04, zorder=5)
# CS =  ax.contour(a1, a2, Z, colors='k',  linewidths=(1) , vmin= Z.min()+0.04, zorder=5)



# df = pd.DataFrame(data=Z_in, columns=a1_in, index=a2_in)
# df2 = pd.DataFrame(np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]),
#                    columns=['a', 'b', 'c'])

# sns.kdeplot(data=df2, x="waiting", y="duration")
# sns.kdeplot(data=df2)

# CS = ax.contour(a1, a2, Z, colors='k',   linewidths=(0.5,))

# CS = ax.contour(a1, a2, Z,10, cmap=plt.cm.gnuplot, extend='both', levels=50)
# CS = ax.contourf(a1, a2, Z,10, colors='k', extend='both', levels=50)
# CS = ax.contourf(a1, a2, Z,10, colors='k')
#
# # CS = ax.contour(tmp1_pos,tmp2_pos, Z_pos,10, cmap=plt.cm.gnuplot, levels=10 )
# # CS = ax.contour(tmp1_pos,tmp2_pos, Z_pos,10, cmap=plt.cm.gnuplot, corner_mask=True)
#
# CS = ax.contour(a1, a2, Z,10,  colors = 'k')
# ax.clabel(CS, inline=True, fontsize=4)

# ax.clabel(CS, inline=True,  fontsize=4)


# 1- ParameterFAMILY CASE :
# manual_locations = [
#                 (1, 1), (-5, -5), (-10, -10 ),(-12.5,-12.5),(-14,-14), (-15,-15),
#                 (5, 5), (10, 10 ),(12.5,12.5), (15,15), (17,17)]


if case == 'c':
    manual_locations = [
                    (1, 1), (0.5, 0.5),(1.5,1.5),(1.75,1.75),(-0.2,0.2), (-0.5,-0.5),
                    (-1,-1),(-1.25,-1.25), (-1.5,-1.5)]
    ax.clabel(CS, inline=True, fontsize=9, colors='black', manual=manual_locations)
if case == 'b':
    manual_locations = [
                    (0.5, 0.5),(-0.2,0.2),(-0.2,-0.2),(-0.5,-0.5), (-0.75,-0.75),
                    (-1,-1),(-1.25,-1.25), (-1.5,-1.5)]
    ax.clabel(CS, inline=True, fontsize=9, colors='black', manual=manual_locations)
if case == 'a':
    manual_locations = [
                  (1.8,1.8),(1.5,1.5) , (1, 1.5),(-1,0.25),(-0.75,-0.75),
                    (-1,-1),(-1.25,-1.25), (-1.5,-1.5)]
    ax.clabel(CS, inline=True, fontsize=9, colors='black', manual=manual_locations)
# else :
#     ax.clabel(CS, inline=True, fontsize=9, colors='black', zorder=5)

# GAMMA = inf
# manual_locations = [
#                 (1, 1), (-5, -5), (-10, -10 ), (-15,-15),
#                 (5, 5), (10, 10 ),(12.5,12.5), (15,15), (17,17)]

#
#
# # GAMMA = 0
# manual_locations = [
#                 (1, 1), (-5, -5), (-10, -10 ), (-15,-15), (-15,7),
#                 (5, 5), (10, 10 ), (15,15), (17,17)]


# ax.clabel(CS, inline=True, fontsize=f_size, colors='black')
# ax.clabel(CS, inline=True, fontsize=9, colors='black', zorder=5)

# cmap = cm.brg(norm(Z))
#
# C_map = cm.inferno(norm(Z))

# ax.imshow(Z, cmap=C_map, extent=[-20, 20, -20, 20], origin='lower', alpha=0.5)

# ax.imshow(norm(Z), extent=[-20, 20, -20, 20], origin='lower',
#                   cmap='bwr', alpha=0.8)

# ax.imshow(norm(Z), extent=[-20, 20, -20, 20],origin='lower', vmin=Z.min(), vmax=Z.max(),
#                   cmap='bwr', alpha=0.6)

# ax.imshow(norm(Z), extent=[-20, 20, -20, 20],origin='lower', norm = norm,
#                   cmap='coolwarm', alpha=0.6)


cmap=mpl.cm.RdBu_r
cmap=mpl.cm.viridis_r
# cmap=mpl.cm.bwr
cmap=mpl.cm.coolwarm
# cmap=mpl.cm.gnuplot
cmap=mpl.cm.magma_r
# cmap=mpl.cm.inferno_r
# cmap=mpl.cm.plasma
# cmap=mpl.cm.plasma_r
# cmap=mpl.cm.cividis_r
# cmap = mpl.colors.LinearSegmentedColormap.from_list("", ["blue","violet","red"])
# cmap = mpl.colors.LinearSegmentedColormap.from_list("", ["blue","orange"])

# cmap = mpl.colors.LinearSegmentedColormap.from_list("", [(0,"red"), (.1,"violet"), (.5, "blue"), (1.0, "green")])
# make a colormap that has land and ocean clearly delineated and of the
# same length (256 + 256)
#
# colors_undersea = plt.cm.terrain(np.linspace(0, 0.17, 256))
# colors_land = plt.cm.terrain(np.linspace(0.25, 1, 256))
# all_colors = np.vstack((colors_undersea, colors_land))
# # cmap = mcolors.LinearSegmentedColormap.from_list(
# #     'terrain_map', all_colors)


# cmap = px.colors.sequential.agsunset
# cmap = plotly.colors.PLOTLY_SCALES["Viridis"]





# cmap = cm.brg(Z)
divnorm=mcolors.TwoSlopeNorm(vmin=Z.min(), vcenter=0., vmax=Z.max())
divnorm=mcolors.TwoSlopeNorm(vmin=Z.min(), vcenter=(Z.max()+Z.min())/2, vmax=Z.max())
# divnorm=mcolors.TwoSlopeNorm(vmin=-500, vcenter=0, vmax=Z.max())

# divnorm=mcolors.TwoSlopeNorm(vmin=-10, vcenter=0. ,vmax=10)
# divnorm=mcolors.TwoSlopeNorm(vmin=-10, vcenter=0., vmax=Z.max())

# divnorm=mcolors.LogNorm(vmin=Z.min(),  vmax=Z.max())  #Test LogNorm

# cmap = cm.brg(divnorm(Z))

# ax.imshow(Z, extent=[-20, 20, -20, 20],origin='lower', norm = norm,
#                   cmap='coolwarm', alpha=0.6)

# ax.imshow(Z, extent=[-20, 20, -20, 20],origin='lower',
#                   cmap='coolwarm', alpha=0.6)
# ax.imshow(Z, extent=[-20, 20, -20, 20],origin='lower',
#                   cmap=cmap, alpha=0.6)

# divnorm=mcolors.TwoSlopeNorm(vmin=Z.min(), vcenter=0., vmax=Z.max())
# plt.imshow(Z, extent=[x.min(), x.max(), y.min(), y.max()],origin='lower',
#                   cmap=cmap, alpha=0.6)

# I = plt.imshow(Z, extent=[x.min(), x.max(), y.min(), y.max()],origin='lower', norm = divnorm,
#                   cmap=cmap, alpha=0.6)
plt.imshow(Z, extent=[x.min(), x.max(), y.min(), y.max()],origin='lower', norm = divnorm,
                  cmap=cmap, alpha=0.9)
# plt.imshow(Z, extent=[x.min(), x.max(), y.min(), y.max()],origin='lower', norm = divnorm,
#                   cmap=cmap, alpha=0.6)


# I = plt.imshow(Z, extent=[x.min(), x.max(), y.min(), y.max()],origin='lower',
                  # cmap=cmap, alpha=0.6)


# I = plt.imshow(Z, extent=[x.min(), x.max(), y.min(), y.max()],origin='lower', norm = mcolors.CenteredNorm(),
#                   cmap=cmap, alpha=0.6)



# COLORBAR :
# cbar = plt.colorbar()
# cbar.ax.tick_params(labelsize=f_size)
# fig.colorbar(I)



##----- ADD RECTANGLE TO COVER QUADRANT :
epsilon = 0.4
epsilon = 0.2
epsilon = 0.005
epsilon = 0.01
# ax.axvspan(0, x.max(), y.min(), 0, alpha=1, color='yellow', zorder=5)#yellow
# ax.fill_between([0, x.max()], y.min(), 0, alpha=0.3, color='yellow', zorder=5)#yellow
# ax.fill_between([x.min(), 0], 0, y.max(), alpha=0.3, color='yellow', zorder=5)#yellow

fillcolor = 'white'
# ax.fill_between([0+epsilon, x.max()], y.min(), 0-epsilon, alpha=0.7, color=fillcolor, zorder=4)#yellow
# ax.fill_between([x.min(), 0-epsilon], 0+epsilon, y.max(), alpha=0.7, color=fillcolor, zorder=4)#yellow
ax.fill_between([0+epsilon, x.max()-epsilon], y.min()+epsilon, 0-epsilon, alpha=1.0, color=fillcolor, zorder=4)#yellow
ax.fill_between([x.min()+epsilon, 0-epsilon], 0+epsilon, y.max()-epsilon, alpha=1.0, color=fillcolor, zorder=4)#yellow

# ax.plot_surface(a1,a2, Z, cmap=cm.coolwarm,
#                        linewidth=0, antialiased=False)



# ax.plot(theta_rho, energy_axial1, 'royalblue', zorder=3, label=r"axialMin1")
# ax.plot(theta_rho, energy_axial2, 'forestgreen', zorder=3, label=r"axialMin2")
# ax.plot(-1.0*alphas, kappas, 'red', zorder=3, )

if case == 'a':
    # PLOT MINIMIZER g*:
    # ax.scatter(g_star[0],g_star[1], s=20, zorder=5)
    # ax.scatter(g_star[0],g_star[1], marker = 'x', s=20, zorder=5)
    # ax.text(g_star[0]+0.1,g_star[1]-0.1, r"$g_*$", color='royalblue', size=15, zorder = 5)
    ax.scatter(g_star[0],g_star[1], marker = 'x', s=40, color='darkslategray', zorder=4)
    ax.text(g_star[0]+0.15,g_star[1]-0.15, r"$g_*$", color='darkslategray', size=13, alpha=1.0, zorder = 5)
if case == 'b':
    # PLOT MINIMIZER g*:
    ax.scatter(g_star[0],g_star[1], s=20, color='royalblue', zorder=5)
    ax.scatter(g_star[0],g_star[1], marker = 'x', s=40, color='darkslategray', zorder=4)
    ax.text(g_star[0]+0.15,g_star[1]-0.15, r"$g_*$", color='darkslategray', size=13, alpha=1.0, zorder = 5)
if case == 'c':
    print('draw Line')
    line = draw_1ParameterLine(q1,q2,q12,q3)


# lg = ax.legend(bbox_to_anchor=(0.0, 0.75), loc='upper left')

if case == 'c':
    # draw MINIMAL VALUES
    print('np.where(Z == Z.min())',np.where(np.round(Z,6) == np.round(Z.min(),6)))
    tmp_x = a1[np.where(np.round(Z,6) == np.round(Z.min(),6))]
    tmp_y = a2[np.where(np.round(Z,6) == np.round(Z.min(),6))]
    # ax.scatter(tmp_x,tmp_y,color='red', zorder=5, s=1)
    # ax.plot(tmp_x,tmp_y,color='forestgreen', zorder=5)

# print('np.where(Z == Z.min())',np.where(np.round(Z,8) == np.round(Z.min(),8)))


if case == 'a':
    ax.scatter(minvec[0],minvec[1],color='royalblue',zorder=5, s=20)


### PLot x and y- Axes
ax.plot(ax.get_xlim(),[0,0],'k--', linewidth=0.5, zorder=4)
ax.plot([0,0],ax.get_ylim(), 'k--', linewidth=0.5, zorder=4)






# ax.set_xlabel(r"$a_1$", fontsize=f_size ,labelpad=0)
# ax.set_ylabel(r"$a_2$", fontsize=f_size ,labelpad=0)
ax.set_xlabel(r"$a_1$", fontsize=f_size )
ax.set_ylabel(r"$a_2$", fontsize=f_size, rotation=0 )

# ax.set_ylabel(r"energy")
ax.tick_params(axis='both', which='major', labelsize=f_size)
ax.tick_params(axis='both', which='minor', labelsize=f_size)



# ax.set_xticks([-np.pi/2, -np.pi/4  ,0,  np.pi/4,  np.pi/2 ])
# labels = ['$0$',r'$\pi/8$', r'$\pi/4$' ,r'$3\pi/8$' , r'$\pi/2$']
# ax.set_yticklabels(labels)


# plt.subplots_adjust(left = 0.25 )


# ax.legend(loc='upper right')

fig.subplots_adjust(left=.18, bottom=.12, right=.90, top=.97)

fig.set_size_inches(width, height)
# fig.set_size_inches(set_size(width, fraction=0.33))
Outputname = 'Energy_ContourG+_case'+ str(case) + '.pdf'
fig.savefig(Outputname)

plt.show()


#
#
#
# # Curve parametrised by \theta_rho = alpha in parameter space
# N=100;
# theta_rho = np.linspace(1, 3, num=N)
# print('theta_rho:', theta_rho)
#
#
# theta_values = []
#
#
# for t in theta_rho:
#
#         s = (1.0/10.0)*t+0.1
#         theta_values.append(s)
#
#
#
#
#
# theta_rho = np.array(theta_rho)
# theta_values = np.array(theta_values)
#
# betas_ = 2.0
#

# alphas, betas, thetas = np.meshgrid(theta_rho, betas_, theta_values, indexing='ij')
#
#
# harmonicMeanVec = np.vectorize(harmonicMean)
# arithmeticMeanVec = np.vectorize(arithmeticMean)
# prestrain_b1Vec = np.vectorize(prestrain_b1)
# prestrain_b2Vec = np.vectorize(prestrain_b2)
#
# GetMuGammaVec = np.vectorize(GetMuGamma)
# muGammas = GetMuGammaVec(betas,thetas,gamma,mu1,rho1,InputFilePath ,OutputFilePath )
#
# q1_vec = harmonicMeanVec(mu1, betas, thetas)
# q2_vec = arithmeticMeanVec(mu1, betas, thetas)
#
# b1_vec = prestrain_b1Vec(rho1, betas, alphas, thetas)
# b2_vec = prestrain_b2Vec(rho1, betas, alphas, thetas)

# special case: q12 == 0!! .. braucht eigentlich nur b1 & b2 ...

# print('type b1_values:', type(b1_values))



# print('size(q1)',q1.shape)
#
#
# energy_axial1 = []
# energy_axial2 = []
#
# # for b1 in b1_values:
# for i in range(len(theta_rho)):
#     print('index i:', i)
#
#     print('theta_rho[i]',theta_rho[i])
#     print('theta_values[i]',theta_values[i])
#
#     q1 = (1.0/6.0)*harmonicMean(mu1, beta, theta_values[i])
#     q2 = (1.0/6.0)*arithmeticMean(mu1, beta, theta_values[i])
#     q12 = 0.0
#     q3 = GetMuGamma(beta, theta_values[i],gamma,mu1,rho1,InputFilePath ,OutputFilePath )
#     b1 = prestrain_b1(rho1,beta, theta_rho[i],theta_values[i] )
#     b2 = prestrain_b2(rho1,beta, theta_rho[i],theta_values[i] )
#
#
#     # q2_vec = arithmeticMean(mu1, betas, thetas)
#     #
#     # b1_vec = prestrain_b1Vec(rho1, betas, alphas, thetas)
#     # b2_vec = prestrain_b2Vec(rho1, betas, alphas, thetas)
#     print('q1[i]',q1)
#     print('q2[i]',q2)
#     print('q3[i]',q3)
#     print('b1[i]',b1)
#     print('b2[i]',b2)
#     # print('q1[i]',q1[0][i])
#     # print('q2[i]',q2[i])
#     # print('b1[i]',b1[i])
#     # print('b2[i]',b2[i])
#     #compute axial energy #1 ...
#
#     a_axial1 = np.array([b1,0])
#     a_axial2 = np.array([0,b2])
#     b = np.array([b1,b2])
#
#     H = np.array([[2*q1, q12+2*q3], [q12+2*q3,2*q2] ])
#     A = np.array([[q1,1/2*q12], [1/2*q12,q2] ])
#
#
#     tmp = H.dot(a_axial1)
#
#     print('H',H)
#     print('A',A)
#     print('b',b)
#     print('a_axial1',a_axial1)
#     print('tmp',tmp)
#
#     tmp = (1/2)*a_axial1.dot(tmp)
#     print('tmp',tmp)
#
#     tmp2 = A.dot(b)
#     print('tmp2',tmp2)
#     tmp2 = 2*a_axial1.dot(tmp2)
#
#     print('tmp2',tmp2)
#     energy_1 = tmp - tmp2
#     print('energy_1',energy_1)
#
#
#     energy_axial1.append(energy_1)
#
#
#     tmp = H.dot(a_axial2)
#
#     print('H',H)
#     print('A',A)
#     print('b',b)
#     print('a_axial2',a_axial2)
#     print('tmp',tmp)
#
#     tmp = (1/2)*a_axial2.dot(tmp)
#     print('tmp',tmp)
#
#     tmp2 = A.dot(b)
#     print('tmp2',tmp2)
#     tmp2 = 2*a_axial2.dot(tmp2)
#
#     print('tmp2',tmp2)
#     energy_2 = tmp - tmp2
#     print('energy_2',energy_2)
#
#
#     energy_axial2.append(energy_2)
#
#
#
#
#
# print('theta_values', theta_values)
#
#
#


#
#
#
#
# kappas = []
# alphas = []
# # G.append(float(s[0]))
#
#
#
#
# for t in T :
#
#     abar_current = sstar*abar+t*abarperp;
#     # print('abar_current', abar_current)
#     abar_current[abar_current < 1e-10] = 0
#     # print('abar_current', abar_current)
#
#     # G = np.array([[2*q1, q12+2*q3], [q12+2*q3,2*q2] ])
#     G = [abar_current[0], abar_current[1] , (2*abar_current[0]*abar_current[1])**0.5 ]
#
#     e = [(abar_current[0]/(abar_current[0]+abar_current[1]))**0.5, (abar_current[1]/(abar_current[0]+abar_current[1]))**0.5]
#     kappa = abar_current[0]+abar_current[1]
#     alpha = math.atan2(e[1], e[0])
#
#     print('angle current:', alpha)
#
#     kappas.append(kappa)
#     alphas.append(alpha)
#
#
#
# alphas = np.array(alphas)
# kappas = np.array(kappas)
#
#
# print('kappas:',kappas)
# print('alphas:',alphas)
# print('min alpha:', min(alphas))
# print('min kappa:', min(kappas))
#
# mpl.rcParams['text.usetex'] = True
# mpl.rcParams["font.family"] = "serif"
# mpl.rcParams["font.size"] = "9"
# width = 6.28 *0.5
# height = width / 1.618
# fig = plt.figure()
# # ax = plt.axes((0.15,0.21 ,0.75,0.75))
# ax = plt.axes((0.15,0.21 ,0.8,0.75))
# ax.tick_params(axis='x',which='major', direction='out',pad=5)
# ax.tick_params(axis='y',which='major', length=3, width=1, direction='out',pad=3)
# # ax.xaxis.set_major_locator(MultipleLocator(0.1))
# # ax.xaxis.set_minor_locator(MultipleLocator(0.05))
# # ax.xaxis.set_major_locator(plt.MultipleLocator(np.pi / 8))
# # ax.xaxis.set_minor_locator(plt.MultipleLocator(np.pi / 16))
# ax.xaxis.set_major_locator(plt.MultipleLocator(np.pi / 2))
# ax.xaxis.set_minor_locator(plt.MultipleLocator(np.pi / 4))
# ax.xaxis.set_major_formatter(plt.FuncFormatter(format_func))
# ax.grid(True,which='major',axis='both',alpha=0.3)
#
#
#
#
# ax.plot(alphas, kappas, 'royalblue', zorder=3, )
# ax.plot(-1.0*alphas, kappas, 'red', zorder=3, )