Skip to content

Commit

Permalink
added cache-friendly sum-computation for permutation
Browse files Browse the repository at this point in the history
  • Loading branch information
@lambday
lambday committed Jun 28, 2016
1 parent 98e24ef commit 84305f0
Show file tree
Hide file tree
Showing 5 changed files with 126 additions and 83 deletions.
3 changes: 2 additions & 1 deletion src/shogun/statistical_testing/MMD.cpp
View file
Expand Up @@ -89,6 +89,7 @@ void CMMD::Self::create_statistic_job()
{
const DataManager& dm = owner.get_data_manager();
auto Bx = dm.blocksize_at(0);
auto By = dm.blocksize_at(1);
switch (statistic_type)
{
case EStatisticType::UNBIASED_FULL:
Expand All @@ -102,7 +103,7 @@ void CMMD::Self::create_statistic_job()
break;
default : break;
};
permutation_job = mmd::WithinBlockPermutation(Bx, statistic_type);
permutation_job = mmd::WithinBlockPermutation(Bx, By, statistic_type);
}

void CMMD::Self::create_variance_job()
Expand Down
3 changes: 2 additions & 1 deletion src/shogun/statistical_testing/QuadraticTimeMMD.cpp
View file
Expand Up @@ -89,6 +89,7 @@ void CQuadraticTimeMMD::Self::create_statistic_job()
SG_SDEBUG("Entering\n");
const DataManager& dm=owner.get_data_manager();
auto Nx=dm.num_samples_at(0);
auto Ny=dm.num_samples_at(1);
switch (owner.get_statistic_type())
{
case EStatisticType::UNBIASED_FULL:
Expand All @@ -102,7 +103,7 @@ void CQuadraticTimeMMD::Self::create_statistic_job()
break;
default : break;
};
permutation_job = mmd::WithinBlockPermutation(Nx, owner.get_statistic_type());
permutation_job = mmd::WithinBlockPermutation(Nx, Ny, owner.get_statistic_type());
SG_SDEBUG("Leaving\n");
}

Expand Down
120 changes: 63 additions & 57 deletions src/shogun/statistical_testing/internals/mmd/WithinBlockPermutation.cpp
View file
Expand Up @@ -16,6 +16,7 @@
* along with this program. If not, see <http:/www.gnu.org/licenses/>.
*/

#include <unordered_map>
#include <shogun/io/SGIO.h>
#include <shogun/lib/SGMatrix.h>
#include <shogun/lib/SGVector.h>
Expand All @@ -32,84 +33,89 @@ using namespace shogun;
using namespace internal;
using namespace mmd;

WithinBlockPermutation::WithinBlockPermutation(index_t n, EStatisticType type)
: n_x(n), stype(type)
WithinBlockPermutation::WithinBlockPermutation(index_t nx, index_t ny, EStatisticType type)
: n_x(nx), n_y(ny), stype(type), terms()
{
SG_SDEBUG("number of samples are %d and %d!\n", n_x, n_y);
std::fill(&terms.term[0], &terms.term[2]+1, 0);
std::fill(&terms.diag[0], &terms.diag[2]+1, 0);
}

void WithinBlockPermutation::add_term(float64_t val, index_t i, index_t j)
{
if (i<n_x && j<n_x)
{
SG_SDEBUG("Adding KernelMatrix(%d,%d)=%f to term 1!\n", i, j, val);
terms.term[0]+=val;
if (i==j)
terms.diag[0]+=val;
}
else if (i>=n_x && j>=n_x)
{
SG_SDEBUG("Adding KernelMatrix(%d,%d)=%f to term 2!\n", i, j, val);
terms.term[1]+=val;
if (i==j)
terms.diag[1]+=val;
}
else if (i>=n_x && j<n_x)
{
SG_SDEBUG("Adding KernelMatrix(%d,%d)=%f to term 3!\n", i, j, val);
terms.term[2]+=val;
if (i-n_x==j)
terms.diag[2]+=val;
}
}

