Line data Source code
1 : // -*- C++ -*-
2 : //# ConvFuncDiskCache.cc: Implementation of the ConvFuncDiskCache class
3 : //# Copyright (C) 1997,1998,1999,2000,2001,2002,2003
4 : //# Associated Universities, Inc. Washington DC, USA.
5 : //#
6 : //# This library is free software; you can redistribute it and/or modify it
7 : //# under the terms of the GNU Library General Public License as published by
8 : //# the Free Software Foundation; either version 2 of the License, or (at your
9 : //# option) any later version.
10 : //#
11 : //# This library is distributed in the hope that it will be useful, but WITHOUT
12 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
14 : //# License for more details.
15 : //#
16 : //# You should have received a copy of the GNU Library General Public License
17 : //# along with this library; if not, write to the Free Software Foundation,
18 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
19 : //#
20 : //# Correspondence concerning AIPS++ should be addressed as follows:
21 : //# Internet email: aips2-request@nrao.edu.
22 : //# Postal address: AIPS++ Project Office
23 : //# National Radio Astronomy Observatory
24 : //# 520 Edgemont Road
25 : //# Charlottesville, VA 22903-2475 USA
26 : //#
27 : //# $Id$
28 : #include <synthesis/MeasurementComponents/ConvFuncDiskCache.h>
29 : #include <synthesis/TransformMachines/Utils.h>
30 : #include <casacore/casa/Exceptions/Error.h>
31 : #include <casacore/lattices/LEL/LatticeExpr.h>
32 : #include <synthesis/TransformMachines/SynthesisError.h>
33 : #include <synthesis/TransformMachines/Utils.h>
34 : #include <casacore/casa/OS/Directory.h>
35 : #include <fstream>
36 :
37 : using namespace casacore;
38 : namespace casa{
39 : //
40 : //-------------------------------------------------------------------------
41 : // Just load the axillary info. if found. The actual functions are loaded
42 : // as and when required using loadConvFunction() method.
43 : //
44 0 : void ConvFuncDiskCache::initCache()
45 : {
46 0 : ostringstream name;
47 0 : String line;
48 0 : Directory dirObj(Dir);
49 :
50 0 : if (Dir.length() == 0)
51 : throw(SynthesisFTMachineError("Got null string for disk cache dir. "
52 0 : "in ConvFuncDiskCache::initCache()"));
53 : //
54 : // If the directory does not exist, create it
55 : //
56 0 : if (!dirObj.exists()) dirObj.create();
57 0 : else if ((!dirObj.isWritable()) || (!dirObj.isReadable()))
58 : {
59 0 : throw(SynthesisFTMachineError(String("Directory \"")+Dir+String("\"")+
60 0 : String(" for convolution function cache"
61 0 : " exists but is unreadable/unwriteable")));
62 : }
63 :
64 : try
65 : {
66 0 : name << Dir << "/" << aux;
67 0 : File file(name);
68 0 : Int Npa=0,Nw=0;
69 0 : ifstream aux;
70 0 : Bool readFromFile=false;
71 0 : if (file.exists() && file.isRegular())
72 : {
73 0 : readFromFile=true;
74 0 : aux.open(name.str().c_str());
75 0 : if (readFromFile && aux.good()) aux >> Npa >> Nw;
76 : else
77 0 : throw(SynthesisFTMachineError(string("Error while reading convolution function cache file ")+name.str( )));
78 : }
79 :
80 0 : if (Npa > 0)
81 : {
82 0 : paList.resize(Npa,true);
83 :
84 0 : IPosition s(3,Nw,1,Npa);
85 0 : XSup.resize(s,true);
86 0 : YSup.resize(s,true);
87 0 : Sampling.resize(Npa,true);
88 0 : for(Int i=0;i<Npa;i++)
89 : {
90 : Float pa, S;
91 : Int XS, YS;
92 0 : s[2]=i;
93 0 : aux >> pa;
94 0 : for(Int iw=0;iw<Nw;iw++)
95 : {
96 0 : s[0]=iw;
97 0 : aux >> XS >> YS;
98 0 : YS = XS;
99 0 : paList[i] = pa*M_PI/180.0;
100 0 : XSup(iw,0,i)=XS;
101 0 : YSup(iw,0,i)=YS;
102 : }
103 0 : aux >> S;
104 0 : Sampling[i]=S;
105 : }
106 : }
107 : }
108 0 : catch(AipsError& x)
109 : {
110 0 : throw(SynthesisFTMachineError(String("Error while initializing CF disk cache: ")
111 0 : +x.getMesg()));
112 : }
113 0 : }
114 0 : ConvFuncDiskCache& ConvFuncDiskCache::operator=(const ConvFuncDiskCache& other)
115 : {
116 0 : paList = other.paList;
117 0 : Sampling = other.Sampling;
118 0 : XSup = other.XSup;
119 0 : YSup = other.YSup;
120 0 : Dir = other.Dir;
121 0 : cfPrefix = other.cfPrefix;
122 0 : aux = other.aux;
123 0 : return *this;
124 : };
125 : //
126 : //-------------------------------------------------------------------------
127 : // Write the conv. functions from the mem. cache to the disk cache.
128 : //
129 0 : void ConvFuncDiskCache::cacheConvFunction(Int which, Float pa,
130 : Array<Complex>& cf,
131 : CoordinateSystem& coords,
132 : CoordinateSystem& ftCoords,
133 : Int &convSize,
134 : Cube<Int> &convSupport,
135 : Float convSampling,
136 : String nameQualifier,
137 : Bool savePA)
138 : {
139 0 : Int N=paList.nelements();
140 0 : if (Dir.length() == 0) return;
141 :
142 : try
143 : {
144 0 : IPosition newConvShape = cf.shape();
145 0 : Int wConvSize = newConvShape(2), directionIndex;
146 0 : for(Int iw=0;iw<wConvSize;iw++)
147 : {
148 0 : IPosition sliceStart(4,0,0,iw,0),
149 0 : sliceLength(4,newConvShape(0),newConvShape(1),1,newConvShape(3));
150 :
151 : // CoordinateSystem ftCoords(coords);
152 0 : ftCoords = coords;
153 0 : directionIndex=ftCoords.findCoordinate(Coordinate::DIRECTION);
154 0 : DirectionCoordinate dc=coords.directionCoordinate(directionIndex);
155 : // AlwaysAssert(directionIndex>=0, AipsError);
156 0 : dc=coords.directionCoordinate(directionIndex);
157 0 : Vector<Bool> axes(2); axes(0)=axes(1)=true;//axes(2)=true;
158 0 : Vector<Int> shape(2); shape(0)=newConvShape(0);shape(1)=newConvShape(1);
159 0 : shape=convSize;
160 0 : Vector<Double>ref(4);
161 0 : ref(0)=ref(1)=ref(2)=ref(3)=0;
162 0 : dc.setReferencePixel(ref);
163 0 : Coordinate* ftdc=dc.makeFourierCoordinate(axes,shape);
164 0 : Vector<Double> refVal;
165 0 : refVal=ftdc->referenceValue();
166 0 : refVal(0)=refVal(1)=0;
167 0 : ftdc->setReferenceValue(refVal);
168 0 : ref(0)=newConvShape(0)/2-1;
169 0 : ref(1)=newConvShape(1)/2-1;
170 0 : ftdc->setReferencePixel(ref);
171 :
172 : // cout << ref << endl << refVal << endl << shape << endl;
173 : // cout << dc.increment() << " " << ftdc->increment() << endl;
174 0 : ftCoords.replaceCoordinate(*ftdc, directionIndex);
175 : {
176 0 : ostringstream name;
177 0 : name << Dir << "/" << cfPrefix << nameQualifier << iw << "_" << which;
178 :
179 0 : IPosition screenShape(4,newConvShape(0),newConvShape(1),newConvShape(3),1);
180 :
181 0 : PagedImage<Complex> thisScreen(screenShape, ftCoords, name);
182 :
183 0 : Array<Complex> buf;
184 0 : buf=((cf(Slicer(sliceStart,sliceLength)).nonDegenerate()));
185 0 : thisScreen.put(buf);
186 : }
187 0 : delete ftdc; ftdc=0;
188 : }
189 0 : if (savePA)
190 : {
191 0 : IPosition s(3,wConvSize,1,N+1);
192 0 : paList.resize(N+1,true);
193 : // XSup.resize(N+1,true);
194 : // YSup.resize(N+1,true);
195 0 : XSup.resize(s,true);
196 0 : YSup.resize(s,true);
197 0 : Sampling.resize(N+1,true);
198 0 : paList[N] = pa;
199 0 : for(Int iw=0;iw<wConvSize;iw++)
200 : {
201 0 : YSup(iw,0,N) = convSupport(iw,0,which);
202 0 : XSup(iw,0,N) = convSupport(iw,0,which);
203 : }
204 0 : Sampling[N]=convSampling;
205 : }
206 : }
207 0 : catch (AipsError& x)
208 : {
209 : throw(SynthesisFTMachineError("Error while caching CF to disk in "
210 : "ConvFuncDiskCache::cacheConvFunction(): "
211 0 : +x.getMesg()));
212 : }
213 : }
214 : //
215 : //-------------------------------------------------------------------------
216 : // Write the weight functions from the mem. cache to the disk cache.
217 : //
218 0 : void ConvFuncDiskCache::cacheWeightsFunction(Int which, Float pa,
219 : Array<Complex>& cfWt,
220 : CoordinateSystem& coords,
221 : Int &convSize,
222 : Cube<Int> &convSupport,
223 : Float convSampling)
224 : {
225 0 : Int N=paList.nelements();
226 0 : if (Dir.length() == 0) return;
227 :
228 : try
229 : {
230 0 : IPosition newConvShape = cfWt.shape();
231 0 : Int wConvSize = newConvShape(2), directionIndex;
232 0 : for(Int iw=0;iw<wConvSize;iw++)
233 : {
234 0 : IPosition sliceStart(4,0,0,iw,0),
235 0 : sliceLength(4,newConvShape(0),newConvShape(1),1,newConvShape(3));
236 :
237 0 : CoordinateSystem ftCoords(coords);
238 0 : directionIndex=ftCoords.findCoordinate(Coordinate::DIRECTION);
239 0 : DirectionCoordinate dc=coords.directionCoordinate(directionIndex);
240 : // AlwaysAssert(directionIndex>=0, AipsError);
241 0 : dc=coords.directionCoordinate(directionIndex);
242 0 : Vector<Bool> axes(2); axes(0)=axes(1)=true;//axes(2)=true;
243 0 : Vector<Int> shape(2); shape(0)=newConvShape(0);shape(1)=newConvShape(1);
244 0 : shape=convSize;
245 0 : Vector<Double>ref(4);
246 0 : ref(0)=ref(1)=ref(2)=ref(3)=0;
247 0 : dc.setReferencePixel(ref);
248 0 : Coordinate* ftdc=dc.makeFourierCoordinate(axes,shape);
249 0 : Vector<Double> refVal;
250 0 : refVal=ftdc->referenceValue();
251 0 : refVal(0)=refVal(1)=0;
252 0 : ftdc->setReferenceValue(refVal);
253 0 : ref(0)=newConvShape(0)/2-1;
254 0 : ref(1)=newConvShape(1)/2-1;
255 0 : ftdc->setReferencePixel(ref);
256 :
257 : // cout << ref << endl << refVal << endl << shape << endl;
258 : // cout << dc.increment() << " " << ftdc->increment() << endl;
259 0 : ftCoords.replaceCoordinate(*ftdc, directionIndex);
260 0 : delete ftdc; ftdc=0;
261 :
262 : {
263 0 : ostringstream name;
264 0 : name << Dir << "/" << cfPrefix << "WT" << iw << "_" << which;
265 :
266 0 : IPosition screenShape(4,newConvShape(0),newConvShape(1),newConvShape(3),1);
267 :
268 0 : PagedImage<Complex> thisScreen(screenShape, ftCoords, name);
269 :
270 0 : Array<Complex> buf;
271 0 : buf=((cfWt(Slicer(sliceStart,sliceLength)).nonDegenerate()));
272 0 : thisScreen.put(buf);
273 : }
274 : }
275 0 : IPosition s(3,wConvSize,1,N+1);
276 0 : paList.resize(N+1,true);
277 : // XSup.resize(N+1,true);
278 : // YSup.resize(N+1,true);
279 0 : XSup.resize(s,true);
280 0 : YSup.resize(s,true);
281 0 : Sampling.resize(N+1,true);
282 0 : paList[N] = pa;
283 0 : for(Int iw=0;iw<wConvSize;iw++)
284 : {
285 0 : YSup(iw,0,N) = convSupport(iw,0,which);
286 0 : XSup(iw,0,N) = convSupport(iw,0,which);
287 : }
288 0 : Sampling[N]=convSampling;
289 : }
290 0 : catch (AipsError& x)
291 : {
292 : throw(SynthesisFTMachineError("Error while caching CFWT to disk in "
293 : "ConvFuncDiskCache::cacheConvFunction(): "
294 0 : +x.getMesg()));
295 : }
296 : }
297 : //
298 : //-------------------------------------------------------------------------
299 : //
300 0 : Bool ConvFuncDiskCache::searchConvFunction(const VisBuffer& vb,
301 : VPSkyJones& vpSJ,
302 : Int& which,
303 : Float &pa)
304 : {
305 0 : Int i,NPA=paList.nelements(); Bool paFound=false;
306 : Float iPA, dPA;
307 0 : dPA = vpSJ.getPAIncrement().getValue("rad");
308 : /*
309 : Vector<Float> antPA = vb.feed_pa(getCurrentTimeStamp(vb));
310 : pa = sum(antPA)/(antPA.nelements()-1);
311 : */
312 0 : pa = getPA(vb);
313 : // cout << dPA*57.295 << " " << pa*57.295 << endl;
314 : // pa = 0;
315 : // cout << "######CFDC::search: " << pa << " " << getPA(vb) << endl;
316 : // if (NPA == 0) return -1;
317 :
318 0 : for(i=0;i<NPA;i++)
319 : {
320 0 : iPA = paList[i];
321 0 : if (fabs(iPA - pa) <= dPA)
322 : {
323 0 : paFound = true;
324 0 : break;
325 : }
326 : }
327 0 : if (paFound) which = i; else which = -i;
328 0 : return paFound;
329 : }
330 : //
331 : //-------------------------------------------------------------------------
332 : //
333 0 : Bool ConvFuncDiskCache::searchConvFunction(const VisBuffer& vb,
334 : ParAngleChangeDetector& vpSJ,
335 : Int& which,
336 : Float &pa)
337 : {
338 0 : if (paList.nelements()==0) initCache();
339 0 : Int i,NPA=paList.nelements(); Bool paFound=false;
340 : Float iPA, dPA;
341 0 : dPA = vpSJ.getParAngleTolerance().getValue("rad");
342 : /*
343 : Vector<Float> antPA = vb.feed_pa(getCurrentTimeStamp(vb));
344 : pa = sum(antPA)/(antPA.nelements()-1);
345 : */
346 0 : pa = getPA(vb);
347 : // pa = 0;
348 : // cout << "######CFDC::search: " << pa << " " << getPA(vb) << endl;
349 : // if (NPA == 0) return -1;
350 :
351 0 : Float paDiff=2*dPA;
352 0 : Int saveNdx=-1;
353 :
354 0 : saveNdx = -1;
355 0 : for(i=0;i<NPA;i++)
356 : {
357 0 : iPA = paList[i];
358 0 : if (fabs(iPA - pa) < paDiff)
359 : {
360 0 : saveNdx = i;
361 0 : paDiff = fabs(iPA-pa);
362 : }
363 : }
364 0 : if (saveNdx > -1)
365 : {
366 0 : iPA = paList[saveNdx];
367 0 : if (fabs(iPA - pa) <= dPA)
368 : {
369 0 : i = saveNdx;
370 0 : paFound=true;
371 : }
372 : }
373 : /*
374 : for(i=0;i<NPA;i++)
375 : {
376 : iPA = paList[i];
377 : if (fabs(iPA - pa) <= dPA)
378 : {
379 : paFound = true;
380 : break;
381 : }
382 : }
383 : */
384 0 : if (paFound)
385 : {
386 : // cout << "SEARCH: " << i << " " << paList[i]*180/M_PI << " " << pa*180/M_PI
387 : // << " " << (paList[i]-pa)*180/M_PI << endl;
388 0 : which = i;
389 : }
390 0 : else which = -i;
391 0 : return paFound;
392 : }
393 : //
394 : //-------------------------------------------------------------------------
395 : //Write the aux. info. also in the disk cache (wonder if this should
396 : //be automatically be called from cacheConvFunction() method).
397 : //
398 0 : void ConvFuncDiskCache::finalize()
399 : {
400 0 : if (Dir.length() == 0) return;
401 0 : ostringstream name;
402 0 : name << Dir << "/aux.dat";
403 : try
404 : {
405 : // cout << "Writing to " << name.str() << endl;
406 0 : IPosition supportShape=XSup.shape();
407 0 : ofstream aux(name.str().c_str());
408 0 : aux << paList.nelements() << " " << supportShape[0] << endl;
409 0 : for(uInt ipa=0;ipa<paList.nelements();ipa++)
410 : {
411 0 : aux << paList[ipa]*180.0/M_PI << " ";
412 0 : for(int iw=0;iw<supportShape[0];iw++)
413 0 : aux << XSup(iw,0,ipa) << " " << YSup(iw,0,ipa) << " ";
414 0 : aux << " " << Sampling[ipa] <<endl;
415 : }
416 : }
417 0 : catch(AipsError &x)
418 : {
419 0 : throw(SynthesisFTMachineError(string("Error while writing ")
420 0 : + name.str( ) + (string) x.getMesg()));
421 : }
422 : }
423 : //
424 : //-------------------------------------------------------------------------
425 : //Along with the aux. info., also save the average PB in the disk cache.
426 : //
427 0 : void ConvFuncDiskCache::finalize(ImageInterface<Float>& avgPB)
428 : {
429 0 : if (Dir.length() == 0) return;
430 0 : finalize();
431 0 : ostringstream Name;
432 0 : Name << Dir <<"/avgPB";
433 : try
434 : {
435 0 : IPosition newShape(avgPB.shape());
436 0 : PagedImage<Float> tmp(newShape, avgPB.coordinates(), Name);
437 :
438 0 : LatticeExpr<Float> le(avgPB);
439 0 : tmp.copyData(le);
440 : }
441 0 : catch(AipsError &x)
442 : {
443 0 : throw(SynthesisFTMachineError(string("Error while writing ")
444 0 : + Name.str( ) + (string) x.getMesg()));
445 : }
446 : }
447 : //
448 : //-------------------------------------------------------------------------
449 : //Load the average PB from the disk cache.
450 : //
451 0 : void ConvFuncDiskCache::loadAvgPB(ImageInterface<Float>& avgPB)
452 : {
453 0 : if (Dir.length() == 0) return;
454 0 : ostringstream name;
455 0 : name << Dir << "/avgPB";
456 : // cout << name.str() << endl;
457 : try
458 : {
459 0 : PagedImage<Float> tmp(name.str().c_str());
460 0 : avgPB.resize(tmp.shape());
461 0 : avgPB.put(tmp.get());
462 : }
463 0 : catch(AipsError& x) // Just rethrowing the exception for now.
464 : // Ultimately, this should be used to make
465 : // the state of this object consistant.
466 : {
467 0 : throw(SynthesisFTMachineError(string("Error while loading \"")+
468 0 : name.str( ) + string("\": ") + (string) x.getMesg()));
469 : }
470 : }
471 : //
472 : //-------------------------------------------------------------------------
473 : //Load a conv. func. from the disk. This is non-optimal due to the
474 : //data structure used for the conv. func. in-memory cache (it's an
475 : //array of pointers where it should really be a List of pointers).
476 : //The conf. func. index, which is also used as a key to located them
477 : //in the mem. cache, are not assured to be contiguous. As a result,
478 : //in the current implementation there can be gaps in the
479 : //convFuncCache array. These gaps are initialized to NULL pointers.
480 : //It's not much of a memory waste, but still non-optimal! Leaving
481 : //it like this for now.
482 : //
483 0 : Bool ConvFuncDiskCache::loadConvFunction(Int where, Int Nw,
484 : PtrBlock < Array<Complex> *> &convFuncCache,
485 : Cube<Int> &convSupport,
486 : Vector<Float>& convSampling,
487 : Double& cfRefFreq, CoordinateSystem& coordSys,
488 : String prefix)
489 : {
490 0 : if (Dir.length() == 0) return false;
491 0 : if (where < (Int)convFuncCache.nelements() && (convFuncCache[where] != NULL)) return false;
492 :
493 0 : Int wConvSize, polInUse=2;
494 0 : Int N=convFuncCache.nelements();
495 :
496 : //
497 : // Re-size the conv. func. memory cache if required, and set the
498 : // new members of the resized cache to NULL. This is used in the
499 : // loop below to make a decision about allocating new memory or
500 : // not.
501 : //
502 0 : convFuncCache.resize(max(where+1,N), true);
503 0 : for(Int i=N;i<=where;i++) convFuncCache[i]=NULL;
504 : //
505 : // Each w-plan is in a separate disk file. Each file contains all
506 : // polarization planes though. Memory cache holds all w-planes and
507 : // poln-planes in a single complex array. The loop below read
508 : // each w-plane image from the disk, and fills in the 3D
509 : // mem. cache for each computed PA.
510 : //
511 0 : wConvSize = Nw;
512 0 : for(Int iw=0;iw<Nw;iw++)
513 : {
514 0 : ostringstream name;
515 : // name << Dir << "/CF" << iw << "_" << where;
516 0 : name << Dir << prefix << iw << "_" << where;
517 : // cout << "CF File name = " << name.str() << endl;
518 : try
519 : {
520 0 : PagedImage<Complex> tmp(name.str().c_str());
521 0 : Int index= tmp.coordinates().findCoordinate(Coordinate::SPECTRAL);
522 0 : coordSys = tmp.coordinates();
523 0 : SpectralCoordinate spCS = coordSys.spectralCoordinate(index);
524 :
525 0 : cfRefFreq=spCS.referenceValue()(0);
526 :
527 0 : polInUse = tmp.shape()(2);
528 0 : IPosition ts=tmp.shape(),ndx(4,0,0,0,0),ts2(4,0,0,0,0);
529 0 : Array<Complex> buf=tmp.get();
530 0 : if (convFuncCache[where] == NULL)
531 0 : convFuncCache[where] = new Array<Complex>(IPosition(4,ts(0),ts(1),
532 0 : wConvSize,polInUse));
533 :
534 0 : ndx(2)=iw; // The w-axis
535 0 : for(ndx(3)=0;ndx(3)<polInUse;ndx(3)++) // The Poln. axis.
536 0 : for(ndx(0)=0;ndx(0)<ts(0);ndx(0)++)
537 0 : for(ndx(1)=0;ndx(1)<ts(1);ndx(1)++)
538 : {
539 0 : ts2(0)=ndx(0);ts2(1)=ndx(1);
540 0 : ts2(2)=ndx(3); // The Poln. axis of the disk-cache
541 0 : ts2(3)=0; // The freq. axis of the disk-cache
542 0 : (*convFuncCache[where])(ndx)=buf(ts2);
543 : }
544 : }
545 0 : catch(AipsError &x)
546 : {
547 0 : throw(SynthesisFTMachineError(string("Error while loading \"")+
548 0 : name.str( ) + string("\": ") + (string) x.getMesg()));
549 : }
550 :
551 : }
552 0 : convSupport.resize(wConvSize,polInUse,where+1,true);
553 0 : for(Int i=0;i<wConvSize;i++)
554 0 : for(Int j=0;j<polInUse;j++)
555 0 : convSupport(i,j,where) = XSup(i,0,where);
556 : // cout << "##### " << convFuncCache.nelements() << endl;
557 0 : convSampling = Sampling;
558 0 : return true;
559 : }
560 : }
|
