close #37 + clean ups

functions/best_fit_synth_cl.py

@@ -48,7 +48,8 @@ def best_fit(err_lst, memb_prob_avrg_sort, completeness, ip_list, bf_params,
             isoch_fit_params = b_f(err_lst, obs_clust, completeness,
                                    ip_list, sc_params, ga_params, cmd_sel)
             # Assign errors as the steps in each parameter.
-            isoch_fit_errors = [ps_params[i + 3][2] for i in range(4)]
+            #param_rs = ip_list[2]
+            isoch_fit_errors = [p_rs[2] for p_rs in ip_list[2]]
 
         elif best_fit_algor == 'genet':
 
@@ -61,9 +62,8 @@ def best_fit(err_lst, memb_prob_avrg_sort, completeness, ip_list, bf_params,
                 obs_clust, completeness, ip_list, sc_params,
                 ga_params, cmd_sel)
 
-        print 'Best fit params obtained ({:g}, {:g}, {:g}, {:g}).'.format(
-            isoch_fit_params[0][0], isoch_fit_params[0][1],
-            isoch_fit_params[0][2], isoch_fit_params[0][3])
+        print ("Best fit params obtained ({:g}, {:g}, {:g}, {:g}, {:g},"
+            " {:g}).".format(*isoch_fit_params[0]))
 
         if best_fit_algor == 'genet':
             if N_b >= 2:
@@ -79,13 +79,13 @@ def best_fit(err_lst, memb_prob_avrg_sort, completeness, ip_list, bf_params,
                     isoch_fit_errors)
             else:
                 print 'Skipping bootstrap process.'
-                isoch_fit_errors = [-1.] * 4
+                isoch_fit_errors = [-1.] * len(isoch_fit_params[0])
 
         # For plotting purposes.
         # Get list of stored isochrones and their parameters.
         isoch_list, param_values = ip_list[0], ip_list[1]
         # Read best fit values for all parameters.
-        m, a, e, d = isoch_fit_params[0]
+        m, a, e, d, mass, binar_f = isoch_fit_params[0]
         # Find indexes for metallicity and age. If indexes are not found due
         # to some difference in the significant figures, use the indices
         # [0, 0] to prevent the code from halting.
@@ -97,13 +97,14 @@ def best_fit(err_lst, memb_prob_avrg_sort, completeness, ip_list, bf_params,
         shift_isoch = move_isoch(cmd_sel, isoch_list[m_i][a_i][:2], e, d)
         # Generate best fit synthetic cluster.
         synth_clst = s_c(err_lst, completeness, sc_params,
-                         isoch_list[m_i][a_i], [-1., -1., e, d], cmd_sel)
+                         isoch_list[m_i][a_i], [-1., -1., e, d, mass, binar_f],
+                         cmd_sel)
 
     else:
         # Pass empty lists to make_plots.
         print 'Skipping parameters fitting process.'
         isoch_fit_params, isoch_fit_errors, shift_isoch, synth_clst = \
-        [[-1., -1., -1., -1.]], [-1., -1., -1., -1.], [], []
+        [[-1., -1., -1., -1., -1.]], [-1., -1., -1., -1., -1.], [], []
 
     bf_return = [isoch_fit_params, isoch_fit_errors, shift_isoch, synth_clst]
 

functions/bootstrap_func.py

@@ -62,8 +62,8 @@ def bootstrap(err_lst, obs_clust0, completeness, ip_list, bf_params,
     # Calculate errors for each parameter.
     isoch_fit_errors = np.std(params_boot, 0)
     # Errors can not be smaller than the steps in each parameter.
-    for i in range(4):
-        isoch_fit_errors[i] = max(ps_params[i + 3][2],
-        isoch_fit_errors[i])
+    param_rs = ip_list[2]
+    for i, p_er in enumerate(isoch_fit_errors):
+        isoch_fit_errors[i] = max(param_rs[i][2], p_er)
 
     return isoch_fit_errors

functions/brute_force_algorithm.py

@@ -19,7 +19,7 @@ def brute_force(err_lst, obs_clust, completeness, ip_list, sc_params,
     isoch_list, param_values, param_rs = ip_list
 
     # Unpack parameters values.
-    m_lst, a_lst, e_lst, d_lst = param_values
+    m_lst, a_lst, e_lst, d_lst, mass_lst, bin_lst = param_values
 
     # Initiate list that will hold the likelihood values telling us how well
     # each isochrone (syhtnetic cluster) fits the observed data.
@@ -41,25 +41,31 @@ def brute_force(err_lst, obs_clust, completeness, ip_list, sc_params,
                 # Iterate through all distance modulus.
                 for d in d_lst:
 
-                    params = [m, a, e, d]
-
-                    # Call likelihood function with m,a,e,d values.
-                    isochrone = isoch_list[m_i][a_i]
-                    # Call likelihood function with m,a,e,d values.
-                    likel_val = i_l(err_lst, obs_clust, completeness, sc_params,
-                                    isochrone, params, cmd_sel)
-                    # Store the likelihood for each synthetic cluster.
-                    score[0].append(likel_val)
-                    score[1].append(params)
-
-                    # Print percentage done.
-                    i += 1
-                    percentage_complete = (100.0 * (i + 1) / tot_sols)
-                    while len(milestones) > 0 and percentage_complete >= \
-                    milestones[0]:
-                        print "  {}% done".format(milestones[0])
-                        # Remove that milestone from the list.
-                        milestones = milestones[1:]
+                    # Iterate through all masses.
+                    for mass in mass_lst:
+
+                        # Iterate through all binary fractions.
+                        for bin_frac in bin_lst:
+
+                            params = [m, a, e, d, mass, bin_frac]
+
+                            # Call likelihood function with m,a,e,d values.
+                            isochrone = isoch_list[m_i][a_i]
+                            # Call likelihood function with m,a,e,d values.
+                            likel_val = i_l(err_lst, obs_clust, completeness,
+                                sc_params, isochrone, params, cmd_sel)
+                            # Store the likelihood for each synthetic cluster.
+                            score[0].append(likel_val)
+                            score[1].append(params)
+
+                            # Print percentage done.
+                            i += 1
+                            percentage_complete = (100.0 * (i + 1) / tot_sols)
+                            while len(milestones) > 0 and \
+                            percentage_complete >= milestones[0]:
+                                print "  {}% done".format(milestones[0])
+                                # Remove that milestone from the list.
+                                milestones = milestones[1:]
 
     # Find index of function with smallest likelihood value.
     # This index thus points to the isochrone that best fits the observed

functions/checker.py

@@ -39,8 +39,8 @@ def check(mypath, cl_files):
         # Read input parameters from params_input.dat file.
         mode, done_dir, gd_params, gh_params, gc_params, cr_params, kp_flag,\
         im_flag, er_params, fr_number, pv_params, da_params, ps_params,\
-        bf_params, sc_params, ga_params, rm_params, pl_params, flag_move_file,\
-        axes_params = gip(mypath)
+        bf_params, sc_params, par_ranges, ga_params, rm_params, pl_params,\
+        flag_move_file, axes_params = gip(mypath)
     except Exception:
         # Halt code.
         print traceback.format_exc()
@@ -118,7 +118,8 @@ def check(mypath, cl_files):
     if bf_params[0]:
 
         # Unpack.
-        iso_path, cmd_select, iso_select, m_rs, a_rs, e_rs, d_rs = ps_params
+        iso_path, cmd_select, iso_select = ps_params
+        m_rs, a_rs, e_rs, d_rs, mass_rs, bin_rs = par_ranges
 
         # Check that CMD is correctly set.
         if cmd_select not in {1, 2, 3, 4, 5, 6, 7}:
@@ -144,10 +145,10 @@ def check(mypath, cl_files):
         # *WE ASUME ALL METALLICITY FILES HAVE THE SAME NUMBER OF AGE VALUES*
         age_vals_all = isochp.get_ages(met_files[0], isoch_format[1])
         # Get parameters ranges stored in params_input.dat file.
-        param_ranges, param_rs = isochp.get_ranges(m_rs, a_rs, e_rs, d_rs)
+        param_ranges, param_rs = isochp.get_ranges(par_ranges)
         # Check that ranges are properly defined.
         p_names = [['metallicity', m_rs], ['age', a_rs], ['extinction', e_rs],
-            ['distance', d_rs]]
+            ['distance', d_rs], ['mass', mass_rs], ['binary', bin_rs]]
         for i, p in enumerate(param_ranges):
             if not p.size:
                 sys.exit("ERROR: No values exist for {} range defined:\n\n"
@@ -179,7 +180,7 @@ def check(mypath, cl_files):
             # isoch_list, isoch_ma, isoch_ed = ip_list
             # Only read files if best fit process is set to run.
             # bf_flag = bf_params[0]
-            ip_list = isochp.ip(ps_params, bf_params[0])
+            ip_list = isochp.ip(ps_params, par_ranges, bf_params[0])
         except:
             print traceback.format_exc()
             sys.exit("ERROR: unknown error reading metallicity files.")

functions/func_caller.py

@@ -43,8 +43,8 @@ def asteca_funcs(mypath, cl_file, ip_list, R_in_place):
     # Read input parameters from params_input.dat file.
     mode, done_dir, gd_params, gh_params, gc_params, cr_params, kp_flag, \
     im_flag, er_params, fr_number, pv_params, da_params, ps_params, bf_params,\
-    sc_params, ga_params, rm_params, pl_params, flag_move_file, axes_params =\
-    gip(mypath)
+    sc_params, par_ranges, ga_params, rm_params, pl_params, flag_move_file, \
+    axes_params = gip(mypath)
 
     # Define system of coordinates used.
     px_deg = gd_params[-1]
@@ -182,7 +182,7 @@ def asteca_funcs(mypath, cl_file, ip_list, R_in_place):
             flag_area_stronger, field_region, flag_pval_test,
             pval_test_params, memb_prob_avrg_sort, lum_func, completeness,
             da_params, bf_params, red_return, err_lst, bf_return, ga_params,
-            er_params, axes_params, ps_params, pl_params)
+            er_params, axes_params, par_ranges, pl_params)
         print 'Plots created.'
 
     # Move file to 'done' dir if flag is set.

functions/genetic_algorithm.py

@@ -66,6 +66,7 @@ def crossover(chromosomes, p_cross, cr_sel):
             # Skip crossover operation.
             cross_chrom.append(chrom_pair[0])
             cross_chrom.append(chrom_pair[1])
+
     return cross_chrom
 
 
@@ -77,6 +78,7 @@ def mutation(cross_chrom, p_mut):
     for i, elem in enumerate(cross_chrom):
         cross_chrom[i] = ''.join(char if random.random() > p_mut else
         str(1 - int(char)) for char in elem)
+
     return cross_chrom
 
 
@@ -258,7 +260,7 @@ def gen_algor(flag_print_perc, err_lst, obs_clust, completeness, ip_list,
     # only depends on the total number of chromosomes n_pop and the fitness
     # differential fdif.
     fitness = [1. / n_pop + fdif * (n_pop + 1. - 2. * (i + 1.)) /
-    (n_pop * (n_pop + 1.)) for i in range(n_pop)]
+        (n_pop * (n_pop + 1.)) for i in range(n_pop)]
 
     ### Initial random population evaluation. ###
     p_lst_r = random_population(param_values, n_pop)
@@ -322,8 +324,8 @@ def gen_algor(flag_print_perc, err_lst, obs_clust, completeness, ip_list,
         # Evaluate each new solution in the objective function and sort
         # according to the best solutions found.
         generation, lkl, isoch_done = evaluation(err_lst, obs_clust,
-        completeness, isoch_list, param_values, p_lst_e, sc_params,
-        isoch_done, cmd_sel)
+            completeness, isoch_list, param_values, p_lst_e, sc_params,
+            isoch_done, cmd_sel)
 
         ### Extinction/Immigration ###
         # If the best solution has remained unchanged for n_ei

functions/get_in_params.py

@@ -98,11 +98,17 @@ def get_in_params(mypath):
                     e_rs = map(float, reader[1:4])
                 elif reader[0] == 'PS_d':
                     d_rs = map(float, reader[1:4])
+
+                # Synthetic cluster parameters.
                 elif reader[0] == 'SC':
                     IMF_name = str(reader[1])
-                    tot_mass = float(reader[2])
-                    f_bin = float(reader[3])
-                    q_bin = float(reader[4])
+                elif reader[0] == 'TM':
+                    mass_rs = map(float, reader[1:4])
+                elif reader[0] == 'FB':
+                    bin_mr = float(reader[1])
+                    bin_rs = map(float, reader[2:5])
+
+                # Genetic algorithm parameters.
                 elif reader[0] == 'GA':
                     n_pop = int(reader[1])
                     n_gen = int(reader[2])
@@ -157,14 +163,16 @@ def get_in_params(mypath):
     axes_params = [m_1, m_2, m_ord, xy_minmax]
 
     # Store photometric system params in lists.
-    ps_params = [iso_path, cmd_select, iso_select, m_rs, a_rs, e_rs, d_rs]
+    ps_params = [iso_path, cmd_select, iso_select]
 
     # Store GA params in lists.
     bf_params = [bf_flag, best_fit_algor, N_b]
-    sc_params = [IMF_name, tot_mass, f_bin, q_bin]
+    sc_params = [IMF_name, bin_mr]
+    par_ranges = [m_rs, a_rs, e_rs, d_rs, mass_rs, bin_rs]
     ga_params = [n_pop, n_gen, fdif, p_cross, cr_sel, p_mut, n_el, n_ei, n_es]
     rm_params = [flag_red_memb, min_prob]
 
     return mode, done_dir, gd_params, gh_params, gc_params, cr_params, kp_flag,\
     im_flag, er_params, fr_number, pv_params, da_params, ps_params, bf_params,\
-    sc_params, ga_params, rm_params, pl_params, flag_move_file, axes_params
+    sc_params, par_ranges, ga_params, rm_params, pl_params, flag_move_file, \
+    axes_params