diff options
author | mcarrier <mcarrier@636b058d-ea47-450e-bf9e-a15bfbe3eedb> | 2018-02-15 16:41:04 +0000 |
---|---|---|
committer | mcarrier <mcarrier@636b058d-ea47-450e-bf9e-a15bfbe3eedb> | 2018-02-15 16:41:04 +0000 |
commit | 69c683e663329d8410ca77c371f877bcc3bef906 (patch) | |
tree | fb18c914cf4056881b2b31875eb6b44e5ce23895 /src/Persistence_representations/include/gudhi/Persistence_weighted_gaussian.h | |
parent | be131d6f74a9264e15a0b1c1e72fa8967c4518bd (diff) | |
parent | 265484997185f3bf900744406206a2d64ca0a20d (diff) |
integrated kernel code in pers_representation branch
git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/kernels@3249 636b058d-ea47-450e-bf9e-a15bfbe3eedb
Former-commit-id: 58e77263a0da3674e4699cef832b6d357dcf12e2
Diffstat (limited to 'src/Persistence_representations/include/gudhi/Persistence_weighted_gaussian.h')
-rw-r--r-- | src/Persistence_representations/include/gudhi/Persistence_weighted_gaussian.h | 143 |
1 files changed, 143 insertions, 0 deletions
diff --git a/src/Persistence_representations/include/gudhi/Persistence_weighted_gaussian.h b/src/Persistence_representations/include/gudhi/Persistence_weighted_gaussian.h new file mode 100644 index 00000000..2884885c --- /dev/null +++ b/src/Persistence_representations/include/gudhi/Persistence_weighted_gaussian.h @@ -0,0 +1,143 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Mathieu Carriere + * + * Copyright (C) 2018 INRIA (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * 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 License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef PERSISTENCE_WEIGHTED_GAUSSIAN_H_ +#define PERSISTENCE_WEIGHTED_GAUSSIAN_H_ + +#ifdef GUDHI_USE_TBB +#include <tbb/parallel_for.h> +#endif + +// gudhi include +#include <gudhi/read_persistence_from_file.h> + +// standard include +#include <cmath> +#include <iostream> +#include <vector> +#include <limits> +#include <fstream> +#include <sstream> +#include <algorithm> +#include <string> +#include <utility> +#include <functional> +#include <boost/math/constants/constants.hpp> + +double pi = boost::math::constants::pi<double>(); +using PD = std::vector<std::pair<double,double> >; + +namespace Gudhi { +namespace Persistence_representations { + +class Persistence_weighted_gaussian{ + + protected: + PD diagram; + + public: + + Persistence_weighted_gaussian(PD _diagram){diagram = _diagram;} + PD get_diagram(){return this->diagram;} + + + // ********************************** + // Utils. + // ********************************** + + + static double pss_weight(std::pair<double,double> P){ + if(P.second > P.first) return 1; + else return -1; + } + + static double arctan_weight(std::pair<double,double> P){ + return atan(P.second - P.first); + } + + template<class Weight = std::function<double (std::pair<double,double>) > > + std::vector<std::pair<double,double> > Fourier_feat(PD D, std::vector<std::pair<double,double> > Z, Weight weight = arctan_weight){ + int m = D.size(); std::vector<std::pair<double,double> > B; int M = Z.size(); + for(int i = 0; i < M; i++){ + double d1 = 0; double d2 = 0; double zx = Z[i].first; double zy = Z[i].second; + for(int j = 0; j < m; j++){ + double x = D[j].first; double y = D[j].second; + d1 += weight(D[j])*cos(x*zx + y*zy); + d2 += weight(D[j])*sin(x*zx + y*zy); + } + B.emplace_back(d1,d2); + } + return B; + } + + std::vector<std::pair<double,double> > random_Fourier(double sigma, int M = 1000){ + std::normal_distribution<double> distrib(0,1); std::vector<std::pair<double,double> > Z; std::random_device rd; + for(int i = 0; i < M; i++){ + std::mt19937 e1(rd()); std::mt19937 e2(rd()); + double zx = distrib(e1); double zy = distrib(e2); + Z.emplace_back(zx/sigma,zy/sigma); + } + return Z; + } + + + + // ********************************** + // Scalar product + distance. + // ********************************** + + + template<class Weight = std::function<double (std::pair<double,double>) > > + double compute_scalar_product(Persistence_weighted_gaussian second, double sigma, Weight weight = arctan_weight, int m = 1000){ + + PD diagram1 = this->diagram; PD diagram2 = second.diagram; + + if(m == -1){ + int num_pts1 = diagram1.size(); int num_pts2 = diagram2.size(); double k = 0; + for(int i = 0; i < num_pts1; i++) + for(int j = 0; j < num_pts2; j++) + k += weight(diagram1[i])*weight(diagram2[j])*exp(-((diagram1[i].first - diagram2[j].first) * (diagram1[i].first - diagram2[j].first) + + (diagram1[i].second - diagram2[j].second) * (diagram1[i].second - diagram2[j].second)) + /(2*sigma*sigma)); + return k; + } + else{ + std::vector<std::pair<double,double> > z = random_Fourier(sigma, m); + std::vector<std::pair<double,double> > b1 = Fourier_feat(diagram1,z,weight); + std::vector<std::pair<double,double> > b2 = Fourier_feat(diagram2,z,weight); + double d = 0; for(int i = 0; i < m; i++) d += b1[i].first*b2[i].first + b1[i].second*b2[i].second; + return d/m; + } + } + + template<class Weight = std::function<double (std::pair<double,double>) > > + double distance(Persistence_weighted_gaussian second, double sigma, Weight weight = arctan_weight, int m = 1000, double power = 1) { + return std::pow(this->compute_scalar_product(*this, sigma, weight, m) + second.compute_scalar_product(second, sigma, weight, m)-2*this->compute_scalar_product(second, sigma, weight, m), power/2.0); + } + + +}; + +} // namespace Persistence_weighted_gaussian +} // namespace Gudhi + +#endif // PERSISTENCE_WEIGHTED_GAUSSIAN_H_ |