Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2023 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include <iostream>
23 : #include <complex>
24 : #include "gridtools/ActionWithGrid.h"
25 : #include "core/ActionRegister.h"
26 : #ifdef __PLUMED_HAS_FFTW
27 : #include <fftw3.h> // FFTW interface
28 : #endif
29 :
30 : namespace PLMD {
31 : namespace fourier {
32 :
33 : //+PLUMEDOC GRIDANALYSIS FOURIER_TRANSFORM
34 : /*
35 : Compute the Discrete Fourier Transform (DFT) by means of FFTW of data stored on a 2D grid.
36 :
37 : This action can operate on any other action that outputs scalar data on a two-dimensional grid.
38 :
39 : Up to now, even if the input data are purely real the action uses a complex DFT.
40 :
41 : Just as a quick reference, given a 1D array \f$\mathbf{X}\f$ of size \f$n\f$, this action computes the vector \f$\mathbf{Y}\f$ given by
42 :
43 : \f[
44 : Y_k = \sum_{j=0}^{n-1} X_j e^{2\pi\, j k \sqrt{-1}/n}.
45 : \f]
46 :
47 : This can be easily extended to more than one dimension. All the other details can be found at http://www.fftw.org/doc/What-FFTW-Really-Computes.html#What-FFTW-Really-Computes.
48 :
49 : The keyword "FOURIER_PARAMETERS" deserves just a note on the usage. This keyword specifies how the Fourier transform will be normalized. The keyword takes two numerical parameters (\f$a,\,b\f$) that define the normalization according to the following expression
50 :
51 : \f[
52 : \frac{1}{n^{(1-a)/2}} \sum_{j=0}^{n-1} X_j e^{2\pi b\, j k \sqrt{-1}/n}
53 : \f]
54 :
55 : The default values of these parameters are: \f$a=1\f$ and \f$b=1\f$.
56 :
57 : \par Examples
58 :
59 : The following example tells Plumed to compute the complex 2D 'backward' Discrete Fourier Transform by taking the data saved on a grid called 'density', and normalizing the output by \f$ \frac{1}{\sqrt{N_x\, N_y}}\f$, where \f$N_x\f$ and \f$N_y\f$ are the number of data on the grid (it can be the case that \f$N_x\neq N_y\f$):
60 :
61 : \plumedfile
62 : FOURIER_TRANSFORM STRIDE=1 GRID=density FT_TYPE=complex FOURIER_PARAMETERS=0,-1
63 : \endplumedfile
64 :
65 : */
66 : //+ENDPLUMEDOC
67 :
68 :
69 : class FourierTransform : public gridtools::ActionWithGrid {
70 : private:
71 : bool firsttime;
72 : std::string output_type;
73 : bool real_output, store_norm;
74 : std::vector<int> fourier_params;
75 : gridtools::GridCoordinatesObject gridcoords;
76 : public:
77 : static void registerKeywords( Keywords& keys );
78 : explicit FourierTransform(const ActionOptions&ao);
79 0 : void setupOnFirstStep( const bool incalc ) override { plumed_error(); }
80 : unsigned getNumberOfDerivatives() override ;
81 : const gridtools::GridCoordinatesObject& getGridCoordinatesObject() const override ;
82 : std::vector<std::string> getGridCoordinateNames() const override ;
83 0 : void performTask( const unsigned& current, MultiValue& myvals ) const override { plumed_error(); }
84 : void calculate() override ;
85 : };
86 :
87 : PLUMED_REGISTER_ACTION(FourierTransform,"FOURIER_TRANSFORM")
88 :
89 3 : void FourierTransform::registerKeywords( Keywords& keys ) {
90 3 : ActionWithGrid::registerKeywords( keys ); keys.use("ARG");
91 6 : keys.add("optional","FT_TYPE","choose what kind of data you want as output on the grid. Possible values are: ABS = compute the complex modulus of Fourier coefficients (DEFAULT); NORM = compute the norm (i.e. ABS^2) of Fourier coefficients; COMPLEX = store the FFTW complex output on the grid (as a vector).");
92 6 : keys.add("compulsory","FOURIER_PARAMETERS","default","what kind of normalization is applied to the output and if the Fourier transform in FORWARD or BACKWARD. This keyword takes the form FOURIER_PARAMETERS=A,B, where A and B can be 0, 1 or -1. The default values are A=1 (no normalization at all) and B=1 (forward FFT). Other possible choices for A are: "
93 : "A=-1: normalize by the number of data, "
94 : "A=0: normalize by the square root of the number of data (one forward and followed by backward FFT recover the original data). ");
95 6 : keys.addOutputComponent("real","FT_TYPE","the real part of the function");
96 6 : keys.addOutputComponent("imag","FT_TYPE","the imaginary part of the function");
97 3 : keys.setValueDescription("the fourier transform of the input grid");
98 3 : }
99 :
100 1 : FourierTransform::FourierTransform(const ActionOptions&ao):
101 : Action(ao),
102 : ActionWithGrid(ao),
103 1 : firsttime(true),
104 1 : real_output(true),
105 1 : store_norm(false),
106 1 : fourier_params(2)
107 : {
108 1 : if( getPntrToArgument(0)->getRank()!=2 ) error("fourier transform currently only works with two dimensional grids");
109 :
110 : // Get the type of FT
111 2 : parse("FT_TYPE",output_type);
112 1 : if (output_type.length()==0) {
113 0 : log<<" keyword FT_TYPE unset. By default output grid will contain REAL Fourier coefficients\n";
114 2 : } else if ( output_type=="ABS" || output_type=="abs") {
115 0 : log << " keyword FT_TYPE is '"<< output_type << "' : will compute the MODULUS of Fourier coefficients\n";
116 2 : } else if ( output_type=="NORM" || output_type=="norm") {
117 0 : log << " keyword FT_TYPE is '"<< output_type << "' : will compute the NORM of Fourier coefficients\n";
118 0 : store_norm=true;
119 2 : } else if ( output_type=="COMPLEX" || output_type=="complex" ) {
120 1 : log<<" keyword FT_TYPE is '"<< output_type <<"' : output grid will contain the COMPLEX Fourier coefficients\n";
121 1 : real_output=false;
122 0 : } else error("keyword FT_TYPE unrecognized!");
123 :
124 : // Normalize output?
125 2 : std::string params_str; parse("FOURIER_PARAMETERS",params_str);
126 1 : if (params_str=="default") {
127 0 : fourier_params.assign( fourier_params.size(), 1 );
128 0 : log.printf(" default values of Fourier parameters A=%i, B=%i : the output will NOT be normalized and BACKWARD Fourier transform is computed \n", fourier_params[0],fourier_params[1]);
129 : } else {
130 1 : std::vector<std::string> fourier_str = Tools::getWords(params_str, "\t\n ,");
131 1 : if (fourier_str.size()>2) error("FOURIER_PARAMETERS can take just two values");
132 3 : for (unsigned i=0; i<fourier_str.size(); ++i) {
133 2 : Tools::convert(fourier_str[i],fourier_params[i]);
134 2 : if (fourier_params[i]>1 || fourier_params[i]<-1) error("values accepted for FOURIER_PARAMETERS are only -1, 1 or 0");
135 : }
136 1 : log.printf(" Fourier parameters are A=%i, B=%i \n", fourier_params[0],fourier_params[1]);
137 1 : }
138 :
139 1 : std::vector<unsigned> shape( getPntrToArgument(0)->getRank() );
140 1 : if (real_output) {
141 0 : addValueWithDerivatives( shape );
142 : } else {
143 1 : addComponentWithDerivatives( "real", shape );
144 2 : addComponentWithDerivatives( "imag", shape );
145 : }
146 :
147 : unsigned dimension = getPntrToArgument(0)->getRank();
148 1 : gridtools::ActionWithGrid* ag=dynamic_cast<gridtools::ActionWithGrid*>( getPntrToArgument(0)->getPntrToAction() );
149 1 : if( !ag ) error("input action should be a grid");
150 1 : const gridtools::GridCoordinatesObject & gcoords( ag->getGridCoordinatesObject() );
151 2 : if( gcoords.getGridType()=="fibonacci" ) error("cannot fourier transform fibonacci grids");
152 3 : std::vector<bool> ipbc( dimension ); for(unsigned i=0; i<dimension; ++i) ipbc[i] = gcoords.isPeriodic(i);
153 2 : gridcoords.setup( "flat", ipbc, 0, 0.0 ); checkRead();
154 : #ifndef __PLUMED_HAS_FFTW
155 : error("this feature is only available if you compile PLUMED with FFTW");
156 : #endif
157 1 : }
158 :
159 4 : unsigned FourierTransform::getNumberOfDerivatives() {
160 4 : return 2;
161 : }
162 :
163 7 : const gridtools::GridCoordinatesObject& FourierTransform::getGridCoordinatesObject() const {
164 7 : return gridcoords;
165 : }
166 :
167 2 : std::vector<std::string> FourierTransform::getGridCoordinateNames() const {
168 2 : gridtools::ActionWithGrid* ag=dynamic_cast<gridtools::ActionWithGrid*>( getPntrToArgument(0)->getPntrToAction() );
169 2 : return ag->getGridCoordinateNames();
170 : }
171 :
172 1 : void FourierTransform::calculate() {
173 1 : if( firsttime ) {
174 1 : gridtools::ActionWithGrid* ag=dynamic_cast<gridtools::ActionWithGrid*>( getPntrToArgument(0)->getPntrToAction() );
175 1 : const gridtools::GridCoordinatesObject & gcoords( ag->getGridCoordinatesObject() );
176 1 : std::vector<double> fspacing; std::vector<unsigned> snbins( getGridCoordinatesObject().getDimension() );
177 1 : std::vector<std::string> smin( gcoords.getDimension() ), smax( gcoords.getDimension() );
178 3 : for(unsigned i=0; i<getGridCoordinatesObject().getDimension(); ++i) {
179 6 : smin[i]=gcoords.getMin()[i]; smax[i]=gcoords.getMax()[i];
180 : // Compute k-grid extents
181 2 : double dmin, dmax; snbins[i]=gcoords.getNbin(false)[i];
182 2 : Tools::convert(smin[i],dmin); Tools::convert(smax[i],dmax);
183 2 : dmax=2.0*pi*snbins[i]/( dmax - dmin ); dmin=0.0;
184 2 : Tools::convert(dmin,smin[i]); Tools::convert(dmax,smax[i]);
185 : }
186 1 : gridcoords.setBounds( smin, smax, snbins, fspacing ); firsttime=false;
187 3 : for(unsigned i=0; i<getNumberOfComponents(); ++i) getPntrToComponent(i)->setShape( gcoords.getNbin(true) );
188 1 : }
189 :
190 : #ifdef __PLUMED_HAS_FFTW
191 : // *** CHECK CORRECT k-GRID BOUNDARIES ***
192 : //log<<"Real grid boundaries: \n"
193 : // <<" min_x: "<<mygrid->getMin()[0]<<" min_y: "<<mygrid->getMin()[1]<<"\n"
194 : // <<" max_x: "<<mygrid->getMax()[0]<<" max_y: "<<mygrid->getMax()[1]<<"\n"
195 : // <<"K-grid boundaries:"<<"\n"
196 : // <<" min_x: "<<ft_min[0]<<" min_y: "<<ft_min[1]<<"\n"
197 : // <<" max_x: "<<ft_max[0]<<" max_y: "<<ft_max[1]<<"\n";
198 :
199 : // Get the size of the input data arrays (to allocate FFT data)
200 1 : std::vector<unsigned> N_input_data( gridcoords.getNbin(true) );
201 3 : size_t fft_dimension=1; for(unsigned i=0; i<N_input_data.size(); ++i) fft_dimension*=static_cast<size_t>( N_input_data[i] );
202 : // FFT arrays
203 1 : std::vector<std::complex<double> > input_data(fft_dimension), fft_data(fft_dimension);
204 :
205 : // Fill real input with the data on the grid
206 : Value* arg=getPntrToArgument(0);
207 1 : unsigned nargs=arg->getNumberOfValues();
208 1 : std::vector<unsigned> ind( arg->getRank() );
209 10202 : for (unsigned i=0; i<arg->getNumberOfValues(); ++i) {
210 : // Get point indices
211 10201 : gridcoords.getIndices(i, ind);
212 : // Fill input data in row-major order
213 10201 : input_data[ind[0]*N_input_data[0]+ind[1]].real( arg->get( i ) );
214 10201 : input_data[ind[0]*N_input_data[0]+ind[1]].imag( 0.0 );
215 : }
216 :
217 : // *** HERE is the only clear limitation: I'm computing explicitly a 2D FT. It should not happen to deal with other than two-dimensional grid ...
218 1 : fftw_plan plan_complex = fftw_plan_dft_2d(N_input_data[0], N_input_data[1], reinterpret_cast<fftw_complex*>(&input_data[0]), reinterpret_cast<fftw_complex*>(&fft_data[0]), fourier_params[1], FFTW_ESTIMATE);
219 :
220 : // Compute FT
221 1 : fftw_execute( plan_complex );
222 :
223 : // Compute the normalization constant
224 : double norm=1.0;
225 3 : for (unsigned i=0; i<N_input_data.size(); ++i) {
226 2 : norm *= pow( N_input_data[i], (1-fourier_params[0])/2 );
227 : }
228 :
229 : // Save FT data to output grid
230 1 : std::vector<unsigned> N_out_data ( getGridCoordinatesObject().getNbin(true) );
231 1 : std::vector<unsigned> out_ind ( getPntrToArgument(0)->getRank() );
232 10202 : for(unsigned i=0; i<getPntrToArgument(0)->getNumberOfValues(); ++i) {
233 10201 : gridcoords.getIndices( i, out_ind );
234 10201 : if (real_output) {
235 : double ft_value;
236 : // Compute abs/norm and fix normalization
237 0 : if (!store_norm) ft_value=std::abs( fft_data[out_ind[0]*N_out_data[0]+out_ind[1]] / norm );
238 0 : else ft_value=std::norm( fft_data[out_ind[0]*N_out_data[0]+out_ind[1]] / norm );
239 : // Set the value
240 0 : getPntrToComponent(0)->set( i, ft_value);
241 : } else {
242 : double ft_value_real, ft_value_imag;
243 10201 : ft_value_real=fft_data[out_ind[0]*N_out_data[0]+out_ind[1]].real() / norm;
244 10201 : ft_value_imag=fft_data[out_ind[0]*N_out_data[0]+out_ind[1]].imag() / norm;
245 : // Set values
246 10201 : getPntrToComponent(0)->set( i, ft_value_real );
247 10201 : getPntrToComponent(1)->set( i, ft_value_imag );
248 : }
249 : }
250 :
251 : // Free FFTW stuff
252 1 : fftw_destroy_plan(plan_complex);
253 : #endif
254 1 : }
255 :
256 : } // end namespace 'gridtools'
257 : } // end namespace 'PLMD'
|
