Onboarding project: Refactoring momentum conservation

As discussed in the group meeting we wanted to onboard another student to the contraction code. @maowerner and I argued that we first need to clean up in order to make a meaningful onboarding process. Today I realized that we have at least one piece of code where one can make a local change without (hopefully) needing to understand much of the code. The following task is a good candidate to get the toolchain (C++, doxygen, git, CMake, tests, GitHub pull request, Travis CI) set up and do a little change without needing to dive into the whole code. At the same time we have integration tests covering this such that one does not need to fear to break anything.

@maowerner and I would provide help along the way as needed. I will formulate the problem with only some hints at first such that we can expand as needed. Please do not hesitate to ask us!

The goal is to have the following feature implemented. It must be send as a pull request with sensible commit message(s) to this repository. All unit tests must successfully complete on Travis CI. Then @maowerner or I will do a code review and then we merge this into master.

In the file src/init_lookup_tables.cpp there is a function build_quantum_numbers_from_correlator_list. This contains one else if-clause for each diagram that we have in the code. Such a (Feynman) diagram is just a thing with vertices (particles going in or out) and propagators between the vertices. The function builds all combination of momenta that the incoming and outgoing particles can have while enforcing momentum conservation. We have an operator list qn_op which contains the operators. These are a struct with a Eigen::Vector3i called momentum. This is a class that represents the three-momentum of the particle at the vertex.

Momentum conservation is easy: The sum of all momenta must be the zero vector. Right now this is implemented in a very redundant fashion. We want to generalize this such that there is general code instead of this long else if-cascade.

This is the code for a diagram that we call “C30” where two particles go in and one particle comes out:

  else if (correlator.type == "C3c" || correlator.type == "C30") {
    std::map<int, int> counter; /** initialized with zero */

    for (const auto &op0 : qn_op[0]) {
      for (const auto &op2 : qn_op[2]) {
        Vector p_so_1 = op0.momentum;
        Vector p_so_2 = op2.momentum;
        Vector p_so = p_so_1 + p_so_2;

        if (desired_total_momentum(p_so, correlator.tot_mom) &&
            momenta_below_cutoff(p_so_1, p_so_2, momentum_cutoff)) {
          for (const auto &op1 : qn_op[1]) {
            Vector p_si = op1.momentum;

            if (desired_total_momentum(p_si, correlator.tot_mom)) {
              if (p_so == -p_si) {
                const int p_tot = p_si.squaredNorm();
                counter[p_tot]++;

                quantum_numbers.emplace_back(std::vector<QuantumNumbers>{op0, op1, op2});
              }
            }
          }
        }
      }
    }
    int total_number_of_combinations = 0;
    for (const auto c : counter) {
      std::cout << "\tCombinations for P = " << c.first << ": " << c.second << std::endl;
      total_number_of_combinations += c.second;
    }
    std::cout << "\tTest finished - Combinations: " << total_number_of_combinations
              << std::endl;
  }

Let us go through this code in steps. We first loop over all operators (and therefore momenta) that are available at the two source vertices. Then we extract the momenta of the two source particle with the following (copied from above):

        Vector p_so_1 = op0.momentum;
        Vector p_so_2 = op2.momentum;
        Vector p_so = p_so_1 + p_so_2;

You see that the operator + is defined for these vectors and just does a componentwise sum of the elements. In p_so we have the total momentum P at the source. Then we only want P² below a certain threshold, so the function desired_total_momentum checks whether we want to keep it. momenta_below_cutoff checks whether the combination of these two momenta are below some other cutoff. This needs only to be done when we have two source vertices. With three vertices we will have to think about this, but for the scope of the issue it can be assumed that there is only one or two source vertices.

We then loop over the operators that are available on the sink. The momentum conservation is done with p_so == -p_si. We can rewrite this as p_si + p_so == Vector{0, 0, 0} and generalize this to the sum of arbitrary many source and sink vertices.

