Committing a working Chemistry Observable with integral generation delegated...

Committing a working Chemistry Observable with integral generation delegated to libint. added some examples showing its use
Signed-off-by: Alex McCaskey <mccaskeyaj@ornl.gov>

compiler/clang/LambdaVisitor.cpp

@@ -6,6 +6,7 @@
 
 #include "qcor.hpp"
 #include "clang/AST/Decl.h"
+#include "clang/AST/Expr.h"
 #include "clang/Basic/IdentifierTable.h"
 #include "clang/Basic/SourceLocation.h"
 #include "clang/Sema/Sema.h"
@@ -66,18 +67,34 @@ bool LambdaVisitor::CallExprToGateInstructionVisitor::VisitIntegerLiteral(
   return true;
 }
 
+bool LambdaVisitor::CallExprToGateInstructionVisitor::VisitUnaryOperator(UnaryOperator* op) {
+    if (op->getOpcode() == UnaryOperator::Opcode::UO_Minus) {
+        addMinus = true;
+    }
+    return true;
+}
+
 bool LambdaVisitor::CallExprToGateInstructionVisitor::VisitFloatingLiteral(
     FloatingLiteral *literal) {
-  InstructionParameter p(literal->getValue().convertToDouble());
+  double value = literal->getValue().convertToDouble();
+  InstructionParameter p(addMinus ? -1.0*value : value);
+  addMinus = false;
   parameters.push_back(p);
   return true;
 }
 
 bool LambdaVisitor::CallExprToGateInstructionVisitor::VisitDeclRefExpr(
     DeclRefExpr *decl) {
+  auto declName = decl->getNameInfo().getAsString();
+  if (addMinus) {
+     declName = "-"+declName;
+  }
   if (dyn_cast<ParmVarDecl>(decl->getDecl())) {
     parameters.push_back(
-        InstructionParameter(decl->getNameInfo().getAsString()));
+        InstructionParameter(declName));
+  } else if (dyn_cast<VarDecl>(decl->getDecl())) {
+      std::cout << "THIS IS A VARDECL: " << declName << "\n";
+      parameters.push_back(InstructionParameter(declName));
   }
   return true;
 }
@@ -305,7 +322,11 @@ bool LambdaVisitor::VisitLambdaExpr(LambdaExpr *LE) {
             {varName, (int)int_value->getValue().signedRoundToDouble()});
         continue;
       } else if (float_value) {
-
+         std::cout << varName << ", THIS DOUBLE VALUE IS KNOWN AT COMPILE TIME: "
+                  << float_value->getValue().convertToDouble()
+                  << "\n";
+        captures.insert(
+            {varName, float_value->getValue().convertToDouble()});
         continue;
       }
 
@@ -322,7 +343,28 @@ bool LambdaVisitor::VisitLambdaExpr(LambdaExpr *LE) {
     visitor.TraverseStmt(LE);
 
     auto function = visitor.getFunction();
-    // std::cout << "\n\nXACC IR:\n" << function->toString() << "\n";
+    for (auto &inst : function->getInstructions()) {
+        if (!inst->isComposite() && inst->nParameters() > 0) {
+            int counter = 0;
+            for (auto& p : inst->getParameters()) {
+                if (p.isVariable()) {
+                    // see if we have a runtime value in the captures map
+                    for (auto& kv : captures) {
+                        if (p.toString() == kv.first  && kv.second.isNumeric()) {
+                            inst->setParameter(counter, kv.second);
+                        } else if (p.toString() == "-"+kv.first && kv.second.which() == 1) {
+                            InstructionParameter pp(-1.0 * mpark::get<double>(kv.second));
+                            inst->setParameter(counter, pp);
+                        }
+                    }
+                }
+                counter++;
+            }
+        }
+    }
+
+
+    std::cout << "\n\nXACC IR:\n" << function->toString() << "\n";
 
     // Check if we have IRGenerators in the tree
     if (function->hasIRGenerators()) {

compiler/clang/LambdaVisitor.hpp

@@ -60,6 +60,7 @@ protected:
     std::vector<InstructionParameter> parameters;
     std::string name;
     std::shared_ptr<IRProvider> provider;
+    bool addMinus = false;
 
   public:
     CallExprToGateInstructionVisitor(const std::string n,
@@ -67,6 +68,7 @@ protected:
         : name(n), provider(p) {}
     std::shared_ptr<Instruction> getInstruction();
     bool VisitIntegerLiteral(IntegerLiteral *il);
+    bool VisitUnaryOperator(UnaryOperator* op);
     bool VisitFloatingLiteral(FloatingLiteral *fl);
     bool VisitDeclRefExpr(DeclRefExpr *expr);
   };

compiler/clang/tests/LambdaVisitorTester.in.cpp

@@ -85,6 +85,29 @@ int main() {
     };
     return 0;
 })param0";
