Tests for evolven3fit_new and lhapdf .info files.

n3fit/src/evolven3fit_new/eko_utils.py

@@ -51,6 +51,7 @@ def construct_eko_cards(
         theory_card_dict = {}
     if op_card_dict is None:
         op_card_dict = {}
+
     # theory_card construction
     theory = Loader().check_theoryID(theoryID).get_description()
     theory.pop("FNS")
@@ -59,19 +60,29 @@ def construct_eko_cards(
         theory["nfref"] = NFREF_DEFAULT
     if "nf0" not in theory:
         theory["nf0"] = NF0_DEFAULT
+
     # The Legacy function is able to construct a theory card for eko starting from an NNPDF theory
     legacy_class = runcards.Legacy(theory, {})
     theory_card = legacy_class.new_theory
+
+    # Prepare the thresholds according to MaxNfPdf
+    thresholds = {"c": theory["kcThr"], "b": theory["kbThr"], "t": theory["kbThr"]}
+    if theory["MaxNfPdf"] < 5:
+        thresholds["b"] = np.inf
+    if theory["MaxNfPdf"] < 6:
+        thresholds["t"] = np.inf
+
     # construct operator card
     q2_grid = utils.generate_q2grid(
         theory["Q0"],
         q_fin,
         q_points,
-        {theory["mb"]: theory["kbThr"], theory["mt"]: theory["ktThr"]},
+        {theory["mb"]: thresholds["b"], theory["mt"]: thresholds["t"]},
     )
     op_card = default_op_card
-    masses = np.array([theory["mc"],theory["mb"],theory["mt"]]) ** 2
-    thresholds_ratios=np.array([theory["kcThr"],theory["kbThr"],theory["ktThr"]]) ** 2
+    masses = np.array([theory["mc"], theory["mb"], theory["mt"]]) ** 2
+    thresholds_ratios = np.array([thresholds["c"], thresholds["b"], thresholds["t"]]) ** 2
+
     atlas = Atlas(
         matching_scales=MatchingScales(masses * thresholds_ratios),
         origin=(theory["Q0"]**2, theory["nf0"])
@@ -100,3 +111,24 @@ def construct_eko_cards(
 
     op_card = runcards.OperatorCard.from_dict(op_card)
     return theory_card, op_card
+
+
+def split_evolgrid(evolgrid):
+    """Split the evolgrid in blocks according to the number of flavors and repeating the last entry of one block in the first entry of the next."""
+    evolgrid_index_list = []
+    evolgrid.sort()
+    starting_nf = evolgrid[0][1]
+    for evo_point in evolgrid:
+        current_nf = evo_point[1]
+        if current_nf != starting_nf:
+            evolgrid_index_list.append(evolgrid.index(evo_point))
+            starting_nf = current_nf
+    start_index = 0
+    evolgrid_list = []
+    for index in evolgrid_index_list:
+        evolgrid_list.append(evolgrid[start_index:index+1])
+        start_index = index
+    evolgrid_list.append(evolgrid[start_index:])
+    return evolgrid_list
+
+

n3fit/src/evolven3fit_new/evolve.py

@@ -77,14 +77,9 @@ def evolve_fit(
 
     usr_path = pathlib.Path(fit_folder)
     initial_PDFs_dict = load_fit(usr_path)
-    x_grid = np.array(
-        initial_PDFs_dict[list(initial_PDFs_dict.keys())[0]]["xgrid"]
-    ).astype(float)
+    x_grid = np.array(initial_PDFs_dict[list(initial_PDFs_dict.keys())[0]]["xgrid"]).astype(float)
     theoryID = utils.get_theoryID_from_runcard(usr_path)
-    theory, op = eko_utils.construct_eko_cards(
-        theoryID, q_fin, q_points, x_grid, op_card_dict, theory_card_dict
-    )
-    qed = theory.order[1] > 0
+
     if eko_path is not None:
         eko_path = pathlib.Path(eko_path)
         _logger.info(f"Loading eko from : {eko_path}")
@@ -94,23 +89,38 @@ def evolve_fit(
             eko_path = (Loader().check_theoryID(theoryID).path) / "eko.tar"
         except FileNotFoundError:
             _logger.info(f"eko not found in theory {theoryID}, we will construct it")
+            theory, op = eko_utils.construct_eko_cards(
+                theoryID, q_fin, q_points, x_grid, op_card_dict, theory_card_dict
+            )
             runner.solve(theory, op, dump_eko)
             eko_path = dump_eko
+
     with eko.EKO.edit(eko_path) as eko_op:
         x_grid_obj = eko.interpolation.XGrid(x_grid)
         eko.io.manipulate.xgrid_reshape(eko_op, targetgrid=x_grid_obj, inputgrid=x_grid_obj)
-    info = info_file.build(theory, op, 1, info_update={})
-    info["NumMembers"] = "REPLACE_NREP"
-    info["ErrorType"] = "replicas"
-    info["XMin"] = float(x_grid[0])
-    info["XMax"] = float(x_grid[-1])
+
     with eko.EKO.read(eko_path) as eko_op:
+        # Read the cards directly from the eko to make sure they are consistent
+        theory = eko_op.theory_card
+        op = eko_op.operator_card
+
+        qed = theory.order[1] > 0
+
+        # Modify the info file with the fit-specific info
+        info = info_file.build(theory, op, 1, info_update={})
+        info["NumMembers"] = "REPLACE_NREP"
+        info["ErrorType"] = "replicas"
+        info["XMin"] = float(x_grid[0])
+        info["XMax"] = float(x_grid[-1])
+        # Save the PIDs in the info file in the same order as in the evolution
+        info["Flavors"] = basis_rotation.flavor_basis_pids
+        info["NumFlavors"] = theory.couplings.max_num_flavs
         dump_info_file(usr_path, info)
-        for replica in initial_PDFs_dict.keys():
-            evolved_block = evolve_exportgrid(initial_PDFs_dict[replica], eko_op, x_grid, qed)
-            dump_evolved_replica(
-                evolved_block, usr_path, int(replica.removeprefix("replica_"))
-            )
+
+        for replica, pdf_data in initial_PDFs_dict.items():
+            evolved_blocks = evolve_exportgrid(pdf_data, eko_op, x_grid, qed)
+            dump_evolved_replica(evolved_blocks, usr_path, int(replica.removeprefix("replica_")))
+
     # remove folder:
     # The function dump_evolved_replica dumps the replica files in a temporary folder
     # We need then to remove it after fixing the position of those replica files
@@ -157,8 +167,8 @@ def evolve_exportgrid(exportgrid, eko, x_grid, qed):
             whether qed is activated or not
     Returns
     -------
-        : np.array
-        evolved block
+        : list(np.array)
+        list of evolved blocks
     """
     # construct LhapdfLike object
     pdf_grid = np.array(exportgrid["pdfgrid"]).transpose()
@@ -170,25 +180,28 @@ def evolve_exportgrid(exportgrid, eko, x_grid, qed):
 
     def ev_pdf(pid, x, Q2):
         return x * evolved_pdf[Q2]["pdfs"][pid][targetgrid.index(x)]
-
-    block = genpdf.generate_block(
-        ev_pdf,
-        xgrid=targetgrid,
-        evolgrid=eko.evolgrid,
-        pids=basis_rotation.flavor_basis_pids,
-    )
-    return block
-
-
-def dump_evolved_replica(evolved_block, usr_path, replica_num):
+    evolgrid_list = eko_utils.split_evolgrid(eko.evolgrid)
+    blocks = []
+    for evgrid in evolgrid_list:
+        block = genpdf.generate_block(
+            ev_pdf,
+            xgrid=targetgrid,
+            evolgrid=evgrid,
+            pids=basis_rotation.flavor_basis_pids,
+        )
+        blocks.append(block)
+    return blocks
+
+
+def dump_evolved_replica(evolved_blocks, usr_path, replica_num):
     """
     Dump the evolved replica given by evolved_block as the replica num "replica_num" in
     the folder usr_path/nnfit/replica_<replica_num>/usr_path.stem.dat
 
     Parameters
     ----------
-        evolved_block: numpy.array
-            block of an evolved PDF
+        evolved_block: list(numpy.array)
+            list of blocks of an evolved PDF
         usr_path: pathlib.Path
             path of the fit folder
         replica_num: int
@@ -199,7 +212,7 @@ def dump_evolved_replica(evolved_block, usr_path, replica_num):
     path_where_dump.mkdir(exist_ok=True)
     to_write_in_head = f"PdfType: replica\nFromMCReplica: {replica_num}\n"
     genpdf.export.dump_blocks(
-        path_where_dump, replica_num, [evolved_block], pdf_type=to_write_in_head
+        path_where_dump, replica_num, evolved_blocks, pdf_type=to_write_in_head
     )
     # fixing_replica_path
     utils.fix_replica_path(usr_path, replica_num)

n3fit/src/n3fit/tests/regressions/hyper-quickcard.yml

@@ -33,7 +33,7 @@ datacuts:
 
 ############################################################
 theory:
-  theoryid: 162        # database id
+  theoryid: 399        # database id
 
 hyperscan_config:
     stopping:
@@ -116,12 +116,11 @@ fitting:
 positivity:
   posdatasets:
     - { dataset: POSF2U,   maxlambda: 1e6 }  # Positivity Lagrange Multiplier
-    - { dataset: POSDYC    , maxlambda: 1e5 }
+    - { dataset: POSDYS    , maxlambda: 1e5 }
 
 integrability:
   integdatasets:
     - {dataset: INTEGXT8, maxlambda: 1e2}
-    - {dataset: INTEGXT3, maxlambda: 1e2}
 
 ############################################################
 closuretest:

n3fit/src/n3fit/tests/regressions/noval-quickcard.yml

@@ -31,7 +31,7 @@ datacuts:
 
 ############################################################
 theory:
-  theoryid: 162        # database id
+  theoryid: 399        # database id
 
 ############################################################
 genrep: True    # on = generate MC replicas, False = use real data
@@ -73,7 +73,6 @@ positivity:
 integrability:
   integdatasets:
     - {dataset: INTEGXT8, maxlambda: 1e2}
-    - {dataset: INTEGXT3, maxlambda: 1e2}
 
 ############################################################
 closuretest:

n3fit/src/n3fit/tests/regressions/pc-quickcard.yml

@@ -31,7 +31,7 @@ datacuts:
 
 ############################################################
 theory:
-  theoryid: 162        # database id
+  theoryid: 399        # database id
 
 ############################################################
 genrep: True    # on = generate MC replicas, False = use real data

n3fit/src/n3fit/tests/regressions/quickcard.yml

@@ -15,7 +15,7 @@ dataset_inputs:
 - { dataset: NMC, frac: 0.5 }
 - { dataset: SLACP, frac: 0.5}
 - { dataset: CMSZDIFF12, frac: 0.5, cfac: ['QCD'], sys: 10 }
-- { dataset: CMSTTBARTOT, frac: 1.0, cfac: ['QCD'] }
+- { dataset: ATLASTTBARTOT8TEV, frac: 1.0, cfac: ['QCD'] }
 
 ############################################################
 datacuts:
@@ -34,7 +34,7 @@ datacuts:
 
 ############################################################
 theory:
-  theoryid: 162        # database id
+  theoryid: 399        # database id
 
 ############################################################
 genrep: True    # on = generate MC replicas, False = use real data
@@ -52,11 +52,11 @@ parameters: # This defines the parameter dictionary that is passed to the Model
     optimizer_name: 'RMSprop'
     learning_rate: 0.00001
     clipnorm: 1.0
-  epochs: 1000
+  epochs: 1100
   positivity:
     multiplier: 1.05
     initial: 1.5
-  stopping_patience: 0.30 # percentage of the number of epochs
+  stopping_patience: 0.10 # percentage of the number of epochs
   layer_type: 'dense'
   dropout: 0.0
   threshold_chi2: 10.0
@@ -77,12 +77,11 @@ fitting:
 positivity:
   posdatasets:
     - { dataset: POSF2U,   maxlambda: 1e6 }  # Positivity Lagrange Multiplier
-    - { dataset: POSDYC, maxlambda: 1e5 }
+    - { dataset: POSDYS, maxlambda: 1e5 }
 
 integrability:
   integdatasets:
     - {dataset: INTEGXT8, maxlambda: 1e2}
-    - {dataset: INTEGXT3, maxlambda: 1e2}
 
 ############################################################
 debug: true