Skip to content

Commit

Permalink
Restore n3fit files from master
Browse files Browse the repository at this point in the history
  • Loading branch information
@achiefa
achiefa committed Oct 17, 2024
1 parent 333e597 commit d75f1b9
Show file tree
Hide file tree
Showing 4 changed files with 9 additions and 193 deletions.
132 changes: 3 additions & 129 deletions n3fit/src/n3fit/layers/DIS.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,6 @@
"""

import numpy as np
from scipy import interpolate as scint

from n3fit.backends import operations as op

Expand All @@ -40,122 +39,6 @@ class DIS(Observable):
while the input pdf is rank 4 of shape (batch_size, replicas, xgrid, flavours)
"""

def __init__(self, fktable_data,
fktable_arr,
dataset_name,
boundary_condition=None,
operation_name="NULL",
nfl=14,
n_replicas=1,
power_corrections=False,
ht_type=None,
exp_kinematics=None,
**kwargs):
super().__init__(fktable_data, fktable_arr, dataset_name, boundary_condition, operation_name, nfl, n_replicas, **kwargs)

self.compute_power_corrections = power_corrections
self.power_corrections = None

# NOTE
# Ratio of SFs are not implemented yet. Work in progress.
if self.compute_power_corrections and exp_kinematics is not None:
self.exp_kinematics = exp_kinematics
if ht_type is None:
self.ht_type = 'ABMP'
raise NotImplementedError("This part should be reimplemented.")
else:
self.ht_type = ht_type

if self.ht_type == 'ABMP':
self.power_corrections = self.compute_abmp_parametrisation()
elif self.ht_type == 'custom':
self.power_corrections = self.compute_custom_parametrisation()
else:
raise Exception(f"HT type {ht_type} is not implemented.")


def compute_abmp_parametrisation(self):
"""
This function is very similar to `compute_ht_parametrisation` in
validphys.theorycovariance.construction.py. However, the latter
accounts for shifts in the 5pt prescription. As of now, this function
is meant to work only for DIS NC data, using the ABMP16 result.
"""
x_knots = [0.0, 0.1, 0.3, 0.5, 0.7, 0.9, 1]
y_h2 = [0.023, -0.032, -0.005, 0.025, 0.051, 0.003, 0.0]
y_ht = [-0.319, -0.134, -0.052, 0.071, 0.030, 0.003, 0.0]
#h2_sigma = [0.019, 0.013, 0.009, 0.006, 0.005, 0.004]
#ht_sigma = [0.126, 0.040, 0.030, 0.025, 0.012, 0.007]
H_2 = scint.CubicSpline(x_knots, y_h2)
H_T = scint.CubicSpline(x_knots, y_ht)

# Reconstruct HL from HT and H2
def H_L(x):
return (H_2(x) - np.power(x, 0.05) * H_T(x))

H_2 = np.vectorize(H_2)
H_L = np.vectorize(H_L)

x = self.exp_kinematics['kin1']
y = self.exp_kinematics['kin3']
Q2 = self.exp_kinematics['kin2']
N2, NL = 1#compute_normalisation_by_experiment(self.dataname, x, y, Q2)

PC_2 = N2 * H_2(x) / Q2
PC_L = NL * H_L(x) / Q2
power_correction = PC_2 + PC_L
power_correction = power_correction.to_numpy()

return power_correction


def compute_custom_parametrisation(self):
"""
This function is very similar to `compute_ht_parametrisation` in
validphys.theorycovariance.construction.py. However, the latter
accounts for shifts in the 5pt prescription. As of now, this function
is meant to work only for DIS NC data, using the ABMP16 result.
"""
# Posteriors from 240812-01-ABMP-large-prior-7k
x_knots = [0.0, 0.1, 0.3, 0.5, 0.7, 0.9, 1]
y_h2_p = [-0.00441, 0.11169, -0.01632, 0.00000, -0.08742, -0.07279, 0.00000]
y_hl_p = [0.00000, -0.06241, -0.08655, -0.03306, 0.00000, -0.05987, 0.0000]
y_h2_d = [-0.04117, 0.00000, 0.03124, -0.01059, 0.04763, 0.00000, 0.00000]
y_hl_d = [0.00316, 0.00469, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000]