+const std::string unary0 = R"unary0(
+int main() {
+    auto l = [&](double t) {
+        X(0);
+        Ry(t, 1);
+        CX(1,0);
+        Rx(-1.57,1);
+    };
+    return 0;
+})unary0";
+
+const std::string rtimeCapture = R"rtimeCapture(
+int main() {
+    double pi = 3.1415;
+    auto l = [&](double t) {
+        X(0);
+        Ry(t, 1);
+        CX(1,0);
+        Rx(pi,1);
+        Rx(-pi,1);
+    };
+    return 0;
+})rtimeCapture";
 
 const std::string hwe0 = R"hwe0(#include <vector>
 int main(int argc, char** argv){
@@ -316,6 +339,61 @@ return "lambda_visitor_tester";
 
     EXPECT_EQ(exp1,src2);
 }
+
+TEST(LambdaVisitorTester, checkUnary) {
+    Rewriter rewriter1, rewriter2;
+    auto action1 = new TestQCORFrontendAction(rewriter1);
+
+    xacc::setOption("qcor-compiled-filename", "lambda_visitor_tester");
+
+    std::vector<std::string> args{"-std=c++11"};
+
+    std::cout << "Source Code:\n" << unary0 << "\n";
+    // first case, I know compile time values, so ahead-of-time compilation
+    EXPECT_TRUE(tooling::runToolOnCodeWithArgs(action1, unary0, args));
+
+    std::ifstream t1(".output.cpp");
+    std::string src2((std::istreambuf_iterator<char>(t1)),
+                     std::istreambuf_iterator<char>());
+    std::remove(".output.cpp");
+
+    std::cout << "HELLO:\n" << src2 <<"\n";
+    const std::string expectedSrc = R"expectedSrc(
+int main() {
+    auto l = [&](){return "lambda_visitor_tester";};
+    return 0;
+})expectedSrc";
+
+  EXPECT_EQ(expectedSrc, src2);
+}
+
+TEST(LambdaVisitorTester, checkRuntimeCapture) {
+    Rewriter rewriter1, rewriter2;
+    auto action1 = new TestQCORFrontendAction(rewriter1);
+
+    xacc::setOption("qcor-compiled-filename", "lambda_visitor_tester");
+
+    std::vector<std::string> args{"-std=c++11"};
+
+    std::cout << "Source Code:\n" << rtimeCapture << "\n";
+    // first case, I know compile time values, so ahead-of-time compilation
+    EXPECT_TRUE(tooling::runToolOnCodeWithArgs(action1, rtimeCapture, args));
+
+    std::ifstream t1(".output.cpp");
+    std::string src2((std::istreambuf_iterator<char>(t1)),
+                     std::istreambuf_iterator<char>());
+    std::remove(".output.cpp");
+
+    std::cout << "HELLO:\n" << src2 <<"\n";
+    const std::string expectedSrc = R"expectedSrc(
+int main() {
+    double pi = 3.1415;
+    auto l = [&](){return "lambda_visitor_tester";};
+    return 0;
+})expectedSrc";
+
+  EXPECT_EQ(expectedSrc, src2);
+}
 int main(int argc, char **argv) {
   qcor::Initialize(argc, argv);
   ::testing::InitGoogleTest(&argc, argv);

examples/deuteron.cpp

@@ -21,6 +21,6 @@ int main(int argc, char **argv) {
   });
 
   auto results = future.get();
-  std::cout << "Results:\n";
-  results->print();
+//   std::cout << "Results:\n";
+//   results->print();
 }

examples/h2_vqe.cpp → examples/h2_hwe_vqe.cpp

@@ -12,20 +12,22 @@ int main(int argc, char **argv) {
 H          0.00000        0.00000        0.00000
 H          0.00000        0.00000        0.7474)geom";
 
