source: trunk/src/MathUtils.cpp@ 2222

Last change on this file since 2222 was 2186, checked in by WataruKawasaki, 13 years ago

New Development: Yes

JIRA Issue: Yes CAS-3149

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: scantable.*sinusoid_baseline() params

Test Programs:

Put in Release Notes: Yes

Module(s):

Description: (1) Implemented an automated sinusoidal fitting functionality

(2) FFT available with scantable.fft()
(3) fixed a bug of parsing 'edge' param used by linefinder.
(4) a function to show progress status for row-based iterations.


  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 9.1 KB
Line 
1//#---------------------------------------------------------------------------
2//# MathUtilities.cc: General math operations
3//#---------------------------------------------------------------------------
4//# Copyright (C) 2004
5//# ATNF
6//#
7//# This program is free software; you can redistribute it and/or modify it
8//# under the terms of the GNU General Public License as published by the Free
9//# Software Foundation; either version 2 of the License, or (at your option)
10//# any later version.
11//#
12//# This program is distributed in the hope that it will be useful, but
13//# WITHOUT ANY WARRANTY; without even the implied warranty of
14//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
15//# Public License for more details.
16//#
17//# You should have received a copy of the GNU General Public License along
18//# with this program; if not, write to the Free Software Foundation, Inc.,
19//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
20//#
21//# Correspondence concerning this software should be addressed as follows:
22//# Internet email: Malte.Marquarding@csiro.au
23//# Postal address: Malte Marquarding,
24//# Australia Telescope National Facility,
25//# P.O. Box 76,
26//# Epping, NSW, 2121,
27//# AUSTRALIA
28//#
29//# $Id:
30//#---------------------------------------------------------------------------
31
32#include <casa/aips.h>
33#include <casa/Arrays/Vector.h>
34#include <casa/Arrays/Slice.h>
35#include <casa/Arrays/MaskedArray.h>
36#include <casa/Arrays/MaskArrMath.h>
37#include <casa/Arrays/VectorSTLIterator.h>
38#include <casa/BasicSL/String.h>
39#include <scimath/Mathematics/MedianSlider.h>
40#include <casa/Exceptions/Error.h>
41
42#include <scimath/Fitting/LinearFit.h>
43#include <scimath/Functionals/Polynomial.h>
44#include <scimath/Mathematics/AutoDiff.h>
45
46
47#include "MathUtils.h"
48
49using namespace casa;
50
51float mathutil::statistics(const String& which,
52 const MaskedArray<Float>& data)
53{
54 String str(which);
55 str.upcase();
56 if (str.matches(String("MIN"))) {
57 return min(data);
58 } else if (str.matches(String("MAX"))) {
59 return max(data);
60 } else if (str.matches(String("SUMSQ"))) {
61 return sumsquares(data);
62 } else if (str.matches(String("SUM"))) {
63 return sum(data);
64 } else if (str.matches(String("MEAN"))) {
65 return mean(data);
66 } else if (str.matches(String("VAR"))) {
67 return variance(data);
68 } else if (str.matches(String("STDDEV"))) {
69 return stddev(data);
70 } else if (str.matches(String("AVDEV"))) {
71 return avdev(data);
72 } else if (str.matches(String("RMS"))) {
73 uInt n = data.nelementsValid();
74 return sqrt(sumsquares(data)/n);
75 } else if (str.matches(String("MEDIAN"))) {
76 return median(data);
77 } else {
78 String msg = str + " is not a valid type of statistics";
79 throw(AipsError(msg));
80 }
81 return 0.0;
82}
83
84IPosition mathutil::minMaxPos(const String& which,
85 const MaskedArray<Float>& data)
86{
87 Float minVal, maxVal;
88 IPosition minPos(data.ndim(), 0), maxPos(data.ndim(), 0);
89 minMax(minVal, maxVal, minPos, maxPos, data);
90 String str(which);
91 str.upcase();
92 if (str.contains(String("MIN"))) {
93 return minPos;
94 } else if (str.contains(String("MAX"))) {
95 return maxPos;
96 } else {
97 String msg = str + " is not a valid type of statistics";
98 throw(AipsError(msg));
99 }
100 //return 0.0;
101}
102
103void mathutil::replaceMaskByZero(Vector<Float>& data, const Vector<Bool>& mask)
104{
105 for (uInt i=0; i<data.nelements(); i++) {
106 if (!mask[i]) data[i] = 0.0;
107 }
108}
109
110
111std::vector<std::string> mathutil::tovectorstring(const Vector<String>& in)
112{
113 std::vector<std::string> out;
114 out.reserve(in.nelements());
115 for (Array<String>::const_iterator it = in.begin(); it != in.end(); ++it) {
116 out.push_back(*it);
117 }
118 return out;
119}
120
121Vector<String> mathutil::toVectorString(const std::vector<std::string>& in)
122{
123 Vector<String> out(in.size());
124 Array<String>::iterator oit = out.begin();
125 for (std::vector<std::string>::const_iterator it=in.begin() ;
126 it != in.end(); ++it,++oit) {
127 *oit = *it;
128 }
129 return out;
130}
131
132void mathutil::hanning(Vector<Float>& out, Vector<Bool>& outmask,
133 const Vector<Float>& in, const Vector<Bool>& mask,
134 Bool relaxed, Bool ignoreOther) {
135 (void) ignoreOther; //suppress unused warning
136 Vector< Vector<Float> > weights(8);
137 Vector<Float> vals(3);
138 vals = 0.0;weights[0] = vals;// FFF
139 vals[0] = 1.0; vals[1] = 0.0; vals[2] = 0.0; weights[1] = vals;// TFF
140 vals[0] = 0.0; vals[1] = 1.0; vals[2] = 0.0; weights[2] = vals;// FTF
141 vals[0] = 1.0/3.0; vals[1] = 2.0/3.0; vals[2] = 0.0; weights[3] = vals;// TTF
142 vals[0] = 0.0; vals[1] = 0.0; vals[2] = 1.0;weights[4] = vals;// FFT
143 vals[0] = 0.5; vals[1] = 0.0; vals[2] = 0.5; weights[5] = vals;// TFT
144 vals[0] = 0.0; vals[1] = 2.0/3.0; vals[2] = 1.0/3.0; weights[6] = vals;// FTT
145 vals[0] = 0.25; vals[1] = 0.5; vals[2] = 0.25; weights[7] = vals;// TTT
146 // Chris' case
147 Vector<Bool> weighted(8);
148 if (relaxed) {
149 weighted = False;
150 weighted[7] = True;
151
152 } else {
153 weighted = True;
154 weighted[0] = False;
155 }
156
157 out.resize(in.nelements());
158 outmask.resize(mask.nelements());
159 // make special case for first and last
160 /// ...here
161 // loop from 1..n-2
162 out.resize(in.nelements());
163 out[0] = in[0];out[out.nelements()-1] = in[in.nelements()-1];
164 outmask.resize(mask.nelements());
165 outmask[0] = mask[0]; outmask[outmask.nelements()-1] = mask[mask.nelements()-1];
166 uInt m;Vector<Float>* w;
167 for (uInt i=1; i < out.nelements()-1;++i) {
168 m = mask[i-1] + 2*mask[i] + 4*mask[i+1];
169 w = &(weights[m]);
170 if (weighted[m]) {
171 out[i] = (*w)[0]*in[i-1] + (*w)[1]*in[i] + (*w)[2]*in[i+1];
172 } else { // mask it
173 out[i] = in[i];//use arbitrary value
174 }
175 outmask[i] = mask[i];
176 }
177}
178
179
180void mathutil::runningMedian(Vector<Float>& out, Vector<Bool>& outflag,
181 const Vector<Float>& in, const Vector<Bool>& flag,
182 float width)
183{
184 uInt hwidth = Int(width+0.5);
185 uInt fwidth = hwidth*2+1;
186 out.resize(in.nelements());
187 outflag.resize(flag.nelements());
188 MedianSlider ms(hwidth);
189 Slice sl(0, fwidth-1);
190 Float medval = ms.add(const_cast<Vector<Float>& >(in)(sl),
191 const_cast<Vector<Bool>& >(flag)(sl));
192 (void) medval;//suppress unused warning
193 uInt n = in.nelements();
194 for (uInt i=hwidth; i<(n-hwidth); ++i) {
195 // add data value
196 out[i] = ms.add(in[i+hwidth], flag[i+hwidth]);
197 outflag[i] = (ms.nval() == 0);
198 }
199 // replicate edge values from first value with full width of values
200 for (uInt i=0;i<hwidth;++i) {
201 out[i] = out[hwidth];
202 outflag[i] = outflag[hwidth];
203 out[n-1-i] = out[n-1-hwidth];
204 outflag[n-1-i] = outflag[n-1-hwidth];
205 }
206}
207
208void mathutil::polyfit(Vector<Float>& out, Vector<Bool>& outmask,
209 const Vector<Float>& in, const Vector<Bool>& mask,
210 float width, int order)
211{
212 uInt hwidth = Int(width+0.5);
213 uInt fwidth = hwidth*2+1;
214 out.resize(in.nelements());
215 outmask.resize(mask.nelements());
216 LinearFit<Float> fitter;
217 Polynomial<Float> poly(order);
218 fitter.setFunction(poly);
219 Vector<Float> sigma(fwidth);
220 sigma = 1.0;
221 Vector<Float> parms;
222 Vector<Float> x(fwidth);
223 indgen(x);
224
225 uInt n = in.nelements();
226
227 for (uInt i=hwidth; i<(n-hwidth); ++i) {
228 // add data value
229 if (mask[i]) {
230 Slice sl(i-hwidth, fwidth);
231 const Vector<Float> &y = const_cast<Vector<Float>& >(in)(sl);
232 const Vector<Bool> &m = const_cast<Vector<Bool>& >(mask)(sl);
233 parms = fitter.fit(x, y, sigma, &m);
234
235 poly.setCoefficients(parms);
236 out[i] = poly(x[hwidth]);//cout << in[i] <<"->"<<out[i]<<endl;
237 } else {
238 out[i] = in[i];
239 }
240 outmask[i] = mask[i];
241 }
242 // replicate edge values from first value with full width of values
243 for (uInt i=0;i<hwidth;++i) {
244 out[i] = out[hwidth];
245 outmask[i] = outmask[hwidth];
246 out[n-1-i] = out[n-1-hwidth];
247 outmask[n-1-i] = outmask[n-1-hwidth];
248 }
249}
250
251void mathutil::doZeroOrderInterpolation(casa::Vector<casa::Float>& data,
252 std::vector<bool>& mask) {
253 int fstart = -1;
254 int fend = -1;
255 for (uInt i = 0; i < mask.size(); ++i) {
256 if (!mask[i]) {
257 fstart = i;
258 while (!mask[i] && i < mask.size()) {
259 fend = i;
260 i++;
261 }
262 }
263
264 // execute interpolation as the following criteria:
265 // (1) for a masked region inside the spectrum, replace the spectral
266 // values with the mean of those at the two channels just outside
267 // the both edges of the masked region.
268 // (2) for a masked region at the spectral edge, replace the values
269 // with the one at the nearest non-masked channel.
270 // (ZOH, but bilateral)
271 Float interp = 0.0;
272 if (fstart-1 > 0) {
273 interp = data[fstart-1];
274 if (fend+1 < Int(data.nelements())) {
275 interp = (interp + data[fend+1]) / 2.0;
276 }
277 } else {
278 interp = data[fend+1];
279 }
280 if (fstart > -1 && fend > -1) {
281 for (int j = fstart; j <= fend; ++j) {
282 data[j] = interp;
283 }
284 }
285
286 fstart = -1;
287 fend = -1;
288 }
289}
Note: See TracBrowser for help on using the repository browser.