Add first version of a fix for #361

Author: cschwan

pineappl/src/grid.rs

@@ -628,6 +628,128 @@ impl Grid {
             .for_each(|subgrid| subgrid.scale(factor));
     }
 
+    /// TODO
+    ///
+    /// # Errors
+    ///
+    /// TODO
+    pub fn reint_node_opt(&mut self) -> Result<bool> {
+        let mut modified = false;
+
+        let self_interpolations = self.interpolations().to_vec();
+        let grid_node_values: Vec<_> = self_interpolations
+            .iter()
+            .map(|interps| interps.node_values())
+            .collect();
+
+        for subgrid in &mut self.subgrids {
+            if let SubgridEnum::ImportSubgridV1(subgrid) = subgrid {
+                let node_values = subgrid.node_values();
+                let mut affected_indices = vec![None; node_values.len()];
+
+                for ((values, grid_values), affected) in node_values
+                    .iter()
+                    .zip(&grid_node_values)
+                    .zip(&mut affected_indices)
+                {
+                    if let Some(index) = values.iter().copied().position(|value| {
+                        !grid_values
+                            .iter()
+                            .find(|&&grid_value| subgrid::node_value_eq(value, grid_value))
+                            .is_some()
+                    }) {
+                        if affected.is_some() {
+                            // TODO: return an error
+                            todo!();
+                        } else {
+                            *affected = Some(index);
+                        }
+                    }
+                }
+
+                if affected_indices.iter().any(Option::is_some) {
+                    // filter out affected node values
+                    let new_node_values: Vec<_> = affected_indices
+                        .iter()
+                        .zip(&node_values)
+                        .map(|(affected, values)| {
+                            if let &Some(index) = affected {
+                                values[..index]
+                                    .iter()
+                                    .copied()
+                                    .chain(values[index + 1..].iter().copied())
+                                    .collect()
+                            } else {
+                                values.to_vec()
+                            }
+                        })
+                        .collect();
+                    let mut array = super::packed_array::PackedArray::new(
+                        new_node_values.iter().map(|values| values.len()).collect(),
+                    );
+
+                    let mut filtered = Vec::new();
+
+                    for (index, value) in subgrid.indexed_iter() {
+                        let Some(new_index): Option<Vec<_>> = index
+                            .iter()
+                            .copied()
+                            .zip(affected_indices.iter().copied())
+                            .map(|(idx, affected)| {
+                                if let Some(filter) = affected {
+                                    if idx < filter {
+                                        Some(idx)
+                                    } else if idx > filter {
+                                        Some(idx - 1)
+                                    } else {
+                                        None
+                                    }
+                                } else {
+                                    Some(idx)
+                                }
+                            })
+                            .collect()
+                        else {
+                            filtered.push(index);
+                            continue;
+                        };
+
+                        array[new_index.as_slice()] = value;
+                    }
+
+                    let mut subg = InterpSubgridV1::new(&self_interpolations);
+                    let weight = subgrid
+                        .indexed_iter()
+                        .filter_map(|(index, value)| {
+                            filtered
+                                .iter()
+                                .find(|&filter| *filter == index)
+                                .map(|_| value)
+                        })
+                        .sum();
+                    dbg!(weight);
+                    let ntuple: Vec<_> = affected_indices
+                        .iter()
+                        .zip(&node_values)
+                        .map(|(affected, node_value)| affected.map(|index| node_value[index]))
+                        .collect::<Option<_>>()
+                        // TODO: generalize this to include interpolated dimensions
+                        .unwrap();
+                    dbg!(&ntuple);
+                    subg.fill(&self_interpolations, &ntuple, weight);
+
+                    let mut new_subgrid = ImportSubgridV1::new(array, new_node_values);
+                    new_subgrid.merge_impl(&subg.into(), None);
+
+                    *subgrid = new_subgrid;
+                    modified = true;
+                }
+            }
+        }
+
+        Ok(modified)
+    }
+
     /// Repair the grid if it was written by bugged versions to disk.
     ///
     /// Returns `true` if this operations did anything. Currently, this scans for these problems:

