How do we implement custom global symmetries (generalized parity)?

Question

I'm trying to write a generalized parity for a double-wire model of Fermions that conserves global number and the parity of the number of fermions in a given wire. Namely the model proposed in https://arxiv.org/pdf/1705.01786.pdf (equations 1-3) with the g term set to 0. The model has the usual nearest-neighbor tunneling, a nearest-neighbor potential energy term, and a cooper pair exchance term between two coupled quanutm wires. To implement this, I've mapped the two 1-d chains into a single 1-d chain with odd numbered sites corresponding to chain a and even numbered sites corresponding to chain b. The parity that is conserved is the parity of the number of fermions on one of the chains, say the odd number of sites. So I need to keep track of that parity as well as enforce global number conservation. The number conservation works fine, but the global parity does not.

I've attempted to adapt the Fermion.h file to include this generalized parity. The adapting the fermion.h file is based on some of the comments here, http://www.itensor.org/support/2030/3-component-hubbard-model?show=2030#q2030, for a 3-component Fermi-hubbard model.

Here is my attempt to generalize the fermion.h file names have been changed to avoid an error that was being thrown about redefining things like Fermion and FermionSite.
#pragma once

 #include "itensor/mps/siteset.h"
 #include "itensor/util/str.h"

 namespace itensor {
 class FermionSiteAlt;
 using FermionOddSitePar = BasicSiteSet<FermionSiteAlt>;
class FermionSiteAlt
{
Index s;
public:

FermionSiteAlt(Index I) : s(I) { }
FermionSiteAlt(Args const& args = Args::global())
    {
    auto ts = TagSet("Site,Fermion");
    auto n = 1;
    if(args.defined("SiteNumber"))
      {
      n = args.getInt("SiteNumber");
      ts.addTags("n="+str(n));
      }
    auto conserveQNs = args.getBool("ConserveQNs",true);
    auto conserve_Nf = args.getBool("ConserveNf",conserveQNs);
    auto conserve_Pfodd = args.getBool("ConservePfodd",false);
    auto oddevenupdown = args.getBool("OddEvenUpDown",false);
    if(not conserveQNs)
        {
        s = Index(2,ts);
        }
    else if(not oddevenupdown) //usual case
        {
          if(conserve_Nf) //usual case
          {
            if(not conserve_Pfodd) //usual case 
                {
                s = Index(QN({"Nf",0,-1}),1,
                         QN({"Nf",1,-1}),1,Out,ts);
                } 
            else // conserve_Pfodd not usual case
                { 
                if( n%2==1 ) { 
                    s = Index(QN({"Nf",0,-1}),1,QN({"Pfodd",0,-2}),1,
                              QN({"Nf",1,-1}),1,QN({"Pfodd",1,-2}),1,
                                Out,ts);
                    }
                if( n%2==0 ) {  
                    s = Index(QN({"Nf",0,-1}),1,QN({"Pfodd",0,-2}),1,
                              QN({"Nf",1,-1}),1,QN({"Pfodd",0,-2}),1,
                                Out,ts);
                    }     
                }
          }
          else
            {
            s = Index(QN({"Pf",0,-2}),1,
                      QN({"Pf",1,-2}),1,Out,ts);
            }
        }
    else
        {
        auto q_emp = QN({"Sz",0},{"Nf",0,-1});
        QN q_occ;           
        if(n%2==1) q_occ = QN({"Sz",+1},{"Nf",1,-1});  
        else       q_occ = QN({"Sz",-1},{"Nf",1,-1});  

        s = Index(q_emp,1,
                  q_occ,1,Out,ts); 
        }
    }

Index
index() const { return s; }

IndexVal
state(std::string const& state)
    {
    if(state == "Emp" || state == "0") 
        {
        return s(1);
        }
    else 
    if(state == "Occ" || state == "1") 
        {
        return s(2);
        }
    else
        {
        throw ITError("State " + state + " not recognized");
        }
    return IndexVal{};
    }

ITensor
op(std::string const& opname,
   Args const& args) const
    {
    auto sP = prime(s);

    auto Emp  = s(1);
    auto EmpP = sP(1);
    auto Occ  = s(2);
    auto OccP = sP(2);

    auto Op = ITensor(dag(s),sP);

    if(opname == "N" || opname == "n")
        {
        Op.set(Occ,OccP,1);
        }
    else
    if(opname == "C")
        {
        Op.set(Occ,EmpP,1);
        }
    else
    if(opname == "Cdag")
        {
        Op.set(Emp,OccP,1);
        }
    else
    if(opname == "A")
        {
        Op.set(Occ,EmpP,1);
        }
    else
    if(opname == "Adag")
        {
        Op.set(Emp,OccP,1);
        }
    else
    if(opname == "F" || opname == "FermiPhase")
        {
        Op.set(Emp,EmpP,1);
        Op.set(Occ,OccP,-1);
        }
    else
    if(opname == "projEmp")
        {
        Op.set(Emp,EmpP,1);
        }
    else
    if(opname == "projOcc")
        {
        Op.set(Occ,OccP,1); 
        }
    else
        {
        throw ITError("Operator \"" + opname + "\" name not recognized");
        }

    return Op;
    }

//
// Deprecated, for backwards compatibility
//

FermionSiteAlt(int n, Args const& args = Args::global())
    {
    *this = FermionSiteAlt({args,"SiteNumber=",n});
    }

};

    //
    // Deprecated, for backwards compatability
    //


// Commented out
 //// Rebels don't care about backwards compatability. Burn it down! 
    // using SpinlessSite = FermionSiteAlt;
    // using Spinless = BasicSiteSet<SpinlessSite>;

    } //namespace itensor

