diff --git a/cpp/daal/src/algorithms/dtrees/gbt/regression/gbt_regression_predict_dense_default_batch_impl.i b/cpp/daal/src/algorithms/dtrees/gbt/regression/gbt_regression_predict_dense_default_batch_impl.i index 8bc60f19307..6e5858d19ed 100644 --- a/cpp/daal/src/algorithms/dtrees/gbt/regression/gbt_regression_predict_dense_default_batch_impl.i +++ b/cpp/daal/src/algorithms/dtrees/gbt/regression/gbt_regression_predict_dense_default_batch_impl.i @@ -199,36 +199,36 @@ protected: template inline services::Status predictContributionInteractions(size_t iTree, size_t nTrees, size_t nRowsData, const algorithmFPType * x, - algorithmFPType * res, const DimType & dim); + const algorithmFPType * nominal, algorithmFPType * res, const DimType & dim); template inline services::Status predictContributionInteractions(size_t iTree, size_t nTrees, size_t nRowsData, const algorithmFPType * x, - algorithmFPType * res, const DimType & dim) + const algorithmFPType * nominal, algorithmFPType * res, const DimType & dim) { if (_featHelper.hasUnorderedFeatures()) { - return predictContributionInteractions(iTree, nTrees, nRowsData, x, res, dim); + return predictContributionInteractions(iTree, nTrees, nRowsData, x, nominal, res, dim); } else { - return predictContributionInteractions(iTree, nTrees, nRowsData, x, res, dim); + return predictContributionInteractions(iTree, nTrees, nRowsData, x, nominal, res, dim); } } // TODO: Add vectorBlockSize templates, similar to predict // template inline services::Status predictContributionInteractions(size_t iTree, size_t nTrees, size_t nRowsData, const algorithmFPType * x, - algorithmFPType * res, const DimType & dim) + const algorithmFPType * nominal, algorithmFPType * res, const DimType & dim) { const size_t nColumnsData = dim.nCols; const bool hasAnyMissing = checkForMissing(x, nTrees, nRowsData, nColumnsData); if (hasAnyMissing) { - return predictContributionInteractions(iTree, nTrees, nRowsData, x, res, dim); + return predictContributionInteractions(iTree, nTrees, nRowsData, x, nominal, res, dim); } else { - return predictContributionInteractions(iTree, nTrees, nRowsData, x, res, dim); + return predictContributionInteractions(iTree, nTrees, nRowsData, x, nominal, res, dim); } } @@ -352,15 +352,17 @@ void PredictRegressionTask::predictContributions(size_t iT const gbt::internal::GbtDecisionTree * currentTree = _aTree[currentTreeIndex]; const void * endAddr = static_cast(&(*uniquePathData.end())); - printf("\n\n\n--> depth = %d | uniquePathData.end() = %p\n\n", depth, endAddr); gbt::treeshap::treeShap(currentTree, currentX, phi, nColumnsData, &_featHelper, uniquePathData.data(), condition, conditionFeature); - printf("treeShap is done\n"); } - for (int iFeature = 0; iFeature < nColumnsData; ++iFeature) + if (condition == 0) { - phi[biasTermIndex] -= phi[iFeature]; + // find bias term by leveraging bias = nominal - sum_i phi_i + for (int iFeature = 0; iFeature < nColumnsData; ++iFeature) + { + phi[biasTermIndex] -= phi[iFeature]; + } } } } @@ -379,12 +381,16 @@ void PredictRegressionTask::predictContributions(size_t iT template template services::Status PredictRegressionTask::predictContributionInteractions(size_t iTree, size_t nTrees, size_t nRowsData, - const algorithmFPType * x, algorithmFPType * res, + const algorithmFPType * x, + const algorithmFPType * nominal, algorithmFPType * res, const DimType & dim) { Status st; - const size_t nColumnsData = dim.nCols; - const size_t nColumnsPhi = nColumnsData + 1; + const size_t nColumnsData = dim.nCols; + const size_t nColumnsPhi = nColumnsData + 1; + const size_t biasTermIndex = nColumnsPhi - 1; + + const size_t interactionMatrixSize = nColumnsPhi * nColumnsPhi; // Allocate buffer for 3 matrices for algorithmFPType of size (nRowsData, nColumnsData) const size_t elementsInMatrix = nRowsData * nColumnsPhi; @@ -396,33 +402,47 @@ services::Status PredictRegressionTask::predictContributio return st; } - // Initialize buffer - service_memset_seq(buffer, algorithmFPType(0), 3 * elementsInMatrix); - // Get pointers into the buffer for our three matrices algorithmFPType * contribsDiag = buffer + 0 * elementsInMatrix; algorithmFPType * contribsOff = buffer + 1 * elementsInMatrix; algorithmFPType * contribsOn = buffer + 2 * elementsInMatrix; + // Initialize nominal buffer + service_memset_seq(contribsDiag, algorithmFPType(0), elementsInMatrix); + + // Copy nominal values (for bias term) to the condition == 0 buffer + PRAGMA_IVDEP + PRAGMA_VECTOR_ALWAYS + for (size_t i = 0; i < nRowsData; ++i) + { + contribsDiag[i * nColumnsPhi + biasTermIndex] = nominal[i]; + } + predictContributions(iTree, nTrees, nRowsData, x, contribsDiag, 0, 0, dim); - for (size_t i = 0; i < nColumnsPhi + 1; ++i) + for (size_t i = 0; i < nColumnsPhi; ++i) { + // initialize/reset the on/off buffers + service_memset_seq(contribsOff, algorithmFPType(0), 2 * elementsInMatrix); + predictContributions(iTree, nTrees, nRowsData, x, contribsOff, -1, i, dim); predictContributions(iTree, nTrees, nRowsData, x, contribsOn, 1, i, dim); for (size_t j = 0; j < nRowsData; ++j) { - for (size_t k = 0; k < nColumnsPhi + 1; ++k) + const unsigned o_offset = j * interactionMatrixSize + i * nColumnsPhi; + const unsigned c_offset = j * nColumnsPhi; + res[o_offset + i] = 0; + for (size_t k = 0; k < nColumnsPhi; ++k) { // fill in the diagonal with additive effects, and off-diagonal with the interactions if (k == i) { - res[i] += contribsDiag[k]; + res[o_offset + i] += contribsDiag[c_offset + k]; } else { - res[k] = (contribsOn[k] - contribsOff[k]) / 2.0; - res[i] -= res[k]; + res[o_offset + k] = (contribsOn[c_offset + k] - contribsOff[c_offset + k]) / 2.0f; + res[o_offset + i] -= res[o_offset + k]; } } } @@ -483,11 +503,20 @@ services::Status PredictRegressionTask::runInternal(servic WriteOnlyRows resRow(result, iStartRow, nRowsToProcess); DAAL_CHECK_BLOCK_STATUS_THR(resRow); - // TODO: Might need the nominal prediction to account for bias terms - // predict(iTree, nTreesToUse, nRowsToProcess, xBD.get(), resRow.get(), dim, resultNColumns); + // nominal values are required to calculate the correct bias term + algorithmFPType * nominal = static_cast(daal_malloc(nRowsToProcess * sizeof(algorithmFPType))); + if (!nominal) + { + safeStat.add(ErrorMemoryAllocationFailed); + return; + } + service_memset_seq(nominal, algorithmFPType(0), nRowsToProcess); + predict(iTree, nTreesToUse, nRowsToProcess, xBD.get(), nominal, dim, 1); // TODO: support tree weights - safeStat |= predictContributionInteractions(iTree, nTreesToUse, nRowsToProcess, xBD.get(), resRow.get(), dim); + safeStat |= predictContributionInteractions(iTree, nTreesToUse, nRowsToProcess, xBD.get(), nominal, resRow.get(), dim); + + daal_free(nominal); } else { diff --git a/cpp/daal/src/algorithms/dtrees/gbt/treeshap.cpp b/cpp/daal/src/algorithms/dtrees/gbt/treeshap.cpp index 6c5fe008d56..85cc4c175f4 100644 --- a/cpp/daal/src/algorithms/dtrees/gbt/treeshap.cpp +++ b/cpp/daal/src/algorithms/dtrees/gbt/treeshap.cpp @@ -12,39 +12,30 @@ namespace internal { // extend our decision path with a fraction of one and zero extensions -void treeShapExtendPath(PathElement * uniquePath, size_t uniqueDepth, float zeroFraction, float oneFraction, FeatureIndexType featureIndex) +void extendPath(PathElement * uniquePath, size_t uniqueDepth, float zeroFraction, float oneFraction, int featureIndex) { uniquePath[uniqueDepth].featureIndex = featureIndex; uniquePath[uniqueDepth].zeroFraction = zeroFraction; uniquePath[uniqueDepth].oneFraction = oneFraction; uniquePath[uniqueDepth].partialWeight = (uniqueDepth == 0 ? 1.0f : 0.0f); + const float constant = 1.0f / static_cast(uniqueDepth + 1); for (int i = uniqueDepth - 1; i >= 0; i--) { - uniquePath[i + 1].partialWeight += oneFraction * uniquePath[i].partialWeight * (i + 1) / static_cast(uniqueDepth + 1); - uniquePath[i].partialWeight = zeroFraction * uniquePath[i].partialWeight * (uniqueDepth - i) / static_cast(uniqueDepth + 1); + uniquePath[i + 1].partialWeight += oneFraction * uniquePath[i].partialWeight * (i + 1) * constant; + uniquePath[i].partialWeight = zeroFraction * uniquePath[i].partialWeight * (uniqueDepth - i) * constant; } } // undo a previous extension of the decision path -void treeShapUnwindPath(PathElement * uniquePath, size_t uniqueDepth, size_t pathIndex) +void unwindPath(PathElement * uniquePath, size_t uniqueDepth, size_t pathIndex) { - printf("treeShapUnwindPath: Going through path elements\n"); - printf("uniquePath = %p\n", uniquePath); - printf("uniqueDepth = %lu\n", uniqueDepth); - printf("pathIndex = %lu\n", pathIndex); - printf("---- start\n"); - printf("%p\n", uniquePath + pathIndex); - const float oneFraction = uniquePath[pathIndex].oneFraction; const float zeroFraction = uniquePath[pathIndex].zeroFraction; - - printf("%p\n", uniquePath + uniqueDepth); - float nextOnePortion = uniquePath[uniqueDepth].partialWeight; + float nextOnePortion = uniquePath[uniqueDepth].partialWeight; for (int i = uniqueDepth - 1; i >= 0; --i) { - printf("%p\n", uniquePath + i); if (oneFraction != 0) { const float tmp = uniquePath[i].partialWeight; @@ -59,7 +50,6 @@ void treeShapUnwindPath(PathElement * uniquePath, size_t uniqueDepth, size_t pat for (size_t i = pathIndex; i < uniqueDepth; ++i) { - printf("%p <- %p\n", uniquePath + i, uniquePath + i + 1); uniquePath[i].featureIndex = uniquePath[i + 1].featureIndex; uniquePath[i].zeroFraction = uniquePath[i + 1].zeroFraction; uniquePath[i].oneFraction = uniquePath[i + 1].oneFraction; @@ -67,19 +57,11 @@ void treeShapUnwindPath(PathElement * uniquePath, size_t uniqueDepth, size_t pat } // determine what the total permutation weight would be if we unwound a previous extension in the decision path -float treeShapUnwoundPathSum(const PathElement * uniquePath, size_t uniqueDepth, size_t pathIndex) +float unwoundPathSum(const PathElement * uniquePath, size_t uniqueDepth, size_t pathIndex) { - printf("treeShapUnwoundPathSum: Going through path elements\n"); - printf("uniquePath = %p\n", uniquePath); - printf("uniqueDepth = %lu\n", uniqueDepth); - printf("pathIndex = %lu\n", pathIndex); - printf("---- start\n"); - printf("%p\n", uniquePath + pathIndex); - const float oneFraction = uniquePath[pathIndex].oneFraction; const float zeroFraction = uniquePath[pathIndex].zeroFraction; - printf("%p\n", uniquePath + uniqueDepth); float nextOnePortion = uniquePath[uniqueDepth].partialWeight; float total = 0; // if (oneFraction != 0) @@ -109,8 +91,6 @@ float treeShapUnwoundPathSum(const PathElement * uniquePath, size_t uniqueDepth, for (int i = uniqueDepth - 1; i >= 0; --i) { - printf("%p\n", uniquePath + i); - if (oneFraction != 0) { const float tmp = nextOnePortion * (uniqueDepth + 1) / static_cast((i + 1) * oneFraction); diff --git a/cpp/daal/src/algorithms/dtrees/gbt/treeshap.h b/cpp/daal/src/algorithms/dtrees/gbt/treeshap.h index bba13c0dbf5..0235aaa6859 100644 --- a/cpp/daal/src/algorithms/dtrees/gbt/treeshap.h +++ b/cpp/daal/src/algorithms/dtrees/gbt/treeshap.h @@ -53,8 +53,7 @@ using FeatureTypes = algorithms::dtrees::internal::FeatureTypes; // the partialWeight of the i'th path element is the permutation weight of paths with i-1 ones in them struct PathElement { - // featureIndex -1 underflows, we use it as a reserved value for the initial node - FeatureIndexType featureIndex = -1; + int featureIndex = 0xaaaaaaaa; float zeroFraction = 0; float oneFraction = 0; float partialWeight = 0; @@ -65,9 +64,9 @@ struct PathElement namespace internal { -void treeShapExtendPath(PathElement * uniquePath, size_t uniqueDepth, float zeroFraction, float oneFraction, FeatureIndexType featureIndex); -void treeShapUnwindPath(PathElement * uniquePath, size_t uniqueDepth, size_t pathIndex); -float treeShapUnwoundPathSum(const PathElement * uniquePath, size_t uniqueDepth, size_t pathIndex); +void extendPath(PathElement * uniquePath, size_t uniqueDepth, float zeroFraction, float oneFraction, int featureIndex); +void unwindPath(PathElement * uniquePath, size_t uniqueDepth, size_t pathIndex); +float unwoundPathSum(const PathElement * uniquePath, size_t uniqueDepth, size_t pathIndex); /** Extension of * \param nodeIndex the index of the current node in the tree, counted from 1 @@ -81,8 +80,8 @@ float treeShapUnwoundPathSum(const PathElement * uniquePath, size_t uniqueDepth, template void treeShap(const gbt::internal::GbtDecisionTree * tree, const algorithmFPType * x, algorithmFPType * phi, size_t nFeatures, FeatureTypes * featureHelper, size_t nodeIndex, size_t depth, size_t uniqueDepth, PathElement * parentUniquePath, - float parentZeroFraction, float parentOneFraction, FeatureIndexType parentFeatureIndex, int condition, - FeatureIndexType conditionFeature, float conditionFraction) + float parentZeroFraction, float parentOneFraction, int parentFeatureIndex, int condition, FeatureIndexType conditionFeature, + float conditionFraction) { DAAL_ASSERT(parentUniquePath); @@ -100,21 +99,15 @@ void treeShap(const gbt::internal::GbtDecisionTree * tree, const algorithmFPType const int * const defaultLeft = tree->getDefaultLeftForSplit() - 1; PathElement * uniquePath = parentUniquePath + uniqueDepth + 1; - const size_t nBytes = uniqueDepth * sizeof(PathElement); + const size_t nBytes = (uniqueDepth + 1) * sizeof(PathElement); const int copy_status = daal::services::internal::daal_memcpy_s(uniquePath, nBytes, parentUniquePath, nBytes); DAAL_ASSERT(copy_status == 0); - if (condition == 0 || conditionFeature != parentFeatureIndex) + if (condition == 0 || conditionFeature != static_cast(parentFeatureIndex)) { - treeShapExtendPath(uniquePath, uniqueDepth, parentZeroFraction, parentOneFraction, parentFeatureIndex); + extendPath(uniquePath, uniqueDepth, parentZeroFraction, parentOneFraction, parentFeatureIndex); } - printf("--------------------------------------------------------------------------------\n"); - printf("depth = %lu\n", depth); - printf("uniqueDepth = %lu\n", uniqueDepth); - printf("uniquePath = %p\n", uniquePath); - printf("uniquePath + uniqueDepth = %p\n", uniquePath + uniqueDepth); - const bool isLeaf = gbt::internal::ModelImpl::nodeIsLeaf(nodeIndex, *tree, depth); // leaf node @@ -122,7 +115,7 @@ void treeShap(const gbt::internal::GbtDecisionTree * tree, const algorithmFPType { for (size_t i = 1; i <= uniqueDepth; ++i) { - const float w = treeShapUnwoundPathSum(uniquePath, uniqueDepth, i); + const float w = unwoundPathSum(uniquePath, uniqueDepth, i); const PathElement & el = uniquePath[i]; phi[el.featureIndex] += w * (el.oneFraction - el.zeroFraction) * splitValues[nodeIndex] * conditionFraction; } @@ -130,11 +123,11 @@ void treeShap(const gbt::internal::GbtDecisionTree * tree, const algorithmFPType return; } - const FeatureIndexType splitFeature = fIndexes[nodeIndex]; - const auto dataValue = x[splitFeature]; + const FeatureIndexType splitIndex = fIndexes[nodeIndex]; + const algorithmFPType dataValue = x[splitIndex]; gbt::prediction::internal::PredictDispatcher dispatcher; - FeatureIndexType hotIndex = updateIndex(nodeIndex, dataValue, splitValues, defaultLeft, *featureHelper, splitFeature, dispatcher); + FeatureIndexType hotIndex = updateIndex(nodeIndex, x[splitIndex], splitValues, defaultLeft, *featureHelper, splitIndex, dispatcher); const FeatureIndexType coldIndex = 2 * nodeIndex + (hotIndex == (2 * nodeIndex)); const float w = nodeCoverValues[nodeIndex]; @@ -152,7 +145,12 @@ void treeShap(const gbt::internal::GbtDecisionTree * tree, const algorithmFPType size_t previousSplitPathIndex = 0; for (; previousSplitPathIndex <= uniqueDepth; ++previousSplitPathIndex) { - if (uniquePath[previousSplitPathIndex].featureIndex == splitFeature) + const FeatureIndexType castIndex = static_cast(uniquePath[previousSplitPathIndex].featureIndex); + + // It cannot be that a feature that is ignored is in the uniquePath + DAAL_ASSERT((condition == 0) || (castIndex != conditionFeature)); + + if (castIndex == splitIndex) { break; } @@ -161,20 +159,19 @@ void treeShap(const gbt::internal::GbtDecisionTree * tree, const algorithmFPType { incomingZeroFraction = uniquePath[previousSplitPathIndex].zeroFraction; incomingOneFraction = uniquePath[previousSplitPathIndex].oneFraction; - treeShapUnwindPath(uniquePath, uniqueDepth, previousSplitPathIndex); - printf("previousSplitPathIndex != uniqueDepth + 1 : %lu != %lu\n", previousSplitPathIndex, uniqueDepth + 1); + unwindPath(uniquePath, uniqueDepth, previousSplitPathIndex); uniqueDepth -= 1; } // divide up the conditionFraction among the recursive calls float hotConditionFraction = conditionFraction; float coldConditionFraction = conditionFraction; - if (condition > 0 && splitFeature == conditionFeature) + if (condition > 0 && splitIndex == conditionFeature) { coldConditionFraction = 0; uniqueDepth -= 1; } - else if (condition < 0 && splitFeature == conditionFeature) + else if (condition < 0 && splitIndex == conditionFeature) { hotConditionFraction *= hotZeroFraction; coldConditionFraction *= coldZeroFraction; @@ -183,9 +180,9 @@ void treeShap(const gbt::internal::GbtDecisionTree * tree, const algorithmFPType treeShap(tree, x, phi, nFeatures, featureHelper, hotIndex, depth + 1, uniqueDepth + 1, uniquePath, hotZeroFraction * incomingZeroFraction, incomingOneFraction, - splitFeature, condition, conditionFeature, hotConditionFraction); + splitIndex, condition, conditionFeature, hotConditionFraction); treeShap(tree, x, phi, nFeatures, featureHelper, coldIndex, depth + 1, uniqueDepth + 1, - uniquePath, coldZeroFraction * incomingZeroFraction, 0, splitFeature, condition, + uniquePath, coldZeroFraction * incomingZeroFraction, 0, splitIndex, condition, conditionFeature, coldConditionFraction); } } // namespace internal @@ -209,7 +206,6 @@ void treeShap(const gbt::internal::GbtDecisionTree * tree, const algorithmFPType DAAL_ASSERT(phi); DAAL_ASSERT(featureHelper); - // parentFeatureIndex -1 underflows, we use it as a reserved value for the initial node treeshap::internal::treeShap(tree, x, phi, nFeatures, featureHelper, 1, 0, 0, parentUniquePath, 1, 1, -1, condition, conditionFeature, 1); }