Line data Source code
1 : // -*- C++ -*-
2 : //# CFCache.cc: Implementation of the CFCache 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: casa-feedback@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/TransformMachines/SynthesisError.h>
29 : #include <synthesis/TransformMachines2/CFCache.h>
30 : #include <synthesis/TransformMachines2/Utils.h>
31 : #include <imageanalysis/Utilities/SpectralImageUtil.h>
32 : #include <casacore/lattices/LEL/LatticeExpr.h>
33 : #include <casacore/casa/System/ProgressMeter.h>
34 : #include <casacore/casa/Exceptions/Error.h>
35 : #include <casacore/casa/Utilities/Regex.h>
36 : #include <fstream>
37 : #include <algorithm>
38 : // #include <tables/Tables/TableDesc.h>
39 : // #include <tables/Tables/SetupNewTab.h>
40 : // #include <tables/Tables/Table.h>
41 :
42 : using namespace casacore;
43 : namespace casa{
44 : using namespace refim;
45 100 : CFCache::~CFCache()
46 : {
47 : //cerr << "#################" << "~CFCache() called" << endl;
48 100 : }
49 : //
50 : //-------------------------------------------------------------------------
51 : // Load just the axillary info. if found. The convolution functions
52 : // are loaded on-demand.
53 : //
54 0 : void CFCache::initCache()
55 : {
56 0 : LogOrigin logOrigin("CFCache2", "initCache");
57 0 : LogIO log_l(logOrigin);
58 :
59 0 : ostringstream name;
60 0 : String line;
61 0 : Directory dirObj(Dir);
62 :
63 0 : if (Dir.length() == 0)
64 0 : throw(SynthesisFTMachineError(LogMessage("Got null string for disk cache dir. ",
65 0 : logOrigin).message()));
66 :
67 0 : Bool ReadOnly = SynthesisUtils::getenv("CFCache.READONLY",0);
68 :
69 0 : if (!dirObj.exists()) dirObj.create();
70 0 : else if ((!dirObj.isReadable()))
71 : {
72 0 : throw(SynthesisFTMachineError(String("Directory \"")+Dir+String("\"")+
73 0 : String(" for convolution function cache"
74 0 : " exists but is unreadable")));
75 : }
76 0 : else if ((!dirObj.isWritable()) && (!ReadOnly))
77 : {
78 0 : throw(SynthesisFTMachineError(String("Directory \"")+Dir+String("\"")+
79 0 : String(" for convolution function cache"
80 0 : " exists but is unwriteable")));
81 : }
82 :
83 :
84 0 : if(ReadOnly)
85 : {
86 0 : log_l << "CF Cache in read only mode." << LogIO::POST;
87 : }
88 :
89 :
90 : try
91 : {
92 0 : name << Dir << "/" << aux;
93 0 : File file(name);
94 0 : Int Npa=0,Nw=0;
95 0 : ifstream aux;
96 0 : Bool readFromFile=false;
97 0 : if (file.exists() && file.isRegular())
98 : {
99 0 : readFromFile=true;
100 0 : aux.open(name.str().c_str());
101 0 : if (readFromFile && aux.good()) aux >> Npa >> Nw;
102 : else
103 0 : throw(SynthesisFTMachineError(string("Error while reading convolution "
104 0 : "function cache file ") + name.str( )));
105 : }
106 :
107 0 : if (Npa > 0)
108 : {
109 0 : paList.resize(Npa,true);
110 :
111 0 : IPosition s(2,Nw,Npa);
112 0 : XSup.resize(s,true);
113 0 : YSup.resize(s,true);
114 0 : Sampling.resize(Npa,true);
115 0 : for(Int i=0;i<Npa;i++)
116 : {
117 : Float pa, S;
118 : Int XS, YS;
119 0 : s[2]=i;
120 0 : aux >> pa;
121 0 : for(Int iw=0;iw<Nw;iw++)
122 : {
123 0 : s[0]=iw;
124 0 : aux >> XS >> YS;
125 0 : YS = XS;
126 0 : paList[i] = pa*M_PI/180.0;
127 0 : XSup(iw,i)=XS;
128 0 : YSup(iw,i)=YS;
129 : }
130 0 : aux >> S;
131 0 : Sampling[i]=S;
132 : }
133 0 : }
134 0 : }
135 0 : catch(AipsError& x)
136 : {
137 0 : throw(SynthesisFTMachineError(String("Error while initializing CF disk cache: ")
138 0 : +x.getMesg()));
139 0 : }
140 0 : }
141 : //
142 : //-----------------------------------------------------------------------
143 : //
144 79 : void CFCache::initPolMaps(PolMapType& polMap, PolMapType& conjPolMap)
145 : {
146 79 : if (OTODone()==false)
147 : {
148 158 : for(Int i=0;i<(Int)memCache2_p.nelements();i++)
149 79 : memCache2_p[i].initPolMaps(polMap, conjPolMap);
150 158 : for(Int i=0;i<(Int)memCacheWt2_p.nelements();i++)
151 79 : memCacheWt2_p[i].initPolMaps(polMap, conjPolMap);
152 79 : OTODone_p=true;
153 : }
154 79 : }
155 : //
156 : //-----------------------------------------------------------------------
157 : //
158 0 : void CFCache::summarize(CFStoreCacheType2& memStore, const String& message, const Bool cfsInfo)
159 : {
160 0 : LogOrigin logOrigin("CFCache2", "summarize");
161 0 : LogIO log_l(logOrigin);
162 :
163 0 : IPosition cfsShp=memStore[0].getShape();
164 0 : Int ipol=0;
165 :
166 0 : if (cfsInfo)
167 : {
168 0 : log_l << "PA: ";
169 0 : for (Int iBL=0; iBL<cfsShp(1); iBL++)
170 0 : for (Int iPA=0; iPA<cfsShp(0); iPA++)
171 : {
172 0 : Quantity pa; Int ant1, ant2;
173 0 : memStore[0].getParams(pa, ant1, ant2, iPA, iBL);
174 0 : log_l << pa.getValue("deg") << " ";
175 0 : }
176 0 : log_l << LogIO::NORMAL1;
177 : }
178 0 : log_l << message << LogIO::NORMAL1;
179 0 : for(Int iBL=0; iBL<cfsShp(1); iBL++)
180 0 : for(Int iPA=0; iPA<cfsShp(0); iPA++)
181 : {
182 0 : CFBuffer& cfb=memStore[0](iPA,iBL);
183 0 : IPosition cfbShp=cfb.getShape();
184 0 : for (Int iw=0; iw<cfbShp[1]; iw++)
185 : {
186 0 : log_l << "Support Size (w:"<< iw << ", PA:" << iPA << ", BL:" << iBL << ", C:*): ";
187 : {
188 0 : for (Int inu=0; inu<cfbShp[0]; inu++)
189 : {
190 0 : CFCell& cc=cfb(inu,iw,ipol);
191 0 : if (!cc.storage_p.null())
192 0 : log_l << cfb(inu, iw, ipol).xSupport_p << " ";
193 : }
194 0 : log_l << LogIO::NORMAL1;
195 : }
196 : }
197 0 : }
198 0 : }
199 : //
200 : //-----------------------------------------------------------------------
201 : //
202 : // By default (i.e., when called without any agruments), load all
203 : // the CFs found in the CF disk cache.
204 7 : void CFCache::initCacheFromList2(const String& path,
205 : const Vector<String>& cfFileNames,
206 : const Vector<String>& cfWtFileNames,
207 : Float selectedPA, Float dPA,
208 : const Int verbose)
209 : {
210 7 : Vector<String> cf, wtcf;
211 7 : cf=cfFileNames;
212 7 : wtcf=cfWtFileNames;
213 : // for (int i = 0; i < cfFileNames.nelements(); i++)
214 : // cf[i] = path+"/"+cfFileNames[i];
215 : // for (int i = 0; i < cfWtFileNames.nelements(); i++)
216 : // wtcf[i] = path+"/"+cfWtFileNames[i];
217 7 : fillCFListFromDisk(cf, path, memCache2_p, true, selectedPA, dPA,verbose);
218 7 : fillCFListFromDisk(wtcf, path, memCacheWt2_p, false, selectedPA, dPA, verbose);
219 7 : memCache2_p[0].primeTheCFB();
220 7 : memCacheWt2_p[0].primeTheCFB();
221 7 : if (verbose > 0) summarize(memCache2_p, "CFS", true);
222 : //summarize(memCacheWt2_p, "WTCFS", false);
223 7 : }
224 :
225 100 : void CFCache::setLazyFill(const Bool& lazyFill)
226 : {
227 100 : loadPixBuf_p=!lazyFill;
228 100 : if (lazyFill)
229 : {
230 162 : LogIO os( LogOrigin("CFCache","setLazyFill",WHERE));
231 81 : os << "Lazy fill is On" << LogIO::POST;
232 81 : }
233 100 : }
234 :
235 100 : void CFCache::initCache2(Bool verbose, Float selectedPA, Float dPA,casacore::String prefix)
236 : {
237 200 : LogOrigin logOrigin("CFCache2", "initCache2");
238 100 : LogIO log_l(logOrigin);
239 100 : Directory dirObj(Dir);
240 :
241 100 : if (Dir.length() == 0)
242 0 : throw(SynthesisFTMachineError(LogMessage("Got null string for disk cache dir. ",
243 0 : logOrigin).message()));
244 :
245 100 : Bool ReadOnly = SynthesisUtils::getenv("CFCache.READONLY",0);
246 :
247 100 : if (!dirObj.exists()) dirObj.create();
248 86 : else if ((!dirObj.isReadable()))
249 : {
250 0 : throw(SynthesisFTMachineError(String("Directory \"")+Dir+String("\"")+
251 0 : String(" for convolution function cache"
252 0 : " exists but is unreadable")));
253 : }
254 86 : else if ((!dirObj.isWritable()) && (!ReadOnly))
255 : {
256 0 : throw(SynthesisFTMachineError(String("Directory \"")+Dir+String("\"")+
257 0 : String(" for convolution function cache"
258 0 : " exists but is unwriteable")));
259 : }
260 :
261 :
262 100 : if(ReadOnly)
263 : {
264 0 : log_l << "CF Cache in read only mode." << LogIO::POST;
265 : }
266 :
267 100 : String memUnit="B";
268 100 : Double memUsed0=0, memUsed1=0;
269 : //
270 : // Lambda function to fill CFs with a given prefix.
271 : //
272 200 : auto fillCF_l = [&](CFStoreCacheType2& memCache2_l, casacore::String& cfprefix)
273 : {
274 200 : fillCFSFromDisk(dirObj,cfprefix, memCache2_l, true, selectedPA, dPA, verbose);
275 200 : memCache2_l[0].primeTheCFB();
276 200 : if (verbose > 0)
277 0 : summarize(memCache2_l, cfprefix, true);
278 200 : return 0.0;//memCache2_l[0].memUsage();
279 100 : };
280 :
281 100 : if (prefix == "CFS*") memUsed0 = fillCF_l(memCache2_p, prefix);
282 100 : else if (prefix == "WTCFS*") memUsed1 = fillCF_l(memCacheWt2_p, prefix);
283 : else // Load both, CFS* and WTCFS* if prefix is blank or of unexpected value
284 : {
285 100 : casacore::String s="CFS*";
286 100 : memUsed0 = fillCF_l(memCache2_p, s);
287 100 : s="WTCFS*";
288 100 : memUsed1 = fillCF_l(memCacheWt2_p, s);
289 100 : }
290 :
291 100 : double tt=memUsed0 + memUsed1;
292 100 : if (tt > 1024.0)
293 : {
294 0 : tt /= 1024.0;
295 0 : memUnit="KB";
296 : }
297 100 : if (tt > 0)
298 0 : log_l << "Total CF Cache memory footprint: " << tt << " (" << memUsed0 << "," << memUsed1 << ") " << memUnit << LogIO::POST;
299 100 : }
300 :
301 :
302 : // void CFCache::initCache2(Bool verbose, Float selectedPA, Float dPA,casacore::String /*prefix*/)
303 : // {
304 : // LogOrigin logOrigin("CFCache2", "initCache2");
305 : // LogIO log_l(logOrigin);
306 :
307 : // Directory dirObj(Dir);
308 :
309 : // if (Dir.length() == 0)
310 : // throw(SynthesisFTMachineError(LogMessage("Got null string for disk cache dir. ",
311 : // logOrigin).message()));
312 :
313 : // Bool ReadOnly = SynthesisUtils::getenv("CFCache.READONLY",0);
314 :
315 : // if (!dirObj.exists()) dirObj.create();
316 : // else if ((!dirObj.isReadable()))
317 : // {
318 : // throw(SynthesisFTMachineError(String("Directory \"")+Dir+String("\"")+
319 : // String(" for convolution function cache"
320 : // " exists but is unreadable")));
321 : // }
322 : // else if ((!dirObj.isWritable()) && (!ReadOnly))
323 : // {
324 : // throw(SynthesisFTMachineError(String("Directory \"")+Dir+String("\"")+
325 : // String(" for convolution function cache"
326 : // " exists but is unwriteable")));
327 : // }
328 :
329 :
330 : // if(ReadOnly)
331 : // {
332 : // log_l << "CF Cache in read only mode." << LogIO::POST;
333 : // }
334 :
335 : // fillCFSFromDisk(dirObj,"CFS*", memCache2_p, true, selectedPA, dPA, verbose);
336 : // fillCFSFromDisk(dirObj,"WTCFS*", memCacheWt2_p, false, selectedPA, dPA, verbose);
337 : // // memCache2_p[0].show("Re-load CFS",cerr);
338 : // // memCacheWt2_p[0].show("Re-load WTCFS",cerr);
339 : // memCache2_p[0].primeTheCFB();
340 : // memCacheWt2_p[0].primeTheCFB();
341 :
342 : // // memCache2_p[0].show("CF Cache: ");
343 : // // memCacheWt2_p[0].show("WTCF Cache: ");
344 :
345 : // Double memUsed0,memUsed1;
346 : // String memUnit="B";
347 : // memUsed0=memCache2_p[0].memUsage();
348 : // memUsed1=memCacheWt2_p[0].memUsage();
349 : // if (memUsed0 > 1024.0)
350 : // {
351 : // memUsed0 /= 1024.0;
352 : // memUsed1 /= 1024.0;
353 : // memUnit="KB";
354 : // }
355 :
356 : // if (verbose > 0)
357 : // {
358 : // summarize(memCache2_p, "CFS", True);
359 : // summarize(memCacheWt2_p, "WTCFS", False);
360 : // }
361 :
362 : // if (memUsed0+memUsed1 > 0)
363 : // log_l << "Total CF Cache memory footprint: " << (memUsed0+memUsed1) << " (" << memUsed0 << "," << memUsed1 << ") " << memUnit << LogIO::POST;
364 :
365 : // // memCache2_p[0].makePersistent("./junk.cf");
366 : // // memCacheWt2_p[0].makePersistent("./junk.cf","","WT");
367 : // }
368 :
369 :
370 : //
371 : //-----------------------------------------------------------------------
372 : //
373 214 : void CFCache::fillCFListFromDisk(const Vector<String>& fileNames,
374 : const String& CFCDir, CFStoreCacheType2& memStore,
375 : Bool showInfo, Float selectPAVal, Float dPA,
376 : const Int verbose)
377 : {
378 428 : LogOrigin logOrigin("CFCache2", "fillCFListFromDisk");
379 214 : LogIO log_l(logOrigin);
380 : (void)showInfo;
381 214 : Bool selectPA = (fabs(selectPAVal) <= 360.0);
382 : try
383 : {
384 214 : if (memStore.nelements() == 0) memStore.resize(1,true);
385 214 : memStore[0].setLazyFill(!loadPixBuf_p);
386 214 : memStore[0].setCFCacheDir(getCacheDir());
387 214 : CFCacheTableType cfCacheTable_l;
388 : // Regex regex(Regex::fromPattern(pattern));
389 : // Vector<String> fileNames(dirObj.find(regex));
390 :
391 214 : if (fileNames.nelements() > 0)
392 : {
393 : // String CFCDir=dirObj.path().absoluteName();
394 : // if (showInfo)
395 : // log_l << "No. of " << pattern << " found in "
396 : // << dirObj.path().originalName() << ": "
397 : // << fileNames.nelements() << LogIO::POST;
398 :
399 : //
400 : // Gather the list of PA values
401 : //
402 : {
403 186 : ProgressMeter pm(1.0, Double(fileNames.nelements()),
404 558 : "Reading CFCache aux. info.", "","","",true);
405 1302 : for (uInt i=0; i < fileNames.nelements(); i++)
406 : {
407 2232 : PagedImage<Complex> thisCF(CFCDir+'/'+fileNames[i]);
408 1116 : TableRecord miscinfo = thisCF.miscInfo();
409 : // miscinfo.print(cerr);
410 : Double paVal;
411 : //UNUSED: Double wVal; Int mVal;
412 1116 : miscinfo.get("ParallacticAngle",paVal);
413 1116 : paList_p.push_back(paVal);
414 1116 : pm.update(Double(i));
415 1116 : }
416 186 : }
417 : //
418 : // Make the PA-value list unique
419 : //
420 186 : sort( paList_p.begin(), paList_p.end() );
421 186 : paList_p.erase( unique( paList_p.begin(), paList_p.end() ), paList_p.end() );
422 186 : cfCacheTable_l.resize(paList_p.size());
423 :
424 : //
425 : // For each CF, load the PA, Muelller element, WValue and
426 : // the Ref. Freq. Insert these values in the lists in the
427 : // cfCacheTable
428 : //
429 186 : Array<Complex> pixBuf;
430 :
431 : // Int cfCount=0;
432 : // for (uInt i=0; i<fileNames.nelements(); i++)
433 : // {
434 : // Double paVal, wVal, fVal, sampling; Int mVal, xSupport, ySupport;
435 : // CoordinateSystem coordSys;
436 :
437 : // getCFParams(fileNames[i], pixBuf, coordSys, sampling, paVal,
438 : // xSupport, ySupport, fVal, wVal, mVal,false);
439 : // Bool pickThisCF=true;
440 : // if (selectPA) pickThisCF = (fabs(paVal - selectPAVal) <= dPA);
441 : // cerr << fileNames[i] << " " << paVal << " " << selectPAVal << " " << dPA << " " << pickThisCF << endl;
442 : // if (pickThisCF) cfCount++;
443 : // }
444 : // cerr << "Will load " << cfCount << "CFs." << endl;
445 186 : log_l << "Loading misc info from CFs" << LogIO::POST;
446 186 : TableRecord miscInfo;
447 : {
448 186 : ProgressMeter pm(1.0, Double(fileNames.nelements()),
449 558 : "Loading CFs", "","","",true);
450 1302 : for (uInt i=0; i < fileNames.nelements(); i++)
451 : {
452 : Double paVal, wVal, fVal, sampling, conjFreq; Int mVal, xSupport, ySupport, conjPoln;
453 1116 : CoordinateSystem coordSys;
454 :
455 1116 : IPosition cfShape;
456 2232 : miscInfo = SynthesisUtils::getCFParams(Dir, fileNames[i], cfShape, pixBuf, coordSys, sampling, paVal,
457 1116 : xSupport, ySupport, fVal, wVal, mVal,conjFreq, conjPoln,false);
458 :
459 1116 : Bool pickThisCF=true;
460 1116 : if (selectPA) pickThisCF = (fabs(paVal - selectPAVal) <= dPA);
461 1116 : if (pickThisCF)
462 : {
463 1116 : Int ipos; SynthesisUtils::stdNearestValue(paList_p, (Float)paVal,ipos);
464 1116 : uInt paPos=ipos;
465 :
466 1116 : if (paPos < paList_p.size())
467 : {
468 1116 : cfCacheTable_l[paPos].freqList.push_back(fVal);
469 1116 : cfCacheTable_l[paPos].wList.push_back(wVal);
470 1116 : cfCacheTable_l[paPos].muellerList.push_back(mVal);
471 1116 : cfCacheTable_l[paPos].cfNameList.push_back(fileNames[i]);
472 : //cerr << paPos << " " << fileNames[i] << endl;
473 : }
474 : }
475 1116 : pm.update(Double(fileNames.nelements()));
476 1116 : }
477 186 : }
478 372 : for (uInt ipa=0; ipa < cfCacheTable_l.size(); ipa++)
479 : {
480 : //
481 : // Resize the CFStore (poorly named private variable
482 : // memCache2_p) to add CFBuffer for the each entry in
483 : // the paList.
484 : //
485 186 : vector<String> fileNames(cfCacheTable_l[ipa].cfNameList);
486 :
487 186 : Quantity paQuant(paList_p[ipa],"deg"), dPA(1.0,"deg");
488 186 : memStore[0].resize(paQuant, dPA, 0,0);
489 186 : CountedPtr<CFBuffer> cfb=memStore[0].getCFBuffer(paQuant, dPA, 0, 0);
490 :
491 : //
492 : // Get the list of f, w, mVals from cfCacheTable_l for
493 : // the current ipa index. Sort them. And convert
494 : // them into a list of unique entires.
495 : //
496 186 : vector<Double> fList(cfCacheTable_l[ipa].freqList),
497 186 : wList(cfCacheTable_l[ipa].wList);
498 186 : vector<Int> mList(cfCacheTable_l[ipa].muellerList);
499 186 : sort( fList.begin(), fList.end() );
500 186 : sort( wList.begin(), wList.end() );
501 186 : sort( mList.begin(), mList.end() );
502 186 : fList.erase(SynthesisUtils::Unique(fList.begin(), fList.end()), fList.end());
503 186 : wList.erase(SynthesisUtils::Unique(wList.begin(), wList.end()), wList.end());
504 186 : mList.erase(SynthesisUtils::Unique(mList.begin(), mList.end()), mList.end());
505 186 : PolMapType muellerElements;
506 186 : Int npol=mList.size();
507 186 : muellerElements.resize(npol);
508 558 : for (Int ii=0;ii<npol;ii++)
509 : {
510 372 : muellerElements[ii].resize(1);
511 372 : muellerElements[ii][0]=mList[ii];
512 : }
513 186 : Double wIncr; miscInfo.get("WIncr", wIncr);
514 186 : Vector<Double> const wListV(wList);
515 186 : Vector<Double> const fListV(fList);
516 186 : cfb->resize(wIncr,0.0,wListV,fListV,
517 : muellerElements,muellerElements,muellerElements,muellerElements);
518 186 : cfb->setPA(paList_p[ipa]);
519 186 : cfb->setDir(Dir);
520 : //
521 : // Now go over the list of fileNames corresponding to
522 : // the current PA value and them the current CFBuffer.
523 : //
524 1302 : for (uInt nf=0; nf<fileNames.size(); nf++)
525 : {
526 : Double paVal, wVal, fVal, sampling, conjFreq;
527 : Int mVal, xSupport, ySupport, conjPoln;
528 1116 : CoordinateSystem coordSys;
529 : //
530 : // Get the parameters from the CF file
531 : //
532 1116 : IPosition cfShape;
533 2232 : TableRecord miscInfo = SynthesisUtils::getCFParams(Dir,fileNames[nf], cfShape, pixBuf, coordSys, sampling, paVal,
534 1116 : xSupport, ySupport, fVal, wVal, mVal, conjFreq, conjPoln,loadPixBuf_p,True);
535 : //
536 : // Get the storage buffer from the CFBuffer and
537 : // fill it in what we got from the getCFParams
538 : // call above.
539 : //
540 1116 : if (loadPixBuf_p)
541 : {
542 228 : Array<Complex> &cfBuf=(*(cfb->getCFCellPtr(fVal, wVal,mVal)->storage_p));
543 : //
544 : // Fill the cfBuf with the pixel array from the
545 : // disk file. Add it, along with the extracted CF
546 : // parameters to the CFBuffer.
547 : //
548 228 : cfBuf.assign(pixBuf);
549 : }
550 :
551 : //cfb->addCF(&cfBuf,coordSys,fsampling,xSupport,ySupport,fVal,wVal,mVal);
552 : Int fndx,wndx, mndx;
553 1116 : SynthesisUtils::stdNearestValue(fList, fVal, fndx);
554 1116 : SynthesisUtils::stdNearestValue(wList, wVal, wndx);
555 1116 : SynthesisUtils::stdNearestValue(mList, mVal, mndx);
556 : //Float fsampling=sampling;
557 : //
558 : // The coordSys, sampling, xSupport, ySuport
559 : // params are set for the CFCell at the location
560 : // determined by fndx and wndx. The mndx is
561 : // determined inside using mVal (why this
562 : // treatment for mndx, please don't ask. Not just
563 : // yet (SB)).
564 : // String telescopeName;miscInfo.get("TelescopeName",telescopeName);
565 : //Float diameter; miscInfo.get("Diameter",diameter);
566 : //cfb->setParams(fndx, wndx, 0,0, coordSys, fsampling, xSupport, ySupport,
567 : // fVal, wVal, mVal,fileNames[nf],conjFreq, conjPoln,
568 : // telescopeName, diameter);
569 :
570 : // cfb->setParams(fndx, wndx, mVal, miscInfo);
571 1116 : auto ndx=cfb->setParams(fndx, wndx, 0,0, fVal, wVal, mVal, coordSys,miscInfo);
572 1116 : (cfb->getCFCellPtr(ndx(0), ndx(1), ndx(2)))->shape_p=cfShape;
573 :
574 1116 : if (verbose > 0) log_l << cfCacheTable_l[ipa].cfNameList[nf]
575 : << "[" << fndx << "," << wndx << "," << mndx << "] "
576 0 : << paList_p[ipa] << " " << xSupport << LogIO::POST;
577 1116 : }
578 : //cfb->show("cfb: ");
579 186 : }
580 :
581 186 : }
582 214 : }
583 0 : catch(AipsError& x)
584 : {
585 0 : throw(SynthesisFTMachineError(String("Error while initializing CF disk cache: ")
586 0 : +x.getMesg()));
587 0 : }
588 :
589 : //memStore[0].getCFBuffer(0,0)->show("CFB0: ");
590 :
591 214 : }
592 : //
593 : //-----------------------------------------------------------------------
594 : //
595 200 : void CFCache::fillCFSFromDisk(const Directory dirObj, const String& pattern,
596 : CFStoreCacheType2& memStore,
597 : Bool showInfo, Float selectPAVal, Float dPA,
598 : const Int verbose)
599 : {
600 400 : LogOrigin logOrigin("CFCache2", "fillCFSFromDisk");
601 200 : LogIO log_l(logOrigin);
602 : try
603 : {
604 200 : Regex regex(Regex::fromPattern(pattern));
605 200 : Vector<String> fileNames(dirObj.find(regex));
606 200 : String CFCDir=dirObj.path().absoluteName();
607 200 : if (showInfo)
608 : log_l << "No. of " << pattern << " found in "
609 400 : << dirObj.path().originalName() << ": "
610 200 : << fileNames.nelements() << LogIO::POST;
611 :
612 200 : fillCFListFromDisk(fileNames, CFCDir, memStore, showInfo, selectPAVal, dPA, verbose);
613 200 : }
614 0 : catch(AipsError& x)
615 : {
616 0 : throw(SynthesisFTMachineError(String("Error while initializing CF disk cache: ")
617 0 : +x.getMesg()));
618 0 : }
619 200 : }
620 : //
621 : //-----------------------------------------------------------------------
622 : //
623 0 : TableRecord CFCache::getCFParams(const String& fileName,
624 : Array<Complex>& pixelBuffer,
625 : CoordinateSystem& coordSys,
626 : Double& sampling,
627 : Double& paVal,
628 : Int& xSupport, Int& ySupport,
629 : Double& fVal, Double& wVal, Int& mVal,
630 : Double& conjFreq, Int& conjPoln,
631 : Bool loadPixels)
632 : {
633 : try
634 : {
635 0 : PagedImage<Complex> thisCF(Dir+'/'+fileName);
636 0 : TableRecord miscinfo = thisCF.miscInfo();
637 :
638 0 : if (loadPixels) pixelBuffer.assign(thisCF.get());
639 0 : miscinfo.get("ParallacticAngle", paVal);
640 0 : miscinfo.get("MuellerElement", mVal);
641 0 : miscinfo.get("WValue", wVal);
642 0 : miscinfo.get("Xsupport", xSupport);
643 0 : miscinfo.get("Ysupport", ySupport);
644 0 : miscinfo.get("Sampling", sampling);
645 0 : miscinfo.get("ConjFreq", conjFreq);
646 0 : miscinfo.get("ConjPoln", conjPoln);
647 0 : Int index= thisCF.coordinates().findCoordinate(Coordinate::SPECTRAL);
648 0 : coordSys = thisCF.coordinates();
649 0 : SpectralCoordinate spCS = coordSys.spectralCoordinate(index);
650 0 : fVal=static_cast<casacore::Float>(spCS.referenceValue()(0));
651 0 : return miscinfo;
652 0 : }
653 0 : catch(AipsError& x)
654 : {
655 0 : throw(SynthesisFTMachineError(String("Error in CFCache::getCFParams(): ")
656 0 : +x.getMesg()));
657 0 : }
658 : }
659 : //
660 : //-----------------------------------------------------------------------
661 : //
662 0 : CFCache& CFCache::operator=(const CFCache& other)
663 : {
664 : // if (this != other)
665 : {
666 0 : paList = other.paList;
667 0 : Sampling = other.Sampling;
668 0 : XSup = other.XSup;
669 0 : YSup = other.YSup;
670 0 : Dir = other.Dir;
671 0 : cfPrefix = other.cfPrefix;
672 0 : WtImagePrefix = other.WtImagePrefix;
673 0 : aux = other.aux;
674 0 : paCD_p = other.paCD_p;
675 0 : memCache_p = other.memCache_p;
676 0 : memCacheWt_p = other.memCacheWt_p;
677 0 : cfCacheTable_p = other.cfCacheTable_p;
678 0 : OTODone_p = other.OTODone_p;
679 0 : loadPixBuf_p=other.loadPixBuf_p;
680 : }
681 0 : return *this;
682 : };
683 : //
684 : //-----------------------------------------------------------------------
685 : //
686 0 : Long CFCache::size()
687 : {
688 0 : Long s=0;
689 0 : for(uInt i=0;i<memCache_p.nelements();i++)
690 0 : s+=memCache_p[0].data->size();
691 0 : for(uInt i=0;i<memCacheWt_p.nelements();i++)
692 0 : s+=memCacheWt_p[0].data->size();
693 :
694 0 : return s*sizeof(Complex);
695 : }
696 : //
697 : //-----------------------------------------------------------------------
698 : //
699 0 : void CFCache::makeFTCoordSys(const CoordinateSystem& coords,
700 : const Int& convSize,
701 : const Vector<Double>& ftRef,
702 : CoordinateSystem& ftCoords)
703 : {
704 : Int directionIndex;
705 :
706 0 : ftCoords = coords;
707 0 : directionIndex=ftCoords.findCoordinate(Coordinate::DIRECTION);
708 : // The following line follows the (lame) logic that if a
709 : // DIRECTION axis was not found, the coordinate system must be of
710 : // the FT domain already
711 0 : if (directionIndex == -1) return;
712 :
713 0 : DirectionCoordinate dc;//=coords.directionCoordinate(directionIndex);
714 : // AlwaysAssert(directionIndex>=0, AipsError);
715 0 : dc=coords.directionCoordinate(directionIndex);
716 0 : Vector<Bool> axes(2); axes(0)=axes(1)=true;//axes(2)=true;
717 0 : Vector<Int> shape(2,convSize);
718 :
719 : //cerr << "CFC: " << shape << endl;
720 :
721 0 : Vector<Double>ref(4);
722 0 : ref(0)=ref(1)=ref(2)=ref(3)=0;
723 0 : dc.setReferencePixel(ref);
724 0 : Coordinate* ftdc=dc.makeFourierCoordinate(axes,shape);
725 0 : Vector<Double> refVal;
726 0 : refVal=ftdc->referenceValue();
727 : // refVal(0)=refVal(1)=0;
728 : // ftdc->setReferenceValue(refVal);
729 0 : ref(0)=ftRef(0);
730 0 : ref(1)=ftRef(1);
731 0 : ftdc->setReferencePixel(ref);
732 :
733 0 : ftCoords.replaceCoordinate(*ftdc, directionIndex);
734 0 : delete ftdc; ftdc=0;
735 0 : }
736 : //
737 : //-------------------------------------------------------------------------
738 : //
739 0 : Int CFCache::addToMemCache(CFStoreCacheType& memCache_l,
740 : Float pa, CFType* cf,
741 : CoordinateSystem& coords,
742 : Vector<Int>& xConvSupport,
743 : Vector<Int>& yConvSupport,
744 : Float convSampling)
745 : {
746 0 : Float dPA=paCD_p.getParAngleTolerance().getValue("rad");
747 :
748 0 : Int where=-1, wConvSize = cf->shape()(CFDefs::NWPOS);
749 0 : Bool found=searchConvFunction(where, pa, dPA);
750 : //
751 : // If the PA value was not found, the return value in "where" is
752 : // the negative of the location in which the PA value should be
753 : // found. Convert it to positive value to be used for resizing
754 : // etc.
755 : //
756 0 : where=abs(where);
757 : //
758 : // Resize the arrays if the CF for the relevant PA was not in the
759 : // MEM cache. Note that if the arrays are already of size
760 : // where+1, Array<>::resize() is a no-op.
761 : //
762 0 : Int N=memCache_l.nelements();
763 :
764 0 : memCache_l.resize(max(N,where+1), true);
765 0 : if ((Int)paList.nelements() <= where)
766 : {
767 0 : IPosition s(2,wConvSize,where+1);
768 0 : paList.resize(where+1,true);
769 0 : XSup.resize(s,true); YSup.resize(s,true);
770 0 : Sampling.resize(where+1,true);
771 0 : }
772 : //
773 : // If the PA was not found, enter the aux. values in the internal
774 : // arrays.
775 : //
776 0 : if (!found)
777 : {
778 0 : paList[where] = pa;
779 0 : for(Int iw=0;iw<wConvSize;iw++)
780 : {
781 0 : YSup(iw,where) = xConvSupport(iw);
782 0 : XSup(iw,where) = yConvSupport(iw);
783 : }
784 0 : Sampling[where]=convSampling;
785 : }
786 : //
787 : // If the CF was not in the mem. cache, add it.
788 : //
789 0 : if (memCache_l[where].null())
790 : {
791 0 : Vector<Float> sampling(1);sampling[0]=convSampling;
792 :
793 0 : Int maxXSup=max(xConvSupport), maxYSup=max(yConvSupport);
794 0 : memCache_l[where] = CFStore(cf,coords,sampling,
795 : xConvSupport,yConvSupport,
796 0 : maxXSup,maxYSup,Quantity(pa,"rad"),
797 0 : 0);
798 0 : }
799 :
800 0 : return where;
801 : }
802 : //-------------------------------------------------------------------------
803 : // Write the conv. functions from the mem. cache to the disk cache.
804 : //
805 0 : Int CFCache::cacheConvFunction(Int which, const Float& pa, CFType& cf,
806 : CoordinateSystem& coords,
807 : CoordinateSystem& ftCoords,
808 : Int &convSize,
809 : Vector<Int> &xConvSupport,
810 : Vector<Int> &yConvSupport,
811 : Float convSampling,
812 : String nameQualifier,
813 : Bool savePA)
814 : {
815 0 : LogIO log_l(LogOrigin("CFCache2","cacheConvFunction"));
816 0 : Int whereCached_l=-1;
817 0 : if (Dir.length() == 0) return whereCached_l;
818 0 : if (which < 0)
819 : {
820 : Int i;
821 0 : searchConvFunction(i,pa,paCD_p.getParAngleTolerance().get("rad"));
822 : // which = paList.nelements();
823 0 : which = abs(i);
824 : }
825 :
826 : try
827 : {
828 0 : IPosition newConvShape = cf.shape();
829 0 : Int wConvSize = newConvShape(CFDefs::NWPOS);
830 0 : for(Int iw=0;iw<wConvSize;iw++)
831 : {
832 0 : IPosition sliceStart(4,0,0,iw,0),
833 0 : sliceLength(4,newConvShape(CFDefs::NXPOS),
834 0 : newConvShape(CFDefs::NYPOS),
835 : 1,
836 0 : newConvShape(CFDefs::NPOLPOS));
837 :
838 0 : Vector<Double> ftRef(2);
839 0 : ftRef(0)=newConvShape(CFDefs::NXPOS)/2-1;
840 0 : ftRef(1)=newConvShape(CFDefs::NYPOS)/2-1;
841 0 : makeFTCoordSys(coords, convSize, ftRef, ftCoords);
842 0 : ostringstream name;
843 0 : name << Dir << "/" << cfPrefix << nameQualifier << iw << "_" << which;
844 0 : const CFType tmpArr=cf(Slicer(sliceStart,sliceLength));
845 :
846 : // storeArrayAsImage(name, ftCoords,tmpArr);
847 :
848 0 : IPosition screenShape(4,newConvShape(CFDefs::NXPOS),
849 0 : newConvShape(CFDefs::NYPOS),
850 0 : newConvShape(CFDefs::NPOLPOS),
851 0 : 1);
852 0 : Record miscinfo;
853 0 : miscinfo.define("Xsupport",xConvSupport);
854 0 : miscinfo.define("Ysupport",yConvSupport);
855 0 : miscinfo.define("sampling", convSampling);
856 0 : miscinfo.define("ParallacticAngle",pa);
857 0 : PagedImage<Complex> thisScreen(screenShape, ftCoords, name);
858 0 : thisScreen.setMiscInfo(miscinfo);
859 0 : Array<Complex> buf;
860 0 : buf=((cf(Slicer(sliceStart,sliceLength)).nonDegenerate()));
861 0 : thisScreen.put(buf);
862 0 : }
863 0 : if (savePA)
864 : {
865 0 : CFStoreCacheType& memCacheObj = getMEMCacheObj(nameQualifier);
866 0 : whereCached_l=addToMemCache(memCacheObj, pa,&cf, ftCoords,
867 : xConvSupport, yConvSupport, convSampling);
868 : }
869 0 : }
870 0 : catch (AipsError& x)
871 : {
872 0 : throw(SynthesisFTMachineError("Error while caching CF to disk in "+x.getMesg()));
873 0 : }
874 0 : return whereCached_l;
875 0 : }
876 : //
877 : //-------------------------------------------------------------------------
878 : //
879 0 : void CFCache::cacheConvFunction(const Float pa, CFStore& cfs,
880 : String nameQualifier,Bool savePA)
881 : {
882 0 : if (cfs.data.null())
883 0 : throw(SynthesisError(LogMessage("Won't cache a NULL CFStore",
884 0 : LogOrigin("CFCache::cacheConvFunction")).message()));
885 0 : CoordinateSystem ftcoords;
886 0 : Int which=-1, whereCached=-1;
887 0 : Int convSize=(Int)cfs.data->shape()(0);
888 :
889 0 : whereCached = cacheConvFunction(which, pa, *(cfs.data), cfs.coordSys, ftcoords, convSize,
890 0 : cfs.xSupport, cfs.ySupport, cfs.sampling[0],
891 : nameQualifier,savePA);
892 0 : cfs.coordSys = ftcoords;
893 0 : if (whereCached > -1)
894 : {
895 0 : CFStoreCacheType& memCacheObj=getMEMCacheObj(nameQualifier);
896 0 : cfs=memCacheObj[whereCached];
897 : }
898 0 : }
899 : //
900 : //-------------------------------------------------------------------------
901 : //
902 0 : Bool CFCache::searchConvFunction(Int& which,
903 : const Float pa, const Float dPA)
904 : {
905 0 : if (paList.nelements()==0) initCache();
906 0 : Int i,NPA=paList.nelements(); Bool paFound=false;
907 : Float iPA;
908 :
909 0 : Float paDiff=2*dPA;
910 0 : Int saveNdx=-1;
911 :
912 0 : saveNdx = -1;
913 0 : for(i=0;i<NPA;i++)
914 : {
915 0 : iPA = paList[i];
916 0 : if (fabs(iPA - pa) < paDiff)
917 : {
918 0 : saveNdx = i;
919 0 : paDiff = fabs(iPA-pa);
920 : }
921 : }
922 0 : if (saveNdx > -1)
923 : {
924 0 : iPA = paList[saveNdx];
925 0 : if (fabs(iPA - pa) <= dPA)
926 : {
927 0 : i = saveNdx;
928 0 : paFound=true;
929 : }
930 : }
931 0 : if (paFound) which = i;
932 0 : else which = -i;
933 0 : return paFound;
934 : }
935 : //
936 : //-------------------------------------------------------------------------
937 : //Write the aux. info. also in the disk cache (wonder if this should
938 : //be automatically be called from cacheCFtion() method).
939 : //
940 169 : void CFCache::flush()
941 : {
942 : // If WtImagePrefix is set, no need to save avgPB in the CFCache
943 169 : if (WtImagePrefix != "") return;
944 :
945 0 : LogIO log_l(LogOrigin("CFCache2", "flush"));
946 :
947 0 : if (Dir.length() == 0) return;
948 0 : ostringstream name;
949 :
950 0 : name << Dir << "/aux.dat";
951 0 : Int n=memCache_p.nelements(),nw;
952 0 : if (n>0)
953 : try
954 : {
955 0 : nw=memCache_p[0].xSupport.nelements();
956 0 : ofstream aux(name.str().c_str());
957 0 : aux << n << " " << nw << endl;
958 0 : for(Int ipa=0;ipa<n;ipa++)
959 : {
960 0 : aux << paList[ipa]*180.0/M_PI << " ";
961 0 : for(int iw=0;iw<nw;iw++)
962 0 : aux << memCache_p[ipa].xSupport(iw) << " " << memCache_p[ipa].ySupport(iw) << " ";
963 0 : aux << " " << Sampling[ipa] <<endl;
964 : }
965 0 : }
966 0 : catch(AipsError &x)
967 : {
968 0 : throw(SynthesisFTMachineError(string("Error while writing ")
969 0 : + name.str( ) + (string) x.getMesg()));
970 0 : }
971 0 : }
972 : //
973 : //-------------------------------------------------------------------------
974 : //Along with the aux. info., also save the average PB in the disk cache.
975 : //
976 52 : void CFCache::flush(ImageInterface<Float>& avgPB, String qualifier)
977 : {
978 : // If WtImagePrefix is set, no need to save avgPB in the CFCache
979 52 : if (WtImagePrefix != "") return;
980 :
981 0 : LogIO log_l(LogOrigin("CFCache2", "flush"));
982 :
983 0 : if (Dir.length() == 0) return;
984 0 : flush();
985 0 : ostringstream Name;
986 0 : Name << Dir <<"/avgPB" << qualifier;
987 : try
988 : {
989 0 : storeImg(Name, avgPB);
990 0 : avgPBReady_p=true;
991 0 : avgPBReadyQualifier_p = qualifier;
992 : }
993 0 : catch(AipsError &x)
994 : {
995 0 : throw(SynthesisFTMachineError(string("Error while writing ")
996 0 : + Name.str( ) + (string) x.getMesg()));
997 0 : }
998 0 : }
999 : //
1000 : //-------------------------------------------------------------------------
1001 : //Load the average PB from the disk cache.
1002 : //
1003 164 : Int CFCache::loadWtImage(ImageInterface<Float>& avgPB, String qualifier, std::tuple<int, double> cubeinfo)
1004 : {
1005 328 : LogIO log_l(LogOrigin("CFCache2", "loadWtImage"));
1006 164 : ostringstream name, sumWtName;
1007 164 : name << WtImagePrefix << ".weight" << qualifier;
1008 164 : if (qualifier != "") sumWtName << WtImagePrefix << ".sumwt.tt0";// << qualifier;
1009 52 : else sumWtName << WtImagePrefix << ".sumwt";
1010 :
1011 : try
1012 : {
1013 : // First try to load .weight image. If this fails, AipsError
1014 : // will be caught and a NOTCACHED returned.
1015 512 : PagedImage<Float> tmp(name.str().c_str());
1016 : //cerr << "BBBBBBBBBBBBBB " << tmp.shape() << endl;
1017 : // Now try to load .sumwt. If .sumwt is not found, this is a
1018 : // fatal error (inconsistancy on the disk). So this time
1019 : // throw a SEVER exception.
1020 72 : Float sumwt=1.0;
1021 : try
1022 : {
1023 144 : PagedImage<Float> sumWtTmp(sumWtName.str().c_str());
1024 72 : sumwt=max(sumWtTmp.get());
1025 72 : }
1026 0 : catch (AipsError& x)
1027 : {
1028 0 : log_l << "Sum-of-weights not found " << x.getMesg() << LogIO::SEVERE;
1029 0 : }
1030 :
1031 72 : int nchan=std::get<0> (cubeinfo);
1032 72 : Double freqofBegChan=std::get<1>(cubeinfo);
1033 72 : if(freqofBegChan > 0.0 && (tmp.shape()(3) > nchan)){
1034 :
1035 : //get freq of first chan of chunk
1036 0 : CoordinateSystem cs=tmp.coordinates();
1037 0 : SpectralCoordinate fsys=cs.spectralCoordinate(cs.findCoordinate(Coordinate::SPECTRAL));
1038 : Double startchan;
1039 0 : fsys.toPixel(startchan, freqofBegChan);
1040 0 : Int endchan=startchan+nchan-1;
1041 0 : ImageInterface<Float>* subim=SpectralImageUtil::getChannel(tmp,Int(startchan), endchan);
1042 0 : avgPB.resize(subim->shape());
1043 0 : LatticeExpr<Float> myexpr( (*subim)*sumwt);
1044 0 : avgPB.copyData(myexpr);
1045 0 : avgPB.setCoordinateInfo(subim->coordinates());
1046 0 : delete subim;
1047 0 : }
1048 : else{
1049 72 : avgPB.resize(tmp.shape());
1050 144 : LatticeExpr<Float> myexpr( tmp*sumwt);
1051 72 : avgPB.copyData(myexpr);
1052 72 : avgPB.setCoordinateInfo(tmp.coordinates());
1053 72 : }
1054 : //avgPB.resize(tmp.shape());
1055 : //avgPB.put(tmp.get()*sumwt);
1056 : //avgPB.setCoordinateInfo(tmp.coordinates());
1057 : //cerr << "peak = " << max(tmp.get()*sumwt) << endl;
1058 72 : }
1059 92 : catch(AipsError& x) // Just rethrowing the exception for now.
1060 : // Ultimately, this should be used to make
1061 : // the state of this object consistant.
1062 : {
1063 92 : return NOTCACHED;
1064 92 : }
1065 72 : log_l << "Loaded \"" << name.str() << "\"" << LogIO::POST;
1066 72 : avgPBReady_p=true;
1067 72 : return DISKCACHE;
1068 :
1069 164 : }
1070 : //
1071 : //-------------------------------------------------------------------------
1072 : //Load the average PB from the disk cache.
1073 : //
1074 164 : Int CFCache::loadAvgPB(ImageInterface<Float>& avgPB, String qualifier, std::tuple<int, double>cubeinfo)
1075 : {
1076 328 : LogIO log_l(LogOrigin("CFCache2", "loadAvgPB"));
1077 :
1078 164 : if (WtImagePrefix != "")
1079 : {
1080 164 : return loadWtImage(avgPB, qualifier, cubeinfo);
1081 : }
1082 :
1083 0 : if (Dir.length() == 0)
1084 0 : throw(SynthesisFTMachineError("Cache dir. name null"));
1085 :
1086 :
1087 0 : ostringstream name;
1088 0 : name << Dir << "/avgPB" << qualifier;
1089 : // cout << name.str() << endl;
1090 : try
1091 : {
1092 0 : PagedImage<Float> tmp(name.str().c_str());
1093 0 : int nchan=std::get<0> (cubeinfo);
1094 0 : Double freqofBegChan=std::get<1>(cubeinfo);
1095 0 : if(freqofBegChan > 0.0 && (tmp.shape()(3) > nchan)){
1096 :
1097 : //get freq of first chan of chunk
1098 0 : CoordinateSystem cs=tmp.coordinates();
1099 0 : SpectralCoordinate fsys=cs.spectralCoordinate(cs.findCoordinate(Coordinate::SPECTRAL));
1100 : Double startchan;
1101 0 : fsys.toPixel(startchan, freqofBegChan);
1102 0 : Int endchan=startchan+nchan-1;
1103 0 : ImageInterface<Float>* subim=SpectralImageUtil::getChannel(tmp,Int(startchan), endchan);
1104 0 : avgPB.resize(subim->shape());
1105 0 : avgPB.copyData(*subim);
1106 0 : avgPB.setCoordinateInfo(subim->coordinates());
1107 0 : delete subim;
1108 0 : }
1109 : else{
1110 0 : avgPB.resize(tmp.shape());
1111 0 : avgPB.copyData(tmp);
1112 0 : avgPB.setCoordinateInfo(tmp.coordinates());
1113 : }
1114 :
1115 0 : }
1116 0 : catch(AipsError& x) // Just rethrowing the exception for now.
1117 : // Ultimately, this should be used to make
1118 : // the state of this object consistant.
1119 : {
1120 0 : return NOTCACHED;
1121 0 : }
1122 0 : log_l << "Loaded \"" << name.str() << "\"" << LogIO::POST;
1123 0 : avgPBReady_p=true;
1124 0 : return DISKCACHE;
1125 164 : }
1126 : //
1127 : //-------------------------------------------------------------------------
1128 : //Load a conv. func. from the disk. This is non-optimal due to the
1129 : //data structure used for the conv. func. in-memory cache (it's an
1130 : //array of pointers where it should really be a List of pointers).
1131 : //The conf. func. index, which is also used as a key to located them
1132 : //in the mem. cache, are not assured to be contiguous. As a result,
1133 : //in the current implementation there can be gaps in the
1134 : //convFuncCache array. These gaps are initialized to NULL pointers.
1135 : //It's not much of a memory waste, but still non-optimal! Leaving
1136 : //it like this for now.
1137 : //
1138 : // Return TRUE if loaded from disk and FLASE if found in the mem. cache.
1139 : //
1140 0 : Int CFCache::loadFromDisk(Int where, Float /*pa*/, Float /*dPA*/,
1141 : Int Nw, CFStoreCacheType &convFuncCache,
1142 : CFStore& cfs, String nameQualifier)
1143 : {
1144 0 : LogIO log_l(LogOrigin("CFCache2", "loadFromDisk"));
1145 :
1146 0 : Vector<Int> xconvSupport,yconvSupport;;
1147 0 : Vector<Float> convSampling;
1148 : //Double cfRefFreq;
1149 0 : CoordinateSystem coordSys;
1150 0 : Array<Complex> cfBuf;
1151 : Float samplingFromMisc, paFromMisc;
1152 0 : if (Dir.length() == 0)
1153 0 : throw(SynthesisFTMachineError("Cache dir. name not set"));
1154 :
1155 0 : if (where < (Int)convFuncCache.nelements() && (!convFuncCache[where].data.null()))
1156 0 : return MEMCACHE;
1157 :
1158 0 : Int wConvSize, polInUse=2;
1159 0 : Int N=convFuncCache.nelements();
1160 :
1161 : //
1162 : // Re-size the conv. func. memory cache if required, and set the
1163 : // new members of the resized cache to NULL. This is used in the
1164 : // loop below to make a decision about allocating new memory or
1165 : // not.
1166 : //
1167 0 : convFuncCache.resize(max(where+1,N), true);
1168 : // for(Int i=N;i<=where;i++) convFuncCache[i].data=NULL;
1169 : //
1170 : // Each w-plan is in a separate disk file. Each file contains all
1171 : // polarization planes. Memory cache holds all w-planes and
1172 : // poln-planes in a single complex array. The loop below read
1173 : // each w-plane image from the disk, and fills in the 3D
1174 : // mem. cache for each computed PA.
1175 : //
1176 0 : wConvSize = Nw;
1177 0 : for(Int iw=0;iw<Nw;iw++)
1178 : {
1179 0 : ostringstream name;
1180 0 : name << Dir << "/" << cfPrefix << nameQualifier << iw << "_" << where;
1181 : try
1182 : {
1183 0 : PagedImage<Complex> tmp(name.str().c_str());
1184 0 : Record miscInfo;
1185 0 : miscInfo = tmp.miscInfo();
1186 :
1187 0 : miscInfo.get("Xsupport", xconvSupport);
1188 0 : miscInfo.get("Ysupport", yconvSupport);
1189 0 : miscInfo.get("sampling", samplingFromMisc);
1190 0 : miscInfo.get("ParallacticAngle", paFromMisc);
1191 0 : convSampling = samplingFromMisc;
1192 :
1193 :
1194 0 : Int index= tmp.coordinates().findCoordinate(Coordinate::SPECTRAL);
1195 0 : coordSys = tmp.coordinates();
1196 0 : SpectralCoordinate spCS = coordSys.spectralCoordinate(index);
1197 :
1198 : //cfRefFreq=spCS.referenceValue()(0);
1199 :
1200 0 : polInUse = tmp.shape()(2);
1201 0 : IPosition ts=tmp.shape(),ndx(4,0,0,0,0),ts2(4,0,0,0,0);
1202 0 : Array<Complex> imBuf=tmp.get();
1203 0 : if (convFuncCache[where].data.null())
1204 0 : cfBuf.resize(IPosition(4,ts(0),ts(1), wConvSize,polInUse));
1205 : // cfBuf = new CFType(IPosition(4,ts(0),ts(1), wConvSize,polInUse));
1206 :
1207 0 : ndx(CFDefs::NWPOS)=iw; // The w-axis
1208 0 : for(ndx(CFDefs::NPOLPOS)=0;ndx(CFDefs::NPOLPOS)<polInUse;ndx(CFDefs::NPOLPOS)++) // The Poln. axis.
1209 0 : for(ndx(CFDefs::NXPOS)=0;ndx(CFDefs::NXPOS)<ts(CFDefs::NXPOS);ndx(CFDefs::NXPOS)++)
1210 0 : for(ndx(CFDefs::NYPOS)=0;ndx(CFDefs::NYPOS)<ts(CFDefs::NYPOS);ndx(CFDefs::NYPOS)++)
1211 : {
1212 0 : ts2(CFDefs::NXPOS)=ndx(CFDefs::NXPOS);
1213 0 : ts2(CFDefs::NYPOS)=ndx(CFDefs::NYPOS);
1214 0 : ts2(2)=ndx(CFDefs::NPOLPOS); // The Poln. axis of the disk-cache. The LHS index is different!
1215 0 : ts2(3)=0; // The freq. axis of the disk-cache
1216 0 : (cfBuf)(ndx)=imBuf(ts2);
1217 : }
1218 0 : }
1219 0 : catch(AipsError &x)
1220 : {
1221 0 : throw(SynthesisFTMachineError(string("Error while loading \"")+
1222 0 : name.str( ) + string("\": ") + (string) x.getMesg()));
1223 0 : }
1224 0 : }
1225 : // xconvSupport.resize(wConvSize,true);
1226 : // yconvSupport.resize(wConvSize,true);
1227 : // for(Int i=0;i<wConvSize;i++)
1228 : // {
1229 : // xconvSupport(i) = XSup(i,where);
1230 : // yconvSupport(i) = YSup(i,where);
1231 : // }
1232 :
1233 : // convSampling = Sampling[where];
1234 :
1235 : // where=addToMemCache(convFuncCache, pa, &cfBuf, coordSys,
1236 : // xconvSupport, yconvSupport, Sampling[where]);
1237 0 : where=addToMemCache(convFuncCache, paFromMisc, &cfBuf, coordSys,
1238 : xconvSupport, yconvSupport, samplingFromMisc);
1239 0 : cfs=convFuncCache[where];
1240 : // convFuncCache[where].show("loadFromDisk: ");
1241 :
1242 0 : return DISKCACHE;
1243 0 : };
1244 : //
1245 : //-----------------------------------------------------------------------
1246 : //
1247 0 : Int CFCache::locateConvFunction(CFStore& cfs, CFStore& cfwts,
1248 : const Int Nw, const Float pa, const Float dPA,
1249 : const Int mosXPos, const Int mosYPos)
1250 : {
1251 : Int retVal;
1252 : // This assumes that the return state of locateConvFunction() for
1253 : // "CF" and "CFWT" will be the same. I.e. if "CF" is found in the
1254 : // cache, "CFWT" will be found. If "CF" is not found "CFWT" won't
1255 : // be found either.
1256 0 : retVal=locateConvFunction(cfs, Nw, pa, dPA,"",mosXPos,mosYPos);
1257 0 : locateConvFunction(cfwts, Nw, pa, dPA,"WT",mosXPos,mosYPos);
1258 :
1259 0 : return retVal;
1260 : }
1261 : //
1262 : //-----------------------------------------------------------------------
1263 : //
1264 : // Locate a convlution function in either mem. or disk cache.
1265 : // Return CFCache::DISKCACHE (=1) if found in the disk cache.
1266 : // CFCache::MEMCACHE (=2) if found in the mem. cache.
1267 : // <0 if not found in either cache. In this case, absolute of
1268 : // the return value corresponds to the index in the list of
1269 : // conv. funcs. where this conv. func. should be filled
1270 : //
1271 0 : Int CFCache::locateConvFunction(CFStore& cfs,
1272 : const Int Nw, const Float pa, const Float dPA,
1273 : const String& nameQualifier,
1274 : const Int /*mosXPos*/, const Int /*mosYPos*/)
1275 : {
1276 0 : LogIO log_l(LogOrigin("CFCache2", "locatedConvFunction"));
1277 :
1278 0 : Int paKey,retVal=NOTCACHED;
1279 : Bool found;
1280 : // Search for the PA corresponding to the supplied VB to find a
1281 : // paKey in memCache_p which has a Conv. Func. within dPA (dPA is
1282 : // given by paCD). If found, return the key in paKey.
1283 0 : found = searchConvFunction(paKey,pa, dPA);
1284 :
1285 0 : if (found)
1286 : {
1287 0 : CFStoreCacheType &memCacheObj=getMEMCacheObj(nameQualifier);
1288 0 : retVal=loadFromDisk(paKey,pa,dPA,Nw,memCacheObj,cfs,nameQualifier);
1289 :
1290 0 : switch (retVal)
1291 : {
1292 0 : case DISKCACHE:
1293 : {
1294 0 : if (paKey < (Int)memCacheObj.nelements())
1295 : log_l << "Loaded from disk cache: Conv. func. # "
1296 0 : << paKey << LogIO::POST;
1297 : else
1298 0 : throw(SynthesisFTMachineError("Internal error: paKey out of range"));
1299 0 : break;
1300 : }
1301 0 : case MEMCACHE:
1302 : {
1303 0 : cfs=(memCacheObj[paKey]);
1304 0 : break;
1305 : }
1306 0 : case NOTCACHED: {break;}
1307 : };
1308 : }
1309 0 : return retVal;
1310 0 : }
1311 : //
1312 : //-----------------------------------------------------------------------
1313 : //
1314 0 : CFCache::CFStoreCacheType& CFCache::getMEMCacheObj(const String& nameQualifier)
1315 : {
1316 0 : LogIO log_l(LogOrigin("CFCache::getMEMCacheObj"));
1317 0 : if (nameQualifier == "") return memCache_p;
1318 0 : else if (nameQualifier == "WT") return memCacheWt_p;
1319 : else
1320 0 : log_l << "Internal error. Unknown name qualifier '"+nameQualifier+"'."
1321 0 : << LogIO::EXCEPTION << endl;
1322 : //
1323 : // Return to suppress compiler warning. Control will never reach
1324 : // the following line.
1325 : //
1326 0 : return memCache_p;
1327 0 : }
1328 : //
1329 : //-----------------------------------------------------------------------
1330 : //
1331 : // void CFCache::constructTable_p(CFCache::CFCacheTable& tab)
1332 : // {
1333 : // TableDesc td("CFCache2","0.0",TableDesc::Scratch);
1334 :
1335 : // add PA and baseline type info.
1336 :
1337 : // td.addColumn(ScalarColumnDesc<Double>("Ref. Frequency"));
1338 : // td.addColumn(ScalarColumnDesc<Double>("W Value"));
1339 : // td.addColumn(ScalarColumnDesc<uInt>("Mueller Index"));
1340 : // td.addColumn(ScalarColumnDesc<String>("CF Disk File Name"));
1341 : // SetupNewTable newTab("CFCache.dat",td,Table::Scratch);
1342 : // tab = Table(newTab);
1343 : // }
1344 :
1345 : } // end casa namespace
|