Continued in first comment below - character limit prevented posting both sets of code.

Comments

Here is the code to actually compute the MPSs:

     #include "itensor/all.h"
     #include "itensor/FermionSiteOddSitesParity.h"
     #include <string>
     #include <iostream>
     #include <sstream>
     #include <iomanip>
     #include <map>
     #include <random>  
     #include <chrono>
     #include "itensor/util/print_macro.h"

     using namespace itensor;



     int main()
     {

    //
    // Set up model parameters
    //
     
     int PhysicalLatticeLength = 5 ;
     int numberOfParticles = 5;
   
   
     int NumberFermionicSpecies = 2;  
     int L = PhysicalLatticeLength*NumberFermionicSpecies;
   
     int Lp1by2m1 = (L+1)/2 - 1 ;


     println("L is \n");
     std::cout << L << " \n ";

     auto minW = 1.0 ;
     auto maxW = 1.0 ;
     auto deltaW = 0.5 ;    
     int     NumberWIndices = floor((maxW - minW)/ deltaW + 1.01);
 
     int NumberOfStates = 9; // 5
   
   
     double EnergySpectrum[NumberWIndices][NumberOfStates];
     double randphysicalonebodypot[PhysicalLatticeLength];
     double randonebodypotential[L];
     double overlapMatrix[NumberWIndices][NumberOfStates][NumberOfStates];
   
     double globalParityArray[NumberWIndices][NumberOfStates];   
   
     int   globalNumberTemp;
     std::vector<int> globalNumberVector;
   
     // Build the space, conserving Fermionic particle number     
     auto psites = FermionOddSitePar(L,{"ConserveNf",true, "ConservePfodd", true});
     // auto psites = Fermion(L,{"ConserveNf",true});
   
     // Parameters for Iemini et al. Wire Model   

     double t = 1.0 ; // Tunneling sets energy & time scale
     double U  = 0.5;
     double G = 0.0;
     // W will be swept over
 
 
   
    //
    // Build Hamiltonian and initial MPS state to stick into DMRG
    //
   
   
     for( auto WIndex: range(NumberWIndices) )   
     {    // start Delta Inex loop
   
    auto W = minW + deltaW*WIndex;
   
    println("About to begin the calculation for W =   \n" );
    std::cout << W << " \n ";
   
     // Create the operators for building the Hamiltonian
     auto ampo = AutoMPO(psites);
   
    // Nearest Neighbor Couplings
     for( auto ii:range1(PhysicalLatticeLength-1))
    {
    // Tunneling for all for both wires species  
    ampo += -t, "Cdag", NumberFermionicSpecies*ii-1, "C", NumberFermionicSpecies*ii+1; // Wire a
    ampo += -t, "Cdag", NumberFermionicSpecies*ii+1, "C", NumberFermionicSpecies*ii-1;      // Wire a
    ampo += -t, "Cdag", NumberFermionicSpecies*ii-0, "C", NumberFermionicSpecies*ii+2; // Wire b
    ampo += -t, "Cdag", NumberFermionicSpecies*ii+2, "C", NumberFermionicSpecies*ii-0; // Wire b
   
    // Potential Term of adjascent sites
    ampo += -2*U, "N", NumberFermionicSpecies*ii-1, "N", NumberFermionicSpecies*ii+1;
    ampo += -2*U, "N", NumberFermionicSpecies*ii-0, "N", NumberFermionicSpecies*ii+2;      
        
    // Nearest-Neighbor Cooper Pair exchange between wires
    ampo += W, "Cdag", NumberFermionicSpecies*ii-0, "Cdag", NumberFermionicSpecies*ii+2, "C", NumberFermionicSpecies*ii+1, "C", NumberFermionicSpecies*ii-1;
         ampo += W, "Cdag", NumberFermionicSpecies*ii-1, "Cdag", NumberFermionicSpecies*ii+1, "C", NumberFermionicSpecies*ii+2, "C", NumberFermionicSpecies*ii-0;
   
   
    }

    // Single-site Terms
     for( auto ii:range1(PhysicalLatticeLength))
    {
   
    // One-body term to set zero
    ampo += 2*U, "N", NumberFermionicSpecies*ii-1; // na(j)
    ampo += 2*U, "N", NumberFermionicSpecies*ii-0; // nb(j)  
   
   
    }

     // Generate the Hamiltonian
     auto H = toMPO(ampo);
   
     // Generating an initial state over the space specified above   
   
     auto state = InitState(psites);
     // Set the initial state to be an MPS over the set of sites call psites
     //// With defined QNs from FermionSiteOddSitesParity
   
     // Set QNs for total particle number and fixed-parity
     for( int ll : range1(L))
     {
       state.set(ll,"0");    
     }
     for( int ll : range1(numberOfParticles) )
     {
       state.set(2*ll-1,"1");    
     }

     // Generate the initial MPS to actually use
     auto psi0 = state;
     

    //
    // Set DMRG properties
    //
   
      for( auto ii : range1(L) ) if( ii%2==1 ) state.set(ii,"Occ");
    // Set the number of DMRG sweeps
     auto sweeps = Sweeps(8);
     // Set the maximum bond dimension per sweep
     // Setting it low initially allows the computation to be much faster
     sweeps.maxdim() = 80, 80, 100, 150, 150,  250, 350, 400;  
     // Set the cutoff for the bond dimension
     sweeps.cutoff() = 1E-6, 1E-6, 1E-8, 1E-8, 1E-10, 1E-10, 1E-10, 1E-11;  
     // Set the number of Davidson iterations per sweep
     sweeps.niter() =  4, 4, 4, 3, 3, 3, 2, 2;
     // Set the noise to add to the MPS itself each sweep   
     sweeps.noise() = 1E-5, 1E-5, 1E-6, 1E-7, 1E-8, 1E-10, 1E-10, 1E-11;
   
    //
    //  Run the DMRG
    //
    int nstates = NumberOfStates;
    std::vector<MPS> states;
    std::vector<double> energies;
   
        //     int nn = 0 ;
    auto [energy,psi] = dmrg(H,states,psi0,sweeps,"Silent");
        //// Calculate the energy and state
        //// Well, that was short
        states.push_back(psi);
        energies.push_back(energy);
   
    // Can save ITensor objects to file, but not generic C data arrays - not great?
       writeToFile("psig", psi);
        
    EnergySpectrum[WIndex][0] = energy;
               
    for(int nn = 1; nn < nstates; nn++){
        auto [energye,psie] = dmrg(H,states,psi0,sweeps,{"Silent=", false, "Weight=", 5.0});

        states.push_back(psie);
        energies.push_back(energye);
        EnergySpectrum[WIndex][nn] = energye;
        }
   
        //
        //  Post-process the MPS for parity
        //
   
     for( int nn : range(NumberOfStates) )   
     {
    psi = states[nn];
    psi.noPrime();  
   
    // 'gauge' the MPS to site 1
          //  // any 'position' between 1 and L, inclusive, would work here
    psi.position(1);
    psidag = dag(psi);
       for( int kk =1; kk < L; ++kk){
        auto kkSiteInd = siteIndex(psi, kk) ;
        if(kk%2==1) psidag.prime(kkSiteInd) ;
    }
      psidag.prime("Link");
   
      auto lii_1 = leftLinkIndex(psi,1);
    psi.prime(lii_1);
    auto Par = psi(1);
         Par *= op(psites,"F",1) ; // op(psites,"Id",1) - 2*op(psites,"N", 1);
         Par *= psidag(1);
   
    for(int kk = 2; kk < L; ++kk)
    {   
        if( kk%2==1 ){ // odd sites, calculate 1-2N
            // Need to count the number of Fermions between sites ii and jj
            //// including phase for Jordan-Wigner string
            Par *= psi(kk);
            Par *= op(psites,"F",kk); // this is the operator we need!
                // op(psites,"Id",kk) - 2*op(psites,"N", kk);
              // J-W string - need to add the "F" argument in between
              //// to get proper phase with Fermions
            Par *= psidag(kk);
            std::cout << " Computed parity operator contribution for an odd site ";
        } else { // even sites, use an identity operator
            // Need to count the number of Fermions between sites ii and jj
            //// including phase for Jordan-Wigner string
            Par *= psi(kk);
            // Par *= op(psites,"F",kk); // Dont need for commuting operators?
              // J-W string - need to add the "F" argument in between
              //// to get proper phase with Fermions
            Par *= psidag(kk);
            std::cout << " Computed parity operator contribution for an even site ";
        } // else an even site
        std::cout << " Computed parity operator contribution for site " << kk  << " \n " << " \n ";
    } // over site kk

    // Final site   
    // index linking jj to jj+1 and use an identity operator
    auto lL = rightLinkIndex(psi,L);
    Par *= prime(psi(L),lL);
        // Like above for site ii, prime again at site jj (this time on the RHS)
     // Par *= op(psites,"F",L);
       Par *= psidag(L);
    std::cout << " Computed parity operator contribution for site " << L  << " \n ";
   
    std::cout << " About to print the operator Par" << " \n ";
    Print(Par);
    std::cout << " \n " << " Printed the operator Par" << " \n ";
   
    double Parityresult = elt(Par); // use eltC if expecting complex
    globalParityArray[WIndex][nn] = Parityresult;
   
     } // nn state of interest loop   
   
   
    for( int nn : range(NumberOfStates) )   
    {
    auto psinn=states[nn];
    for(int mm : range(NumberOfStates))
    {
        auto psimm=states[mm] ;
        overlapMatrix[WIndex][nn][mm]=inner(psinn,psimm);
    } // over mm
    } // over nn
   
   
    println("Energies at this parameter set are \n" );
    for( int nn : range(NumberOfStates) )   
    {
        std::cout << EnergySpectrum[WIndex][nn]  << " \n ";
    }
   
    println("Global Parities of Eigenstates at this parameter set are \n" );
    for( int nn : range(NumberOfStates) )   
    {
        std::cout << globalParityArray[WIndex][nn]  << " \n ";
    }
   
     } // WIndex ending loop    
   
   

     string parityFileName = "ParityDoubleWire_L" + str(L) + "_N" + str(numberOfParticles) + ".txt" ;
     std::ofstream paritySpectrumFile(parityFileName); // ("EnergySpectrumFile.txt");
     paritySpectrumFile << NumberWIndices << NumberOfStates << "\n";
     for (int i = 0; i < NumberWIndices; i++)
       {   
    for (int n = 0; n < NumberOfStates; n++)
         {
            paritySpectrumFile << globalParityArray[i][n] << " ";
        }
            paritySpectrumFile << "\n";
       }           
   
   
      string energyFileName = "EnergiesDoubleWire_L" + str(L) + "_N" + str(numberOfParticles) + ".txt" ;
     std::ofstream EnergySpectrumFile(energyFileName); // ("EnergySpectrumFile.txt");
     EnergySpectrumFile << NumberWIndices << NumberOfStates << "\n";
     for (int i = 0; i < NumberWIndices; i++)
       {   
    for (int n = 0; n < NumberOfStates; n++)
          {
            EnergySpectrumFile << EnergySpectrum[i][n] << " ";
        }
            EnergySpectrumFile << "\n";
         }           
   
     
     string overlapFileName = "overlapMatrixDoubleWire_L" + str(L) + "_N" + str(numberOfParticles) + ".txt" ;
     std::ofstream overlapFile(overlapFileName);  
     overlapFile << NumberWIndices*NumberOfStates << " " << L << "\n";
     for (int i = 0; i < NumberWIndices; i++)
       {   
    for (int nn = 0; nn < NumberOfStates; nn++)
      {
          for (int mm = 0; mm < NumberOfStates; mm++)
        {
            overlapFile << overlapMatrix[i][nn][mm] << " ";
        } // j
          overlapFile << "\n";
         } // n
       }    // i
   
   
   
     return 0;
    }