H_2p = scint.CubicSpline(x_knots, y_h2_p)
H_lp = scint.CubicSpline(x_knots, y_hl_p)
H_2d = scint.CubicSpline(x_knots, y_h2_d)
H_ld = scint.CubicSpline(x_knots, y_hl_d)

H_2p = np.vectorize(H_2p)
H_lp = np.vectorize(H_lp)
H_2d = np.vectorize(H_2d)
H_ld = np.vectorize(H_ld)

x = self.exp_kinematics['kin1']
y = self.exp_kinematics['kin3']
Q2 = self.exp_kinematics['kin2']
N2, NL = compute_normalisation_by_experiment(self.dataname, x, y, Q2)

if "_P_" in self.dataname or "HERA" in self.dataname:
PC_2 = N2 * H_2p(x) / Q2
PC_L = NL * H_lp(x) / Q2
elif "_D_" in self.dataname:
PC_2 = N2 * H_2d(x) / Q2
PC_L = NL * H_ld(x) / Q2
else:
# TODO
# Need to implement this
PC_2 = 0 / Q2 #N2 * H_2d(x) / Q2
PC_L = 0 / Q2 #NL * H_ld(x) / Q2

power_correction = PC_2 + PC_L
power_correction = power_correction.to_numpy()

return power_correction


