This Python file performs a complex analysis of star formation histories by loading results from previously executed fits, defining conditions for valid combinations of parameters, and applying a series of functions to filter, optimize, and calculate statistics like mean chi-squared values. It includes importing external libraries, defining custom lambda functions for filtering based on stellar age, distance modulus, and color excess, and employs a bolometric correction to calculate apparent magnitudes. The script concludes by saving filtered and optimized results considering predefined conditions into a serialized (pickled) file, incorporating a datetime format to label the output file distinctly.

from itertools import combinations
import numpy as np
import pickle as pkl
from fidanka.isochrone.isochrone import interp_isochrone_age
from tqdm import tqdm
with open("./fitResults.denseAlpha.fixedLowmass.pkl", 'rb') as f:
    fitResults = pkl.load(f)
 
bfr = fitResults['r']
 
 
# Function to calculate mean chi2
def mean_chi2(optimization_results):
    return np.mean([opt.fun for opt in optimization_results])
 
 
# Function to check if A and E satisfy conditions
def satisfy_conditions(A, E, cond):
    if not cond['age'](A['opt'].x[0], E['opt'].x[0]):
        return False
    if not cond['mu'](A['opt'].x[1], E['opt'].x[1]):
        return False
    if not cond['Av'](A['opt'].x[2], E['opt'].x[2]):
        return False
conditions = {
    'age': lambda a, e: a < e,
    'mu': lambda a_mu, e_mu: 0.95 * e_mu <= a_mu <= 1.05 * e_mu,
    'Av': lambda a_av, e_av: 0.95 * e_av <= a_av <= 1.05 * e_av
}
# conditions = {
#     # Age condition: A is younger than E and within 1 percent of E's age
#     'age': lambda a_age, e_age: a_age < e_age and (e_age * 0.99) <= a_age,
 
#     # Distance modulus condition: Agree within 0.1 percent
#     'mu': lambda a_mu, e_mu: (e_mu * 0.999) <= a_mu <= (e_mu * 1.001),
 
#     # Color excess condition: Agree within 0.1 percent
#     'Av': lambda a_av, e_av: (e_av * 0.999) <= a_av <= (e_av * 1.001)
# }
 
 
# Function to apply bolometric corrector and get apparent mags
def apply_bc(iso, bc, dm, excess):
    return bc.apparent_mags(iso[:, 1], iso[:, 2], iso[:, 3], mu=dm, Av=excess)
 
bfr_list = [
    (population, Y, alphaML, fitResults)
    for population, popFitResults in bfr.items()
    for Y, YFitResults in popFitResults.items()
    for alphaML, fitResults in YFitResults.items()
]
 
 
# Generate combinations
combos = combinations(bfr_list, 2)
 
# Filter to keep only valid combinations and calculate mean chi2
valid_combos = [(A, E) for A, E in combos if satisfy_conditions(A[3], E[3], conditions)]
 
 
# Create the result list
result_list = []
for combo in tqdm(valid_combos):
    (popA, A_y, A_a, A), (popE, E_y, E_a, E) = combo
    A_age = A['opt'].x[0]
    E_age = E['opt'].x[0]
 
    A_dm = A['opt'].x[1]
    E_dm = E['opt'].x[1]
 
    A_excess = A['opt'].x[2]
    E_excess = E['opt'].x[2]
 
    A_iso = interp_isochrone_age(A['iso'], A_age)
    E_iso = interp_isochrone_age(E['iso'], E_age)
 
    A_apparent_mags = apply_bc(A_iso, A['bc'], A_dm, A_excess)
    E_apparent_mags = apply_bc(E_iso, E['bc'], E_dm, E_excess)
 
    mean_chi2_val = mean_chi2([A['opt'], E['opt']])
    result_list.append((mean_chi2_val, [(popA, A_y, A_a, A['opt']), (popE, E_y, E_a, E['opt'])], [A_apparent_mags, E_apparent_mags]))
 
 
# Filter out those with mean goodness of fit greater than 1
filtered_chi2_results = [result for result in result_list if result[0] <= 1]
 
# Sort the combos from best fit to worst
filtered_chi2_results.sort(key=lambda x: x[0])
print(filtered_chi2_results[0])
 
import datetime
currentDate = datetime.datetime.now().strftime("%m-%d-%Y")
with open(f"ReducedResults.denseAlpha.fixedLowmass_{currentDate}.pkl", 'wb') as f:
    pkl.dump(filtered_chi2_results, f)