I am doing the following fiiting of parameters of some matrices and the code seems to have some problem

[shukla-fits]

import numpy as np
import math
from scipy.linalg import eigh
from iminuit import Minuit

Define the loop function h1 (same as f1)

def h1(x):
return -1 / (16 * np.pi2) * (np.log(x2) - 1)

def h2(x):
return -1 / (4 * 16 * np.pi2) * (2.0 * np.log(x2) - 1)

def f2(t):
return -1 / (4 * 16 * np.pi2) * (2.0 * np.log(t2) - 1)

def funcMN(mm, rr):
if mm == rr:
return h1(mm)

numerator = mm**2 * np.log(mm**2) - rr**2 * np.log(rr**2)
denominator = mm**2 - rr**2
return -1 / (16 * np.pi**2) * ((numerator / denominator) - 1)

def funcloop(c, t):
if c == 1:
raise ValueError("Input c must not be 1 to avoid division by zero.")

term1 = 2 * c**2 * np.log(c)
term2 = (c - 1) * (1 - 3 * c + 2 * (c - 1) * np.log(t**2))
numerator = term1 + term2
denominator = 4 * (c - 1)**2
return -1 / (16 * np.pi**2) * (numerator / denominator)

Generate real symmetric matrix from flattened parameters

def generate_real_symmetric_matrix_from_params(params, shape):
assert shape[0] == shape[1], "Matrix must be square."
matrix = np.zeros(shape)
indices = np.triu_indices(shape[0])
matrix[indices] = params
matrix[(indices[1], indices[0])] = params # Make symmetric
return matrix

Generate complex matrix from flattened parameters

def generate_complex_matrix_from_params(params, shape):
half = len(params) // 2
real_part = np.reshape(params[:half], shape)
imaginary_part = np.reshape(params[half:], shape)
return real_part + 1j * imaginary_part

Define the chi-square function

def chi_square(params, Md_expected, Mu_expected, Me_expected):
# Split parameters into matrices
y1_params = params[:3]
y2_params = params[3:21]
y3_params = params[21:24]
MT_param = params[24]
MN_param = params[25]
MV_param = 2.0

# f1 is the same as h1
f1 = h1

# Assuming g is defined
g = 0.6  # Replace with the actual value of g

# Generate matrices from parameters
y1 = np.diag(y1_params)
y2 = generate_complex_matrix_from_params(y2_params, (3, 3))
y3 = y3_params[:, np.newaxis]  # Make y3 a column matrix

# Calculate delyd
delyd = 2 * g**2 * y2 * h1(MV_param) + (y1 @ np.conjugate(y1) @ y2) * f1(MT_param) + y2 @ np.conjugate(y3) @ y3.T * funcMN(MT_param, MN_param)

# Calculate delye
delye = 6 * g**2 * y2.T * h1(MV_param) + 3 * y2.T @ np.conjugate(y1) @ y1 * f1(MT_param)

# Calculate delyu
delyu = 4 * g**2 * y1 * h1(MV_param) + (y1 @ np.conjugate(y2) @ y2.T + y2 @ np.conjugate(y2.T) @ y1.T) * f1(MT_param)

identity_matrix = np.eye(3)

# Calculate kql
kql = -3 * g**2 * h2(MV_param) * identity_matrix + (0.5 * np.conjugate(y1) @ y1 + np.conjugate(y2) @ y2.T) * f2(MT_param)

# Calculate kuc
kuc = -4 * g**2 * h2(MV_param) * identity_matrix + (np.conjugate(y1) @ y1.T + 2 * np.conjugate(y2) @ y2.T) * f2(MT_param)

# Calculate kdc_
kdc = -2 * g**2 * h2(MV_param) * identity_matrix + (2 * np.conjugate(y2).T @ y2 * f2(MT_param)) + np.conjugate(y3) @ y3.T * funcloop(((MN_param**2) / (MT_param**2)), MT_param)

# Calculate kll
kll = -3 * g**2 * h2(MV_param) * identity_matrix + 3 * np.conjugate(y2).T @ y2 * f2(MT_param)

# Calculate kec
kec = -6 * g**2 * h2(MV_param) * identity_matrix + 3 * np.conjugate(y1).T @ y1 * f2(MT_param)

yu = y1 + delyu - 0.5 * (kql.T @ y1 + y1 @ kuc)
yd = y2 + delyd - 0.5 * (kql.T @ y2 + y2 @ kdc)
ye = y2.T + delye - 0.5 * (kll.T @ y2.T + y2.T @ kec)

# Perform eigenvalue decomposition
evalyd, evecyd = eigh(yd @ np.conjugate(yd.T))
evalyu, evecyu = eigh(yu @ np.conjugate(yu.T))  # Eigenvalues of yu
evalye, evecye = eigh(ye @ np.conjugate(ye.T))