-  auto op = qcor::getObservable("chemistry", {{"basis","sto-3g"}, {"geometry", geom}});
+  auto op = qcor::getObservable("chemistry",
+                                {{"basis", "sto-3g"}, {"geometry", geom}});
   int nq = op->nBits();
 
-  std::vector<std::pair<int,int>> coupling{{0,1},{1,2},{2,3}};
+  std::vector<std::pair<int, int>> coupling{{0, 1}, {1, 2}, {2, 3}};
 
   auto future = qcor::submit([&](qcor::qpu_handler &qh) {
     qh.vqe(
         [&](std::vector<double> x) {
-          hwe(x, {{"n-qubits",nq},{"layers",1}, {"coupling",coupling}});
+          hwe(x, {{"n-qubits", nq}, {"layers", 1}, {"coupling", coupling}});
         },
         op, optimizer);
   });
 
   auto results = future.get();
-  std::cout << "Results:\n";
-//   results->print();
+  auto energy = mpark::get<double>(results->getInformation("opt-val"));
+  std::cout << "Results: " << energy << "\n";
+  //   results->print();
 }

examples/h2_ucc1_vqe.cpp

+#include "qcor.hpp"
+
+int main(int argc, char **argv) {
+
+  qcor::Initialize(argc, argv);
+
+  auto optimizer = qcor::getOptimizer(
+      "nlopt", {{"nlopt-optimizer", "cobyla"}, {"nlopt-maxeval", 1000}});
+
+  auto geom = R"geom(2
+
+H          0.00000        0.00000        0.00000
+H          0.00000        0.00000        0.7474)geom";
+
+  auto op = qcor::getObservable("chemistry",
+                                {{"basis", "sto-3g"}, {"geometry", geom}});
+  auto future = qcor::submit([&](qcor::qpu_handler &qh) {
+    qh.vqe(
+        [&](double x) {
+          X(0);
+          X(2);
+          Rx(1.57,0);
+          H(1);
+          H(2);
+          H(3);
+          CX(0,1);
+          CX(1,2);
+          CX(2,3);
+          Rz(x, 3);
+          CX(2,3);
+          CX(1,2);
+          CX(0,1);
+          Rx(-1.57,0);
+          H(1);
+          H(2);
+          H(3);
+        },
+        op, optimizer);
+  });
+
+  auto results = future.get();
+  auto energy = mpark::get<double>(results->getInformation("opt-val"));
+  std::cout << "Results: " << energy << "\n";
+  //   results->print();
+}

examples/hwe_deuteron.cpp

@@ -11,6 +11,7 @@ int main(int argc, char **argv) {
       "pauli", "5.907 - 2.1433 X0X1 - 2.1433 Y0Y1 + .21829 Z0 - 6.125 Z1");
 
   int nq = op->nBits();