The energies printed are noticably out of order, and only some have a defined odd-site (single wire) parity.

       Energies at this parameter set are:
     -2.2028
     -2.22071
     -1.46725
     -2.21413
     -2.20282
     -1.95795
     -1.95795
     -1.62416
     -1.69117

       Global Parity Measurements of Eigenstates are
     -1
     1
     1
     -0.997684
     0.997674
     0.999984
     -0.999984
     0.859696
     -0.678753

Of course, the energies can be tweaked by playing with the term that is proportional to the overlap between the states, but it doesn't seem to be conserving the parity. Any insight would be greatly appreciated.

---

After moving      
 for( auto ii : range1(L) ) if( ii%2==1 ) state.set(ii,"Occ");
to just before assigning
  auto psi0 = state; ,
the issue seems to be that the states simply aren't keeping track of the parity.

Without the noise term, the energies are wildly off, but the variant of parity is conserved. With the noise term on, the energies seem correct but the parity is not conserved. Using the randomMPS then redoing the symmetrization of the initial MPS labeled "psi0" doesn't do anything to help. The issue is I think in the modified fermion.h file (called FermionSiteOddSitesParity.h here, a not great filename but it's there).

The main problem is constructing a parity operator that only acts on half of the lattice non-trivially:
// conserves total particle number and the parity of the number of particles on an odd-indexed sites
                {
                if( n%2==1 ) {
                    s = Index(QN({"Nf",0,-1}),1,QN({"Pfodd",0,-2}),1,
                              QN({"Nf",1,-1}),1,QN({"Pfodd",1,-2}),1,
                                Out,ts);
                    }
                if( n%2==0 ) {  
                    s = Index(QN({"Nf",0,-1}),1,QN({"Pfodd",0,-2}),1,
                              QN({"Nf",1,-1}),1,QN({"Pfodd",0,-2}),1,
                                Out,ts);
                    }     
                }