float64_t WithinBlockPermutation::operator()(SGMatrix<float64_t> km)
{
SG_SDEBUG("Entering!\n");
inds.resize(km.num_rows);
std::iota(inds.data(), inds.data()+inds.size(), 0);
SGVector<index_t> permuted_inds(inds.data(), inds.size(), false);
SGVector<index_t> permuted_inds(n_x+n_y);
std::iota(permuted_inds.vector, permuted_inds.vector+permuted_inds.vlen, 0);
CMath::permute(permuted_inds);

const index_t n_y=km.num_rows-n_x;
SG_SDEBUG("number of samples are %d and %d!\n", n_x, n_y);
std::unordered_map<index_t, index_t> inds;
for (int i=0; i<permuted_inds.vlen; ++i)
inds.insert(std::make_pair(permuted_inds[i], i));

auto term_1=0.0;
for (auto i=0; i<n_x; ++i)
{
for (auto j=0; j<n_x; ++j)
{
if (i>j)
term_1+=km(permuted_inds[i], permuted_inds[j]);
}
}
term_1*=2;
SG_SDEBUG("term_1 sum (without diagonal) = %f!\n", term_1);
if (stype==EStatisticType::BIASED_FULL)
{
for (auto i=0; i<n_x; ++i)
term_1+=km(permuted_inds[i], permuted_inds[i]);
SG_SDEBUG("term_1 sum (with diagonal) = %f!\n", term_1);
term_1/=n_x*n_x;
}
else
term_1/=n_x*(n_x-1);
SG_SDEBUG("term_1 (normalized) = %f!\n", term_1);
std::fill(&terms.term[0], &terms.term[2]+1, 0);
std::fill(&terms.diag[0], &terms.diag[2]+1, 0);

auto term_2=0.0;
for (auto i=n_x; i<n_x+n_y; ++i)
for (auto j=0; j<n_x+n_y; ++j)
{
for (auto j=n_x; j<n_x+n_y; ++j)
{
if (i>j)
term_2+=km(permuted_inds[i], permuted_inds[j]);
}
for (auto i=0; i<n_x+n_y; ++i)
add_term(km(i, j), inds.find(i)->second, inds.find(j)->second);
}
term_2*=2.0;
SG_SDEBUG("term_2 sum (without diagonal) = %f!\n", term_2);
if (stype==EStatisticType::BIASED_FULL)

SG_SDEBUG("term_1 sum (with diagonal) = %f!\n", terms.term[0]);
SG_SDEBUG("term_2 sum (with diagonal) = %f!\n", terms.term[1]);
if (stype!=EStatisticType::BIASED_FULL)
{
for (auto i=n_x; i<n_x+n_y; ++i)
term_2+=km(permuted_inds[i], permuted_inds[i]);
SG_SDEBUG("term_2 sum (with diagonal) = %f!\n", term_2);
term_2/=n_y*n_y;
terms.term[0]-=terms.diag[0];
terms.term[1]-=terms.diag[1];
SG_SDEBUG("term_1 sum (without diagonal) = %f!\n", terms.term[0]);
SG_SDEBUG("term_2 sum (without diagonal) = %f!\n", terms.term[1]);
terms.term[0]/=n_x*(n_x-1);
terms.term[1]/=n_y*(n_y-1);
}
else
term_2/=n_y*(n_y-1);
SG_SDEBUG("term_2 (normalized) = %f!\n", term_2);

auto term_3=0.0;
for (auto i=n_x; i<n_x+n_y; ++i)
{
for (auto j=0; j<n_x; ++j)
term_3+=km(permuted_inds[i], permuted_inds[j]);
terms.term[0]/=n_x*n_x;
terms.term[1]/=n_y*n_y;
}
SG_SDEBUG("term_3 sum (with diagonal) = %f!\n", term_3);
SG_SDEBUG("term_1 (normalized) = %f!\n", terms.term[0]);
SG_SDEBUG("term_2 (normalized) = %f!\n", terms.term[1]);

SG_SDEBUG("term_3 sum (with diagonal) = %f!\n", terms.term[2]);
if (stype==EStatisticType::UNBIASED_INCOMPLETE)
{
for (auto i=0; i<n_x; ++i)
term_3-=km(permuted_inds[i+n_x], permuted_inds[i]);
SG_SDEBUG("term_3 sum (without diagonal) = %f!\n", term_3);
term_3/=n_x*(n_x-1);
terms.term[2]-=terms.diag[2];
SG_SDEBUG("term_3 sum (without diagonal) = %f!\n", terms.term[2]);
terms.term[2]/=n_x*(n_x-1);
}
else
term_3/=n_x*n_y;
SG_SDEBUG("term_3 (normalized) = %f!\n", term_3);
terms.term[2]/=n_x*n_y;
SG_SDEBUG("term_3 (normalized) = %f!\n", terms.term[2]);

SG_SDEBUG("Leaving!\n");
return term_1+term_2-2*term_3;
return terms.term[0]+terms.term[1]-2*terms.term[2];
}
15 changes: 12 additions & 3 deletions src/shogun/statistical_testing/internals/mmd/WithinBlockPermutation.h
View file
Expand Up @@ -34,15 +34,24 @@ namespace internal
namespace mmd
{

struct WithinBlockPermutation
class WithinBlockPermutation
{
using return_type = float64_t;
WithinBlockPermutation(index_t, EStatisticType);
public:
WithinBlockPermutation(index_t, index_t, EStatisticType);
return_type operator()(SGMatrix<float64_t> kernel_matrix);
// return_type operator()(CGPUMatrix<float64_t> kernel_matrix);
private:
const index_t n_x;
const index_t n_y;
const EStatisticType stype;
std::vector<index_t> inds;
void add_term(float64_t, index_t, index_t);
struct terms_t
{
float64_t term[3];
float64_t diag[3];
};
terms_t terms;
};

}
Expand Down
68 changes: 47 additions & 21 deletions tests/unit/statistical_testing/internals/WithinBlockPermutation_unittest.cc
View file
Expand Up @@ -51,16 +51,25 @@ using namespace Eigen;
TEST(WithinBlockPermutation, biased_full)
{
const index_t dim=2;
const index_t n=6;
const index_t m=4;
const index_t n=13;
const index_t m=7;

using operation=std::function<float64_t(SGMatrix<float64_t>)>;

SGMatrix<float64_t> data(dim, n+m);
std::iota(data.matrix, data.matrix+dim*(n+m), 1);
SGMatrix<float64_t> data_p(dim, n);
std::iota(data_p.matrix, data_p.matrix+dim*n, 1);
std::for_each(data_p.matrix, data_p.matrix+dim*n, [&n](float64_t& val) { val/=n; });

auto feats=new CDenseFeatures<float64_t>(data);
SGMatrix<float64_t> data_q(dim, m);
std::iota(data_q.matrix, data_q.matrix+dim*m, n+1);
std::for_each(data_q.matrix, data_q.matrix+dim*m, [&m](float64_t& val) { val/=2*m; });

auto feats_p=new CDenseFeatures<float64_t>(data_p);
auto feats_q=new CDenseFeatures<float64_t>(data_q);
auto feats=feats_p->create_merged_copy(feats_q);
SG_REF(feats);
SG_UNREF(feats_p);
SG_UNREF(feats_q);

auto kernel=some<CGaussianKernel>();
kernel->set_width(2.0);
Expand All @@ -69,7 +78,7 @@ TEST(WithinBlockPermutation, biased_full)
auto mat=kernel->get_kernel_matrix();

// compute using within-block-permutation functor
operation compute=shogun::internal::mmd::WithinBlockPermutation(n, EStatisticType::BIASED_FULL);
operation compute=shogun::internal::mmd::WithinBlockPermutation(n, m, EStatisticType::BIASED_FULL);
sg_rand->set_seed(12345);
auto result_1=compute(mat);

Expand Down Expand Up @@ -101,22 +110,31 @@ TEST(WithinBlockPermutation, biased_full)
EXPECT_NEAR(result_1, result_2, 1E-15);
EXPECT_NEAR(result_1, result_3, 1E-15);

kernel->remove_lhs_and_rhs();
SG_UNREF(feats);
}

TEST(WithinBlockPermutation, unbiased_full)
{
const index_t dim=2;
const index_t n=3;
const index_t m=4;
const index_t n=13;
const index_t m=7;

using operation=std::function<float64_t(SGMatrix<float64_t>)>;

SGMatrix<float64_t> data(dim, n+m);
std::iota(data.matrix, data.matrix+dim*(n+m), 1);
SGMatrix<float64_t> data_p(dim, n);
std::iota(data_p.matrix, data_p.matrix+dim*n, 1);
std::for_each(data_p.matrix, data_p.matrix+dim*n, [&n](float64_t& val) { val/=n; });

SGMatrix<float64_t> data_q(dim, m);
std::iota(data_q.matrix, data_q.matrix+dim*m, n+1);
std::for_each(data_q.matrix, data_q.matrix+dim*m, [&m](float64_t& val) { val/=2*m; });

auto feats=new CDenseFeatures<float64_t>(data);
auto feats_p=new CDenseFeatures<float64_t>(data_p);
auto feats_q=new CDenseFeatures<float64_t>(data_q);
auto feats=feats_p->create_merged_copy(feats_q);
SG_REF(feats);
SG_UNREF(feats_p);
SG_UNREF(feats_q);

auto kernel=some<CGaussianKernel>();
kernel->set_width(2.0);
Expand All @@ -125,7 +143,7 @@ TEST(WithinBlockPermutation, unbiased_full)
auto mat=kernel->get_kernel_matrix();

// compute using within-block-permutation functor
operation compute=shogun::internal::mmd::WithinBlockPermutation(n, EStatisticType::UNBIASED_FULL);
operation compute=shogun::internal::mmd::WithinBlockPermutation(n, m, EStatisticType::UNBIASED_FULL);
sg_rand->set_seed(12345);
auto result_1=compute(mat);

Expand Down Expand Up @@ -157,22 +175,30 @@ TEST(WithinBlockPermutation, unbiased_full)
EXPECT_NEAR(result_1, result_2, 1E-15);
EXPECT_NEAR(result_1, result_3, 1E-15);

kernel->remove_lhs_and_rhs();
SG_UNREF(feats);
}

TEST(WithinBlockPermutation, unbiased_incomplete)
{
const index_t dim=2;
const index_t n=5;
const index_t m=5;
const index_t n=10;

using operation=std::function<float64_t(SGMatrix<float64_t>)>;

SGMatrix<float64_t> data(dim, n+m);
std::iota(data.matrix, data.matrix+dim*(n+m), 1);
SGMatrix<float64_t> data_p(dim, n);
std::iota(data_p.matrix, data_p.matrix+dim*n, 1);
std::for_each(data_p.matrix, data_p.matrix+dim*n, [&n](float64_t& val) { val/=n; });

SGMatrix<float64_t> data_q(dim, n);
std::iota(data_q.matrix, data_q.matrix+dim*n, n+1);
std::for_each(data_q.matrix, data_q.matrix+dim*n, [&n](float64_t& val) { val/=2*n; });

auto feats=new CDenseFeatures<float64_t>(data);
auto feats_p=new CDenseFeatures<float64_t>(data_p);
auto feats_q=new CDenseFeatures<float64_t>(data_q);
auto feats=feats_p->create_merged_copy(feats_q);
SG_REF(feats);
SG_UNREF(feats_p);
SG_UNREF(feats_q);

auto kernel=some<CGaussianKernel>();
kernel->set_width(2.0);
Expand All @@ -181,7 +207,7 @@ TEST(WithinBlockPermutation, unbiased_incomplete)
auto mat=kernel->get_kernel_matrix();

// compute using within-block-permutation functor
operation compute=shogun::internal::mmd::WithinBlockPermutation(n, EStatisticType::UNBIASED_INCOMPLETE);
operation compute=shogun::internal::mmd::WithinBlockPermutation(n, n, EStatisticType::UNBIASED_INCOMPLETE);
sg_rand->set_seed(12345);
auto result_1=compute(mat);

Expand Down Expand Up @@ -213,5 +239,5 @@ TEST(WithinBlockPermutation, unbiased_incomplete)
EXPECT_NEAR(result_1, result_2, 1E-15);
EXPECT_NEAR(result_1, result_3, 1E-15);

kernel->remove_lhs_and_rhs();
SG_UNREF(feats);
}

0 comments on commit 84305f0

Please sign in to comment.