When you look at the function in the file you will notice that there are blocks that are very similar to each other. They differ in the operators that they use. In this C30 diagram you have seen that qn_op[0] and qn_op[2] are source operators and qn_op[1] is the sink operator. We can basically reduce this to the statement that we have vertices (0, 2) at the source and (1) at the sink. Other diagrams then just have different numbers there. The goal is to generalize this such that one can specify a diagram just like this:

std::pair<std::vector<int>, std::vector<int>> diagram_C30{{0, 2}, {1}};

This way we can eliminate some 35 lines of code per diagram and replace it one or two lines. In the future we want to add more diagrams (like three pions to three pions, eta to three pions); therefore adding more diagrams shall be simple.

For this task you need to create a data structure to hold the information about the diagram. We think that std::map<std::string, std::pair<std::vector<int>, std::vector<int>>> is what should do the job, but that it does not necessarily need to be this exactly. Then one has to generalize the code such that it does not depend on hard-coded numbers (the 0, 2, 1) but retrieves this from the data structure.

These are the tasks that should give some guidance. First setting everything up:

[ ] Create an account on GitHub.com, tell me your username such that I can add you to out HISKP-LQCD organization
[ ] Set up a keypair for SSH (ssh-keygen on the command line) if you do not already have one. Register the public key with GitHub
[ ] Set up git to know your name and email address.
[ ] If you use your personal computer, make sure that the dependencies are installed, see the README.
[ ] Clone this repository via SSH.
[ ] Create a build directory and build the code (in Debug mode).
[ ] Run all the tests with ctest. (You can use ctest --output-on-failure later on to reduce the output.)

Now we have the code working on your machine and we can start developing.

[ ] Create a new branch in git for this pull request named refactor-momentum-conservation with git checkout -b refactor-momentum-conservation.
[ ] Open the file src/init_lookup_tables.cpp and find the function build_quantum_numbers_from_correlator_list.
[ ] Add the data structure to hold the vertex numbers and fill in the numbers for one of the diagrams. Most if else constructs handle multiple related diagrams, you need to create an entry for both C3c and C30 with the same numbers. This is a bit redundant, but we can resolve that some time later.
[ ] Convert the else if block of the C3c and C30 diagrams to use the new data structure. All other diagrams shall still use the old structure such that the code still works. We do not want to break the code, we want to keep it working the whole time. This way we do not build up large amounts of errors to fix at the end.
[ ] Compile the code. If this fails, you need to fix the compilation errors.
[ ] Run all the tests. If the tests succeed it means that everything works as before. In case the tests fail you need to go back to the code and see what you need to change.
[ ] Commit your changes into the repository with a sensible commit message, e.g. git commit -am 'Refactor momentum conservation of C30 and C3c'. Commit messages must start with a capital letter and be in present tense.
[ ] Push your changes to GitHub such that the other members of the group can see your code. The command for this is git push origin refactor-momentum-conservation.
[ ] Convert the next else if block to use the new data structure. It is likely that you need to work a bit on the generalization to support the second case. Perhaps some other abstraction is needed to do this. This part likely is the hardest, and if you need help. please ask!
[ ] Compile, run the tests, commit the change to git.
[ ] Convert the remaining blocks to the new data structure and abstraction (change code, compile, run tests, commit).

Once we are done with all those items, we have a refactored code which has the exact same functionality as before but is cleaner and easier to extend. Now we need to get these changes approved and into the main branch (master) of the project.

[ ] Push your changes to GitHub (git push origin refactor-momentum-conservation).
[ ] Create a pull request. After you have pushed, this should be offered directly on https://github.com/HISKP-LQCD/sLapH-contractions. Otherwise go to https://github.com/HISKP-LQCD/sLapH-contractions/pulls and create it manually.

Now the pull request is online and @maowerner or I will do a review and then merge it.

HISKP-LQCD / sLapH-contractions

Onboarding project: Refactoring momentum conservation #76