Remove single level implementatioN

pynbody · Dec 2, 2022 · d6c50a3 · d6c50a3
1 parent d5e3322
commit d6c50a3
Showing 1 changed file with 2 additions and 83 deletions.
diff --git a/genetIC/src/simulation/modifications/splice.hpp b/genetIC/src/simulation/modifications/splice.hpp
@@ -38,88 +38,6 @@ namespace modifications {
     return mask;
   }
 
-  //! Return the field f which satisfies f = a in flagged region while minimising (f-b).C^-1.(f-b) elsewhere
-  template<typename DataType, typename T=tools::datatypes::strip_complex<DataType>>
-  fields::Field<DataType,T> spliceOneLevel(fields::Field<DataType,T> & a,
-                                           fields::Field<DataType,T> & b,
-                                           const fields::Field<DataType,T> & cov) {
-
-      // To understand the implementation below, first read Appendix A of Cadiou et al (2021),
-      // and/or look at the 1D toy implementation (in tools/toy_implementation/gene_splicing.ipynb) which
-      // contains a similar derivation and near-identical implementation.
-
-      assert (&a.getGrid() == &b.getGrid());
-      assert (&a.getGrid() == &cov.getGrid());
-      auto mask = generateMaskFromFlags(a.getGrid());
-      auto maskCompl = generateMaskComplementFromFlags(a.getGrid());
-
-      a.toFourier();
-      b.toFourier();
-
-      // The preconditioner should be almost equal to the covariance.
-      // We however set the fundamental of the power spectrum to a non-null value,
-      // otherwise, the mean value in the spliced region is unconstrained.
-      fields::Field<DataType,T> preconditioner(cov);
-      preconditioner.setFourierCoefficient(0, 0, 0, 1);
-
-      fields::Field<DataType,T> delta_diff = b-a;
-      delta_diff.applyTransferFunction(preconditioner, 0.5);
-      delta_diff.toReal();
-
-      fields::Field<DataType,T> z(delta_diff);
-      z*=mask;
-      z.toFourier();
-      z.applyTransferFunction(preconditioner, -1.0);
-      z.toReal();
-      z*=maskCompl;
-      z.toFourier();
-      z.applyTransferFunction(preconditioner, 0.5);
-      z.toReal();
-
-      auto X = [&preconditioner, &maskCompl](const fields::Field<DataType,T> & input) -> fields::Field<DataType,T>
-      {
-        fields::Field<DataType,T> v(input);
-        assert (!input.isFourier());
-        assert (!v.isFourier());
-        v.toFourier();
-        v.applyTransferFunction(preconditioner, 0.5);
-        v.toReal();
-        v*=maskCompl;
-        v.toFourier();
-        v.applyTransferFunction(preconditioner, -1.0);
-        v.toReal();
-        v*=maskCompl;
-        v.toFourier();
-        v.applyTransferFunction(preconditioner, 0.5);
-        v.toReal();
-        return v;
-      };
-
-      fields::Field<DataType,T> alpha = tools::numerics::conjugateGradient<DataType>(X, z);
-
-      alpha.toFourier();
-      alpha.applyTransferFunction(preconditioner, 0.5);
-      alpha.toReal();
-
-      fields::Field<DataType,T> bInDeltaBasis(b);
-      bInDeltaBasis.toFourier();
-      bInDeltaBasis.applyTransferFunction(preconditioner, 0.5);
-      bInDeltaBasis.toReal();
-
-      alpha*=maskCompl;
-      alpha+=bInDeltaBasis;
-
-      delta_diff*=mask;
-      alpha-=delta_diff;
-
-      assert(!alpha.isFourier());
-      alpha.toFourier();
-      alpha.applyTransferFunction(preconditioner, -0.5);
-      alpha.toReal();
-
-      return alpha;
-  }
-
   template<typename DataType, typename T=tools::datatypes::strip_complex<DataType>>
   fields::OutputField<DataType> spliceMultiLevel(
     const fields::OutputField<DataType> &noiseB,
@@ -331,7 +249,8 @@ namespace modifications {
     rhs.getFieldForLevel(0).dumpGridData("rhs-0.npz");
     if (Nlevel > 1) rhs.getFieldForLevel(1).dumpGridData("rhs-1.npz");
     std::cout << "rhs[0]=" << rhs.getFieldForLevel(0).getDataVector()[0] << std::endl;
-    auto n_alpha = tools::numerics::minresYolo<DataType>(Q, rhs, 1e-5);
+    // auto n_alpha = tools::numerics::minresYolo<DataType>(Q, rhs, 1e-6);
+    auto n_alpha = tools::numerics::minres<DataType>(Q, rhs, 1e-6);
     // auto n_alpha = tools::numerics::bicgstab<DataType>(Q, rhs);
     preconditioner(n_alpha);