1 /***************************************************************************
2 * datac.h is part of Math Graphic Library
3 * Copyright (C) 2007-2016 Alexey Balakin <mathgl.abalakin@gmail.ru> *
4 * *
5 * This program is free software; you can redistribute it and/or modify *
6 * it under the terms of the GNU Lesser General Public License as *
7 * published by the Free Software Foundation; either version 3 of the *
8 * License, or (at your option) any later version. *
9 * *
10 * This program is distributed in the hope that it will be useful, *
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
13 * GNU General Public License for more details. *
14 * *
15 * You should have received a copy of the GNU Lesser General Public *
16 * License along with this program; if not, write to the *
17 * Free Software Foundation, Inc., *
18 * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
19 ***************************************************************************/
20 #ifndef _MGL_DATAC_H_
21 #define _MGL_DATAC_H_
22
23 #include "mgl2/data.h"
24 #include "mgl2/datac_cf.h"
25
26 #if MGL_HAVE_ARMA
27 #include <armadillo>
28 #endif
29 //-----------------------------------------------------------------------------
30 #include <vector>
31 #include <string>
32 //-----------------------------------------------------------------------------
33 #ifndef SWIG
34 dual MGL_EXPORT mglLinearC(const dual *a, long nx, long ny, long nz, mreal x, mreal y, mreal z);
35 dual MGL_EXPORT mglSpline3C(const dual *a, long nx, long ny, long nz, mreal x, mreal y, mreal z,dual *dx=0, dual *dy=0, dual *dz=0);
36 dual MGL_EXPORT mglSpline3Cs(const dual *a, long nx, long ny, long nz, mreal x, mreal y, mreal z);
37 //-----------------------------------------------------------------------------
38 /// Class for working with complex data array
39 class MGL_EXPORT mglDataC : public mglDataA
40 {
41 public:
42 using mglDataA::Momentum;
43 using mglDataA::Last;
44 long nx; ///< number of points in 1st dimensions ('x' dimension)
45 long ny; ///< number of points in 2nd dimensions ('y' dimension)
46 long nz; ///< number of points in 3d dimensions ('z' dimension)
47 dual *a; ///< data array
48 bool link; ///< use external data (i.e. don't free it)
49
50 /// Initiate by other mglDataC variable
mglDataC(const mglDataC & d)51 mglDataC(const mglDataC &d) { a=0; mgl_datac_set(this,&d); } // NOTE: must be constructor for mglDataC& to exclude copy one
mglDataC(const mglDataA & d)52 mglDataC(const mglDataA &d) { a=0; mgl_datac_set(this,&d); }
53 #if MGL_HAVE_RVAL
mglDataC(mglDataC && d)54 mglDataC(mglDataC &&d):nx(d.nx),ny(d.ny),nz(d.nz),a(d.a),link(d.link)
55 { s=d.s; temp=d.temp; func=d.func; o=d.o; id=d.id; d.a=0; d.func=0; }
56 #endif
mglDataC(const mglDataA & re,const mglDataA & im)57 mglDataC(const mglDataA &re, const mglDataA &im) { a=0; mgl_datac_set_ri(this,&re,&im); }
mglDataC(HCDT d)58 mglDataC(HCDT d) { a=0; mgl_datac_set(this, d); }
mglDataC(HCDT re,HCDT im)59 mglDataC(HCDT re, HCDT im) { a=0; mgl_datac_set_ri(this, re, im); }
mglDataC(bool,mglDataC * d)60 mglDataC(bool, mglDataC *d) // NOTE: Variable d will be deleted!!!
61 { if(d)
62 { nx=d->nx; ny=d->ny; nz=d->nz; a=d->a; d->a=0;
63 temp=d->temp; func=d->func; o=d->o; s=d->s;
64 id=d->id; link=d->link; delete d; }
65 else { a=0; Create(1); } }
66 /// Initiate by flat array
mglDataC(int size,const dual * d)67 mglDataC(int size, const dual *d) { a=0; Set(d,size); }
mglDataC(int rows,int cols,const dual * d)68 mglDataC(int rows, int cols, const dual *d) { a=0; Set(d,cols,rows); }
mglDataC(int size,const double * d)69 mglDataC(int size, const double *d) { a=0; Set(d,size); }
mglDataC(int rows,int cols,const double * d)70 mglDataC(int rows, int cols, const double *d) { a=0; Set(d,cols,rows); }
mglDataC(int size,const float * d)71 mglDataC(int size, const float *d) { a=0; Set(d,size); }
mglDataC(int rows,int cols,const float * d)72 mglDataC(int rows, int cols, const float *d) { a=0; Set(d,cols,rows); }
mglDataC(const dual * d,int size)73 mglDataC(const dual *d, int size) { a=0; Set(d,size); }
mglDataC(const dual * d,int rows,int cols)74 mglDataC(const dual *d, int rows, int cols) { a=0; Set(d,cols,rows); }
mglDataC(const double * d,int size)75 mglDataC(const double *d, int size) { a=0; Set(d,size); }
mglDataC(const double * d,int rows,int cols)76 mglDataC(const double *d, int rows, int cols) { a=0; Set(d,cols,rows); }
mglDataC(const float * d,int size)77 mglDataC(const float *d, int size) { a=0; Set(d,size); }
mglDataC(const float * d,int rows,int cols)78 mglDataC(const float *d, int rows, int cols) { a=0; Set(d,cols,rows); }
79 /// Allocate memory and copy data from std::vector<T>
mglDataC(const std::vector<int> & d)80 mglDataC(const std::vector<int> &d) { a=0; Set(d); }
mglDataC(const std::vector<float> & d)81 mglDataC(const std::vector<float> &d) { a=0; Set(d); }
mglDataC(const std::vector<double> & d)82 mglDataC(const std::vector<double> &d) { a=0; Set(d); }
mglDataC(const std::vector<std::complex<double>> & d)83 mglDataC(const std::vector<std::complex<double> > &d) { a=0; Set(d); }
mglDataC(const std::vector<std::complex<float>> & d)84 mglDataC(const std::vector<std::complex<float> > &d) { a=0; Set(d); }
85 /// Read data from file
mglDataC(const char * fname)86 mglDataC(const char *fname) { a=0; Read(fname); }
87 /// Allocate the memory for data array and initialize it zero
88 mglDataC(long xx=1,long yy=1,long zz=1) { a=0; Create(xx,yy,zz); }
89 /// Delete the array
~mglDataC()90 virtual ~mglDataC() { if(!link && a) delete []a; }
91
92 /// Move all data from variable d, and delete this variable.
Move(mglDataC * d)93 inline void Move(mglDataC *d) // NOTE: Variable d will be deleted!!!
94 { if(d && d->GetNN()>1)
95 { bool l=link; dual *b=a;
96 nx=d->nx; ny=d->ny; nz=d->nz; a=d->a; d->a=b;
97 temp=d->temp; func=d->func; o=d->o; s=d->s;
98 id=d->id; link=d->link; d->link=l; delete d; }
99 else if(d) { *this = d->a[0]; delete d; }
100 }
101
102 inline dual GetVal(long i, long j=0, long k=0) const
103 { return mgl_datac_get_value(this,i,j,k);}
104 inline void SetVal(dual f, long i, long j=0, long k=0)
105 { mgl_datac_set_value(this,f,i,j,k); }
106 /// Get sizes
GetNx()107 long GetNx() const { return nx; }
GetNy()108 long GetNy() const { return ny; }
GetNz()109 long GetNz() const { return nz; }
110
111 /// Link external data array (don't delete it at exit)
112 inline void Link(dual *A, long NX, long NY=1, long NZ=1)
113 { mgl_datac_link(this,reinterpret_cast<mdual*>(A),NX,NY,NZ); }
Link(mglDataC & d)114 inline void Link(mglDataC &d) { Link(d.a,d.nx,d.ny,d.nz); }
115 /// Allocate memory and copy the data from the gsl_vector
Set(gsl_vector * m)116 inline void Set(gsl_vector *m) { mgl_datac_set_vector(this,m); }
117 /// Allocate memory and copy the data from the gsl_matrix
Set(gsl_matrix * m)118 inline void Set(gsl_matrix *m) { mgl_datac_set_matrix(this,m); }
119
120 /// Allocate memory and copy the data from the (float *) array
121 inline void Set(const float *A,long NX,long NY=1,long NZ=1)
122 { mgl_datac_set_float(this,A,NX,NY,NZ); }
123 /// Allocate memory and copy the data from the (double *) array
124 inline void Set(const double *A,long NX,long NY=1,long NZ=1)
125 { mgl_datac_set_double(this,A,NX,NY,NZ); }
126 /// Allocate memory and copy the data from the (complex *) array
127 inline void Set(const dual *A,long NX,long NY=1,long NZ=1)
128 { mgl_datac_set_complex(this,reinterpret_cast<const mdual*>(A),NX,NY,NZ); }
129 /// Allocate memory and scanf the data from the string
130 inline void Set(const char *str,long NX,long NY=1,long NZ=1)
131 { mgl_datac_set_values(this,str,NX,NY,NZ); }
132 /// Import data from abstract type
Set(HCDT dat)133 inline void Set(HCDT dat) { mgl_datac_set(this, dat); }
Set(const mglDataA & dat)134 inline void Set(const mglDataA &dat) { mgl_datac_set(this, &dat); }
Set(const mglDataA & re,const mglDataA & im)135 inline void Set(const mglDataA &re, const mglDataA &im) { mgl_datac_set_ri(this, &re, &im); }
Set(HCDT re,HCDT im)136 inline void Set(HCDT re, HCDT im) { mgl_datac_set_ri(this, re, im); }
SetAmpl(const mglDataA & ampl,const mglDataA & phase)137 inline void SetAmpl(const mglDataA &l, const mglDataA &phase)
138 { mgl_datac_set_ap(this, &l, &phase); }
139 /// Allocate memory and copy data from std::vector<T>
Set(const std::vector<int> & d)140 inline void Set(const std::vector<int> &d)
141 { if(d.size()>0) { Create(d.size()); for(long i=0;i<nx;i++) a[i] = d[i]; }
142 else Create(1); }
Set(const std::vector<float> & d)143 inline void Set(const std::vector<float> &d)
144 { if(d.size()>0) Set(&(a[0]),d.size()); else Create(1); }
Set(const std::vector<double> & d)145 inline void Set(const std::vector<double> &d)
146 { if(d.size()>0) Set(&(a[0]),d.size()); else Create(1); }
Set(const std::vector<std::complex<double>> & d)147 inline void Set(const std::vector<std::complex<double> > &d)
148 { if(d.size()>0) { Create(d.size()); for(long i=0;i<nx;i++) a[i] = d[i]; }
149 else Create(1); }
Set(const std::vector<std::complex<float>> & d)150 inline void Set(const std::vector<std::complex<float> > &d)
151 { if(d.size()>0) { Create(d.size()); for(long i=0;i<nx;i++) a[i] = d[i]; }
152 else Create(1); }
153
154 #if MGL_HAVE_ARMA
Set(const arma::cx_mat & d)155 inline void Set(const arma::cx_mat &d)
156 { Create(d.n_rows,d.n_cols); for(long i=0;i<nx*ny;i++) a[i] = d[i]; }
Set(const arma::mat & d)157 inline void Set(const arma::mat &d)
158 { Create(d.n_rows,d.n_cols); for(long i=0;i<nx*ny;i++) a[i] = d[i]; }
Set(const arma::cx_cube & d)159 inline void Set(const arma::cx_cube &d)
160 { Create(d.n_rows,d.n_cols,d.n_slices); for(long i=0;i<nx*ny*nz;i++) a[i] = d[i]; }
Set(const arma::cube & d)161 inline void Set(const arma::cube &d)
162 { Create(d.n_rows,d.n_cols,d.n_slices); for(long i=0;i<nx*ny*nz;i++) a[i] = d[i]; }
Set(const arma::cx_vec & d)163 inline void Set(const arma::cx_vec &d)
164 { Create(d.n_elem); for(long i=0;i<nx;i++) a[i] = d[i]; }
Set(const arma::vec & d)165 inline void Set(const arma::vec &d)
166 { Create(d.n_elem); for(long i=0;i<nx;i++) a[i] = d[i]; }
167 inline arma::cx_mat arma_mat(long k=0) const { return arma::cx_mat(a+k*nx*ny,ny,nx); }
arma_cube()168 inline arma::cx_cube arma_cube() const {return arma::cx_cube(a,nx,ny,nz);}
169 #endif
170
171 /// Create or recreate the array with specified size and fill it by zero
172 inline void Create(long mx,long my=1,long mz=1)
173 { mgl_datac_create(this,mx,my,mz); }
174 /// Rearange data dimensions
175 inline void Rearrange(long mx, long my=0, long mz=0)
176 { mgl_datac_rearrange(this,mx,my,mz); }
177 /// Transpose dimensions of the data (generalization of Transpose)
178 inline void Transpose(const char *dim="yx")
179 { mgl_datac_transpose(this,dim); }
180 /// Extend data dimensions
181 inline void Extend(long n1, long n2=0)
182 { mgl_datac_extend(this,n1,n2); }
183 /// Reduce size of the data
184 inline void Squeeze(long rx,long ry=1,long rz=1,bool smooth=false)
185 { mgl_datac_squeeze(this,rx,ry,rz,smooth); }
186 /// Crop the data
187 inline void Crop(long n1, long n2,char dir='x')
188 { mgl_datac_crop(this,n1,n2,dir); }
189 /// Crop the data to be most optimal for FFT (i.e. to closest value of 2^n*3^m*5^l)
190 inline void Crop(const char *how="235x")
191 { mgl_datac_crop_opt(this, how); }
192 /// Insert data
193 inline void Insert(char dir, long at=0, long num=1)
194 { mgl_datac_insert(this,dir,at,num); }
195 /// Delete data
196 inline void Delete(char dir, long at=0, long num=1)
197 { mgl_datac_delete(this,dir,at,num); }
198 /// Join with another data array
Join(const mglDataA & d)199 inline void Join(const mglDataA &d)
200 { mgl_datac_join(this,&d); }
201
202 /// Modify the data by specified formula
203 inline void Modify(const char *eq,long dim=0)
204 { mgl_datac_modify(this, eq, dim); }
205 /// Modify the data by specified formula
Modify(const char * eq,const mglDataA & vdat,const mglDataA & wdat)206 inline void Modify(const char *eq,const mglDataA &vdat, const mglDataA &wdat)
207 { mgl_datac_modify_vw(this,eq,&vdat,&wdat); }
208 /// Modify the data by specified formula
Modify(const char * eq,const mglDataA & vdat)209 inline void Modify(const char *eq,const mglDataA &vdat)
210 { mgl_datac_modify_vw(this,eq,&vdat,0); }
211 /// Modify the data by specified formula assuming x,y,z in range [r1,r2]
212 inline void Fill(mglBase *gr, const char *eq, const char *opt="")
213 { mgl_datac_fill_eq(gr,this,eq,0,0,opt); }
214 inline void Fill(mglBase *gr, const char *eq, const mglDataA &vdat, const char *opt="")
215 { mgl_datac_fill_eq(gr,this,eq,&vdat,0,opt); }
216 inline void Fill(mglBase *gr, const char *eq, const mglDataA &vdat, const mglDataA &wdat,const char *opt="")
217 { mgl_datac_fill_eq(gr,this,eq,&vdat,&wdat,opt); }
218 /// Equidistantly fill the data to range [x1,x2] in direction dir
219 inline void Fill(dual x1,dual x2=mglNaN,char dir='x')
220 { mgl_datac_fill(this,x1,x2,dir); }
221
222 /// Fill the data by interpolated values of vdat parametrically depended on xdat,ydat,zdat for x,y,z in range [p1,p2] using global spline
223 inline void RefillGS(const mglDataA &xdat, const mglDataA &vdat, mreal x1, mreal x2,long sl=-1)
224 { mgl_datac_refill_gs(this,&xdat,&vdat,x1,x2,sl); }
225 /// Fill the data by interpolated values of vdat parametrically depended on xdat,ydat,zdat for x,y,z in range [p1,p2]
226 inline void Refill(const mglDataA &xdat, const mglDataA &vdat, mreal x1, mreal x2,long sl=-1)
227 { mgl_datac_refill_x(this,&xdat,&vdat,x1,x2,sl); }
228 inline void Refill(const mglDataA &xdat, const mglDataA &vdat, mglPoint p1, mglPoint p2,long sl=-1)
229 { mgl_datac_refill_x(this,&xdat,&vdat,p1.x,p2.x,sl); }
230 inline void Refill(const mglDataA &xdat, const mglDataA &ydat, const mglDataA &vdat, mglPoint p1, mglPoint p2,long sl=-1)
231 { mgl_datac_refill_xy(this,&xdat,&ydat,&vdat,p1.x,p2.x,p1.y,p2.y,sl); }
Refill(const mglDataA & xdat,const mglDataA & ydat,const mglDataA & zdat,const mglDataA & vdat,mglPoint p1,mglPoint p2)232 inline void Refill(const mglDataA &xdat, const mglDataA &ydat, const mglDataA &zdat, const mglDataA &vdat, mglPoint p1, mglPoint p2)
233 { mgl_datac_refill_xyz(this,&xdat,&ydat,&zdat,&vdat,p1.x,p2.x,p1.y,p2.y,p1.z,p2.z); }
234 /// Fill the data by interpolated values of vdat parametrically depended on xdat,ydat,zdat for x,y,z in axis range of gr
235 inline void Refill(HMGL gr, const mglDataA &xdat, const mglDataA &vdat, long sl=-1, const char *opt="")
236 { mgl_datac_refill_gr(gr,this,&xdat,0,0,&vdat,sl,opt); }
237 inline void Refill(HMGL gr, const mglDataA &xdat, const mglDataA &ydat, const mglDataA &vdat, long sl=-1, const char *opt="")
238 { mgl_datac_refill_gr(gr,this,&xdat,&ydat,0,&vdat,sl,opt); }
239 inline void Refill(HMGL gr, const mglDataA &xdat, const mglDataA &ydat, const mglDataA &zdat, const mglDataA &vdat, const char *opt="")
240 { mgl_datac_refill_gr(gr,this,&xdat,&ydat,&zdat,&vdat,-1,opt); }
241
242
243 /// Put value to data element(s)
244 inline void Put(dual val, long i=-1, long j=-1, long k=-1)
245 { mgl_datac_put_val(this,val,i,j,k); }
246 /// Put array to data element(s)
247 inline void Put(const mglDataA &dat, long i=-1, long j=-1, long k=-1)
248 { mgl_datac_put_dat(this,&dat,i,j,k); }
249
250 /// Read data from tab-separated text file with auto determining size
Read(const char * fname)251 inline bool Read(const char *fname)
252 { return mgl_datac_read(this,fname); }
253 /// Read data from text file with specifeid size
254 inline bool Read(const char *fname,long mx,long my=1,long mz=1)
255 { return mgl_datac_read_dim(this,fname,mx,my,mz); }
256 /// Save whole data array (for ns=-1) or only ns-th slice to text file
257 void Save(const char *fname,long ns=-1) const
258 { mgl_datac_save(this,fname,ns); }
259 /// Get whole data array (for ns=-1) or only ns-th slice to string
260 std::string Get(long ns=-1) const
261 { return mgl_datac_to_string(this,ns); }
262
263 /// Read data from tab-separated text files with auto determining size which filenames are result of sprintf(fname,templ,t) where t=from:step:to
264 inline bool ReadRange(const char *templ, double from, double to, double step=1, bool as_slice=false)
265 { return mgl_datac_read_range(this,templ,from,to,step,as_slice); }
266 /// Read data from tab-separated text files with auto determining size which filenames are satisfied to template (like "t_*.dat")
267 inline bool ReadAll(const char *templ, bool as_slice=false)
268 { return mgl_datac_read_all(this, templ, as_slice); }
269 /// Read data from text file with size specified at beginning of the file
270 inline bool ReadMat(const char *fname, long dim=2)
271 { return mgl_datac_read_mat(this,fname,dim); }
272
273 /// Read data array from HDF file (parse HDF4 and HDF5 files)
ReadHDF(const char * fname,const char * data)274 inline int ReadHDF(const char *fname,const char *data)
275 { return mgl_datac_read_hdf(this,fname,data); }
276 /// Save data to HDF file
277 void SaveHDF(const char *fname,const char *data,bool rewrite=false) const
278 { mgl_datac_save_hdf(this,fname,data,rewrite); }
279
280 /// Get real part of data values
Real()281 inline mglData Real() const
282 { return mglData(true,mgl_datac_real(this)); }
283 /// Get imaginary part of data values
Imag()284 inline mglData Imag() const
285 { return mglData(true,mgl_datac_imag(this)); }
286 /// Get absolute value of data values, i.e. |u|
Abs()287 inline mglData Abs() const
288 { return mglData(true,mgl_datac_abs(this)); }
289 /// Get square of absolute value of data values, i.e. |u|^2
Norm()290 inline mglData Norm() const
291 { return mglData(true,mgl_datac_norm(this)); }
292 /// Get argument of data values
Arg()293 inline mglData Arg() const
294 { return mglData(true,mgl_datac_arg(this)); }
295
296 /// Get column (or slice) of the data filled by formulas of named columns
Column(const char * eq)297 inline mglDataC Column(const char *eq) const
298 { return mglDataC(true,mgl_datac_column(this,eq)); }
299 /// Get momentum (1D-array) of data along direction 'dir'. String looks like "x1" for median in x-direction, "x2" for width in x-dir and so on.
Momentum(char dir,const char * how)300 inline mglDataC Momentum(char dir, const char *how) const
301 { return mglDataC(true,mgl_datac_momentum(this,dir,how)); }
302 /// Get sub-array of the data with given fixed indexes
303 inline mglDataC SubData(long xx,long yy=-1,long zz=-1) const
304 { return mglDataC(true,mgl_datac_subdata(this,xx,yy,zz)); }
SubData(const mglDataA & xx,const mglDataA & yy,const mglDataA & zz)305 inline mglDataC SubData(const mglDataA &xx, const mglDataA &yy, const mglDataA &zz) const
306 { return mglDataC(true,mgl_datac_subdata_ext(this,&xx,&yy,&zz)); }
SubData(const mglDataA & xx,const mglDataA & yy)307 inline mglDataC SubData(const mglDataA &xx, const mglDataA &yy) const
308 { return mglDataC(true,mgl_datac_subdata_ext(this,&xx,&yy,0)); }
SubData(const mglDataA & xx)309 inline mglDataC SubData(const mglDataA &xx) const
310 { return mglDataC(true,mgl_datac_subdata_ext(this,&xx,0,0)); }
311 /// Get data from sections ids, separated by value val along specified direction.
312 /** If section id is negative then reverse order is used (i.e. -1 give last section). */
313 inline mglDataC Section(const mglDataA &ids, char dir='y', mreal val=NAN) const
314 { return mglDataC(true,mgl_datac_section(this,&ids,dir,val)); }
315 inline mglDataC Section(long id, char dir='y', mreal val=NAN) const
316 { return mglDataC(true,mgl_datac_section_val(this,id,dir,val)); }
317
318 /// Get trace of the data array
Trace()319 inline mglDataC Trace() const
320 { return mglDataC(true,mgl_datac_trace(this)); }
321 /// Get array which is result of summation in given direction or directions
Sum(const char * dir)322 inline mglDataC Sum(const char *dir) const
323 { return mglDataC(true,mgl_datac_sum(this,dir)); }
324 /// Get the data which is direct multiplication (like, d[i,j] = this[i]*a[j] and so on)
Combine(const mglDataA & dat)325 inline mglDataC Combine(const mglDataA &dat) const
326 { return mglDataC(true,mgl_datac_combine(this,&dat)); }
327 /// Resize the data to new size of box [x1,x2]*[y1,y2]*[z1,z2]
328 inline mglDataC Resize(long mx,long my=1,long mz=1, mreal x1=0,mreal x2=1, mreal y1=0,mreal y2=1, mreal z1=0,mreal z2=1) const
329 { return mglDataC(true,mgl_datac_resize_box(this,mx,my,mz,x1,x2,y1,y2,z1,z2)); }
330 /// Get array which values is result of interpolation this for coordinates from other arrays
331 inline mglDataC Evaluate(const mglData &idat, bool norm=true) const
332 { return mglDataC(true,mgl_datac_evaluate(this,&idat,0,0,norm)); }
333 inline mglDataC Evaluate(const mglData &idat, const mglData &jdat, bool norm=true) const
334 { return mglDataC(true,mgl_datac_evaluate(this,&idat,&jdat,0,norm)); }
335 inline mglDataC Evaluate(const mglData &idat, const mglData &jdat, const mglData &kdat, bool norm=true) const
336 { return mglDataC(true,mgl_datac_evaluate(this,&idat,&jdat,&kdat,norm)); }
337
338 /// Find correlation with another data arrays
Correl(const mglData & dat,const char * dir)339 inline mglDataC Correl(const mglData &dat, const char *dir) const
340 { return mglDataC(true,mgl_datac_correl(this,&dat,dir)); }
341 /// Find auto correlation function
AutoCorrel(const char * dir)342 inline mglDataC AutoCorrel(const char *dir) const
343 { return mglDataC(true,mgl_datac_correl(this,this,dir)); }
344
345 /// Create n-th points distribution of this data values in range [v1, v2]
346 inline mglData Hist(long n,mreal v1=0,mreal v2=1, long nsub=0) const
347 { return mglData(true,mgl_data_hist(this,n,v1,v2,nsub)); }
348 /// Create n-th points distribution of this data values in range [v1, v2] with weight w
349 inline mglData Hist(const mglDataA &w, long n,mreal v1=0,mreal v2=1, long nsub=0) const
350 { return mglData(true,mgl_data_hist_w(this,&w,n,v1,v2,nsub)); }
351 /// Get array of positions of first value, which abs() is large val
First(const char * dir,mreal val)352 inline mglData First(const char *dir, mreal val) const
353 { return mglData(true,mgl_data_first_dir(this,dir,val)); }
354 /// Get array of positions of last value, which abs() is large val
Last(const char * dir,mreal val)355 inline mglData Last(const char *dir, mreal val) const
356 { return mglData(true,mgl_data_last_dir(this,dir,val)); }
357
358 /// Get array which is result of maximal values in given direction or directions
Max(const char * dir)359 inline mglData Max(const char *dir) const
360 { return mglData(true,mgl_data_max_dir(this,dir)); }
361 /// Get array which is result of minimal values in given direction or directions
Min(const char * dir)362 inline mglData Min(const char *dir) const
363 { return mglData(true,mgl_data_min_dir(this,dir)); }
364 /// Find roots for set of nonlinear equations defined by textual formula
MultiRoots(const char * eq,const char * vars)365 inline mglDataC MultiRoots(const char *eq, const char *vars) const
366 { return mglDataC(true,mgl_find_roots_txt_c(eq, vars, this)); }
367
368 /// Cumulative summation the data in given direction or directions
CumSum(const char * dir)369 inline void CumSum(const char *dir) { mgl_datac_cumsum(this,dir); }
370 /// Integrate (cumulative summation) the data in given direction or directions
Integral(const char * dir)371 inline void Integral(const char *dir) { mgl_datac_integral(this,dir); }
372 /// Differentiate the data in given direction or directions
Diff(const char * dir)373 inline void Diff(const char *dir) { mgl_datac_diff(this,dir); }
374 /// Differentiate the parametrically specified data along direction v1
Diff(const mglDataA & v1)375 inline void Diff(const mglDataA &v1)
376 { mgl_datac_diff_par(this,&v1,0,0); }
377 /// Differentiate the parametrically specified data along direction v1 with v2=const
Diff(const mglDataA & v1,const mglDataA & v2)378 inline void Diff(const mglDataA &v1, const mglDataA &v2)
379 { mgl_datac_diff_par(this,&v1,&v2,0); }
380 /// Differentiate the parametrically specified data along direction v1 with v2,v3=const
Diff(const mglDataA & v1,const mglDataA & v2,const mglDataA & v3)381 inline void Diff(const mglDataA &v1, const mglDataA &v2, const mglDataA &v3)
382 { mgl_datac_diff_par(this,&v1,&v2,&v3); }
383 /// Double-differentiate (like laplace operator) the data in given direction
Diff2(const char * dir)384 inline void Diff2(const char *dir) { mgl_datac_diff2(this,dir); }
385
386 /// Swap left and right part of the data in given direction (useful for fourier spectrums)
Swap(const char * dir)387 inline void Swap(const char *dir) { mgl_datac_swap(this,dir); }
388 /// Roll data along direction dir by num slices
Roll(char dir,long num)389 inline void Roll(char dir, long num) { mgl_datac_roll(this,dir,num); }
390 /// Mirror the data in given direction (useful for fourier spectrums)
Mirror(const char * dir)391 inline void Mirror(const char *dir) { mgl_datac_mirror(this,dir); }
392 /// Smooth the data on specified direction or directions
393 /** String \a dir may contain:
394 * ‘x’, ‘y’, ‘z’ for 1st, 2nd or 3d dimension;
395 * ‘dN’ for linear averaging over N points;
396 * ‘3’ for linear averaging over 3 points;
397 * ‘5’ for linear averaging over 5 points.
398 * By default quadratic averaging over 5 points is used. */
399 inline void Smooth(const char *dirs="xyz",mreal delta=0)
400 { mgl_datac_smooth(this,dirs,delta); }
401 /// Limit the data to be inside [-v,v], keeping the original phase
Limit(mreal v)402 inline void Limit(mreal v)
403 { mgl_datac_limit(this, v); }
404 /// Keep the data phase/value along line i and j in given direction.
405 /** Parameter "how" may contain: 'x','y' or 'z' for direction (default is 'y'); 'a' for keeping amplitude instead of phase.*/
406 inline void Keep(const char *how, long i, long j=0)
407 { mgl_datac_keep(this, how, i, j); }
408
409 /// Set as the data envelop
410 inline void Envelop(char dir='x') { mgl_datac_envelop(this,dir); }
411 /// Hankel transform
Hankel(const char * dir)412 inline void Hankel(const char *dir) { mgl_datac_hankel(this,dir); }
413 /// Apply Sin-Fourier transform
SinFFT(const char * dir)414 inline void SinFFT(const char *dir) { mgl_datac_sinfft(this,dir); }
415 /// Apply Cos-Fourier transform
CosFFT(const char * dir)416 inline void CosFFT(const char *dir) { mgl_datac_cosfft(this,dir); }
417 /// Fourier transform
FFT(const char * dir)418 inline void FFT(const char *dir) { mgl_datac_fft(this,dir); }
419 /// Calculate one step of diffraction by finite-difference method with parameter q
420 /** Parameter \a how may contain:
421 * ‘x‘,‘y‘,‘z‘ or ‘r‘ for directions or axial along x,
422 * ‘e‘,‘g‘,‘0‘,‘1‘,‘2‘,‘3‘ for boundary conditions: exponential, Gaussian, zero, constant, linear, quadratic.*/
Diffraction(const char * how,mreal q)423 inline void Diffraction(const char *how, mreal q) { mgl_datac_diffr(this,how,q); }
424 /// Apply wavelet transform
425 /** Parameter \a dir may contain:
426 * ‘x‘,‘y‘,‘z‘ for directions,
427 * ‘d‘ for daubechies, ‘D‘ for centered daubechies,
428 * ‘h‘ for haar, ‘H‘ for centered haar,
429 * ‘b‘ for bspline, ‘B‘ for centered bspline,
430 * ‘i‘ for applying inverse transform. */
Wavelet(const char * how,int k)431 inline void Wavelet(const char *how, int k) { mgl_datac_wavelet(this, how, k); }
432
433 /// Interpolate by cubic spline the data to given point x=[0...nx-1], y=[0...ny-1], z=[0...nz-1]
434 inline dual Spline(mreal x,mreal y=0,mreal z=0) const
435 { return mgl_datac_spline(this, x,y,z); }
436 /// Interpolate by cubic spline the data to given point x,\a y,\a z which normalized in range [0, 1]
437 inline dual Spline1(mreal x,mreal y=0,mreal z=0) const
438 { return mgl_datac_spline(this, x*(nx-1),y*(ny-1),z*(nz-1)); }
439 /// Interpolate by linear function the data to given point x=[0...nx-1], y=[0...ny-1], z=[0...nz-1]
440 inline dual Linear(mreal x,mreal y=0,mreal z=0) const
441 { return mgl_datac_linear_ext(this,x,y,z,0,0,0); }
442 /// Interpolate by line the data to given point x,\a y,\a z which normalized in range [0, 1]
443 inline dual Linear1(mreal x,mreal y=0,mreal z=0) const
444 { return mgl_datac_linear_ext(this,x*(nx-1),y*(ny-1),z*(nz-1),0,0,0); }
445 /// Interpolate by linear function the data and return its derivatives at given point x=[0...nx-1], y=[0...ny-1], z=[0...nz-1]
446 inline dual Linear(mglPoint &dif, mreal x,mreal y=0,mreal z=0) const
447 {
448 mdual val,dx,dy,dz;
449 val = mgl_datac_linear_ext(this,x,y,z, &dx, &dy, &dz);
450 dif.Set(dx.real(),dy.real(),dz.real()); return val;
451 }
452 /// Interpolate by line the data and return its derivatives at given point x,\a y,\a z which normalized in range [0, 1]
453 inline dual Linear1(mglPoint &dif, mreal x,mreal y=0,mreal z=0) const
454 {
455 mdual val,dx,dy,dz;
456 val = mgl_datac_linear_ext(this,x,y,z, &dx, &dy, &dz);
457 dif.Set(dx.real(),dy.real(),dz.real());
458 dif.x/=nx>1?nx-1:1; dif.y/=ny>1?ny-1:1; dif.z/=nz>1?nz-1:1;
459 return val;
460 }
461 /// Return an approximated x-value (root) when dat(x) = val
462 inline mreal Solve(mreal val, bool use_spline=true, long i0=0) const
463 { return mgl_data_solve_1d(this, val, use_spline, i0); }
464 /// Return an approximated value (root) when dat(x) = val
465 inline mglData Solve(mreal val, char dir, bool norm=true) const
466 { return mglData(true,mgl_data_solve(this, val, dir, 0, norm)); }
467 inline mglData Solve(mreal val, char dir, const mglData &i0, bool norm=true) const
468 { return mglData(true,mgl_data_solve(this, val, dir, &i0, norm)); }
469
470 /// Copy data from other mglDataA variable
471 inline const mglDataA &operator=(const mglDataA &d)
472 { if(this!=&d) Set(&d); return d; }
473 inline const mglDataC &operator=(const mglDataC &d)
474 { if(this!=&d) Set(&d); return d; }
475 inline dual operator=(dual val)
476 {
477 #pragma omp parallel for
478 for(long i=0;i<nx*ny*nz;i++) a[i]=val; return val; }
479 inline dual operator=(mreal val)
480 {
481 #pragma omp parallel for
482 for(long i=0;i<nx*ny*nz;i++) a[i]=val; return val; }
483 /// Multiply the data by other one for each element
484 inline void operator*=(const mglDataA &d) { mgl_datac_mul_dat(this,&d); }
485 /// Divide the data by other one for each element
486 inline void operator/=(const mglDataA &d) { mgl_datac_div_dat(this,&d); }
487 /// Add the other data
488 inline void operator+=(const mglDataA &d) { mgl_datac_add_dat(this,&d); }
489 /// Subtract the other data
490 inline void operator-=(const mglDataA &d) { mgl_datac_sub_dat(this,&d); }
491 /// Multiply each element by the number
492 inline void operator*=(dual d) { mgl_datac_mul_num(this,d); }
493 /// Divide each element by the number
494 inline void operator/=(dual d) { mgl_datac_div_num(this,d); }
495 /// Add the number
496 inline void operator+=(dual d) { mgl_datac_add_num(this,d); }
497 /// Subtract the number
498 inline void operator-=(dual d) { mgl_datac_sub_num(this,d); }
499 #ifndef SWIG
500 /// Direct access to the data cell
501 inline dual operator[](long i) const { return a[i]; }
502 inline dual &operator[](long i) { return a[i]; }
503 #endif
504
505
506 #ifndef DEBUG
507 /// Get the value in given cell of the data
508 mreal v(long i,long j=0,long k=0) const { return abs(a[i+nx*(j+ny*k)]); }
509 /// Set the value in given cell of the data
510 void set_v(mreal val, long i,long j=0,long k=0) { a[i+nx*(j+ny*k)]=val; }
511 #else
512 /// Get the value in given cell of the data with border checking
513 mreal v(long i,long j=0,long k=0) const { return abs(mgl_datac_get_value(this,i,j,k)); }
514 /// Set the value in given cell of the data
515 void set_v(mreal val, long i,long j=0,long k=0) { mgl_datac_set_value(this,val,i,j,k); }
516 #endif
517 /// Get the complex value in given cell of the data
518 dual vc(long i,long j=0,long k=0) const { return a[i+nx*(j+ny*k)]; }
vcthr(long i)519 dual vcthr(long i) const { return a[i]; }
520 /// Get the interpolated value and its derivatives in given data cell without border checking
521 mreal valueD(mreal x,mreal y=0,mreal z=0,mreal *dx=0,mreal *dy=0,mreal *dz=0) const
522 { dual aa,ax,ay,az; mreal res;
523 aa = mglSpline3C(a,nx,ny,nz,x,y,z,&ax,&ay,&az); res = abs(aa);
524 if(dx) *dx = res?(real(aa)*real(ax)+imag(aa)*imag(ax))/res:0;
525 if(dy) *dy = res?(real(aa)*real(ay)+imag(aa)*imag(ay))/res:0;
526 if(dz) *dz = res?(real(aa)*real(az)+imag(aa)*imag(az))/res:0;
527 return res; }
528 /// Get the interpolated value in given data cell without border checking
529 mreal value(mreal x,mreal y=0,mreal z=0) const
530 { return abs(mglSpline3Cs(a,nx,ny,nz,x,y,z)); }
vthr(long i)531 mreal vthr(long i) const { return abs(a[i]); }
532 // add for speeding up !!!
533 mreal dvx(long i,long j=0,long k=0) const
534 { long i0 = size_t(i)<size_t(nx-1) ? i+nx*(j+ny*k):nx*(1+j+ny*k)-2; return abs(a[i0+1]-a[i0]); }
535 // { long i0=i+nx*(j+ny*k);
536 // return i>0? abs(i<nx-1? (a[i0+1]-a[i0-1])/mreal(2):a[i0]-a[i0-1]) : abs(a[i0+1]-a[i0]); }
537 mreal dvy(long i,long j=0,long k=0) const
538 { long i0 = size_t(j)<size_t(ny-1) ? i+nx*(j+ny*k):i+nx*(ny*(k+1)-2); return abs(a[i0+nx]-a[i0]); }
539 // { long i0=i+nx*(j+ny*k);
540 // return j>0? abs(j<ny-1? (a[i0+nx]-a[i0-nx])/mreal(2):a[i0]-a[i0-nx]) : abs(a[i0+nx]-a[i0]);}
541 mreal dvz(long i,long j=0,long k=0) const
542 { long n=nx*ny, i0 = size_t(k)<size_t(nz-1) ? i+nx*(j+ny*k):i+nx*(j+ny*(nz-2)); return abs(a[i0+n]-a[i0]); }
543 // { long i0=i+nx*(j+ny*k), n=nx*ny;
544 // return k>0? abs(k<nz-1? (a[i0+n]-a[i0-n])/mreal(2):a[i0]-a[i0-n]) : abs(a[i0+n]-a[i0]); }
545 };
546 //-----------------------------------------------------------------------------
547 /// Saves result of PDE solving (|u|^2) for "Hamiltonian" ham with initial conditions ini
548 inline mglDataC mglPDEc(mglBase *gr, const char *ham, const mglDataA &ini_re, const mglDataA &ini_im, mreal dz=0.1, mreal k0=100,const char *opt="")
549 { return mglDataC(true, mgl_pde_solve_c(gr,ham, &ini_re, &ini_im, dz, k0,opt)); }
550 /// Saves result of PDE solving for "Hamiltonian" ham with initial conditions ini along a curve ray (must have nx>=7 - x,y,z,px,py,pz,tau or nx=5 - x,y,px,py,tau)
551 inline mglDataC mglQO2dc(const char *ham, const mglDataA &ini_re, const mglDataA &ini_im, const mglDataA &ray, mreal r=1, mreal k0=100)
552 { return mglDataC(true, mgl_qo2d_solve_c(ham, &ini_re, &ini_im, &ray, r, k0, 0, 0)); }
553 inline mglDataC mglQO2dc(const char *ham, const mglDataA &ini_re, const mglDataA &ini_im, const mglDataA &ray, mglData &xx, mglData &yy, mreal r=1, mreal k0=100)
554 { return mglDataC(true, mgl_qo2d_solve_c(ham, &ini_re, &ini_im, &ray, r, k0, &xx, &yy)); }
555 /// Saves result of PDE solving for "Hamiltonian" ham with initial conditions ini along a curve ray (must have nx>=7 - x,y,z,px,py,pz,tau or nx=5 - x,y,px,py,tau)
556 inline mglDataC mglQO3dc(const char *ham, const mglDataA &ini_re, const mglDataA &ini_im, const mglDataA &ray, mreal r=1, mreal k0=100)
557 { return mglDataC(true, mgl_qo3d_solve_c(ham, &ini_re, &ini_im, &ray, r, k0, 0, 0, 0)); }
558 inline mglDataC mglQO3dc(const char *ham, const mglDataA &ini_re, const mglDataA &ini_im, const mglDataA &ray, mglData &xx, mglData &yy, mglData &zz, mreal r=1, mreal k0=100)
559 { return mglDataC(true, mgl_qo3d_solve_c(ham, &ini_re, &ini_im, &ray, r, k0, &xx, &yy, &zz)); }
560 /// Saves result of ODE solving for var complex variables with right part func (separated by ';') and initial conditions x0 over time interval [0,tmax] with time step dt
561 inline mglDataC mglODEc(const char *func, const char *var, const mglDataA &ini, mreal dt=0.1, mreal tmax=10)
562 { return mglDataC(true, mgl_ode_solve_str_c(func,var, &ini, dt, tmax)); }
563 /// Saves result of ODE solving for var complex variables with right part func (separated by ';') and initial conditions x0 of size n*m over time interval [0,tmax] with time step dt
564 inline mglDataC mglODEcs(const char *func, const char *var, char brd, const mglDataA &ini, mreal dt=0.1, mreal tmax=10)
565 { return mglDataC(true, mgl_ode_solve_set_c(func,var, brd, &ini, dt, tmax)); }
566 //-----------------------------------------------------------------------------
567 /// Get array as solution of tridiagonal system of equations a[i]*x[i-1]+b[i]*x[i]+c[i]*x[i+1]=d[i]
568 /** String \a how may contain:
569 * 'x', 'y', 'z' for solving along x-,y-,z-directions, or
570 * 'h' for solving along hexagonal direction at x-y plain (need nx=ny),
571 * 'c' for using periodical boundary conditions,
572 * 'd' for diffraction/diffuse calculation. */
mglTridMatC(const mglDataA & A,const mglDataA & B,const mglDataA & C,const mglDataC & D,const char * how)573 inline mglDataC mglTridMatC(const mglDataA &A, const mglDataA &B, const mglDataA &C, const mglDataC &D, const char *how)
574 { return mglDataC(true, mgl_datac_tridmat(&A, &B, &C, &D, how)); }
575 //-----------------------------------------------------------------------------
576 /// Get sub-array of the data with given fixed indexes
577 inline mglDataC mglSubDataC(const mglDataA &dat, long xx, long yy=-1, long zz=-1)
578 { return mglDataC(true,mgl_datac_subdata(&dat,xx,yy,zz)); }
mglSubDataC(const mglDataA & dat,const mglDataA & xx,const mglDataA & yy,const mglDataA & zz)579 inline mglDataC mglSubDataC(const mglDataA &dat, const mglDataA &xx, const mglDataA &yy, const mglDataA &zz)
580 { return mglDataC(true,mgl_datac_subdata_ext(&dat,&xx,&yy,&zz)); }
mglSubDataC(const mglDataA & dat,const mglDataA & xx,const mglDataA & yy)581 inline mglDataC mglSubDataC(const mglDataA &dat, const mglDataA &xx, const mglDataA &yy)
582 { return mglDataC(true,mgl_datac_subdata_ext(&dat,&xx,&yy,0)); }
mglSubDataC(const mglDataA & dat,const mglDataA & xx)583 inline mglDataC mglSubDataC(const mglDataA &dat, const mglDataA &xx)
584 { return mglDataC(true,mgl_datac_subdata_ext(&dat,&xx,0,0)); }
585 //-----------------------------------------------------------------------------
586 /// Prepare coefficients for global spline interpolation
mglGSplineCInit(const mglDataA & xdat,const mglDataA & ydat)587 inline mglDataC mglGSplineCInit(const mglDataA &xdat, const mglDataA &ydat)
588 { return mglDataC(true,mgl_gsplinec_init(&xdat, &ydat)); }
589 /// Evaluate global spline (and its derivatives d1, d2 if not NULL) using prepared coefficients \a coef
590 inline dual mglGSplineC(const mglDataA &coef, mreal dx, dual *d1=0, dual *d2=0)
591 { return mgl_gsplinec(&coef, dx, reinterpret_cast<mdual*>(d1), reinterpret_cast<mdual*>(d2)); }
592 //-----------------------------------------------------------------------------
593 /// Evaluate formula 'str' for given list of variables 'vars'.
594 /** NOTE: you need to delete returned data array!*/
595 HADT MGL_EXPORT mglFormulaCalcC(const char *str, const std::vector<mglDataA*> &vars);
596 //-----------------------------------------------------------------------------
597 #define _DN_(a) ((mglDataC *)*(a))
598 #define _DC_ ((mglDataC *)*d)
599 //-----------------------------------------------------------------------------
600 #ifndef SWIG
601 /// Wrapper class for complex expression evaluating
602 class MGL_EXPORT mglExprC
603 {
604 HAEX ex;
mglExprC(const mglExprC &)605 mglExprC(const mglExprC &){ex=0;} // copying is not allowed
606 const mglExprC &operator=(const mglExprC &t){ex=0; return t;} // copying is not allowed
607 public:
mglExprC(const char * expr)608 mglExprC(const char *expr) { ex = mgl_create_cexpr(expr); }
~mglExprC()609 ~mglExprC() { mgl_delete_cexpr(ex); }
610 /// Return value of expression for given x,y,z variables
611 inline dual Eval(dual x, dual y=0, dual z=0)
612 { return mgl_cexpr_eval(ex,x,y,z); }
613 /// Return value of expression for given x,y,z,u,v,w variables
Eval(dual x,dual y,dual z,dual u,dual v,dual w)614 inline dual Eval(dual x, dual y, dual z, dual u, dual v, dual w)
615 {
616 mdual var[26];
617 var['x'-'a']=x; var['y'-'a']=y; var['z'-'a']=z;
618 var['u'-'a']=u; var['v'-'a']=v; var['w'-'a']=w;
619 return mgl_cexpr_eval_v(ex,var); }
620 /// Return value of expression for given variables
Eval(dual var[26])621 inline dual Eval(dual var[26])
622 { return mgl_cexpr_eval_v(ex,reinterpret_cast<mdual*>(var)); }
623 };
624 #endif
625 //-----------------------------------------------------------------------------
626 #endif
627 #endif
628