wproj3ddat_regression.py

Go to the documentation of this file.

import os
import time


# Enable benchmarking?
benchmarking = True
usedasync = False

# 
# Set up some useful variables
#
# Get path to CASA home directory by stipping name from '$CASAPATH'
pathname=os.environ.get('CASAPATH').split()[0]

# This is where the NGC5921 UVFITS data will be
fitsdata=pathname+'/data/demo/3DDAT.fits'

os.system('rm -rf n3ddat*')

# Start benchmarking
if benchmarking:
    startTime = time.time()
    startProc = time.clock()

#
#=====================================================================
#
# Import the data from FITS to MS
#
print '--Import--'

# Safest to start from task defaults
default('importuvfits')

prefix='n3ddat'
# Set up the MS filename and save as new global variable
msfile = prefix + '.ms'

# Use task importuvfits
fitsfile = fitsdata
vis = msfile
antnamescheme = 'new'
importuvfits()

# Note that there will be a ngc5921.ms.flagversions
# containing the initial flags as backup for the main ms flags.

# Record import time
if benchmarking:
    importtime = time.time()

tflagdata(vis=msfile, mode='manual', antenna='VA10', spw='*', field='0')

print '--Widefield--'
default('widefield')

# Pick up our split source data
vis = msfile

# Make an image root file name
imname = prefix + '.clean'
imagename = imname

# Set up the output image cube
mode = 'mfs'


# This is a single-source MS with one spw
field = '0'
spw = '*'

selectdata=True
#uvrange='0.15~10000 klambda'
# Standard gain factor 0.1
gain = 0.1

# Set the output image size and cell size (arcsec)
imsize = [3000,3000]

# Do a simple Clark clean
psfmode = 'clark'

cell = ['30arcsec','30arcsec']

# Fix maximum number of iterations
niter = 6000

# Also set flux residual threshold (in mJy)
threshold=0.0

# Set up the weighting
weighting = 'natural'

#ftmachine
ftmachine='wproject'
wprojplanes=200
cyclefactor=3

# No clean mask or cleanbox for now
mask = ''
#cleanbox = []

# But if you had a cleanbox saved in a file, e.g. "regionfile.txt"
# you could use it:
#cleanbox='regionfile.txt'
#
# and if you wanted to use interactive clean
#cleanbox='interactive'

widefield()

#
# Following is should be the equivalent clean task setup
#
#clean(vis="n3ddat.ms",imagename="ttt",outlierfile="",field="0",spw="*",selectdata=False,timerange="",uvrange="",antenna="",scan="",mode="mfs",gridmode="widefield",wprojplanes=200,facets=1,cfcache="cfcache.dir",painc=360.0,epjtable="",interpolation="nearest",niter=6000,gain=0.1,threshold=0.0,psfmode="clark",imagermode="csclean",ftmachine="mosaic",mosweight=False,scaletype="SAULT",multiscale=[],negcomponent=-1,smallscalebias=0.6,interactive=False,mask=[],nchan=1,start=0,width=1,imsize=[3000, 3000],cell=['30arcsec', '30arcsec'],phasecenter="",restfreq="",stokes="I",weighting="natural",robust=0.0,uvtaper=False,outertaper=[],innertaper=[],modelimage="",restoringbeam=[''],pbcor=False,minpb=0.1,calready=False,noise="1.0Jy",npixels=0,npercycle=100,cyclefactor=3,cyclespeedup=-1,nterms=1,reffreq="")



# Should find stuff in the logger like:
#
# Fitted beam used in restoration: 51.5643 by 45.6021 (arcsec)
#     at pa 14.5411 (deg)
#
# It will have made the images:
# -----------------------------
# ngc5921.clean.image
# ngc5921.clean.model
# ngc5921.clean.residual

clnimage = imname+'.image'

# Record clean completion time
if benchmarking:
    cleantime = time.time()

print '--Imstat--'
default('imstat')

imagename = clnimage

cubestats =imstat()

#
#=====================================================================
#
# Can do some image statistics if you wish
# Treat this like a regression script
# WARNING: currently requires toolkit
#
print ""
print ' 3DDAT results '
print ' =============== '

print ''
print ' --Regression Tests--'
print ''
# Pull the max from the cubestats dictionary
# created above using imstat
thistest_immax=cubestats['max'][0]
oldtest_immax = 7.50
print ' Clean image max should be ',oldtest_immax
print ' Found : Image Max = ',thistest_immax
diff_immax = abs((oldtest_immax-thistest_immax)/oldtest_immax)
print ' Difference (fractional) = ',diff_immax

print ''
# Pull the rms from the cubestats dictionary
thistest_imrms=cubestats['rms'][0]
oldtest_imrms = 0.134
print ' Clean image rms should be ',oldtest_imrms
print ' Found : Image rms = ',thistest_imrms
diff_imrms = abs((oldtest_imrms-thistest_imrms)/oldtest_imrms)
print ' Difference (fractional) = ',diff_imrms



# Record processing completion time
if benchmarking:
    endProc = time.clock()
    endTime = time.time()

#
#=====================================================================
#
if not benchmarking:
    print ''
    print '--- Done ---'
else:
    import datetime
    datestring=datetime.datetime.isoformat(datetime.datetime.today())
    outfile='wproj3ddat.'+datestring+'.log'
    logfile=open(outfile,'w')

    print >>logfile,''
    print >>logfile,''


    if (diff_immax < 0.05):
        passimmax = True
        print >>logfile,'* Passed image max test'
    else:
        passimmax = False
        print >>logfile,'* FAILED image max test'
    print >>logfile,'*  Image max '+str(thistest_immax)

    if (diff_imrms < 0.05):
        passimrms = True
        print >>logfile,'* Passed image rms test'
    else:
        passimrms = False
        print >>logfile,'* FAILED image rms test'
    print >>logfile,'*  Image rms '+str(thistest_imrms)

    if (passimmax and passimrms ): 
        regstate=True
        print >>logfile,'---'
        print >>logfile,'Passed Regression test for wproj3ddat'
        print >>logfile,'---'
    else: 
        regstate=False
        print >>logfile,'---'
        print >>logfile,'FAILED Regression test for wproj3ddat'
        print >>logfile,'---'
    print >>logfile,'*********************************'

    print >>logfile,''
    print >>logfile,''
    print >>logfile,'********* Benchmarking *****************'
    print >>logfile,'*                                      *'
    print >>logfile,'Total wall clock time was: ', endTime - startTime
    print >>logfile,'Total CPU        time was: ', endProc - startProc
    print >>logfile,'Processing rate MB/s  was: ', 35.1/(endTime - startTime)
    print >>logfile,'* Breakdown:                           *'
    print >>logfile,'*   import       time was: '+str(importtime-startTime)
    print >>logfile,'*   clean        time was: '+str(cleantime-importtime)
    print >>logfile,'*****************************************'
    print >>logfile,'basho (test cpu) time was: ?? seconds'

    logfile.close()

    print ""
    print "Done!"
        
#exit()
#
#=====================================================================