Line data Source code
1 : //# CTTimeInterp1.cc: Implementation of CTTimeInterp1.h
2 : //# Copyright (C) 2012
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be adressed as follows:
20 : //# Internet email: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //#
27 :
28 : #include <synthesis/CalTables/CTTimeInterp1.h>
29 : #include <synthesis/CalTables/CTMainColumns.h>
30 : #include <synthesis/CalTables/RIorAParray.h>
31 : #include <casacore/casa/aips.h>
32 : #include <casacore/casa/BasicSL/Constants.h>
33 : #include <casacore/casa/Arrays/Array.h>
34 : #include <casacore/casa/IO/ArrayIO.h>
35 : #include <casacore/casa/Arrays/ArrayMath.h>
36 : #include <casacore/casa/Logging/LogMessage.h>
37 : #include <casacore/casa/Logging/LogSink.h>
38 : #include <casacore/scimath/Functionals/Interpolate1D.h>
39 : #include <casacore/scimath/Functionals/ScalarSampledFunctional.h>
40 : #include <casacore/scimath/Functionals/ArraySampledFunctional.h>
41 :
42 : #define CTTIMEINTERPVERB1 false
43 :
44 : using namespace casacore;
45 : namespace casa { //# NAMESPACE CASA - BEGIN
46 :
47 : // Null ctor
48 : //CTTimeInterp1::CTTimeInterp1() {};
49 :
50 : // From NewCalTable
51 10286 : CTTimeInterp1::CTTimeInterp1(NewCalTable& ct,
52 : const String& timetype,
53 : Array<Float>& result,
54 10286 : Array<Bool>& rflag) :
55 10286 : ct_(ct),
56 10286 : mcols_p(NULL),
57 10286 : timeType_(timetype),
58 10286 : currTime_(-999.0),
59 10286 : currIdx_(-1),
60 10286 : lastWasExact_(false),
61 10286 : timeRef_(0.0),
62 10286 : timelist_(),
63 10286 : domain_(2,0.0),
64 10286 : flaglist_(),
65 10286 : tInterpolator_p(NULL),
66 10286 : cfreq_(-1.0),
67 10286 : cycles_(),
68 10286 : result_(),
69 20572 : rflag_()
70 : {
71 :
72 : if (CTTIMEINTERPVERB1) cout << "CTTimeInterp1::ctor()" << endl;
73 :
74 : // Access to main columns (move to NewCalTable?)
75 10286 : mcols_p = new ROCTMainColumns(ct_);
76 :
77 : // Flags
78 10286 : mcols_p->flag().getColumn(flaglist_);
79 :
80 : // Record referenced timelist_
81 : // TBD: handle flagged times
82 10286 : Vector<Double> dtime;
83 10286 : mcols_p->time().getColumn(dtime);
84 10286 : domain_(0)=min(dtime);
85 10286 : domain_(1)=max(dtime);
86 :
87 : // NB: time is rendered as Vector<Float> for processing
88 : // because Interpolate1D doesn't work if it is V<Double>...
89 : //timeRef_=86400.0*floor(dtime(0)/86400.0);
90 : // 2016Aug22 (gmoellen): use ~current time (not prior midnight)
91 : // as timeRef_ to ~minimize loss of precision for faster sampling
92 10286 : timeRef_=floor(dtime(0)-1.0f);
93 10286 : dtime-=timeRef_; // Relative to reference time
94 10286 : timelist_.resize(dtime.shape());
95 10286 : convertArray(timelist_,dtime); // store referenced times as Floats
96 :
97 : // Create the _time_ interpolator
98 : // TBD: f(par) (because flags may be different for each par...)
99 10286 : tInterpolator_p = new Interpolate1D<Float,Array<Float> >();
100 10287 : setInterpType(timeType_);
101 :
102 : // Get fiducial frequency
103 :
104 : // Get cycles, if nec.
105 10285 : if (timetype.contains("PD")) {
106 0 : Int spw=mcols_p->spwId()(0);
107 0 : MSSpWindowColumns spwcol(ct_.spectralWindow());
108 0 : Int nChan=spwcol.numChan()(spw);
109 0 : cfreq_=Vector<Double>(spwcol.chanFreq()(spw))(nChan/2);
110 : //cout << "cfreq_ = " << cfreq_ << endl;
111 0 : mcols_p->cycles(cycles_);
112 : //cout << "ant = " << mcols_p->antenna1()(0) << ": cycles_ = " << cycles_ << endl;
113 0 : }
114 :
115 : // Reference supplied arrays for results
116 : // Elsewhere, must always use strict (non-shape-changing) assignment for these!
117 : // TBD: verify shapes..
118 10285 : result_.reference(result);
119 10285 : rflag_.reference(rflag);
120 :
121 10294 : }
122 :
123 : // Destructor
124 20250 : CTTimeInterp1::~CTTimeInterp1() {
125 10285 : if (tInterpolator_p)
126 10285 : delete tInterpolator_p;
127 10285 : if (mcols_p)
128 10285 : delete mcols_p;
129 20250 : };
130 :
131 2826041 : Bool CTTimeInterp1::interpolate(Double newtime) {
132 :
133 : if (CTTIMEINTERPVERB1) {cout.precision(12);cout << "CTTimeInterp1::interpolate("<<newtime<<"):" << endl;}
134 :
135 : // Don't work unnecessarily
136 2826041 : if (newtime==currTime_)
137 0 : return false; // no change
138 :
139 : if (CTTIMEINTERPVERB1) cout << " newtime is new..." << endl;
140 :
141 : // A new time is specified, so some work may be required
142 :
143 : // Convert supplied time value to Float (referenced to timeRef_)
144 2826041 : Float fnewtime(newtime-timeRef_);
145 :
146 : // Establish registration in time
147 2826041 : Bool exact(false);
148 2826041 : Int newIdx(currIdx_);
149 2826041 : Bool newReg=findTimeRegistration(newIdx,exact,fnewtime);
150 :
151 : if (CTTIMEINTERPVERB1)
152 : cout <<boolalpha<< " newReg="<<newReg<< " newIdx="<<newIdx<< " exact="<<exact
153 : << " lastWasExact_=" << lastWasExact_ << endl;
154 :
155 : // Update interpolator coeffs if new registr. and not nearest
156 2826041 : if (newReg || (!exact && lastWasExact_)) {
157 : // Only bother if not 'nearest' nor exact...
158 601736 : if (!timeType().contains("nearest") && !exact) {
159 24171 : if (timeType().contains("linear")) {
160 21811 : ScalarSampledFunctional<Float> xf(timelist_(Slice(newIdx,2)));
161 21811 : Vector<uInt> rows(2); indgen(rows); rows+=uInt(newIdx);
162 21811 : Array<Float> ya(mcols_p->fparamArray("",rows));
163 21811 : ArraySampledFunctional<Array<Float> > yf(ya);
164 21811 : tInterpolator_p->setData(xf,yf,true);
165 24171 : } else if (timeType().contains("cubic")) { // Added for CAS-10787 (16/2/2018 WK)
166 2360 : Int newIdxCubic(newIdx-1);
167 2360 : if (newIdxCubic < 0) {
168 40 : newIdxCubic = 0;
169 2320 : } else if (newIdxCubic > (Int)timelist_.nelements()-4) {
170 0 : newIdxCubic = timelist_.nelements()-4;
171 : }
172 : //cout << "{newIdxCubic = " << newIdxCubic << " / " << timelist_.nelements() << "}" << flush;
173 2360 : ScalarSampledFunctional<Float> xf(timelist_(Slice(newIdxCubic,4)));
174 2360 : Vector<uInt> rows(4); indgen(rows); rows+=uInt(newIdxCubic);
175 2360 : Array<Float> ya(mcols_p->fparamArray("",rows));
176 2360 : ArraySampledFunctional<Array<Float> > yf(ya);
177 2360 : tInterpolator_p->setData(xf,yf,true);
178 2360 : }
179 : }
180 601736 : }
181 : else
182 : // Escape if registration unchanged and 'nearest' or exact
183 2224305 : if (timeType().contains("nearest") || exact) return false; // no change
184 :
185 : // Now calculate the interpolation result
186 :
187 834116 : if (timeType().contains("nearest") || exact) {
188 : if (CTTIMEINTERPVERB1) cout << " nearest or exact" << endl;
189 1155130 : Cube<Float> r(mcols_p->fparamArray("",Vector<uInt>(1,newIdx)));
190 577565 : result_=r.xyPlane(0);
191 577565 : rflag_=flaglist_.xyPlane(newIdx);
192 577565 : }
193 : else {
194 : if (CTTIMEINTERPVERB1) cout << " non-trivial non-nearest" << endl;
195 : // Delegate to the interpolator
196 : // The next line is needed to restore the given interpolation type, which has been overwritten to linear in setData() above - CAS-10787 (16/2/2018 WK)
197 256551 : setInterpType(timeType());
198 256551 : result_=(*tInterpolator_p)(fnewtime);
199 256551 : rflag_=(flaglist_.xyPlane(newIdx) || flaglist_.xyPlane(newIdx+1));
200 : }
201 :
202 : // Now remember for next round
203 834116 : currTime_=newtime;
204 834116 : currIdx_=newIdx;
205 834116 : lastWasExact_=exact;
206 :
207 834116 : return true;
208 :
209 : }
210 :
211 2715539 : Bool CTTimeInterp1::interpolate(Double newtime, Double freq) {
212 :
213 2715539 : Bool newcal=this->interpolate(newtime);
214 :
215 2715539 : if (newcal && timeType().contains("PD"))
216 0 : applyPhaseDelay(freq);
217 :
218 2715539 : return newcal;
219 :
220 : }
221 :
222 0 : void CTTimeInterp1::state(Bool verbose) {
223 :
224 0 : cout << endl << "-state---------" << endl;
225 0 : cout.precision(12);
226 0 : cout << " timeType_= " << timeType_ << endl;
227 0 : cout << " ntime = " << timelist_.nelements() << endl;
228 0 : cout << " currTime_= " << currTime_ << " (" << Float(currTime_-timeRef_) << ")" << endl;
229 0 : cout << " currIdx_ = " << currIdx_ << endl;
230 0 : cout << " timeRef_ = " << timeRef_ << endl;
231 0 : cout << " domain_ = " << domain_ << endl;
232 0 : if (verbose) {
233 0 : cout << " result_ = " << result_ << endl;
234 0 : cout << boolalpha;
235 0 : cout << " rflag_ = " << rflag_ << endl;
236 : }
237 0 : cout << "---------------" << endl;
238 0 : }
239 :
240 266837 : void CTTimeInterp1::setInterpType(String strtype) {
241 266837 : timeType_=strtype;
242 266837 : currTime_=-999.0; // ensure force refreshed calculation
243 266837 : if (strtype.contains("nearest")) {
244 1529 : tInterpolator_p->setMethod(Interpolate1D<Float,Array<Float> >::nearestNeighbour);
245 1529 : return;
246 : }
247 265308 : if (strtype.contains("linear")) {
248 262747 : tInterpolator_p->setMethod(Interpolate1D<Float,Array<Float> >::linear);
249 262747 : return;
250 : }
251 2561 : if (strtype.contains("cubic")) {
252 2560 : tInterpolator_p->setMethod(Interpolate1D<Float,Array<Float> >::cubic);
253 2560 : return;
254 : }
255 1 : if (strtype.contains("spline")) {
256 0 : tInterpolator_p->setMethod(Interpolate1D<Float,Array<Float> >::spline);
257 0 : return;
258 : }
259 1 : throw(AipsError("Unknown interp type: '"+strtype+"'!!"));
260 : }
261 :
262 2826041 : Bool CTTimeInterp1::findTimeRegistration(Int& idx,Bool& exact,Float newtime) {
263 :
264 : if (CTTIMEINTERPVERB1) cout << " CTTimeInterp1::findTimeRegistration()" << endl;
265 :
266 2826041 : Int ntime=timelist_.nelements();
267 :
268 : // If only one time in timelist, that's it, and its exact (behave as nearest)
269 2826041 : if (ntime==1) {
270 1943054 : idx=0;
271 1943054 : exact=true;
272 : }
273 : else {
274 : // More than one element in timelist, find the right one:
275 :
276 : // Behave as nearest outside absolute range of available calibration
277 : // to avoid wild extrapolation, else search for the correct soln slot
278 882987 : if (newtime<timelist_(0)) {
279 : // Before first timestamp, use first slot as nearest one
280 31125 : idx=0;
281 31125 : exact=true;
282 : }
283 851862 : else if (newtime>timelist_(ntime-1)) {
284 : // After last timestamp, use last slot as nearest one
285 20892 : idx=ntime-1;
286 20892 : exact=true;
287 : }
288 : else
289 : // Find index in timelist where time would be:
290 : // TBD: can we use last result to help this?
291 830970 : idx=binarySearch(exact,timelist_,newtime,ntime,0);
292 :
293 : // If not (yet) an exact match...
294 882987 : if ( !exact ) {
295 :
296 : // identify this time via index just prior
297 258015 : if (idx>0) idx--;
298 :
299 : // If nearest, fine-tune slot to actual nearest:
300 258015 : if ( timeType().contains("nearest") ) {
301 : // exact=true; // Nearest behaves like exact match
302 2928 : if (idx!=(ntime-1) &&
303 1464 : (timelist_(idx+1)-newtime) < (newtime-timelist_(idx)) )
304 688 : idx++;
305 : } else {
306 : // linear modes require one later slot
307 256551 : if (idx==ntime-1) idx--;
308 : }
309 : }
310 :
311 : }
312 : // Return if new
313 2826041 : return (idx!=currIdx_);
314 :
315 : }
316 :
317 0 : void CTTimeInterp1::applyPhaseDelay(Double freq) {
318 :
319 0 : if (freq>0.0) {
320 0 : IPosition blc(2,1,0),trc(result_.shape()-1),stp(2,2,1);
321 0 : Matrix<Float> rph(result_(blc,trc,stp));
322 0 : if (cfreq_>0.0) {
323 0 : rph+=cycles_.xyPlane(currIdx_);
324 0 : rph*=Float(freq/cfreq_);
325 : }
326 0 : }
327 : else
328 0 : throw(AipsError("CTTimeInterp1::applyPhaseDelay: invalid frequency."));
329 :
330 0 : }
331 :
332 : } //# NAMESPACE CASA - END
|