sqrtyd = np.sqrt(evalyd)
sqrtyu = np.sqrt(evalyu)
sqrtye = np.sqrt(evalye)

# Calculate chi-square for eigenvalues
eigenvalue_chisq = np.sum(((sqrtyd - Md_expected) / (0.1 * Md_expected))**2) + \
                   np.sum(((sqrtyu - Mu_expected) / (0.1 * Mu_expected))**2) + \
                   np.sum(((sqrtye - Me_expected) / (0.1 * Me_expected))**2)

# Calculate the product of unitary matrices
product = np.conjugate(evecyu.T) @ evecyd  # Since evecyu is identity matrix for diagonal yu
st12 = np.abs(product[0, 1]) / np.sqrt(np.abs(product[0, 0])**2 + np.abs(product[0, 1])**2)
pt12q = np.arcsin(st12)
pt23q = np.arctan(np.abs(product[1, 2]) / np.abs(product[2, 2]))
st23 = np.sin(pt23q)

pt13q = np.arccos(np.abs(product[2, 2]) / np.cos(pt23q))
st13 = np.sin(pt13q)
JCPq = np.imag(product[0, 0] * product[2, 2] *
               np.conj(product[0, 2]) * np.conj(product[2, 0]))
sndltq = JCPq / (np.cos(pt12q) * np.sin(pt12q) * np.cos(pt23q) *
                 np.sin(pt23q) * np.sin(pt13q) * (np.cos(pt13q))**2.0)

product_chisq = ((st12 - st12exp) / (0.1 * st12exp))**2 + \
                ((st23 - st23exp) / (0.1 * st23exp))**2 + \
                ((st13 - st13exp) / (0.1 * st13exp))**2 + \
                ((sndltq - 0.78) / (0.1 * 0.78))**2

return eigenvalue_chisq + product_chisq

Example usage and minimization

initial_params_y1_real = np.random.uniform(low=-math.sqrt(4math.pi), high=math.sqrt(4math.pi), size=(3,))
initial_params_y2_real = np.random.uniform(low=-math.sqrt(4math.pi), high=math.sqrt(4math.pi), size=(9,))
initial_params_y2_imag = np.random.uniform(low=-math.sqrt(4math.pi), high=math.sqrt(4math.pi), size=(9,))
initial_params_y3_real = np.random.uniform(low=-math.sqrt(4math.pi), high=math.sqrt(4math.pi), size=(3,))
initial_params_MT = np.random.uniform(low=1e-5, high=1e3, size=1)[0]
initial_params_MN = np.random.uniform(low=1e-5, high=1e3, size=1)[0]

initial_params = np.hstack((initial_params_y1_real, initial_params_y2_real, initial_params_y2_imag, initial_params_y3_real, initial_params_MT, initial_params_MN))

Mu_expected = np.array([2.81e-6, 0.00142, 0.427])
Md_expected = np.array([6.14e-6, 1.25e-4, 5.80e-3])
Me_expected = np.array([2.75e-6, 5.72e-4, 9.68e-3])
st12exp = 0.2286
st23exp = 0.04579
st13exp = 0.00424

User-defined errors for each element

y1_real_errors = [1e-7, 1e-4, 1e-2]
y2_real_errors = [1e-7, 1e-7, 1e-7, 1e-4, 1e-4, 1e-4, 1e-2, 1e-2, 1e-2]
y2_imag_errors = [1e-7, 1e-7, 1e-7, 1e-4, 1e-4, 1e-4, 1e-2, 1e-2, 1e-2]
y3_real_errors = [1e-4, 1e-1, 1e-1]
MT_error = [1e1]
MN_error = [1e-1]
error_array = np.array(y1_real_errors + y2_real_errors + y2_imag_errors + y3_real_errors + MT_error + MN_error)

def minimize_chi_square():
def chi_square_wrapper(*params):
return chi_square(np.array(params), Md_expected, Mu_expected, Me_expected)

m = Minuit(chi_square_wrapper, *initial_params)
print(minuit.params)
for _ in range(5):
    m.migrad()
    m.errors = [1e0] * len(initial_params)
    final_chisq = m.fval
    m.limits = [(-np.sqrt(4*np.pi), np.sqrt(4*np.pi))] * (len(initial_params) - 2) + [(1e-5, 1e3), (1e-5, 1e3)]
    print(final_chisq)

return m, final_chisq

minuit, final_chisq = minimize_chi_square()

Print the optimized parameters

print("Optimized Parameters:")
print(minuit.values)

Print final chi-square value

final_chisq = minuit.fval
print("Final Chi-square:", final_chisq)