[Template Fitting] Is it possible to apply this constraint?

[zjkjsd]

Aloha Dr. Dembinski,

Thank you for always being available and helpful! I wonder if you could give me some hints on implementing the following constraint in my template fitting. I have 6 components, say, A, B, C, D, E, F, and a constraint:
$$\frac{yield(A)}{yield(B)}=\frac{yield(C)}{yield(D)}*Constant_{known}$$

Mahalo,
Boyang

[HDembinski]

When you have an equality constraint like this, you would normally express one of the degrees of freedom in your cost function by the others.

I assume that yield(A), yield(B) etc are numbers, and free fit parameters (if we ignore the constraint).

def original_cost_function(yield_a, yield_b, yield_c, yield_d):
     ...

def modified_cost_function(yield_b, yield_c, yield_d):
    yield_a = <compute from yield_b, yield_c, yield_d using the equality constraint>
    return original_cost_function(yield_a, yield_b, yield_c, yield_d)

I assume you are doing a template fit, so in your case, it would work something like this.

template_cost = Template(...)  # insert your templates and data

def modified_cost_function(yield_b, yield_c, yield_d):
    yield_a = <compute from yield_b, yield_c, yield_d using the equality constraint>
    return template_cost(yield_a, yield_b, yield_c, yield_d)

m = Minuit(modified_cost_function, *initial_yields)
m.migrad()

An alternative is to use the Scipy minimizer and use a NonlinearConstraint. In your case that's

scipy.optimize.NonlinearConstraint(lambda yields: yields[0]/yields[1] - yields[2]/yields[3] * constant, 0, 0)

I don't have experience how well the fit works with equality constraints. So far, I always used the first solution to this problem.

[zjkjsd]

Thank you so much for your help!! Sorry for the late response, but I was sick for 2 weeks and didn't work on the fitting until now.

For solution 1, is it correct that the modified_cost_function needs to take as many arguments as the fitting parameters:

def modified_cost_function(yield_a, yield_b, yield_c, yield_d):
    same as above

Some errors are raised if yield_a is not there.

Both solutions gave some output, however, with INVALID errors too. (Scipy solution seems to be less problematic).
(there are 6 fitting parameters in the example below) Could you share some insights on how to fix the problem? Thank you very much!

# solution 1: migrad minimizer

plt.figure(figsize=(16,8))
plt.title('migrad minimizer')
c4 = Template(data, edges, temp, method="da") # 'jsc', 'asy'

def modified_c4(a,b,c,d,e,f):
    b = a*e/d/1.5
    return c4(a,b,c,d,e,f)

m4 = Minuit(modified_c4, *[3000]*6, name=samples_name)
m4.limits = (0, None)
m4.simplex().migrad()
m4.hesse()

# Solution 2: scipy minimizer

plt.figure(figsize=(16,8))
plt.title('scipy minimizer')
from scipy.optimize import NonlinearConstraint
c3 = Template(data, edges, temp, method="da") # 'jsc', 'asy'
m3 = Minuit(c3, *init_yields, name=samples_name)
m3.limits = (None, None)
m3.scipy(constraints=NonlinearConstraint(lambda *yields: yields[0]/yields[1] - yields[3]/yields[4] * 1.5, 0, 0))
m3.hesse()

[zjkjsd]

Here are my templates and data files. It would be greatly appreciated if you have time to take a look.
templates_data.zip

# load template
file_path = 'templates.json'
with open(file_path, encoding="utf-8") as file:
    ws = json.load(file)
electron_channel = ws['channels'][0]
sample_names = [s['name'] for s in electron_channel['samples']]
templates = np.array([s['data'] for s in electron_channel['samples']])

# load test_data
file_path2 = 'test.json'
with open(file_path2, encoding="utf-8") as file:
    ws2 = json.load(file)
electron_channel2 = ws2['channels'][0]
data = np.array([s['data'] for s in electron_channel2['samples']])
counts=np.sum(data,axis=0).reshape(41,47)

# remove empty bins in the template and test_data
indices_non_zero = np.where(np.sum(templates,axis=0).reshape(41,47) != 0) # mask for removing empty bins in the template
temp = [i.reshape(41,47)[indices_non_zero] for i in templates]
test_data = counts[indices_non_zero]
edges = np.linspace(0, len(test_data), len(test_data)+1)
init_yields = [3000]*6

[HDembinski]

def modified_c4(a,b,c,d,e,f):
    b = a*e/d/1.5
    return c4(a,b,c,d,e,f)

This cannot work, you need this:

def modified_c4(a,c,d,e,f):
    b = a*e/d/1.5
    return c4(a,b,c,d,e,f)

or you need to fix the parameter "b" from the point of view of Minuit. If you write a function which effectively does not use the value that you pass to the cost function, then derivatives with respect to that parameter which Minuit computes are always zero, including the second derivatives. That's why HESSE fails, you cannot invert a Hessian matrix which contains a row and column of zeros, it does not have full rank.

[zjkjsd]

Thank you very much! This makes sense and works well!