GFWPowerArray.cxx

#include "GFWPowerArray.h"
int GFWPowerArray::getHighestHarmonic(const HarSet &inhar) {
  //Highest possible harmonic: sum of same-sign harmonics
  int maxPos=0, maxNeg=0;
  for(int val:inhar) if(val>0) maxPos+=val; else maxNeg+=abs(val);
  return maxPos>maxNeg?maxPos:maxNeg;
};
HarSet GFWPowerArray::TrimVec(HarSet hars, int ind) {
  HarSet retVec = hars;
  retVec.erase(retVec.begin()+ind);
  return retVec;
};
HarSet GFWPowerArray::AddConstant(HarSet hars, int offset) {
  HarSet retVec = hars;
  for(int &val : retVec) val+=offset;
  return retVec;
};
void GFWPowerArray::FlushVectorToMaster(HarSet &masterVector, HarSet &comVec, const int &MaxPower) {
  int nPartLoc = MaxPower-comVec.size()+1;
  for(auto &val: comVec) {
    int absVal = abs(val);
    if(masterVector.at(absVal) < nPartLoc) {
      masterVector.at(absVal) = nPartLoc;
    };
  };
};
void GFWPowerArray::RecursiveFunction(HarSet &masterVector, HarSet hars, int offset, const int &MaxPower) {
  HarSet compVec = AddConstant(hars,offset);
  FlushVectorToMaster(masterVector, compVec, MaxPower);
  for(int i=0;i<hars.size();i++) RecursiveFunction(masterVector,TrimVec(hars,i),offset+hars.at(i),MaxPower);;
};
void GFWPowerArray::PrintVector(const HarSet &singleSet) {
  int vcSize = (int)singleSet.size();
  if(!vcSize) printf("Vector is empty!\n");
  printf("{%i",singleSet[0]);
  for(int i=1;i<vcSize;i++) printf(", %i",singleSet[i]);
  printf("}\n");
}
HarSet GFWPowerArray::GetPowerArray(vector<HarSet> inHarmonics) {
  //First, find maximum number of particle correlations ( = max power) and maximum (sum of) harmonics
  int MaxHar=0;
  int nMaxPart=0;
  for(HarSet singleSet: inHarmonics) {
    int harSum = getHighestHarmonic(singleSet);
    MaxHar=harSum>MaxHar?harSum:MaxHar;
  };
  //Make a vector with MaxHar+1 entries (entry 0 for sum=0)
  HarSet retVec = HarSet(MaxHar+1);
  //Then loop over all combinations and calculate max powers
  for(HarSet singleSet: inHarmonics) {
    int lNPart = (int)singleSet.size(); //Total number of particles correlated
    RecursiveFunction(retVec,singleSet,0,lNPart);
    //Harmonic sum = 0 is a special case. In principle all 0 cases with non-zero harmonics are captured by the function above, but to calculate normalization, we set all harmonics to 0. This means that sum=0 power is the max number of harmonics/particles being correlated
    nMaxPart=(lNPart>nMaxPart)?lNPart:nMaxPart;
  };
  //Override the sum=0 power with the number of correlated particles
  if(retVec[0]<nMaxPart) retVec[0]=nMaxPart;
  //Need an extra power ( = 0) for all non-zero powers
  for(int &val:retVec) if(val!=0) val++;
  return retVec;
};
void GFWPowerArray::PowerArrayTest() {
  vector<HarSet> AllHars = {
    HarSet{2},
    HarSet{3},
    HarSet{2, 2},
    HarSet{3, 3}
  };
  printf("Input harmonics are:\n");
  for(HarSet inSet:AllHars) PrintVector(inSet);
  printf("The configuration of powers must then be:\n");
  auto vc = GetPowerArray(AllHars);
  PrintVector(vc);
};