source: branches/alma/src/MathUtils.cpp@ 1535

Last change on this file since 1535 was 1530, checked in by Kana Sugimoto, 16 years ago

New Development: No

JIRA Issue: Yes (CAS-1079)

Ready to Release: Yes

Interface Changes: Yes for scantable.stats()

What Interface Changed:
Return value is now a list of abcissa values
when stat='min_abc' and 'max_abc'.

Test Programs:

run scantable.stats() with the parameter stat='min_abc' or 'max_abc'.

Put in Release Notes: No

Module(s): sdstat

Description:

scantable.stats() are modified so that the tool returns abscissa
chan/freq/velocity of min or max when stat='min_abc' or 'max_abc',
respectively.
Minor changes in IPosition mathutil::minMaxPos() to accept both
min (max) and min_abc (max_abc) as a parameter.


  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 6.4 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/BasicSL/String.h>
38#include <scimath/Mathematics/MedianSlider.h>
39#include <casa/Exceptions/Error.h>
40
41#include "MathUtils.h"
42
43using namespace casa;
44
45float mathutil::statistics(const String& which,
46 const MaskedArray<Float>& data)
47{
48 String str(which);
49 str.upcase();
50 if (str.matches(String("MIN"))) {
51 return min(data);
52 } else if (str.matches(String("MAX"))) {
53 return max(data);
54 } else if (str.matches(String("SUMSQ"))) {
55 return sumsquares(data);
56 } else if (str.matches(String("SUM"))) {
57 return sum(data);
58 } else if (str.matches(String("MEAN"))) {
59 return mean(data);
60 } else if (str.matches(String("VAR"))) {
61 return variance(data);
62 } else if (str.matches(String("STDDEV"))) {
63 return stddev(data);
64 } else if (str.matches(String("AVDEV"))) {
65 return avdev(data);
66 } else if (str.matches(String("RMS"))) {
67 uInt n = data.nelementsValid();
68 return sqrt(sumsquares(data)/n);
69 } else if (str.matches(String("MEDIAN"))) {
70 return median(data);
71 } else {
72 String msg = str + " is not a valid type of statistics";
73 throw(AipsError(msg));
74 }
75 return 0.0;
76}
77
78IPosition mathutil::minMaxPos(const String& which,
79 const MaskedArray<Float>& data)
80{
81 Float minVal, maxVal;
82 IPosition minPos(data.ndim(), 0), maxPos(data.ndim(), 0);
83 minMax(minVal, maxVal, minPos, maxPos, data);
84 String str(which);
85 str.upcase();
86 if (str.contains(String("MIN"))) {
87 return minPos;
88 } else if (str.contains(String("MAX"))) {
89 return maxPos;
90 } else {
91 String msg = str + " is not a valid type of statistics";
92 throw(AipsError(msg));
93 }
94 //return 0.0;
95}
96
97void mathutil::replaceMaskByZero(Vector<Float>& data, const Vector<Bool>& mask)
98{
99 for (uInt i=0; i<data.nelements(); i++) {
100 if (!mask[i]) data[i] = 0.0;
101 }
102}
103
104
105std::vector<std::string> mathutil::tovectorstring(const Vector<String>& in)
106{
107 std::vector<std::string> out;
108 Vector<String>::const_iterator it = in.begin();
109 for (uInt i=0; it != in.end(); ++it,++i) {
110 out.push_back(*it);
111 }
112 return out;
113}
114
115Vector<String> mathutil::toVectorString(const std::vector<std::string>& in)
116{
117 Vector<String> out(in.size());
118 uInt i=0;
119 std::vector<std::string>::const_iterator it;
120 for (it=in.begin();it != in.end();++it,++i) {
121 out[i] = casa::String(*it);
122 }
123 return out;
124}
125
126void mathutil::hanning(Vector<Float>& out, Vector<Bool>& outmask,
127 const Vector<Float>& in, const Vector<Bool>& mask,
128 Bool relaxed, Bool ignoreOther) {
129 Vector< Vector<Float> > weights(8);
130 Vector<Float> vals(3);
131 vals = 0.0;weights[0] = vals;// FFF
132 vals[0] = 1.0; vals[1] = 0.0; vals[2] = 0.0; weights[1] = vals;// TFF
133 vals[0] = 0.0; vals[1] = 1.0; vals[2] = 0.0; weights[2] = vals;// FTF
134 vals[0] = 1.0/3.0; vals[1] = 2.0/3.0; vals[2] = 0.0; weights[3] = vals;// TTF
135 vals[0] = 0.0; vals[1] = 0.0; vals[2] = 1.0;weights[4] = vals;// FFT
136 vals[0] = 0.5; vals[1] = 0.0; vals[2] = 0.5; weights[5] = vals;// TFT
137 vals[0] = 0.0; vals[1] = 2.0/3.0; vals[2] = 1.0/3.0; weights[6] = vals;// FTT
138 vals[0] = 0.25; vals[1] = 0.5; vals[2] = 0.25; weights[7] = vals;// TTT
139 // Chris' case
140 Vector<Bool> weighted(8);
141 if (relaxed) {
142 weighted = False;
143 weighted[7] = True;
144
145 } else {
146 weighted = True;
147 weighted[0] = False;
148 }
149
150 out.resize(in.nelements());
151 outmask.resize(mask.nelements());
152 // make special case for first and last
153 /// ...here
154 // loop from 1..n-2
155 out.resize(in.nelements());
156 out[0] = in[0];out[out.nelements()-1] = in[in.nelements()-1];
157 outmask.resize(mask.nelements());
158 outmask = False;
159 uInt m;Vector<Float>* w;
160 for (uInt i=1; i < out.nelements()-1;++i) {
161 m = mask[i-1] + 2*mask[i] + 4*mask[i+1];
162 w = &(weights[m]);
163 if (weighted[m]) {
164 out[i] = (*w)[0]*in[i-1] + (*w)[1]*in[i] + (*w)[2]*in[i+1];
165 outmask[i] = True;
166 } else { // mask it
167 out[i] = in[i];//use arbitrary value
168 outmask[i] = False;
169 }
170 }
171}
172
173
174void mathutil::runningMedian(Vector<Float>& out, Vector<Bool>& outflag,
175 const Vector<Float>& in, const Vector<Bool>& flag,
176 float width)
177{
178 Int hwidth = Int(width+0.5);
179 Int fwidth = hwidth*2+1;
180 out.resize(in.nelements());
181 outflag.resize(flag.nelements());
182 MedianSlider ms(hwidth);
183 Slice sl(0, fwidth-1);
184 Float medval = ms.add(const_cast<Vector<Float>& >(in)(sl),
185 const_cast<Vector<Bool>& >(flag)(sl));
186 uInt n = in.nelements();
187 for (uInt i=hwidth; i<(n-hwidth); ++i) {
188 // add data value
189 out[i] = ms.add(in[i+hwidth], flag[i+hwidth]);
190 outflag[i] = (ms.nval() == 0);
191 }
192 // replicate edge values from fisrt value with full width of values
193 for (uInt i=0;i<hwidth;++i) {
194 out[i] = out[hwidth];
195 outflag[i] = outflag[hwidth];
196 out[n-1-i] = out[n-1-hwidth];
197 outflag[n-1-i] = outflag[n-1-hwidth];
198 }
199}
Note: See TracBrowser for help on using the repository browser.