---

I attempted to make an asantz from the states that are of the form (1,0,1,0,...,1,0,0,...,0), (0,0,1,0,...,1,0,0,0,0), ... (0,..., 0,1,0,1,0) with particles only on odd sites with desired total number of particles with starting positions that go until the final odd site is filled. When I print this MPS:
- Sites not on the edges have things like this:
ITensor ord=3:
(dim=1|id=41|"l=13,Link") <Out>
  1: 1 QN({"Nf",0,-1})
(dim=4|id=450|"n=13,Site,Fermion") <Out>
  1: 1 QN({"Nf",0,-1})
  2: 1 QN({"Pfodd",0,-2})
  3: 1 QN({"Nf",1,-1})
  4: 1 QN({"Pfodd",1,-2})
(dim=2|id=577|"l=12,Link") <In>
  1: 2 QN({"Nf",0,-1})
{norm=1.41 (QDense Real)}

depending on the site index, the dimension of the in/out of the first and last on that list vary in dimension, but always of the form 1: i QN({"Nf",0,-1}) .

So it seems like the issue is in the definition of the quantum numbers of my altered fermion.h file (called FermionSiteOddSitesParity.h here). Will try messing with that some.

---

Editing this final comment. I commented out some code that I *thought* was just a redundancy in forming the initial ansatz (different from the version posted above), and it seems to have worked. My code is slightly different from above. I'll organize all this Kitaev Wire stuff into a document once I get it working smoothly (and tested some).

---

Hi Jared,
Thanks for being patient with my slow reply on this. So from your last comment above, did you find a solution now?

Looking at your custom site type file, it seems your solution is the one I would suggest: defining the same "Nf" quantum number for both types of sites, but for the "a" sites also defining a "Pf" parity quantum number. (Technically you don't have to define or include it at all for the "b" sites, because QN objects automatically "merge" their QN values and any that are missing are assumed to just have the value zero. You can try that out by adding two individual QN objects yourself if you want.)

My only suggestion looking at your code is that you can probably remove or simplify a lot of the if-else cases implementing user-defined options because I'm guessing you only want one of those cases for your system. The reason the library code is some complex is that we wanted to offer a lot of different options for people.

Miles

---

Hi Miles,

Thanks for taking the time to look at this. I think the issue with the alternative site definition was that I was attempting to define two independent quantum numbers, but since the phase and number operators are both tied to the same variable of whether or not that site is occupied. Once I combined them and made it so the onsite index was binary, it seemed to fix it.

Removing the extra if-thens is a good idea for this system, I don't need this to get bigger and uglier with each new problem.

It does seem that the code only works when I set up the ansatz using the state language then generate an MPS from the state, which from the documentation I thought was antiquated (the InitState documentation gives a warning about the comments being for old code: http://itensor.org/docs.cgi?page=classes/initstate&vers=cppv3).  

There was one more thing that was somewhat difficult for me to understand. From the fermion.h file, it seems that setting either to "occupied" or to "1" should give the same effect, but when I used both (as a redundancy), the code didn't preserve the abstracted parity, but I think when I did this for the Kitaev wire model, it worked fine with the redundancy.

---

I see, so that page is only saying that a certain interface to InitState is antiquated. (Perhaps I should remove that page now too, because very few people may be even using an older version.) But the InitState system is necessary to use for systems with QN conservation (using its more up-to-date interface) because MPS are required to be initialized with well defined quantum numbers, so we can't just construct a 'generic' MPS for you just from the site indices without more information. Here is an example of a code using the newest interface for InitState:
http://itensor.org/docs.cgi?vers=cppv3&page=formulas/ladder

(in the Julia version we've been able to just replace an InitState with an array of strings or numbers, and might pursue that simpler design in the future for the C++ version).

Regarding your other question, when making an InitState for Fermion sites, yes the strings "Occ" and "1" should have exactly the same meaning and effect. The string "occupied" is not one of the recognized inputs, though. I just updated the docs to reflect this: http://itensor.org/docs.cgi?vers=cppv3&page=classes/fermion

Miles

---

Oh, thank you for clarifying that about INitState. I don't know if that page needs to be removed (or if anyone uses a version of ITensor before version 2 these days), but if you added a note saying what you did in the first paragraph, it would remove the confusion.

I think I did use "Occ" or at least I did in the sample code above, but I don't think there's any need to write the same thing twice like I had before since it just invites errors.

Thanks again, Miles.

Jared

---

Good suggestions about the docs. I think what's happening there is that even though those pages aren't linked to or reachable by following links on the ITensor site, Google is indexing them anyway so they are showing up in search results and confusing people. I just updated that InitState one and removed some others that were out of date (and again aren't actually linked to on the site but can still be accessed directly if you find them in a search).

I agree it might be time to take down the ITensor v2 docs!

---

Excellent. Thanks for keeping things updated!

---

Also, if you don't mind, I'm collecting feedback from users here and there about their choice of language (C++ versus Julia). What are your current reasons for using the C++ version of ITensor versus the Julia one? Just more familiarity with C++ and/or already using that from before or another reason? Thanks -

---

I don't mind at all. I'm not very familiar at all with Julia, but more or less assumed C++ would simply compute faster (of course, even if that's true at a fundamental level, user-written code is likely the bigger factor). There two other minor reasons, like C(++) is a more general purpose language, and since I'm uncertain on the probability of a post-doc, having some background in C is a more profitable bet. A more unique-to-me reason is that I work part-time with an experimental rock physics group rather than TA in my department, and use OpenFOAM some for that - OpenFOAM uses C with C++ comments,  so using it for ITensor allows me to consolidate my syntax frustrations.

---

Thanks for that feedback! Yes it definitely makes sense given all those reasons. I'm thinking of adding a small FAQ to the ITensor site to help people weigh this choice.

Just for your info, Julia is just as fast of a language as C++, even for very low-level code. This is because Julia is a compiled language too, and all the types of variables in Julia can be known to the compiler, so just like in C++. There are even BLAS libraries now being written in pure Julia that rival Fortran BLAS! That being said, because of its greater flexibility as a language, it is easier to accidentally write slow code in Julia (slower than the comparable C++ I mean), though one can catch these slowdown by profiling the code.

Answer 1

(Marking this question as answered - see discussion above.)