Commit 7931c344 authored by BRAMAS Berenger's avatar BRAMAS Berenger

test hmat clustering and rank

parent 00bfe742
......@@ -222,7 +222,9 @@ public:
}
int getRank() const{
return rank;
}
};
#endif // FSVDBLOCK_HPP
......
......@@ -118,7 +118,7 @@ public:
free(croot);
}
void fillClusterTree(FClusterTree<FReal>* ctree){
void fillClusterTree(FClusterTree<FReal>* ctree) const {
int* permsOrigToNew = new int[dim];
int* permsNewToOrig = new int[dim];
......
......@@ -29,6 +29,19 @@ protected:
std::pair<int, int>* gclusters;
public:
static FReal GetDefaultRadius(const int inDim, const FReal inDistMat[]){
FReal avg = 0;
FReal min = std::numeric_limits<FReal>::max();
const FReal dimSquare = FReal(inDim*inDim);
for(int idxCol = 0 ; idxCol < inDim ; ++idxCol){
for(int idxRow = 0 ; idxRow < inDim ; ++idxRow){
avg += (inDistMat[idxCol*inDim + idxRow]/dimSquare);
min = (min < inDistMat[idxCol*inDim + idxRow]? min:inDistMat[idxCol*inDim + idxRow]);
}
}
return (avg+min) / FReal(2);
}
FGraphThreshold(const int inDim, const FReal inDistMat[], const FReal inTreshold)
: dim(inDim), treshold(inTreshold),
permutations(new int[dim]),
......@@ -211,7 +224,7 @@ public:
delete[] gclusters;
}
void fillClusterTree(FClusterTree<FReal>* ctree){
void fillClusterTree(FClusterTree<FReal>* ctree) const {
int* permsOrigToNew = new int[dim];
int* permsNewToOrig = new int[dim];
......
......@@ -127,7 +127,7 @@ public:
delete[] gclusters;
}
void fillClusterTree(FClusterTree<FReal>* ctree){
void fillClusterTree(FClusterTree<FReal>* ctree) const {
int* permsOrigToNew = new int[dim];
int* permsNewToOrig = new int[dim];
......
......@@ -12,18 +12,18 @@
#include "../Utils/FMatrixIO.hpp"
template <class FReal>
class FMatDense {
class FMatDensePerm {
protected:
int matDim;
FReal* values;
int* permOrigToNew;
FMatDense(const FMatDense&) = delete;
FMatDense& operator=(const FMatDense&) = delete;
FMatDensePerm(const FMatDensePerm&) = delete;
FMatDensePerm& operator=(const FMatDensePerm&) = delete;
public:
using BlockDescriptor = FDenseBlockPermWrapper<FReal>;
explicit FMatDense(const int inDim, const FReal* inValues = nullptr, const int* inPermOrigToNew = nullptr)
explicit FMatDensePerm(const int inDim, const FReal* inValues = nullptr, const int* inPermOrigToNew = nullptr)
: matDim(inDim), values(nullptr){
values = new FReal[matDim*matDim];
if(inValues){
......@@ -46,7 +46,7 @@ public:
}
}
explicit FMatDense(const char inFilename[], const int* inPermOrigToNew = nullptr)
explicit FMatDensePerm(const char inFilename[], const int* inPermOrigToNew = nullptr)
: matDim(0), values(nullptr){
int readNbRows = 0;
int readNbCols = 0;
......@@ -65,7 +65,7 @@ public:
}
}
~FMatDense(){
~FMatDensePerm(){
delete[] values;
delete[] permOrigToNew;
}
......@@ -96,7 +96,7 @@ public:
FAssertLF(0 < nbRows);
FAssertLF(0 < nbCols);
FAssertLF(rowIdx + nbRows <= matDim);
FAssertLF(colIdx + nbRows <= matDim);
FAssertLF(colIdx + nbCols <= matDim);
return FDenseBlockPermWrapper<FReal>(&values[colIdx*matDim+rowIdx], nbRows, nbCols, matDim, permOrigToNew);
}
};
......
// ===================================================================================
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Berenger Bramas, Matthias Messner
// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// This software is a computer program whose purpose is to compute the FMM.
//
// This software is governed by the CeCILL-C and LGPL licenses and
// abiding by the rules of distribution of free software.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public and CeCILL-C Licenses for more details.
// "http://www.cecill.info".
// "http://www.gnu.org/licenses".
// ===================================================================================
// @SCALFMM_PRIVATE
#include "../Src/Containers/FStaticDiagonalBisection.hpp"
#include "../Src/Viewers/FMatDensePerm.hpp"
#include "../Src/Blocks/FDenseBlock.hpp"
#include "../Src/Blocks/FSVDBlock.hpp"
#include "../Src/Clustering/FMaxDistCut.hpp"
#include "../Src/Clustering/FCCLTreeCluster.hpp"
#include "../Src/Clustering/FGraphThreshold.hpp"
#include "../Src/Utils/FSvgRect.hpp"
#include "Utils/FParameters.hpp"
#include "Utils/FParameterNames.hpp"
#include "Utils/FTic.hpp"
#include <memory>
template <class FReal, class MatrixClass>
void CheckRank(const FClusterTree<double>& ctree, MatrixClass& matrix,
const char outputdir[], const char configName[], const int height){
typedef FDenseBlock<FReal> LeafClass;
typedef FSVDBlock<FReal,7> CellClass;
typedef FStaticDiagonalBisection<FReal, LeafClass, CellClass> GridClass;
const int dim = matrix.getDim();
std::unique_ptr<int[]> permutationOrgToNew(new int[dim]);
for( int idxLevel = 1 ; idxLevel < height ; ++idxLevel){
const int nbPartitions = FMath::pow2(height-1);
std::unique_ptr<int[]> partitions(new int[nbPartitions]);
ctree.getPartitions(height, nbPartitions, partitions.get());
ctree.fillPermutations(permutationOrgToNew.get());
matrix.setPermutOrigToNew(permutationOrgToNew.get());
std::cout << "\tLevel " << idxLevel << " build blocks\n";
FTic timer;
GridClass grid(dim, height, partitions.get(), nbPartitions);
grid.fillBlocks(matrix);
std::cout << "\tdone in " << timer.tacAndElapsed() << "s\n";
{
char svgName[1024];
sprintf(svgName, "%s/%s-%d.svg", outputdir, configName, idxLevel);
std::cout << "\tSave svg to " << svgName << "\n";
FSvgRect output(svgName, dim);
grid.forAllCellBlocks([&](const FBlockDescriptor& info, const CellClass& cell){
output.addRectWithLegend(info.col, info.row, info.nbCols, info.nbRows, info.level, cell.getRank());
});
grid.forAllLeafBlocks([&](const FBlockDescriptor& info, const LeafClass& /*leaf*/){
output.addRectWithLegend(info.col, info.row, info.nbCols, info.nbRows, info.level);
});
}
}
}
int main(int argc, char** argv){
static const FParameterNames SvgOutParam = {
{"-fout", "--out", "-out"} ,
"Svg output directory."
};
FHelpDescribeAndExit(argc, argv, "Test the rank for different clustering.",
FParameterDefinitions::InputFileOne,
FParameterDefinitions::InputFileTwow,
FParameterDefinitions::OctreeHeight,
SvgOutParam);
////////////////////////////////////////////////////////////////////
const char* outputdir = FParameters::getStr(argc, argv, SvgOutParam.options, "/tmp/");
const char* distanceFilename = FParameters::getStr(argc, argv, FParameterDefinitions::InputFileOne.options, "../Addons/HMat/Data/unitCube1000_ONE_OVER_R.bin");
const char* matrixFilename = FParameters::getStr(argc, argv, FParameterDefinitions::InputFileTwow.options, "../Addons/HMat/Data/unitCube1000_ONE_OVER_R.bin");
const int height = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeHeight.options, 4);
std::cout << "Check until height = " << height << "\n";
typedef double FReal;
typedef FMatDensePerm<FReal> MatrixClass;
////////////////////////////////////////////////////////////////////
std::cout << "Load distances file " << distanceFilename << "\n";
std::unique_ptr<FReal[]> distanceValues;
int distanceValuesDim = 0;
{
int readNbRows = 0;
int readNbCols = 0;
FReal* distanceValuesPtr = nullptr;
FAssertLF(FMatrixIO::read(distanceFilename, &distanceValuesPtr, &readNbRows, &readNbCols));
FAssertLF(readNbRows == readNbCols);
distanceValuesDim = readNbRows;
distanceValues.reset(distanceValuesPtr);
}
////////////////////////////////////////////////////////////////////
std::cout << "Load matrix file " << matrixFilename << "\n";
MatrixClass matrix(matrixFilename);
const int matrixDim = matrix.getDim();
FAssertLF(distanceValuesDim == matrixDim);
std::cout << "Matrices dim = " << matrixDim << "\n";
////////////////////////////////////////////////////////////////////
{
std::cout << "Test FMaxDistCut\n";
FMaxDistCut<FReal> partitioner(matrixDim, distanceValues.get());
FClusterTree<double> tclusters;
partitioner.fillClusterTree(&tclusters);
tclusters.checkData();
CheckRank<FReal, MatrixClass>(tclusters, matrix, outputdir, "FMaxDistCut", height);
}
////////////////////////////////////////////////////////////////////
{
std::cout << "Test FGraphThreshold\n";
FGraphThreshold<FReal> partitioner(matrixDim, distanceValues.get(), FGraphThreshold<FReal>::GetDefaultRadius(matrixDim, distanceValues.get()));
FClusterTree<double> tclusters;
partitioner.fillClusterTree(&tclusters);
tclusters.checkData();
CheckRank<FReal, MatrixClass>(tclusters, matrix, outputdir, "FGraphThreshold", height);
}
////////////////////////////////////////////////////////////////////
{
std::cout << "Test FCCLTreeCluster\n";
FCCLTreeCluster<FReal> partitioner(matrixDim, distanceValues.get(), CCL::CCL_TM_MAXIMUM);
FClusterTree<double> tclusters;
partitioner.fillClusterTree(&tclusters);
tclusters.checkData();
CheckRank<FReal, MatrixClass>(tclusters, matrix, outputdir, "FCCLTreeCluster", height);
}
return 0;
}
......@@ -42,7 +42,7 @@ int main(int argc, char** argv){
FAssertLF(readNbRows == readNbCols);
const int dim = readNbRows;
FGraphThreshold<double> partitioner(dim, distances, 1.0);
FGraphThreshold<double> partitioner(dim, distances, FGraphThreshold<double>::GetDefaultRadius(dim, distances));
FClusterTree<double> tclusters;
partitioner.fillClusterTree(&tclusters);
......
......@@ -21,60 +21,23 @@
// FUSE_SCOTCH
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."
};
FHelpDescribeAndExit(argc, argv,"Test the bisection.",SvgOutParam,DimParam,FParameterDefinitions::OctreeHeight,
FParameterDefinitions::NbParticles);
const int nbParticles = (FParameters::getValue(argc, argv, FParameterDefinitions::NbParticles.options, 3000) & ~1);
const int height = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeHeight.options, 4);
const char* outputdir = FParameters::getStr(argc, argv, SvgOutParam.options, "/tmp/");
typedef double FReal;
const FReal radius = 0.5;
const FReal distanceBetweenSpheres = 0.4;
const int nbPointsSphere = nbParticles/2;
std::unique_ptr<FReal[]> spherePoints(new FReal[nbParticles*4]);
unifRandonPointsOnSphere(nbPointsSphere, radius, spherePoints.get()) ;
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;
}
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));
}
}
//FGraphThreshold<FReal> partitioner(dim, distances.get(), radius/4);
//FCCLTreeCluster<FReal> partitioner(dim, distances.get(), CCL::CCL_TM_MAXIMUM);
FMaxDistCut<FReal> partitioner(dim, distances.get());
template <class FReal, class PartitionerClass>
void TestPartitions(const PartitionerClass& partitioner, const int height,
const char outputdir[], const int dim,
const FReal coords[], const char config[]){
FClusterTree<double> tclusters;
partitioner.fillClusterTree(&tclusters);
tclusters.checkData();
tclusters.saveToXml(outputdir, "scotch.xml");
tclusters.saveToDot(outputdir, "gscotch.dot");
{
char filenameBuffer[1024];
sprintf(filenameBuffer, "%s/%s.xml", outputdir, config);
std::cout << "\t save xml to " << filenameBuffer << "\n";
tclusters.saveToXml(filenameBuffer);
sprintf(filenameBuffer, "%s/%s.dot", outputdir, config);
std::cout << "\t save cdot to " << filenameBuffer << "\n";
tclusters.saveToDot(filenameBuffer);
}
std::unique_ptr<int[]> permutations(new int[dim]);
std::unique_ptr<int[]> invpermutations(new int[dim]);
......@@ -87,7 +50,8 @@ int main(int argc, char** argv){
{
char coordfilename[1024];
sprintf(coordfilename, "%s/coord-%d.csv", outputdir, idxLevel);
sprintf(coordfilename, "%s/%s-coord-%d.csv", outputdir, config, idxLevel);
std::cout << "\t save coord for level " << idxLevel << " to " << coordfilename << "\n";
FILE* fcoord = fopen(coordfilename, "w");
int offsetParticles = 0;
......@@ -95,9 +59,9 @@ int main(int argc, char** argv){
for(int idxPart = 0 ; idxPart < partitions[idxPartition] ; ++idxPart){
const int idxUnk = invpermutations[idxPart+offsetParticles];
fprintf(fcoord, "%e,%e,%e,%d\n",
spherePoints[idxUnk*4 + 0],
spherePoints[idxUnk*4 + 1],
spherePoints[idxUnk*4 + 2],
coords[idxUnk*4 + 0],
coords[idxUnk*4 + 1],
coords[idxUnk*4 + 2],
idxPartition);
}
offsetParticles += partitions[idxPartition];
......@@ -114,7 +78,8 @@ int main(int argc, char** argv){
GridClass bissection(dim, idxLevel, partitions.get(), nbPartitions);
char svgfilename[1024];
sprintf(svgfilename, "%s/scotch-%d.svg", outputdir, idxLevel);
sprintf(svgfilename, "%s/%s-%d.svg", outputdir, config, idxLevel);
std::cout << "\t save svg for level " << idxLevel << " to " << svgfilename << "\n";
FSvgRect output(svgfilename, dim);
bissection.forAllBlocksDescriptor([&](const FBlockDescriptor& info){
......@@ -122,7 +87,84 @@ int main(int argc, char** argv){
});
}
}
}
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."
};
FHelpDescribeAndExit(argc, argv,"Test the bisection.",SvgOutParam,DimParam,FParameterDefinitions::OctreeHeight,
FParameterDefinitions::NbParticles);
const int nbParticles = (FParameters::getValue(argc, argv, FParameterDefinitions::NbParticles.options, 3000) & ~1);
const int height = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeHeight.options, 6);
const char* outputdir = FParameters::getStr(argc, argv, SvgOutParam.options, "/tmp/");
typedef double FReal;
/////////////////////////////////////////////////////////////////
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[]> spherePoints(new FReal[nbParticles*4]);
unifRandonPointsOnSphere(nbPointsSphere, radius, spherePoints.get()) ;
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::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));
}
}
/////////////////////////////////////////////////////////////////
{
std::cout << "Test FGraphThreshold\n";
FGraphThreshold<FReal> partitioner(dim, distances.get(), radius/4);
TestPartitions<FReal, FGraphThreshold<FReal>>(partitioner,height, outputdir,
dim, spherePoints.get(), "FGraphThreshold");
}
/////////////////////////////////////////////////////////////////
{
std::cout << "Test FCCLTreeCluster\n";
FCCLTreeCluster<FReal> partitioner(dim, distances.get(), CCL::CCL_TM_MAXIMUM);
TestPartitions<FReal, FCCLTreeCluster<FReal>>(partitioner,height, outputdir,
dim, spherePoints.get(), "FCCLTreeCluster");
}
/////////////////////////////////////////////////////////////////
{
std::cout << "Test FMaxDistCut\n";
FMaxDistCut<FReal> partitioner(dim, distances.get());
TestPartitions<FReal, FMaxDistCut<FReal>>(partitioner,height, outputdir,
dim, spherePoints.get(), "FMaxDistCut");
}
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