Commit b5fc3a80 authored by BRAMAS Berenger's avatar BRAMAS Berenger
parents e89b3547 75174454
// @SCALFMM_PRIVATE
#ifndef FCCLKCLUSTER_HPP
#define FCCLKCLUSTER_HPP
#include "./Utils/FGlobal.hpp"
#include "./Utils/FAssert.hpp"
#include "./Utils/FMath.hpp"
#include "../Utils/FHUtils.hpp"
#include <stack>
#include <vector>
#include <functional>
#include <queue>
#include <limits>
#include <algorithm>
extern "C" {
#include <cluster.h>
}
namespace CCL {
enum ClusterCenterMethod {
CCL_CCM_ARITHMETIC_MEAN,
CCL_CCM_MEDIAN
};
inline char ClusterCenterMethodToChar(const ClusterCenterMethod method){
switch (method) {
case CCL_CCM_ARITHMETIC_MEAN:
return 'a';
break;
case CCL_CCM_MEDIAN:
return 'm';
break;
default:
break;
}
return '?';
}
enum Distance {
CCL_DIST_MEAN,
CCL_DIST_MEDIAN,
CCL_DIST_SHORTEST,
CCL_DIST_LONGEST,
CCL_DIST_AVG
};
inline char DistanceToChar(const Distance method){
switch (method) {
case CCL_DIST_MEAN:
return 'a';
break;
case CCL_DIST_MEDIAN:
return 'm';
break;
case CCL_DIST_SHORTEST:
return 's';
break;
case CCL_DIST_LONGEST:
return 'x';
break;
case CCL_DIST_AVG:
return 'v';
break;
default:
break;
}
return '?';
}
}
template <class FReal>
class FCCLKCluster {
protected:
const int nbPartitions;
const int dim;
CCL::ClusterCenterMethod method;
int* partitions;
public:
FCCLKCluster(const int inNbPartitions, const int inDim, const FReal inDistMat[], const CCL::ClusterCenterMethod inMethod)
: nbPartitions(inNbPartitions), dim(inDim), method(inMethod), partitions(nullptr) {
double** distMatPtrs = new double*[dim];
// Build mask, everyone is here
for(int idxRow = 0 ; idxRow < dim ; ++idxRow){
distMatPtrs[idxRow] = new double[idxRow+1];
for(int idxCol = 0 ; idxCol <= idxRow ; ++idxCol){
distMatPtrs[idxRow][idxCol] = double(inDistMat[idxCol*dim + idxRow]);
}
}
// allocate partitions
partitions = new int[dim];
// Number of call to EM algorithm
const int nbPass = 1;
// Errors
double* error = new double[dim];
// Nb of times the optimal clustering was found
int* ifound = new int[dim];
kmedoids (nbPartitions, dim, distMatPtrs, nbPass, partitions, error, ifound);
for(int idxRow = 0 ; idxRow < dim ; ++idxRow){
delete[] distMatPtrs[idxRow];
}
delete[] distMatPtrs;
}
~FCCLKCluster(){
delete[] partitions;
}
int getPartitions(const int inNbPartitions, int inNbIdxInPartitions[]) const{
/// Map partitions to 0.. nbPartitions
FAssertLF(inNbPartitions == nbPartitions);
// Copy partitions
int* sortedPartitions = new int[dim];
for(int idx = 0 ; idx < dim ; ++idx){
sortedPartitions[idx]=partitions[idx];
}
// sort partitions
std::sort(sortedPartitions,sortedPartitions+dim);
// Map partitions to 0..nbPartitions
int counterPartition=0;
int* mapPartitions = new int[inNbPartitions];
int currentPartition=sortedPartitions[0];
for(int idx = 0 ; idx < dim ; ++idx){
mapPartitions[counterPartition]=currentPartition;
if(sortedPartitions[idx+1]!=currentPartition){
currentPartition=sortedPartitions[idx+1];
++counterPartition;
}
}
FAssertLF(counterPartition == inNbPartitions);
/// Count particles in each partition
int totalGiven = 0;
for(int idxPartition = 0 ; idxPartition < inNbPartitions ; ++idxPartition){
inNbIdxInPartitions[idxPartition]=0;
for(int idx = 0 ; idx < dim ; ++idx){
if(partitions[idx]==mapPartitions[idxPartition])
inNbIdxInPartitions[idxPartition]+=1;
}
totalGiven +=inNbIdxInPartitions[idxPartition];
}
FAssertLF(totalGiven == dim);
return 0; // no empty partition in kclusters/kmedoids algorithms
}
};
#endif // FCCLKCLUSTER_HPP
// @SCALFMM_PRIVATE
#include "../Src/Clustering/FCCLKCluster.hpp"
#include "../Src/Utils/FMatrixIO.hpp"
#include "../Src/Containers/FPartitionsMapping.hpp"
#include "../Src/Utils/FSvgRect.hpp"
#include "../Src/Viewers/FDenseBlockWrapper.hpp"
#include "../Src/Blocks/FDenseBlock.hpp"
#include "Utils/FParameters.hpp"
#include "Utils/FParameterNames.hpp"
#include <memory>
int main(int argc, char** argv){
static const FParameterNames SvgOutParam = {
{"-fout", "--out", "-out"} ,
"Svg output directory."
};
static const FParameterNames DimParam = {
{"-N", "-nb", "-dim"} ,
"Dim of the matrix."
};
static const FParameterNames NbPartitionsParam = {
{"-P", "-nbPartitions", "-nbClusters"} ,
"Dim of the matrix."
};
FHelpDescribeAndExit(argc, argv,"Test the k-medoid algorithm from CCL (using 2 methods for cluster center computation).",SvgOutParam,DimParam,NbPartitionsParam,FParameterDefinitions::InputFile);
const char* filename = FParameters::getStr(argc, argv, FParameterDefinitions::InputFile.options, "../Addons/HMat/Data/unitCube1000.bin");
typedef float FReal;
const char* outputdir = FParameters::getStr(argc, argv, SvgOutParam.options, "/tmp/");
int readNbRows = 0;
int readNbCols = 0;
// Read distances
FReal* distances = nullptr;
FAssertLF(FMatrixIO::read(filename, &distances, &readNbRows, &readNbCols));
FAssertLF(readNbRows == readNbCols);
const int dim = readNbRows;
// Define number of partitions
const int nbPartitions = FParameters::getValue(argc, argv, NbPartitionsParam.options, 4);
// Test Cluster Center Method ARITHMETIC_MEAN
{
FCCLKCluster<FReal> partitioner(nbPartitions, dim, distances, CCL::CCL_CCM_ARITHMETIC_MEAN);
std::unique_ptr<int[]> partitions(new int[nbPartitions]);
partitioner.getPartitions(nbPartitions,partitions.get());
{
typedef FDenseBlock<FReal> CellClass;
typedef FPartitionsMapping<FReal, CellClass> GridClass;
GridClass bissection(dim, partitions.get(), nbPartitions);
char svgName[1024];
sprintf(svgName, "%s/%s-%d.svg", outputdir, "CCL_CCM_ARITHMETIC_MEAN", nbPartitions);
FSvgRect output(svgName, dim);
std::cout << "\tSave svg to " << svgName << "\n";
bissection.forAllBlocksDescriptor([&](const FBlockDescriptor& info){
output.addRectWithLegend(info.col, info.row, info.nbCols, info.nbRows, info.level);
});
}
}
// Test Cluster Center Method MEDIAN
{
FCCLKCluster<FReal> partitioner(nbPartitions, dim, distances, CCL::CCL_CCM_MEDIAN);
std::unique_ptr<int[]> partitions(new int[nbPartitions]);
partitioner.getPartitions(nbPartitions,partitions.get());
{
typedef FDenseBlock<FReal> CellClass;
typedef FPartitionsMapping<FReal, CellClass> GridClass;
GridClass bissection(dim, partitions.get(), nbPartitions);
char svgName[1024];
sprintf(svgName, "%s/%s-%d.svg", outputdir, "CCL_CCM_MEDIAN", nbPartitions);
FSvgRect output(svgName, dim);
std::cout << "\tSave svg to " << svgName << "\n";
bissection.forAllBlocksDescriptor([&](const FBlockDescriptor& info){
output.addRectWithLegend(info.col, info.row, info.nbCols, info.nbRows, info.level);
});
}
}
return 0;
}
......@@ -32,6 +32,7 @@ void TestPartitions(const PartitionerClass& partitioner,
std::unique_ptr<int[]> partitions(new int[nbPartitions]);
partitioner.getPartitions(nbPartitions, partitions.get());
if(coords!=nullptr)
{
char coordfilename[1024];
sprintf(coordfilename, "%s/%s-coord-%d.csv", outputdir, config, idFile);
......@@ -81,51 +82,84 @@ int main(int argc, char** argv){
"Dim of the matrix."
};
FHelpDescribeAndExit(argc, argv,"Test the bisection.",SvgOutParam,DimParam,FParameterDefinitions::OctreeHeight,
FHelpDescribeAndExit(argc, argv,"Test the bisection.",SvgOutParam,DimParam,
FParameterDefinitions::InputFileOne,
FParameterDefinitions::InputFileTwow,
FParameterDefinitions::NbParticles);
const int nbParticles = (FParameters::getValue(argc, argv, FParameterDefinitions::NbParticles.options, 3000) & ~1);
const char* outputdir = FParameters::getStr(argc, argv, SvgOutParam.options, "/tmp/");
typedef double FReal;
typedef float FReal;
/////////////////////////////////////////////////////////////////
std::cout << "Create spheres\n";
////////////////////////////////////////////////////////////////////
const FReal radius = 0.5;
const FReal distanceBetweenSpheres = 0.4;
const int nbPointsSphere = nbParticles/2;
const char* distanceFilename = FParameters::getStr(argc, argv, FParameterDefinitions::InputFileOne.options, "../Addons/HMat/Data/unitCube1000.bin");
std::cout << "\tradius = " << radius << "\n";
std::cout << "\tdistanceBetweenSpheres = " << distanceBetweenSpheres << "\n";
std::cout << "\tnbPointsSphere = " << nbPointsSphere << "\n";
////////////////////////////////////////////////////////////////////
std::unique_ptr<FReal[]> spherePoints(new FReal[nbParticles*4]);
unifRandonPointsOnSphere(nbPointsSphere, radius, spherePoints.get()) ;
std::cout << "Load distances file " << distanceFilename << "\n";
for(int idxPoint = 0 ; idxPoint < nbPointsSphere ; ++idxPoint){
spherePoints[(idxPoint+nbPointsSphere)*4 + 0] = spherePoints[idxPoint*4 + 0] + radius*2 + distanceBetweenSpheres;
spherePoints[(idxPoint+nbPointsSphere)*4 + 1] = spherePoints[idxPoint*4 + 1] + radius*2 + distanceBetweenSpheres;
spherePoints[(idxPoint+nbPointsSphere)*4 + 2] = spherePoints[idxPoint*4 + 2] + radius*2 + distanceBetweenSpheres;
std::unique_ptr<FReal[]> particles(nullptr);
std::unique_ptr<FReal[]> distances;
int dim = 0;
{
int readNbRows = 0;
int readNbCols = 0;
FReal* distanceValuesPtr = nullptr;
FAssertLF(FMatrixIO::read(distanceFilename, &distanceValuesPtr, &readNbRows, &readNbCols));
FAssertLF(readNbRows == readNbCols);
dim = readNbRows;
distances.reset(distanceValuesPtr);
}
/////////////////////////////////////////////////////////////////
std::cout << "Compute distance\n";
const int dim = nbParticles;
std::unique_ptr<FReal[]> distances(new FReal[nbParticles*nbParticles]);
for(int idxCol = 0 ; idxCol < nbParticles ; ++idxCol){
for(int idxRow = 0 ; idxRow < nbParticles ; ++idxRow){
const FReal diffx = spherePoints[idxCol*4 + 0] - spherePoints[idxRow*4 + 0];
const FReal diffy = spherePoints[idxCol*4 + 1] - spherePoints[idxRow*4 + 1];
const FReal diffz = spherePoints[idxCol*4 + 2] - spherePoints[idxRow*4 + 2];
distances[idxCol*nbParticles+idxRow] = FMath::Sqrt((diffx*diffx) + (diffy*diffy) + (diffz*diffz));
}
const FSize displaySize=10;
std::cout<<"\nD=["<<std::endl;
for ( FSize i=0; i<displaySize; ++i) {
for ( FSize j=0; j<displaySize; ++j)
std::cout << distances.get()[i*nbParticles+j] << " ";
std::cout<< std::endl;
}
/////////////////////////////////////////////////////////////////
std::cout<<"]"<<std::endl;
// /////////////////////////////////////////////////////////////////
//
// std::cout << "Create spheres\n";
//
// const FReal radius = 0.5;
// const FReal distanceBetweenSpheres = 0.4;
// const int nbPointsSphere = nbParticles/2;
//
// std::cout << "\tradius = " << radius << "\n";
// std::cout << "\tdistanceBetweenSpheres = " << distanceBetweenSpheres << "\n";
// std::cout << "\tnbPointsSphere = " << nbPointsSphere << "\n";
//
// std::unique_ptr<FReal[]> particles(new FReal[nbParticles*4]);
// unifRandonPointsOnSphere(nbPointsSphere, radius, particles.get()) ;
//
// for(int idxPoint = 0 ; idxPoint < nbPointsSphere ; ++idxPoint){
// particles[(idxPoint+nbPointsSphere)*4 + 0] = particles[idxPoint*4 + 0] + radius*2 + distanceBetweenSpheres;
// particles[(idxPoint+nbPointsSphere)*4 + 1] = particles[idxPoint*4 + 1] + radius*2 + distanceBetweenSpheres;
// particles[(idxPoint+nbPointsSphere)*4 + 2] = particles[idxPoint*4 + 2] + radius*2 + distanceBetweenSpheres;
// }
//
// /////////////////////////////////////////////////////////////////
//
// std::cout << "Compute distance\n";
//
// const int dim = nbParticles;
// std::unique_ptr<FReal[]> distances(new FReal[nbParticles*nbParticles]);
// for(int idxCol = 0 ; idxCol < nbParticles ; ++idxCol){
// for(int idxRow = 0 ; idxRow < nbParticles ; ++idxRow){
// const FReal diffx = particles[idxCol*4 + 0] - particles[idxRow*4 + 0];
// const FReal diffy = particles[idxCol*4 + 1] - particles[idxRow*4 + 1];
// const FReal diffz = particles[idxCol*4 + 2] - particles[idxRow*4 + 2];
// distances[idxCol*nbParticles+idxRow] = FMath::Sqrt((diffx*diffx) + (diffy*diffy) + (diffz*diffz));
// }
// }
//
// /////////////////////////////////////////////////////////////////
FReal distMin = std::numeric_limits<FReal>::max();
FReal distMax = std::numeric_limits<FReal>::min();
......@@ -140,14 +174,14 @@ int main(int argc, char** argv){
std::cout << "Dist max " << distMax << "\n";
/////////////////////////////////////////////////////////////////
const FReal stepThresh = (distMax-distMin)/10;
for(FReal thresh = distMin ; thresh <= distMax ; thresh += stepThresh){
const FReal stepThresh = (distMax-distMin)/100;
for(FReal thresh = distMin ; thresh <= distMax/10 ; thresh += stepThresh){
std::cout << "Test FConnexClustering for thresh " << thresh << "\n";
FConnexClustering<FReal> partitioner(dim, distances.get(), thresh);
std::cout << "\tGot " << partitioner.getNbPartitions() << " partitions\n";
TestPartitions<FReal, FConnexClustering<FReal>>(partitioner,outputdir,
int((thresh-distMin)/stepThresh),
dim, spherePoints.get(), "FConnexClustering");
dim, particles.get(), "FConnexClustering");
}
return 0;
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment