Commit 866ecb5f authored by BRAMAS Berenger's avatar BRAMAS Berenger

add a test to visualize the partitions

parent b4603f17
......@@ -62,85 +62,118 @@ public:
nbValuesUnderThreshold += (inDistMat[idxColReal*dim+ idxRowReal] < treshold && idxRowReal != idxColReal? 1 : 0);
}
}
FAssertLF(nbValuesUnderThreshold != 0);
// Base value for all array indexings
const SCOTCH_Num baseval = 0;
// Number of vertices in graph
const SCOTCH_Num vertnbr = sizeInterval;
// Number of arcs in graph. Since edges are represented by both of their ends,
// the number of edge data in the graph is twice the number of graph edges.
const SCOTCH_Num edgenbr = nbValuesUnderThreshold;// it is already symmetric
// Array of start indices in edgetab of vertex adjacency sub-arrays
SCOTCH_Num* verttab = new SCOTCH_Num[vertnbr+1];
SCOTCH_Num* vendtab = &verttab[1];
// edgetab[verttab[i]] to edgetab[vendtab[i] − 1]
SCOTCH_Num* edgetab = new SCOTCH_Num[edgenbr];
// Optional array, of size vertnbr, holding the integer load associated with every vertex.
SCOTCH_Num* velotab = nullptr;
// Optional array, of a size equal at least to (max i (vendtab[i]) − baseval), holding the integer load associated with every arc.
SCOTCH_Num* edlotab = nullptr;
SCOTCH_Num* vlbltab = nullptr;
verttab[0] = 0;
for(int idxCol = 0 ; idxCol < sizeInterval ; ++idxCol){
const int idxColReal = permutations[idxCol+currentCluster.first];
verttab[idxCol+1] = verttab[idxCol];
for(int idxRow = 0 ; idxRow < sizeInterval ; ++idxRow){
const int idxRowReal = permutations[idxRow+currentCluster.first];
if(inDistMat[idxColReal*dim+ idxRowReal] < treshold
&& idxColReal != idxRowReal){
edgetab[verttab[idxCol+1]] = idxRow;
verttab[idxCol+1] += 1;
if(nbValuesUnderThreshold != 0){
// Base value for all array indexings
const SCOTCH_Num baseval = 0;
// Number of vertices in graph
const SCOTCH_Num vertnbr = sizeInterval;
// Number of arcs in graph. Since edges are represented by both of their ends,
// the number of edge data in the graph is twice the number of graph edges.
const SCOTCH_Num edgenbr = nbValuesUnderThreshold;// it is already symmetric
// Array of start indices in edgetab of vertex adjacency sub-arrays
SCOTCH_Num* verttab = new SCOTCH_Num[vertnbr+1];
SCOTCH_Num* vendtab = &verttab[1];
// edgetab[verttab[i]] to edgetab[vendtab[i] − 1]
SCOTCH_Num* edgetab = new SCOTCH_Num[edgenbr];
// Optional array, of size vertnbr, holding the integer load associated with every vertex.
SCOTCH_Num* velotab = nullptr;
// Optional array, of a size equal at least to (max i (vendtab[i]) − baseval), holding the integer load associated with every arc.
SCOTCH_Num* edlotab = nullptr;
SCOTCH_Num* vlbltab = nullptr;
verttab[0] = 0;
for(int idxCol = 0 ; idxCol < sizeInterval ; ++idxCol){
const int idxColReal = permutations[idxCol+currentCluster.first];
verttab[idxCol+1] = verttab[idxCol];
for(int idxRow = 0 ; idxRow < sizeInterval ; ++idxRow){
const int idxRowReal = permutations[idxRow+currentCluster.first];
if(inDistMat[idxColReal*dim+ idxRowReal] < treshold
&& idxColReal != idxRowReal){
edgetab[verttab[idxCol+1]] = idxRow;
verttab[idxCol+1] += 1;
}
}
}
}
FAssertLF(verttab[vertnbr] == edgenbr);
SCOTCH_Graph grafdat;
FAssertLF(SCOTCH_graphInit(&grafdat) == 0);
FAssertLF(SCOTCH_graphBuild(&grafdat, baseval, vertnbr, verttab, vendtab,
velotab, vlbltab, edgenbr, edgetab, edlotab) == 0);
FAssertLF(SCOTCH_graphCheck(&grafdat) == 0);
SCOTCH_Strat straptr;
FAssertLF(SCOTCH_stratInit(&straptr) == 0);
//const SCOTCH_Num pwgtmax = std::numeric_limits<SCOTCH_Num>::max(); //maximum cluster vertex weight
//const double densmin = 0; // the minimum edge density
//const double bbalval = 0.2;// bipartition imbalance ratio
//SCOTCH_stratGraphClusterBuild (&straptr, SCOTCH_STRATDEFAULT, pwgtmax, densmin, bbalval);
const int partnbr = 2;
SCOTCH_Num* parttab = new SCOTCH_Num[vertnbr];
FAssertLF(SCOTCH_graphPart(&grafdat, partnbr, &straptr, parttab) == 0);
{
int idxLeft = 0;
int idxRight = sizeInterval-1;
while(idxLeft <= idxRight){
if(parttab[idxLeft] == 0){
idxLeft += 1;
FAssertLF(verttab[vertnbr] == edgenbr);
SCOTCH_Graph grafdat;
FAssertLF(SCOTCH_graphInit(&grafdat) == 0);
FAssertLF(SCOTCH_graphBuild(&grafdat, baseval, vertnbr, verttab, vendtab,
velotab, vlbltab, edgenbr, edgetab, edlotab) == 0);
FAssertLF(SCOTCH_graphCheck(&grafdat) == 0);
SCOTCH_Strat straptr;
FAssertLF(SCOTCH_stratInit(&straptr) == 0);
//const SCOTCH_Num pwgtmax = std::numeric_limits<SCOTCH_Num>::max(); //maximum cluster vertex weight
//const double densmin = 0; // the minimum edge density
//const double bbalval = 0.2;// bipartition imbalance ratio
//SCOTCH_stratGraphClusterBuild (&straptr, SCOTCH_STRATDEFAULT, pwgtmax, densmin, bbalval);
const int partnbr = 2;
SCOTCH_Num* parttab = new SCOTCH_Num[vertnbr];
FAssertLF(SCOTCH_graphPart(&grafdat, partnbr, &straptr, parttab) == 0);
{
int idxLeft = 0;
int idxRight = sizeInterval-1;
while(idxLeft <= idxRight){
if(parttab[idxLeft] == 0){
idxLeft += 1;
}
else if(parttab[idxRight] == 1){
idxRight -= 1;
}
else{
const int unk = parttab[idxLeft];
parttab[idxLeft] = parttab[idxRight];
parttab[idxRight] = unk;
const int unkPerm = permutations[idxLeft+currentCluster.first];
permutations[idxLeft+currentCluster.first] = permutations[idxRight+currentCluster.first];
permutations[idxRight+currentCluster.first] = unkPerm;
idxLeft += 1;
idxRight -= 1;
}
}
else if(parttab[idxRight] == 1){
idxRight -= 1;
// idxLeft is on the first 1
if(idxLeft == 1){
gclusters[idxCluster].first = permutations[currentCluster.first];
countToUnknowns += 1;
FAssertLF(countToUnknowns <= dim);
}
else{
const int unk = parttab[idxLeft];
parttab[idxLeft] = parttab[idxRight];
parttab[idxRight] = unk;
const int unkPerm = permutations[idxLeft+currentCluster.first];
permutations[idxLeft+currentCluster.first] = permutations[idxRight+currentCluster.first];
permutations[idxRight+currentCluster.first] = unkPerm;
idxLeft += 1;
idxRight -= 1;
else if(idxLeft > 1){
FAssertLF(idxChildCluster >= 0);
gclusters[idxCluster].first = (-idxChildCluster)-1;
idxChildCluster -= 1;
intervals.push(std::pair<int,int>(currentCluster.first, currentCluster.first+idxLeft));
}
if(idxRight == sizeInterval-2){
gclusters[idxCluster].second = permutations[sizeInterval-1+currentCluster.first];
countToUnknowns += 1;
FAssertLF(countToUnknowns <= dim);
}
else if(idxRight < sizeInterval-2){
FAssertLF(idxChildCluster >= 0);
gclusters[idxCluster].second = (-idxChildCluster)-1;
idxChildCluster -= 1;
intervals.push(std::pair<int,int>(currentCluster.first+idxLeft, currentCluster.first+sizeInterval));
}
}
SCOTCH_stratExit(&straptr);
SCOTCH_graphExit(&grafdat);
delete[] parttab;
delete[] verttab;
delete[] edgetab;
}
else{
int idxLeft = sizeInterval/2;
int idxRight = idxLeft-1;
// idxLeft is on the first 1
if(idxLeft == 1){
gclusters[idxCluster].first = permutations[currentCluster.first];
......@@ -164,16 +197,8 @@ public:
idxChildCluster -= 1;
intervals.push(std::pair<int,int>(currentCluster.first+idxLeft, currentCluster.first+sizeInterval));
}
fflush(stdout);
}
SCOTCH_stratExit(&straptr);
SCOTCH_graphExit(&grafdat);
delete[] parttab;
delete[] verttab;
delete[] edgetab;
idxCluster -= 1;
}
FAssertLF(idxCluster == -1);
......
// @SCALFMM_PRIVATE
#include "../Src/Clustering/FCCLTreeCluster.hpp"
#include "../Src/Clustering/FGraphThreshold.hpp"
#include "../Src/Clustering/FCCLTreeCluster.hpp"
#include "../Src/Utils/FMatrixIO.hpp"
#include "../Src/Containers/FStaticDiagonalBisection.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 "Files/FGenerateDistribution.hpp"
#include <memory>
// 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);
//FCCLTreeCluster<FReal> partitioner(dim, distances.get(), CCL::CCL_TM_MAXIMUM);
FClusterTree<double> tclusters;
partitioner.fillClusterTree(&tclusters);
tclusters.checkData();
tclusters.saveToXml(outputdir, "scotch.xml");
tclusters.saveToDot(outputdir, "gscotch.dot");
std::unique_ptr<int[]> permutations(new int[dim]);
tclusters.fillPermutations(permutations.get());
for(int idxLevel = 2 ; idxLevel <= height ; ++idxLevel){
const int nbPartitions = FMath::pow2(idxLevel-1);
std::unique_ptr<int[]> partitions(new int[nbPartitions]);
tclusters.getPartitions(idxLevel, nbPartitions, partitions.get());
{
char coordfilename[1024];
sprintf(coordfilename, "%s/coord-%d.csv", outputdir, idxLevel);
FILE* fcoord = fopen(coordfilename, "w");
int offsetParticles = 0;
for(int idxPartition = 0 ; idxPartition < nbPartitions ; ++idxPartition){
for(int idxPart = 0 ; idxPart < partitions[idxPartition] ; ++idxPart){
const int idxUnk = permutations[idxPart+offsetParticles];
fprintf(fcoord, "%e,%e,%e,%d\n",
spherePoints[idxUnk*4 + 0],
spherePoints[idxUnk*4 + 1],
spherePoints[idxUnk*4 + 2],
idxPartition);
}
offsetParticles += partitions[idxPartition];
}
FAssertLF(offsetParticles == dim);
fclose(fcoord);
}
{
typedef FDenseBlock<FReal> LeafClass;
typedef FDenseBlock<FReal> CellClass;
typedef FStaticDiagonalBisection<FReal, LeafClass, CellClass> GridClass;
GridClass bissection(dim, idxLevel, partitions.get(), nbPartitions);
char svgfilename[1024];
sprintf(svgfilename, "%s/scotch-%d.svg", outputdir, idxLevel);
FSvgRect output(svgfilename, dim);
bissection.forAllBlocksDescriptor([&](const FBlockDescriptor& info){
output.addRectWithLegend(info.col, info.row, info.nbCols, info.nbRows, info.level);
});
}
}
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