+  
   auto future = qcor::submit([&](qcor::qpu_handler &qh) {
     qh.vqe(
         [&](std::vector<double> x) {

ir/observable/chemistry/CMakeLists.txt

+set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_CURRENT_SOURCE_DIR}/libint-wrapper/cmake/Modules)
+
+find_package(Libint)
+
+if (Libint_FOUND)
+add_subdirectory(libint-wrapper)
 
 set(LIBRARY_NAME qcor-observable-chemistry)
 
@@ -9,9 +15,9 @@ usfunctiongeneratebundleinit(TARGET ${LIBRARY_NAME} OUT SRC)
 
 add_library(${LIBRARY_NAME} SHARED ${SRC})
 
-target_include_directories(${LIBRARY_NAME} PUBLIC . ../..)
+target_include_directories(${LIBRARY_NAME} PUBLIC . ../.. libint-wrapper/src)
 
-target_link_libraries(${LIBRARY_NAME} PUBLIC CppMicroServices PRIVATE xacc)
+target_link_libraries(${LIBRARY_NAME} PUBLIC CppMicroServices PRIVATE xacc libint-wrapper)
 
 set(_bundle_name qcor_chemistry_observable)
 set_target_properties(${LIBRARY_NAME}
@@ -34,3 +40,4 @@ if(QCOR_BUILD_TESTS)
 endif()
 
 install(TARGETS ${LIBRARY_NAME} DESTINATION ${CMAKE_INSTALL_PREFIX}/plugins)
+endif()
\ No newline at end of file

ir/observable/chemistry/ChemistryObservable.cpp

 #include "ChemistryObservable.hpp"
 #include "FermionOperator.hpp"
 
+#include "LibIntWrapper.hpp"
+
+#include "xacc_service.hpp"
+#include "ObservableTransform.hpp"
+
 #include <iostream>
 
 namespace qcor {
 namespace observable {
 
-namespace py = pybind11;
-using namespace py::literals;
-
 std::vector<std::shared_ptr<Function>>
-ChemistryObservable::observe(std::shared_ptr<Function> function) {}
+ChemistryObservable::observe(std::shared_ptr<Function> function) {
+   auto jw = xacc::getService<xacc::ObservableTransform>("jw");
+//    std::cout << "Fermion:\n" << fermionOp->toString() << "\n";
+   auto pauli = jw->transform(fermionOp);
+//    std::cout << "Pauli:\n" << pauli->toString() << "\n";
+   return pauli->observe(function);
+}
 
 const std::string ChemistryObservable::toString() { return ""; }
 
@@ -19,338 +27,55 @@ void ChemistryObservable::fromString(const std::string str) {
       "ChemisryObservable.fromString not implemented.");
 }
 
-const int ChemistryObservable::nBits() { return 0; }
-
+const int ChemistryObservable::nBits() { return fermionOp->nBits(); }
 
 void ChemistryObservable::fromOptions(
     std::map<std::string, InstructionParameter> &&options) {
   fromOptions(options);
 }
 
-
 void ChemistryObservable::fromOptions(
-    std::map<std::string, InstructionParameter> &options) {}
+    std::map<std::string, InstructionParameter> &options) {
+
+  std::map<std::string, std::string> opts;
+  for (auto &kv : options) {
+    opts.insert({kv.first, kv.second.toString()});
+  }
+  libintwrapper::LibIntWrapper libint;
+  libint.generate(opts);
+
+  auto hpq = libint.hpq();
+  auto hpqrs = libint.hpqrs();
+  auto enuc = libint.E_nuclear();
+
+  std::stringstream ss;
+  ss << enuc;
+  for (int i = 0; i < hpq.dimension(0); i++) {
+    for (int j = 0; j < hpq.dimension(1); j++) {
+      if (std::fabs(hpq(i, j)) > 1e-12) {
+        auto negOrPlus = hpq(i, j) < 0.0 ? " - " : " + ";
+        ss << negOrPlus << std::fabs(hpq(i, j)) << " " << i << "^ " << j;
+      }
+    }
+  }
+
+  for (int i = 0; i < hpqrs.dimension(0); i++) {
+    for (int j = 0; j < hpqrs.dimension(1); j++) {
+      for (int k = 0; k < hpqrs.dimension(2); k++) {
+        for (int l = 0; l < hpqrs.dimension(3); l++) {
+          if (std::fabs(hpqrs(i, j, k, l)) > 1e-12) {
+            auto negOrPlus = hpqrs(i, j, k, l) < 0.0 ? " - " : " + ";
+            ss << negOrPlus << std::fabs(hpqrs(i, j, k, l)) << " " << i << "^ "
+               << j << "^ " << k << " " << l;
+          }
+        }
+      }
+    }
+  }
+
+  fermionOp =
+      std::make_shared<xacc::quantum::FermionOperator>(ss.str());
+}
 
 } // namespace observable
 } // namespace qcor
-
-
-// void ChemistryObservable::fromOptions(
-//     std::map<std::string, InstructionParameter> &&options) {
-
-//   auto geom = options["geometry"].toString();
-//   std::ofstream tmp(".chem_obs_geom.xyz");
-//   tmp << geom;
-//   tmp.close();
-//   std::ifstream input_file(".chem_obs_geom.xyz");
-
-//   auto atoms = libint2::read_dotxyz(input_file);
-
-//   // count the number of electrons
-//   auto nelectron = 0;
-//   for (auto i = 0; i < atoms.size(); ++i)
-//     nelectron += atoms[i].atomic_number;
-//   const auto ndocc = nelectron / 2;
-//   std::cout << "# of electrons = " << nelectron << std::endl;
-
-//   auto enuc = 0.0;
-//   for (auto i = 0; i < atoms.size(); i++) {
-//     for (auto j = i + 1; j < atoms.size(); j++) {
-//       auto xij = atoms[i].x - atoms[j].x;
-//       auto yij = atoms[i].y - atoms[j].y;
-//       auto zij = atoms[i].z - atoms[j].z;
-//       auto r2 = xij * xij + yij * yij + zij * zij;
-//       auto r = std::sqrt(r2);
-//       enuc += atoms[i].atomic_number * atoms[j].atomic_number / r;
-//     }
-//   }
-//   std::cout << "Nuclear repulsion energy = " << std::setprecision(15) << enuc
-//             << std::endl;
-
-//   libint2::Shell::do_enforce_unit_normalization(false);
-
-//   //   auto molecule = std::make_shared<psi::Molecule>();
-
-//   std::cout << "Atomic Cartesian coordinates (a.u.):" << std::endl;
-//   for (const auto &a : atoms) {
-//     // molecule->add_atom(a.atomic_number, a.x, a.y, a.z);
-
-//     std::cout << a.atomic_number << " " << a.x << " " << a.y << " " << a.z
-//               << std::endl;
-//   }
-
-//   //   using namespace psi;
-
-//   libint2::initialize();
-//   //   std::shared_ptr<BasisSetParser> parser(new Gaussian94BasisSetParser());
-//   //  std::shared_ptr<BasisSet> aoBasis = BasisSet::construct(parser, molecule,
-//   //  options["basis"].toString());
-//   libint2::BasisSet obs(options["basis"].toString(), atoms, true);
-//   std::copy(begin(obs), end(obs),
-//             std::ostream_iterator<libint2::Shell>(std::cout, "\n"));
-
-//   libint2::Engine eri_engine(libint2::Operator::coulomb, obs.max_nprim(),
-//                              obs.max_l(), 0);
-//   //   eri_engine.set(libint2::BraKet::xs_xs);
-
-//   auto shell2bf = obs.shell2bf(); // maps shell index to basis function index
-//                                   // shell2bf[0] = index of the first basis
-//                                   // function in shell 0 shell2bf[1] = index of
-//                                   // the first basis function in shell 1
-//                                   // ...
-//   const auto &buf = eri_engine.results(); // will point to computed shell sets
-
-//   // const auto& is very important!
-//   auto n = obs.nbf();
-//   auto nshells = obs.size();
-//   Matrix result = Matrix::Zero(n, n);
-
-//   //  for (auto s1 = 0l, s12 = 0l; s1 != nshells; ++s1) {
-//   //       auto bf1 = shell2bf[s1];  // first basis function in this shell
-//   //       auto n1 = obs[s1].size();
-
-//   //       for (auto s2 = 0; s2 <= s1; ++s2, ++s12) {
-
-//   //         auto bf2 = shell2bf[s2];
-//   //         auto n2 = obs[s2].size();
-
-//   //         // compute shell pair; return is the pointer to the buffer
-//   //         eri_engine.compute(obs[s1], obs[s2], obs[s1], obs[s2]);
-//   //         if (buf[0] == nullptr)
-//   //           continue; // if all integrals screened out, skip to next shell
-//   //           set
-
-//   //         // "map" buffer to a const Eigen Matrix, and copy it to the
-//   //         // corresponding blocks of the result
-//   //         Eigen::Map<const Matrix> buf_mat(buf[0], n1, n2);
-//   //         result.block(bf1, bf2, n1, n2) = buf_mat;
-//   //         if (s1 != s2)  // if s1 >= s2, copy {s1,s2} to the corresponding
-//   //         {s2,s1}
-//   //                        // block, note the transpose!
-//   //           result.block(bf2, bf1, n2, n1) = buf_mat.transpose();
-//   //       }
-//   //     }
-
-//   Eigen::Tensor<double,1> data_(16);
-
-//   std::vector<double> data;
-// //   std::vector<double> data;
-//   int ii = 0;
-//   //   std::cout << "MATRIX:\n" << result << "\n";
-//   for (auto s1 = 0; s1 != obs.size(); ++s1) {
-//     for (auto s2 = 0; s2 != obs.size(); ++s2) {
-//       for (auto s3 = 0; s3 != obs.size(); ++s3) {
-//         for (auto s4 = 0; s4 != obs.size(); ++s4) {
-
-//           std::cout << "compute shell set {" << s1 << "," << s2 << "," << s3
-//                     << "," << s4 << "} ... ";
-//           eri_engine.compute(obs[s1], obs[s2], obs[s3], obs[s4]);
-//           std::cout << "done. ";
-//           auto ints_shellset = buf[0]; // location of the computed integrals
-//           if (ints_shellset == nullptr)
-//             continue; // nullptr returned if the entire shell-set was screened
-//                       // out
-
-//           auto bf1 = shell2bf[s1];  // first basis function in first shell
-//           auto n1 = obs[s1].size(); // number of basis functions in first shell
-//           auto bf2 = shell2bf[s2];  // first basis function in second shell
-//           auto n2 = obs[s2].size(); // number of basis functions in second shell
-//           auto bf3 = shell2bf[s3];  // first basis function in first shell
-//           auto n3 = obs[s3].size(); // number of basis functions in first shell
-//           auto bf4 = shell2bf[s4];  // first basis function in second shell
-//           auto n4 = obs[s4].size(); // number of basis functions in second shell
-
-//           // integrals are packed into ints_shellset in row-major (C) form
-//           // this iterates over integrals in this order
-//         //   for (int ii = 0;
-//             //    ii < sizeof(ints_shellset) / sizeof(ints_shellset[0]); ii++) {
-//             std::cout << ii << ", " << ints_shellset[0] << "\n";
-//             data.push_back(ints_shellset[0]);
-//             data_(ii) = ints_shellset[0];
-//         //   }
-//         ii++;
-//         }
-//       }
-//     }
-//   }
-
-// //   std::cout << "EIgen:\n" << data << "\n";
-//   std::vector<int> shape{2,2,2,2};
-//   std::vector<double> idvec {1.,0.,0.,1.};
-
-//  auto arr = xt::adapt(data, shape);
-//  xt::xarray<double> darr = arr;
-
-//  auto identity = xt::adapt(idvec, std::vector<int>{1,1,2,2});
-//  xt::xarray<double> did = identity;
-//  auto tmp1 = xt::linalg::kron(did, darr);
-
-//   auto && sh1 = tmp1.shape();
-//   for (auto& el : sh1) {std::cout << el << ", "; }
-
-//   std::cout << "\n";
-// //  auto kron2 = xt::linalg::kron(identity, xt::transpose(kron, {3,2,1,0}));
-// //  std::cout << arr << "\n" << identity << "\n";
-// //  auto && sh = kron2.shape();
-// // for (auto& el : sh) {std::cout << el << ", "; }
-// // std::cout << "\n";
-//   std::array<int, 4> dims{2,2,2,2};
-
-//   Eigen::Tensor<double, 4> data4 = data_.reshape(dims);
-
-
-//  Eigen::Tensor<double, 2> identityT(2,2);
-//  identityT.setZero();
-//  for (int i = 0; i < 2; i++) identityT(i,i) = 1.0;
-
-// Eigen::Tensor<double,4> idt = identityT.reshape(std::array<int,4>{1,1,2,2});
-
-// Eigen::Tensor<double,8> kroned = idt.contract(data4, std::array<Eigen::IndexPair<int>,0>{});
-
-//  for (int i = 0; i < 8; i++) std::cout << kroned.dimension(i) << " ";
-//  std::cout << "\n";
-
-//  std::vector<double> dd(kroned.data(), kroned.data() + kroned.size());
-
-//  auto adapted = xt::adapt(dd, std::vector<int>{1,1,2,2,2,2,2,2});
-
-//  std::vector<std::vector<int>> v{ {1,2,2,2,2,2,2}, {2,2,2,2,2,2}, {2,2,2,4,2}, {2,2,4,4} };
-//  for (int i = 0; i < 4; i++) {
-//      adapted = xt::concatenate(xt::xtuple(adapted, xt::ones<double>({1,1,2,2,2,2,2,2})), 3);
-//      adapted.reshape(v[i]);
-//     auto && sh11 = adapted.shape();
-//   for (auto& el : sh11) {std::cout << el << ", "; }
-
-//   std::cout << "\n";
-//  }
-
-//  auto && sh11 = adapted.shape();
-//   for (auto& el : sh11) {std::cout << el << ", "; }
-
-//   std::cout << "\n";
-
-// //  Eigen::Tensor<double,7> tmp1(1,2,2,2,2,2,2);
-// //  Eigen::Tensor<double,7> tmp2 = kroned.concatenate(tmp1,3);
-// //  Eigen::Tensor<double,6> tmp2 = tmp1.concatenate(tmp1,3);