The Higher Education and Research forge

Home My Page Projects Code Snippets Project Openings Garamon
Summary Activity SCM

SCM Repository

1 // Copyright (c) 2018 by University Paris-Est Marne-la-Vallee
2 // MetaData.cpp
3 // This file is part of the Garamon Generator.
4 // Authors: Stephane Breuils and Vincent Nozick
5 // Contact: vincent.nozick@u-pem.fr
6 //
7 // Licence MIT
8 // A a copy of the MIT License is given along with this program
11 #include <math.h>
12 #include <set>
14 #include <Eigen/Eigenvalues>
16 #include "MetaData.hpp"
17 #include "ConfigParser.hpp"
20 MetaData::MetaData() : dimension(0), inputMetricDiagonal(false), identityMetric(false), inputMetricPermutationOfDiagonal(false), maxDimPrecomputedProducts(256) {}
22 MetaData::~MetaData() {}
24 void MetaData::display() const {
25     std::cout << "MetaData" << std::endl;
26     std::cout << "dimension         : " << dimension << std::endl;
27     std::cout << "namespace         : " << namespaceName << std::endl;
29     std::cout << "refinement        : " << (useEigenRefinement==true?"true":"false") << std::endl;
30     std::cout << "cleanup           : " << (useNumericalCleanUp==true?"true":"false") << std::endl;
31     std::cout << "maxDim prec func  : " << maxDimPrecomputedProducts << std::endl;
32     std::cout << "maxDim accessors  : " << maxDimBasisAccessor << std::endl;
33     std::cout << "epsilon           : " << epsilon << std::endl;
34     std::cout << "basis vector name : ";
35     for(unsigned int i=0; i<basisVectorName.size(); ++i)
36         std::cout << basisVectorName[i] << " ";
37     std::cout << std::endl;
38     std::cout << "metric            : \n" << metric << std::endl;
39     std::cout << "is initialy diag  : " << (inputMetricDiagonal==true?"true":"false") << std::endl;
40     std::cout << "is identity       : " << (identityMetric==true?"true":"false") << std::endl;
41     std::cout << "is full rank      : " << (fullRankMetric==true?"true":"false") << std::endl;
42     std::cout << "is permutation of diagonal matrix : " << (inputMetricPermutationOfDiagonal==true?"true":"false") << std::endl;
43     std::cout << "diag metric       : " << diagonalMetric.transpose() << std::endl;
44     if(!inputMetricDiagonal){
45         std::cout << "Vector transformation matrix : \n" << transformationMatrix << std::endl;
46         std::cout << "Vector inverse Transformation  : \n" << inverseTransformationMatrix<< std::endl;
47     }
48 }
50 bool MetaData::checkConsistency() const {
52     bool consistencyCheck = true;
54     // maxDimBasisAccessor consistency: at least for vectors
55     if(maxDimBasisAccessor == 0){
56         std::cout << "error: at least vector accessors are required, see 'max dimension basis accessor' in the conf file." << std::endl;
57         consistencyCheck = false;
58     }
60     // dimension consistency
61     if(dimension == 0){
62         std::cout << "error: dimension should not be 0." << std::endl;
63         consistencyCheck = false;
64     }
66     // basis vector name dimension consistency
67     if(basisVectorName.size() != dimension){
68         std::cout << "error: 'basis vector name' size is not consistent with 'dimension'." << std::endl;
69         consistencyCheck = false;
70     }
72     // metric empty
73     bool metricDefined = true;
74     if( (metric.rows()==0) || (metric.cols()==0) ) {
75         metricDefined = false;
76         std::cout << "error: the metric matrix is not defined." << std::endl;
77         consistencyCheck = false;
78     }
80     // metric matrix ?
81     if(metricDefined) {
82         // metric : square matrix ?
83         bool squareMatrix = true;
84         if( metric.rows() != metric.cols() ){
85             std::cout << "error: the metric is not a square matrix." << std::endl;
86             consistencyCheck = false;
87             squareMatrix = false;
88         }
90         // metric : dimension consistency
91         if( ((unsigned int)metric.cols() != dimension) || ((unsigned int)metric.rows() != dimension) ) {
92             std::cout << "error: the metric dimension is not consistent with 'dimension'." << std::endl;
93             consistencyCheck = false;
94         }
96         if(squareMatrix) {
97             // check if metric is symetric
98             bool symetric = true;
99             for(unsigned int i=0; i<(unsigned int)metric.rows(); ++i)
100                 for(unsigned int j=i; j<(unsigned int)metric.cols(); ++j)
101                     if(fabs(metric(i,j) - metric(j,i)) > epsilon)
102                         symetric = false;
103             if (!symetric) {
104                 std::cout << "error: the metric is not a symmetric matrix." << std::endl;
105                 consistencyCheck = false;
106             }
107         }
108     }
110     // check if the name of the vector basis are not ambiguous for high dimensions
111     // i.e. in dimension 15: e12 is for twelve or one-two ?
112     std::set<std::string> basis;
114     // represent a k-vector with a binary number (k-st bit to 1 means the k-st basis is used)
115     for(unsigned int i=1; i<=pow(2,dimension); ++i){
117         std::string kvector;
118         for(unsigned int k=0; k<dimension; ++k)
119             if(i & (1 << k))
120                 kvector = kvector + basisVectorName[k];
122         if(basis.count(kvector) != 0){
123             std::cout << "error in the basis vector name: " << kvector << " is ambiguous." << std::endl;
124             consistencyCheck = false;
125         }else{
126             basis.insert(kvector);
127         }
128     }
130     // namespace name compatible with C++
131     if(isalpha(namespaceName[0]) == 0){
132         std::cout << "error in the namespace name: the first character of '" << namespaceName << "' should be an alphabetic letter for C++ compliance." << std::endl;
133         consistencyCheck = false;
134     }
136     return consistencyCheck;
139 MetaData::MetaData(const std::string &filename):inputMetricDiagonal(false), identityMetric(false) {
141     // open the parser
142     ConfigParser parser(filename);
144     // load all components of the meta data
145     if(!parser.readString("namespace", namespaceName)) {
146         std::cerr << "error: failed to find " << "namespace" << std::endl;
147         exit(EXIT_FAILURE);
148     }
150     if(!parser.readUInt("dimension", dimension)){
151         std::cerr << "error: failed to find " << "dimension" << std::endl;
152         exit(EXIT_FAILURE);
153     }
155     if(!parser.readUInt("max dimension precomputed products", maxDimPrecomputedProducts)){
156         std::cerr << "error: failed to find " << "max dimension precomputed products" << std::endl;
157         exit(EXIT_FAILURE);
158     }
160     if(!parser.readUInt("max dimension basis accessor", maxDimBasisAccessor)){
161         std::cerr << "error: failed to find " << "max dimension basis accessor" << std::endl;
162         exit(EXIT_FAILURE);
163     }
165     if(!parser.readStringList("basis vector name", basisVectorName)){
166         std::cerr << "error: failed to find " << "basis vector name" << std::endl;
167         exit(EXIT_FAILURE);
168     }
170     if(!parser.readMatrix("metric", metric)){
171         std::cerr << "error: failed to find " << "metric" << std::endl;
172         exit(EXIT_FAILURE);
173     }
175     if(!parser.readBool("metric decomposition refinement", useEigenRefinement)){
176         std::cerr << "error: failed to find " << "metric decomposition refinement" << std::endl;
177         exit(EXIT_FAILURE);
178     }
180     if(!parser.readBool("metric decomposition numerical cleanup", useNumericalCleanUp)){
181         std::cerr << "error: failed to find " << "metric decomposition numerical cleanup" << std::endl;
182         exit(EXIT_FAILURE);
183     }
185     if(!parser.readDouble("metric decomposition numerical cleanup espilon", epsilon)){
186         std::cerr << "error: failed to find " << "metric decomposition numerical cleanup espilon" << std::endl;
187         exit(EXIT_FAILURE);
188     }
190     // check the data consistency
191     if(!checkConsistency()){
192         std::cerr << "configuration inconsistent ... abord" << std::endl;
193         exit(EXIT_FAILURE);
194     }
196     // metric diagonalization
197     if(!metricDiagonalization()){
198         std::cerr << "metric diagonalization ... failed" << std::endl;
199         std::cerr << "Please try again without the metric decomposition refinement or without the numerical cleanup." << std::endl;
200         exit(EXIT_FAILURE);
201     }
205 bool MetaData::metricDiagonalization() {
207         // compute the metric rank
208         if (getRank(metric) == dimension)
209                 fullRankMetric = true;
210         else fullRankMetric = false;
212     // check if the metric is already diagonal
213     if(isMatrixDiagonal(metric, epsilon)){
214         diagonalMetric = metric.diagonal();
215         inputMetricDiagonal = true;
216         if(isMatrixIdentity(metric,epsilon))
217             identityMetric = true;
218         return true;
219     }
222     // ckeck if the metric is a permutation of a diagonal matrix (for fast dual)
223     inputMetricPermutationOfDiagonal = isMatrixPermutationOfDiagonal(metric, epsilon);
225     // compute the diagonalization
226     Eigen::MatrixXd diagonalMatrix;
227     eigenDecomposition(metric, transformationMatrix, diagonalMatrix);
228 //    std::cout << "metric\n" << metric << std::endl;
229 //    std::cout << "P\n" << transformationMatrix << std::endl;
230 //    std::cout << "D\n" << diagonalMatrix << std::endl;
231 //    std::cout << "metric ??\n" << transformationMatrix * diagonalMatrix * transformationMatrix.transpose() << std::endl;
233     // invert the transformation matrix
234     inverseTransformationMatrix = Eigen::MatrixXd(transformationMatrix.transpose());
236     // convert floating points to nearest integers (when possible)
237     Eigen::MatrixXd scaleMatrix(metric.rows(),metric.cols());
238     if(useEigenRefinement)
239         scaleMatrix = eigenRefinement(transformationMatrix, diagonalMatrix, inverseTransformationMatrix);
240 //    std::cout << "metric\n" << metric << std::endl;
241 //    std::cout << "P\n" << transformationMatrix << std::endl;
242 //    std::cout << "D\n" << diagonalMatrix << std::endl;
243 //    std::cout << "Pinv\n" << inverseTransformationMatrix << std::endl;
244 //    std::cout << "metric ??\n" << transformationMatrix * diagonalMatrix * inverseTransformationMatrix << std::endl;
247     //numerical clean up
248     if(useNumericalCleanUp) {
249         transformationMatrix = numericalCleanUp(transformationMatrix,epsilon);
250         diagonalMatrix = numericalCleanUp(diagonalMatrix, epsilon);
251         inverseTransformationMatrix = numericalCleanUp(inverseTransformationMatrix,epsilon);
252     }
254     // check if the new metric is identity
255     if(isMatrixIdentity(diagonalMatrix,epsilon))
256         identityMetric = true;
258     // check decomposition
259     if(!checkNumericalCleanUp(metric,transformationMatrix,diagonalMatrix,inverseTransformationMatrix, epsilon))
260         return false;
262     // diagonal metric must be changed if we consider only inverse
263     diagonalMatrix = (scaleMatrix*scaleMatrix)*diagonalMatrix;
264     if(useNumericalCleanUp)
265         diagonalMatrix = numericalCleanUp(diagonalMatrix, epsilon);
267     // compact form of the diagonal matrix
268     diagonalMetric = diagonalMatrix.diagonal();
270     return true;