How to calculate the PDF of a multi-lobe BSDF

In PBRT, a BSDF is represented as an array of BxDF lobes. Each lobe sets a bitflag representing what kind of lobe it is. (Specular, Diffuse, Reflection, Transmission, etc.) I'm investigating how PBRT implemented sampling a BSDF for an input direction, wi.

Following standard Monte Carlo sampling, I was expecting the code to randomly sample an appropriate lobe, then set the $$pdf = pdf_{picked}/n$$ where $$n$$ is the number of lobes of the appropriate type. However this isn't what PBRT does. My question is why? To go in futher, this is the PBRT function:

Spectrum BSDF::Sample_f(const Vector3f &woWorld, Vector3f *wiWorld,
const Point2f &u, Float *pdf, BxDFType type,
BxDFType *sampledType) const {
// Choose which _BxDF_ to sample
int matchingComps = NumComponents(type);
if (matchingComps == 0) {
*pdf = 0;
if (sampledType)
*sampledType = BxDFType(0);
return Spectrum(0);
}
int comp = std::min((int)std::floor(u * matchingComps), matchingComps - 1);

// Get _BxDF_ pointer for chosen component
BxDF *bxdf = nullptr;
int count = comp;
for (int i = 0; i < nBxDFs; ++i)
if (bxdfs[i]->MatchesFlags(type) && count-- == 0) {
bxdf = bxdfs[i];
break;
}

// Remap _BxDF_ sample _u_ to $$[0,1)^2$$
Point2f uRemapped(std::min(u * matchingComps - comp, OneMinusEpsilon),
u);

// Sample chosen _BxDF_
Vector3f wi, wo = WorldToLocal(woWorld);
if (wo.z == 0)
return 0.;
*pdf = 0;
if (sampledType) *sampledType = bxdf->type;
Spectrum f = bxdf->Sample_f(wo, &wi, uRemapped, pdf, sampledType);

if (*pdf == 0) {
if (sampledType) *sampledType = BxDFType(0);
return 0;
}
*wiWorld = LocalToWorld(wi);

// Compute overall PDF with all matching _BxDF_s
if (!(bxdf->type & BSDF_SPECULAR) && matchingComps > 1)
for (int i = 0; i < nBxDFs; ++i)
if (bxdfs[i] != bxdf && bxdfs[i]->MatchesFlags(type))
*pdf += bxdfs[i]->Pdf(wo, wi);
if (matchingComps > 1)
*pdf /= matchingComps;

// Compute value of BSDF for sampled direction
if (!(bxdf->type & BSDF_SPECULAR)) {
bool reflect = Dot(*wiWorld, ng) * Dot(woWorld, ng) > 0;
f = 0.;
for (int i = 0; i < nBxDFs; ++i)
if (bxdfs[i]->MatchesFlags(type) &&
((reflect && (bxdfs[i]->type & BSDF_REFLECTION)) ||
(!reflect && (bxdfs[i]->type & BSDF_TRANSMISSION))))
f += bxdfs[i]->f(wo, wi);
}

return f;
}

In pseudo-code, it looks to be doing the following:

// Get the number of lobes that match the ones we asked for
numMatchingLobes = NumComponents(type);

// Early out check
if (numMatchingLobes == 0):
return

// Randomly pick one of the matching lobes
pickedBxDF = Rand()

// Transform wo from world space to local space
// Then sample the chosen BxDF to get f, pdf, and wi
f, pdf, wi = pickedBxDF.sample_f(woWorld.ToLocal())

// Early out check
if (pdf == 0):
return

// Transform wi to world space
wiWorld = wi.ToWorld()

// Add all of the pdfs of the matching lobes if the picked BxDF wasn't specular
if pickedBxDF.type not Specular && numMatchingLobes > 1:
for BxDF in matchingLobes:
pdf += BxDF.pdf(wi, wo)

// Divide the resulting pdf by the number of matching lobes
pdf /= numMatchingLobes;

// Calculate the combined f if the picked BxDF isn't specular
if pickedBxDF.type not Specular:
// Reset f to zero
f = 0

// Add each BxDF's f, iff it matches reflection / transmission
reflect = wo and wi on same hemisphere as geoNormal
for BxDF in matchingLobes:
if reflect and BxDF.type == Reflection:
f += BxDF.f(wo, wi)
elif !reflect and BxDF.type == Transmission:
f += BxDF.f(wo, wi)

The part I'm not understanding is why it adds the pdf's of all the matching lobes on top of the already calculated pdf. This seems like it would double count the pdf of the picked lobe.

For example, if there are two matching lobes, called A and B. A is picked.

We calculate f, pdf, and wi with A. Then with a for loop, we add the pdf for A and B and divide by 2. Aka:

$$\frac{A + A + B}{2}$$

This doesn't seem right. Is PBRT code wrong? Should it be resetting the pdf, similar to the f calculation code below it?