def gen_mask(self, basis):
"""
Receives a list of active flavours and generates a boolean mask tensor
Expand Down Expand Up @@ -202,11 +85,7 @@ def build(self, input_shape):
if self.num_replicas > 1:
self.compute_observable = compute_dis_observable_many_replica
else:
# Currying the function so that the `Observable` does not need
# to get modified
def compute_dis_observable_one_replica_w_pc(pdf, padded_fk):
return compute_dis_observable_one_replica(pdf, padded_fk, power_corrections = self.power_corrections)
self.compute_observable = compute_dis_observable_one_replica_w_pc
self.compute_observable = compute_dis_observable_one_replica


def compute_dis_observable_many_replica(pdf, padded_fk):
Expand All @@ -228,14 +107,9 @@ def compute_dis_observable_many_replica(pdf, padded_fk):
return op.einsum('brxf, nxf -> brn', pdf[0], padded_fk)


def compute_dis_observable_one_replica(pdf, padded_fk, power_corrections = None):
def compute_dis_observable_one_replica(pdf, padded_fk):
"""
Same operations as above but a specialized implementation that is more efficient for 1 replica,
masking the PDF rather than the fk table.
"""
if power_corrections is None:

return op.tensor_product(pdf[0], padded_fk, axes=[(2, 3), (1, 2)])
else:

return op.tensor_product(pdf[0], padded_fk, axes=[(2, 3), (1, 2)]) + power_corrections
return op.tensor_product(pdf[0], padded_fk, axes=[(2, 3), (1, 2)])
46 changes: 2 additions & 44 deletions n3fit/src/n3fit/model_gen.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -138,9 +138,6 @@ def observable_generator(
positivity_initial=1.0,
integrability=False,
n_replicas=1,
exp_data=None,
power_corrections=None,
ht_type=None
): # pylint: disable=too-many-locals
"""
This function generates the observable models for each experiment.
Expand Down Expand Up @@ -182,10 +179,6 @@ def observable_generator(
set the positivity lagrange multiplier for epoch 1
integrability: bool
switch on/off the integrability constraints
power_corrections: bool
whether to include HT in theory predictions
ht_type: str
type of HT parametrisation
Returns
------
Expand All @@ -201,37 +194,16 @@ def observable_generator(
dataset_xsizes = []
model_inputs = []
model_observables = []
kin_by_dict = {}

if exp_data is not None:
included_processes = [
'DEUTERON',
'NMC',
'NUCLEAR',
#'HERACOMB',
]
for process in exp_data:
commondata = process.load_commondata()
for dataset in commondata:
if process.name in included_processes and "_NC_" in dataset.setname:
kin_by_dict[dataset.setname] = dataset.kinematics
else:
kin_by_dict[dataset.setname] = None

# The first step is to compute the observable for each of the datasets
for dataset in spec_dict["datasets"]:
# Get the generic information of the dataset
dataset_name = dataset.name
kinematics = None

# Look at what kind of layer do we need for this dataset
if dataset.hadronic:
Obs_Layer = DY
else:
Obs_Layer = DIS
if power_corrections:
if exp_data is not None and kin_by_dict[dataset_name] is not None:
kinematics = kin_by_dict[dataset_name]

# Set the operation (if any) to be applied to the fktables of this dataset
operation_name = dataset.operation
Expand All @@ -242,29 +214,15 @@ def observable_generator(
# list of validphys.coredata.FKTableData objects
# these will then be used to check how many different pdf inputs are needed
# (and convolutions if given the case)
if dataset.hadronic:
obs_layer = Obs_Layer(
obs_layer = Obs_Layer(
dataset.fktables_data,
dataset.fktables(),
dataset_name,
boundary_condition,
operation_name,
n_replicas=n_replicas,
name=f"dat_{dataset_name}",
)
else:
obs_layer = Obs_Layer(
dataset.fktables_data,
dataset.fktables(),
dataset_name,
boundary_condition,
operation_name,
n_replicas=n_replicas,
name=f"dat_{dataset_name}",
power_corrections=power_corrections,
exp_kinematics=kinematics if power_corrections else None,
ht_type=None if not power_corrections else ht_type
)
)

# If the observable layer found that all input grids are equal, the splitting will be None
# otherwise the different xgrids need to be stored separately
Expand Down
16 changes: 4 additions & 12 deletions n3fit/src/n3fit/model_trainer.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -109,8 +109,6 @@ def __init__(
theoryid=None,
lux_params=None,
replicas=None,
power_corrections=False,
ht_type=None
):
"""
Parameters
Expand Down Expand Up @@ -145,10 +143,6 @@ def __init__(
dictionary containing the params needed from LuxQED
replicas: list
list with the replicas ids to be fitted
power_corrections: bool
whether to include HT in theory predictions
ht_type: str
type of HT parametrisation
"""
# Save all input information
self.exp_info = list(exp_info)
Expand All @@ -166,8 +160,6 @@ def __init__(
self.lux_params = lux_params
self.replicas = replicas
self.experiments_data = experiments_data
self.power_corrections = power_corrections
self.ht_type = ht_type

# Initialise internal variables which define behaviour
if debug:
Expand Down Expand Up @@ -570,9 +562,6 @@ def _generate_observables(
invcovmat_tr=experiment_data["invcovmat"][i],
invcovmat_vl=experiment_data["invcovmat_vl"][i],
n_replicas=len(self.replicas),
exp_data=self.experiments_data,
power_corrections=self.power_corrections,
ht_type=self.ht_type
)

# Save the input(s) corresponding to this experiment
Expand Down Expand Up @@ -947,7 +936,10 @@ def hyperparametrizable(self, params):
)

if photons:
pdf_model.get_layer("add_photon").register_photon(xinput.input.tensor_content)
if self._scaler: # select only the non-scaled input
pdf_model.get_layer("add_photon").register_photon(xinput.input.tensor_content[:,:,1:])
else:
pdf_model.get_layer("add_photon").register_photon(xinput.input.tensor_content)

# Model generation joins all the different observable layers
# together with pdf model generated above
Expand Down
8 changes: 0 additions & 8 deletions n3fit/src/n3fit/performfit.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,8 +47,6 @@ def performfit(
maxcores=None,
double_precision=False,
parallel_models=False,
power_corrections=False,
ht_type=None
):
"""
This action will (upon having read a validcard) process a full PDF fit
Expand Down Expand Up @@ -135,10 +133,6 @@ def performfit(
whether to use double precision
parallel_models: bool
whether to run models in parallel
power_corrections: bool
whether to include HT in theory predictions
ht_type: str
Type of HT parametrisation
"""
from n3fit.backends import set_initial_state

Expand Down Expand Up @@ -208,8 +202,6 @@ def performfit(
theoryid=theoryid,
lux_params=fiatlux,
replicas=replica_idxs,
power_corrections=power_corrections,
ht_type=ht_type
)

# This is just to give a descriptive name to the fit function
Expand Down

0 comments on commit d75f1b9

Please sign in to comment.