pineappl_cli/src/write.rs

@@ -41,6 +41,7 @@ enum OpsArg {
     MulBinNorm(f64),
     Optimize(bool),
     OptimizeFkTable(FkAssumptions),
+    ReintNodeOpt(bool),
     Repair(bool),
     RewriteChannel((usize, Channel)),
     RewriteOrder((usize, Order)),
@@ -86,7 +87,12 @@ impl FromArgMatches for MoreArgs {
                         });
                     }
                 }
-                "merge_channel_factors" | "optimize" | "repair" | "split_channels" | "upgrade" => {
+                "merge_channel_factors"
+                | "optimize"
+                | "reint_node_opt"
+                | "repair"
+                | "split_channels"
+                | "upgrade" => {
                     let arguments: Vec<Vec<_>> = matches
                         .remove_occurrences(&id)
                         .unwrap()
@@ -99,6 +105,7 @@ impl FromArgMatches for MoreArgs {
                         args[index] = Some(match id.as_str() {
                             "merge_channel_factors" => OpsArg::MergeChannelFactors(arg[0]),
                             "optimize" => OpsArg::Optimize(arg[0]),
+                            "reint_node_opt" => OpsArg::ReintNodeOpt(arg[0]),
                             "repair" => OpsArg::Repair(arg[0]),
                             "split_channels" => OpsArg::SplitChannels(arg[0]),
                             "upgrade" => OpsArg::Upgrade(arg[0]),
@@ -396,6 +403,17 @@ impl Args for MoreArgs {
                     .try_map(|s| s.parse::<FkAssumptions>()),
                 ),
         )
+        .arg(
+            Arg::new("reint_node_opt")
+                .action(ArgAction::Append)
+                .default_missing_value("true")
+                .help("TODO")
+                .long("reint-node-opt")
+                .num_args(0..=1)
+                .require_equals(true)
+                .value_name("ENABLE")
+                .value_parser(clap::value_parser!(bool)),
+        )
         .arg(
             Arg::new("repair")
                 .action(ArgAction::Append)
@@ -596,6 +614,9 @@ impl Subcommand for Opts {
                     // UNWRAP: this cannot fail because we only modify the normalizations
                     .unwrap();
                 }
+                OpsArg::ReintNodeOpt(true) => {
+                    grid.reint_node_opt()?;
+                }
                 OpsArg::Repair(true) => {
                     grid.repair();
                 }
@@ -635,6 +656,7 @@ impl Subcommand for Opts {
                 OpsArg::Upgrade(true) => grid.upgrade(),
                 OpsArg::MergeChannelFactors(false)
                 | OpsArg::Repair(false)
+                | OpsArg::ReintNodeOpt(false)
                 | OpsArg::Optimize(false)
                 | OpsArg::SplitChannels(false)
                 | OpsArg::Upgrade(false) => {}

pineappl_cli/tests/write.rs

@@ -24,6 +24,7 @@ Options:
       --mul-bin-norm <NORM>             Multiply all bin normalizations with the given factor
       --optimize[=<ENABLE>]             Optimize internal data structure to minimize memory and disk usage [possible values: true, false]
       --optimize-fk-table <OPTIMI>      Optimize internal data structure of an FkTable to minimize memory and disk usage [possible values: Nf6Ind, Nf6Sym, Nf5Ind, Nf5Sym, Nf4Ind, Nf4Sym, Nf3Ind, Nf3Sym]
+      --reint-node-opt[=<ENABLE>]       TODO [possible values: true, false]
       --repair[=<ENABLE>]               Repair bugs saved in the grid [possible values: true, false]
       --rewrite-channel <IDX> <CHAN>    Rewrite the definition of the channel with index IDX
       --rewrite-order <IDX> <ORDER>     Rewrite the definition of the order with index IDX