#====================================================================
#                                                                    
# Script for imaging 74MHz widefield data
#                                                                    
# Created STM 2009-06-24 from wproj3ddat_regression.py
# Version STM 2009-06-30 toolkit version for faceting
# Version STM 2009-07-02 use coma ms
#
# Script Notes:                                                         
#    o The results are written out as a dictionary in a pickle file     
#         <scriptprefix>.regression.<datestring>.pickle           
#      as well as in a text file                                        
#         out.<scriptprefix>.<datestring>.log                         
#      (these are not auto-deleted at start of script)                  
#    o This version does not regress against previous runs but
#      instead tests against the input component values
#                                                                       
#====================================================================
# IMPORTANT: VERSIONING HERE
#====================================================================
myscriptvers = '2008-07-02 STM'
mydataset = '74MHz Coma Widefield Simulation'

import os
import time
import pickle

# needed for existence check
from os import F_OK

# Enable benchmarking?
benchmarking = True
usedasync = False

# Prefix for output files
prefix='widefield_coma_simtest'
scriptprefix = 'run_widefield_coma_simtest'

# 
# 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 UVFITS data for template dataset will be, e.g.
# /home/ballista/casa/active/data/regression/wproject_regression/coma.ms
datapath=pathname+'/data/regression/wproject_regression/'
msname='coma.ms'
mspath=datapath+msname
# Path to web archive
webms = 'http://casa.nrao.edu/Data/VLA/Coma/'+msname+'.tgz'

# Clean up old files
os.system('rm -rf '+prefix+'*')

#
#=====================================================================
#
# Set up user parameters
# Natural beam is 172" x 147"

myfield  = "0"
myspw    = "0"
myimsize = 2000
mycell   = "30arcsec"
myniter  = 50
mythreshold = 0.0
myftmachine = "wproject"
mywprojplanes = 128
myfacets = 5
myweighting = 'natural'
myrobust = 0.0
mycyclefactor = 4
mycen    = "B1950 12:57:10.6940 +28.13.45.9620"

mygain = 0.1
mypadding = 1.2

print "Will image with imsize %5i and cell %s " % (myimsize,mycell)
print "Will clean "+str(myniter)+" iterations"
if mycen!="":
    print "Will use phasecenter "+mycen
if myftmachine!="":
    print "Will image using ftmachine "+myftmachine
if myftmachine=="wproject":
    print "Will use "+str(mywprojplanes)+" wprojection planes"
print "Will use "+str(myfacets)+" x "+str(myfacets)+" facets"
print "Will use cyclefactor = "+str(mycyclefactor)
print "Will use loop gain = "+str(mygain)
print "Will use facet padding = "+str(mypadding)

# Save the parameters used in clean etc.
params = {}
# Version of script
params['version'] = myscriptvers
# User set params
params['user'] = {}
params['user']['templatems'] = msname
params['user']['myfield'] = myfield
params['user']['myspw'] = myspw
params['user']['myimsize'] = myimsize
params['user']['mycell'] = mycell
params['user']['myniter'] = myniter
params['user']['mythreshold'] = mythreshold
params['user']['myftmachine'] = myftmachine
params['user']['mywprojplanes'] = mywprojplanes
params['user']['myfacets'] = myfacets
params['user']['myweighting'] = myweighting
params['user']['myrobust'] = myrobust
params['user']['mycyclefactor'] = mycyclefactor
params['user']['mycen'] = mycen
params['user']['mygain'] = mygain
params['user']['mypadding'] = mypadding

# Get center into a direction measure
phc = me.direction(mycen.split()[0],mycen.split()[1],mycen.split()[2])
# Reference pixel of image
myref = int((myimsize+0.5)/2.0)
dcell = qa.convert(mycell,'arcsec')['value']
ucell = qa.convert(mycell,'arcsec')['unit']

#
#=====================================================================
# Start benchmarking
#
#=====================================================================
#
# Get a copy of the MS
#
print '--Copy MS--'

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

if os.access(msname,F_OK):
    # ms in current directory
    print "  Using "+msname+" found in current directory"
    os.system("cp -rf "+msname+" "+msfile)
    os.system('chmod -R a+wx '+msfile)
elif os.access(mspath,F_OK):
    print "  Copying data from "+mspath
    os.system("cp -rf "+mspath+" "+msfile)
    os.system('chmod -R a+wx '+msfile)
else:
    if os.access(msname+'.tgz',F_OK):
        # ms tarball in current directory
        print "  Using "+msname+".tgz found in current directory"
    else:
        # try web retrieval
        print '--Retrieve data from '+webms
        # Use curl (for Mac compatibility)
        os.system('curl '+webms+' -f -o '+msname+'.tgz')
        # NOTE: could also use wget
        #os.system('wget '+webms)
        
    print '--Unpacking tarball '
    os.system('tar xzf '+msname+'.tgz')
    if os.access(msname,F_OK):
        # should now be in current directory
        print "  Using "+msname+" found in current directory"
    else:
        raise IOError," ERROR: "+msname+" not found"

    # ms in current directory
    os.system("cp -rf "+msname+" "+msfile)
    os.system('chmod -R a+wx '+msfile)

print "  Copied "+msname+" to "+msfile

# Use listobs to summarize the MS
print '--Listobs--'
listobs(msfile)

#================================================================================
#MeasurementSet Name:
#/home/sandrock3/smyers/Testing3/Patch4/Widefield/widefield_coma_simtest.ms
#MS Version 2
#================================================================================
#   Observer: AK462     Project:   
#Observation: VLA
#Data records: 607235       Total integration time = 905992 seconds
#   Observed from   22-Nov-1998/11:20:00.0   to   02-Dec-1998/22:59:52.5 (TAI)
#
#   ObservationID = 0         ArrayID = 0
#  Date        Timerange (TAI)          Scan  FldId FieldName    nVis   Int(s)   SpwIds
#  22-Nov-1998/11:20:00.0 - 12:44:37.5     1      0 COMA         35750  5.03     [0]
#              12:50:00.0 - 14:14:37.5     2      0 COMA         35649  5.03     [0]
#              14:20:00.0 - 16:02:37.5     3      0 COMA         40385  5.03     [0]
#              16:11:00.0 - 17:35:37.5     4      0 COMA         33316  5.03     [0]
#              17:41:07.5 - 17:54:15.0     5      0 COMA         5922   5.03     [0]
#
#   ObservationID = 0         ArrayID = 1
#  Date        Timerange (TAI)          Scan  FldId FieldName    nVis   Int(s)   SpwIds
#  27-Nov-1998/16:17:15.0 - 16:24:45.0     1      0 COMA         5406   5.03     [0]
#              16:36:15.0 - 16:47:15.0     2      0 COMA         6849   5.03     [0]
#              16:58:15.0 - 16:58:45.0     3      0 COMA         747    5.04     [0]
#              17:11:52.5 - 17:13:15.0     4      0 COMA         1307   5.03     [0]
#              17:24:52.5 - 17:36:15.0     5      0 COMA         7970   5.03     [0]
#              17:47:52.5 - 17:58:45.0     6      0 COMA         7405   5.03     [0]
#
#   ObservationID = 0         ArrayID = 2
#  Date        Timerange (TAI)          Scan  FldId FieldName    nVis   Int(s)   SpwIds
#  02-Dec-1998/12:36:30.0 - 13:34:20.0     1      0 COMA         42032  5.03     [0]
#              13:39:30.0 - 15:03:07.5     2      0 COMA         63174  5.03     [0]
#              15:17:30.0 - 17:51:15.0     3      0 COMA         116412 5.03     [0]
#              17:56:30.0 - 18:30:45.0     4      0 COMA         26007  5.03     [0]
#              18:38:07.5 - 18:43:15.0     5      0 COMA         5384   5.03     [0]
#              18:49:07.5 - 21:50:45.0     6      0 COMA         129111 5.03     [0]
#              21:56:30.0 - 22:59:52.5     7      0 COMA         44409  5.03     [0]
#           (nVis = Total number of time/baseline visibilities per scan) 
#Fields: 1
#  ID   Code Name         RA            Decl           Epoch   nVis   
#  0         COMA         12:57:10.6940 +28.13.45.9620 B1950   607235 
#   (nVis = Total number of time/baseline visibilities per field) 
#Spectral Windows:  (1 unique spectral windows and 1 unique polarization setups)
#  SpwID  #Chans Frame Ch1(MHz)    ChanWid(kHz)TotBW(kHz)  Ref(MHz)    Corrs   
#  0          18 TOPO  73.0675781  85.4492188  1538.086    73.0675781  RR  LL  
#The SOURCE table is empty: see the FIELD table
#Antennas: 28:
#  ID   Name  Station   Diam.    Long.         Lat.         
#  0    1     VLA:W16   25.0 m   -107.37.57.4  +33.53.33.0  
#  1    2     VLA:N16   25.0 m   -107.37.10.9  +33.54.48.0  
#  2    3     VLA:N12   25.0 m   -107.37.09.0  +33.54.30.0  
#  3    4     VLA:W8    25.0 m   -107.37.21.6  +33.53.53.0  
#  4    5     VLA:N2    25.0 m   -107.37.06.2  +33.54.03.5  
#  5    6     VLA:E10   25.0 m   -107.36.40.9  +33.53.52.0  
#  6    7     VLA:E2    25.0 m   -107.37.04.4  +33.54.01.1  
#  7    8     VLA:N8    25.0 m   -107.37.07.5  +33.54.15.8  
#  8    9     VLA:E8    25.0 m   -107.36.48.9  +33.53.55.1  
#  9    10    VLA:W4    25.0 m   -107.37.10.8  +33.53.59.1  
#  10   11    VLA:N6    25.0 m   -107.37.06.9  +33.54.10.3  
#  11   12    VLA:E16   25.0 m   -107.36.09.8  +33.53.40.0  
#  12   13    VLA:OUT   25.0 m   -107.37.06.0  +33.54.01.7  
#  13   14    VLA:E4    25.0 m   -107.37.00.8  +33.53.59.7  
#  14   15    VLA:N14   25.0 m   -107.37.09.9  +33.54.38.5  
#  15   16    VLA:W14   25.0 m   -107.37.46.9  +33.53.38.9  
#  16   17    VLA:W10   25.0 m   -107.37.28.9  +33.53.48.9  
#  17   18    VLA:W18   25.0 m   -107.38.08.9  +33.53.26.5  
#  18   19    VLA:N4    25.0 m   -107.37.06.5  +33.54.06.1  
#  19   20    VLA:E18   25.0 m   -107.35.57.2  +33.53.35.1  
#  20   21    VLA:E12   25.0 m   -107.36.31.7  +33.53.48.5  
#  21   22    VLA:W12   25.0 m   -107.37.37.4  +33.53.44.2  
#  22   23    VLA:W6    25.0 m   -107.37.15.6  +33.53.56.4  
#  23   24    VLA:N1    25.0 m   -107.37.06.0  +33.54.01.8  
#  24   25    VLA:W2    25.0 m   -107.37.07.4  +33.54.00.9  
#  25   26    VLA:E14   25.0 m   -107.36.21.3  +33.53.44.5  
#  26   27    VLA:N18   25.0 m   -107.37.12.0  +33.54.58.3  
#  27   28    VLA:E6    25.0 m   -107.36.55.6  +33.53.57.7  
#

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

#
#=====================================================================
#
# Set up simulation
#
print '--Set up simulation--'

# Because this is a widefield image, we need to first make a
# blank image from the ms to figure out where the components
# will be.  Uses imager tool for this (rather than clean/widefield).

tempimage = prefix+'.tempimg'
im.open(msfile)
im.defineimage(nx=myimsize, ny=myimsize, cellx=mycell, celly=mycell, mode='mfs', stokes='I', phasecenter=mycen)
im.make(tempimage)

# While we are here, make use of the advise function on facets etc.
fov = str(dcell*myimsize)+ucell
myadv = im.advise(takeadvice=F, amplitudeloss=0.01, fieldofview=fov)
print "  Imager advises facets="+str(myadv['facets'])+" imsize="+str(myadv['pixels'])

im.close()
print '  Created temp image '+tempimage

# Set up dictionary with components to simulated
# Pointing center:
# B1950 12:57:10.6940 +28.13.45.9620

# Canonical image is 2000x2000 at 30"
# Comp A is at center of image
complist = {}
complist['A'] = {}
complist['A']['rf'] = 'B1950'
complist['A']['RA'] = '12:57:10.6940'
complist['A']['Dec'] = '+28.13.45.9620' 
complist['A']['flux'] = 1.0

# Comp B is offset approx 6deg near 1594,1733
complist['B'] = {}
complist['B']['rf'] = 'B1950'
complist['B']['RA'] = '12:33:10.6940'
complist['B']['Dec'] = '+34.13.45.9620'
# To center on a pixel
#complist['B']['RA'] = '12:33:10.6940'
#complist['B']['Dec'] = '+34.13.45.9620'
complist['B']['flux'] = 1.0

# Which of these to use
compnames = ['A','B']
ncomp = compnames.__len__()

# Box around each component in image
# Note that pixel locations in outer parts of the image can be off due to
# geometry distortions if you just do simple offsets, so need to use
# image tool to convert world to pixel coordinates.
dbox = 50
ia.open(tempimage)
for comp in compnames:
    #
    offpix = ia.topixel([complist[comp]['RA'],complist[comp]['Dec']])
    xpixr = offpix['numeric'][0]
    ypixr = offpix['numeric'][1]
    # round off
    xpix = int(xpixr + 0.5)
    ypix = int(ypixr + 0.5)
    # store these in dictionary
    complist[comp]['xpixr'] = xpixr
    complist[comp]['ypixr'] = ypixr
    complist[comp]['xpix'] = xpix
    complist[comp]['ypix'] = ypix
    # now make box for clean and stats
    blcx = xpix - dbox
    blcy = ypix - dbox
    trcx = xpix + dbox
    trcy = ypix + dbox
    mybox = str(blcx)+","+ str(blcy)+","+ str(trcx)+","+ str(trcy)
    complist[comp]['box'] = mybox
    myclnbox = [blcx,blcy,trcx,trcy]
    complist[comp]['clnbox'] = myclnbox

ia.close()
# Set up component list
clfile = prefix + "_sim.cl"
print "  Creating componentlist "

for comp in compnames:
    mydirection = complist[comp]['rf']+" "+complist[comp]['RA']+" "+complist[comp]['Dec']
    myflux = complist[comp]['flux']
    cl.addcomponent(dir=mydirection, flux=myflux, freq='74.0MHz', spectrumtype='constant')
    #print "    Added comp "+comp+" "+str(myflux)+" Jy at "+mydirection+" ("+str(complist[comp]['xpix'])+","+str(complist[comp]['ypix'])+")"
    print "    Added comp %s with %12.6f Jy at %s ( %9.3f , %9.3f )" % ( comp, myflux, mydirection, complist[comp]['xpixr'], complist[comp]['ypixr'] )

cl.rename(clfile)
mycl = cl.torecord()
cl.close()

print "  Created componentlist "+clfile

#
#=====================================================================
#
# Use this MS as template for simulation of point sources
#
print '--Simulator--'

# Now simulate the data
sm.openfromms(msfile)
#sm.setvp(dovp=T, usedefaultvp=T, dosquint=F)
# Turn off primary beam for this widefield simulation
sm.setvp(dovp=F)
sm.predict(complist=clfile)
sm.close()

print "  Completed simulation in MS "+msfile

# Record time
if benchmarking:
    sim1time = time.time()

#
#=====================================================================
#
# Now image the data
#
dirtimag = prefix + ".dirty"
dirimage = dirtimag+'.image'
dirmodel = dirtimag+'.model'
dirresid = dirtimag+'.residual'

contimag = prefix + ".clean"
clnimage = contimag+'.image'
modimage = contimag+'.model'
resimage = contimag+'.residual'

print '--Imager (dirty)--'

# Set up imager
im.open(msfile)

# Imaging options
im.defineimage(nx=myimsize, ny=myimsize,
               cellx=mycell, celly=mycell,
               mode='mfs',
               stokes='I',
               phasecenter=mycen,
               facets=myfacets)

# Other options
im.setoptions(ftmachine=myftmachine, wprojplanes=mywprojplanes, padding=mypadding)

# Cycle controls
im.setmfcontrol(cyclefactor=mycyclefactor)

# No primary beam
im.setvp(dovp=F)

# Set up the weighting
if myweighting=='briggs':
    im.weight(myweighting,rmode='norm',robust=myrobust)
else:
    im.weight(myweighting)

print "First make dirty image"

im.clean('wfclark',
         image=dirimage,
         model=dirmodel,
         residual=dirresid,
         niter=0)

print "Output dirty image is "+dirimage

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

# Clean image
print '--Imager (clean)--'

# Set up clean boxes
clnbox = []
for comp in compnames: 
    clnbox.append(complist[comp]['clnbox'])

print "Will be cleaning using box(es) "+str(clnbox)

# Make mask image
maskimage = contimag+'.mymask'
im.make(maskimage)
im.regiontoimagemask(maskimage,boxes=clnbox)

print "Now actually clean the image using niter = "+str(myniter)

im.clean('wfclark',
         image=clnimage,
         mask=maskimage,
         model=modimage,
         residual=resimage,
         niter=myniter,
         gain=mygain,
         threshold=mythreshold)

print "Output clean restored image is "+clnimage
print "Output clean model image is "+modimage
print "Output clean residual image is "+resimage

im.close()

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

#
#=====================================================================
#
# Expected residuals based on number of iterations/comp and clean gain
# Will key tests to this
#
expectedres = (1.0-mygain)**(float(myniter)/float(ncomp))
print "   "
print "   Expected maximum residual = %8.2e " % expectedres
print "   "

#
#=====================================================================
#
# Some image statistics
#
print '--Imstat--'
default('imstat')

# Set up off-src box
# Offset from first comp
doffx = 150
doffy = 150
doffbox = 100
xpix = complist[compnames[0]]['xpix']
ypix = complist[compnames[0]]['ypix']
blcx = xpix + doffx - doffbox
blcy = ypix + doffy - doffbox
trcx = xpix + doffx + doffbox
trcy = ypix + doffy + doffbox
offbox = str(blcx)+","+ str(blcy)+","+ str(trcx)+","+ str(trcy)
print "Off-source stats in box "+offbox

# Do dirty image stats
dirstats = imstat(imagename=dirimage,box = '')
dirtyimg_max = dirstats['max'][0]
dirtyimg_rms = dirstats['sigma'][0]

diroffstats = imstat(imagename=dirimage,box=offbox)
dirtyoff_rms = diroffstats['sigma'][0]

# Do clean image
allstats = imstat(imagename=clnimage,box = '')
imageall_max = allstats['max'][0]
imageall_rms = allstats['sigma'][0]

offstats = imstat(imagename=clnimage,box=offbox)
imageoff_rms = offstats['sigma'][0]

# Do residual image
resstats = imstat(imagename=resimage,box = '')
residall_max = resstats['max'][0]
residall_rms = resstats['sigma'][0]

resoffstats = imstat(imagename=resimage,box=offbox)
residoff_rms = resoffstats['sigma'][0]

# Do whole model
modstats = imstat(imagename=modimage,box = '')
modelall_flux = modstats['sum'][0]

maxflux = 0.0
totalflux = 0.0
# Box around each component in image
for comp in compnames:
    myflux = complist[comp]['flux']
    totalflux = totalflux + myflux
    if myflux > maxflux:
        maxflux = myflux
    mybox = complist[comp]['box']
    print "Component "+comp+" stats in box "+mybox
    # Dirty image
    xstat_dirimg = imstat(imagename=dirimage,box=mybox)
    # Clean image
    xstat_img = imstat(imagename=clnimage,box=mybox)
    # Clean model
    xstat_mod = imstat(imagename=modimage,box=mybox)
    # Residual image
    xstat_res = imstat(imagename=resimage,box=mybox)
    # Save in complist
    complist[comp]['dirtystat'] = xstat_dirimg
    complist[comp]['imstat'] = xstat_img
    complist[comp]['modstat'] = xstat_mod
    complist[comp]['resstat'] = xstat_res

stats1time=time.time()

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

walltime = (endTime - startTime)
cputime = (endProc - startProc)
#
#=====================================================================
# Print results
#=====================================================================
#
import datetime
datestring=datetime.datetime.isoformat(datetime.datetime.today())

myvers = casalog.version()
myuser = os.getenv('USER')
myhost = str( os.getenv('HOST') )
mycwd = os.getcwd()
myos = os.uname()
mypath = os.environ.get('CASAPATH')

# Print to terminal
print 'Running '+myvers+' on host '+myhost
print '  at '+datestring
print '  using '+mypath
#
# Save in logfile
outfile = 'out.'+prefix+'.'+datestring+'.log'
logfile=open(outfile,'w')

print >>logfile,'Running '+myvers+' on host '+myhost
print >>logfile,'  at '+datestring
print >>logfile,'  using '+mypath

#
##########################################################################
# Write some more stuff to logfile
print >>logfile,''
print >>logfile,'---'
print >>logfile,'Script version: '+params['version']
print >>logfile,'User set parameters used in execution:'
print >>logfile,'---'
print ''
print '---'
print 'Script version: '+params['version']
print 'User set parameters used in execution:'
print '---'

for keys in params['user'].keys():
    print >>logfile,'  %s  =  %s ' % ( keys, params['user'][keys] )
    print '  %s  =  %s ' % ( keys, params['user'][keys] )

print >>logfile,''
print ''

#
##########################################################################
#
# Save info in regression dictionary
new_regression = {}

new_regression['date'] = datestring
new_regression['version'] = myvers
new_regression['user'] = myuser
new_regression['host'] = myhost
new_regression['cwd'] = mycwd
new_regression['os'] = myos
new_regression['aipspath'] = mypath

new_regression['dataset'] = mydataset

# Store the parameters used in clean etc.
new_regression['params'] = params

# Store the results of this regression
results = {}

# First the simulation values
op = 'none'
tol = 0.01
results['sim_comp_max'] = {}
results['sim_comp_max']['name'] = 'Sim Model max flux'
results['sim_comp_max']['value'] = maxflux
results['sim_comp_max']['op'] = op
results['sim_comp_max']['tol'] = tol
results['sim_comp_max']['test'] = 'none'
results['sim_comp_max']['testval'] = 0.0

results['sim_comp_flux'] = {}
results['sim_comp_flux']['name'] = 'Sim Model total flux'
results['sim_comp_flux']['value'] = totalflux
results['sim_comp_flux']['op'] = op
results['sim_comp_flux']['tol'] = tol
results['sim_comp_flux']['test'] = 'none'
results['sim_comp_flux']['testval'] = 0.0

# Now the results
op = 'divf'
tol = 0.10
results['clean_dirtyimg_max'] = {}
results['clean_dirtyimg_max']['name'] = 'Dirty image (all) max'
results['clean_dirtyimg_max']['value'] = dirtyimg_max
results['clean_dirtyimg_max']['op'] = op
results['clean_dirtyimg_max']['tol'] = tol
results['clean_dirtyimg_max']['test'] = 'sim_comp_max'
results['clean_dirtyimg_max']['testval'] = results['sim_comp_max']['value']

op = 'none'
tol = 0.0
results['clean_dirtyimg_rms'] = {}
results['clean_dirtyimg_rms']['name'] = 'Dirty image (all) rms'
results['clean_dirtyimg_rms']['value'] = dirtyimg_rms
results['clean_dirtyimg_rms']['op'] = op
results['clean_dirtyimg_rms']['tol'] = tol
results['clean_dirtyimg_rms']['test'] = 'none'
results['clean_dirtyimg_rms']['test'] = 0.0

op = 'logf'
tol = 0.3
results['clean_dirtyoff_rms'] = {}
results['clean_dirtyoff_rms']['name'] = 'Dirty image (off-src) rms'
results['clean_dirtyoff_rms']['value'] = dirtyoff_rms
results['clean_dirtyoff_rms']['op'] = op
results['clean_dirtyoff_rms']['tol'] = tol
results['clean_dirtyoff_rms']['test'] = 'value'
results['clean_dirtyoff_rms']['testval'] = expectedres

op = 'divf'
tol = 0.02
results['clean_imageall_max'] = {}
results['clean_imageall_max']['name'] = 'Clean image (all) max'
results['clean_imageall_max']['value'] = imageall_max
results['clean_imageall_max']['op'] = op
results['clean_imageall_max']['tol'] = tol
results['clean_imageall_max']['test'] = 'sim_comp_max'
results['clean_imageall_max']['testval'] = results['sim_comp_max']['value']

op = 'none'
tol = 0.0
results['clean_imageall_rms'] = {}
results['clean_imageall_rms']['name'] = 'Clean image (all) rms'
results['clean_imageall_rms']['value'] = imageall_rms
results['clean_imageall_rms']['op'] = op
results['clean_imageall_rms']['tol'] = tol
results['clean_imageall_rms']['test'] = 'none'
results['clean_imageall_rms']['testval'] = 0.0

op = 'logf'
tol = 0.3
results['clean_imageoff_rms'] = {}
results['clean_imageoff_rms']['name'] = 'Clean image (off-src) rms'
results['clean_imageoff_rms']['value'] = imageoff_rms
results['clean_imageoff_rms']['op'] = op
results['clean_imageoff_rms']['tol'] = tol
results['clean_imageoff_rms']['test'] = 'value'
results['clean_imageoff_rms']['testval'] = expectedres

results['clean_residall_max'] = {}
results['clean_residall_max']['name'] = 'Residual image (all) max'
results['clean_residall_max']['value'] = residall_max
results['clean_residall_max']['op'] = op
results['clean_residall_max']['tol'] = tol
results['clean_residall_max']['test'] = 'value'
results['clean_residall_max']['testval'] = expectedres

results['clean_residall_rms'] = {}
results['clean_residall_rms']['name'] = 'Residual image (all) rms'
results['clean_residall_rms']['value'] = residall_rms
results['clean_residall_rms']['op'] = op
results['clean_residall_rms']['tol'] = tol
results['clean_residall_rms']['test'] = 'value'
results['clean_residall_rms']['testval'] = expectedres

results['clean_residoff_rms'] = {}
results['clean_residoff_rms']['name'] = 'Residual image (off-src) rms'
results['clean_residoff_rms']['value'] = residoff_rms
results['clean_residoff_rms']['op'] = op
results['clean_residoff_rms']['tol'] = tol
results['clean_residoff_rms']['test'] = 'value'
results['clean_residoff_rms']['testval'] = expectedres

op = 'divf'
tol = expectedres
results['clean_modelall_flux'] = {}
results['clean_modelall_flux']['name'] = 'Model image (all) flux'
results['clean_modelall_flux']['value'] = modelall_flux
results['clean_modelall_flux']['op'] = op
results['clean_modelall_flux']['tol'] = tol
results['clean_modelall_flux']['test'] = 'sim_comp_flux'
results['clean_modelall_flux']['testval'] = results['sim_comp_flux']['value']

# Done filling results
new_regression['results'] = results

# Now the component results
new_regression['compnames'] = compnames
#Store raw statistics
new_regression['compstats'] = complist

compresults = {}
for comp in compnames:
    compresults[comp] = {}
    
    xstat_dirimg = complist[comp]['dirtystat']
    xstat_img = complist[comp]['imstat'] 
    xstat_mod = complist[comp]['modstat']
    xstat_res = complist[comp]['resstat']

    compresults[comp]['box'] = complist[comp]['box']

    # First the simulated input values
    op = 'none'
    tol = 0.0
    compresults[comp]['sim_comp_flux'] = {}
    compresults[comp]['sim_comp_flux']['name'] = 'Simulated comp flux'
    compresults[comp]['sim_comp_flux']['value'] = complist[comp]['flux']
    compresults[comp]['sim_comp_flux']['op'] = op
    compresults[comp]['sim_comp_flux']['tol'] = tol
    compresults[comp]['sim_comp_flux']['test'] = 'none'
    compresults[comp]['sim_comp_flux']['testval'] = 0.0

    compresults[comp]['sim_comp_posx'] = {}
    compresults[comp]['sim_comp_posx']['name'] = 'Simulated comp pos x'
    compresults[comp]['sim_comp_posx']['value'] = complist[comp]['xpix']
    compresults[comp]['sim_comp_posx']['op'] = op
    compresults[comp]['sim_comp_posx']['tol'] = tol
    compresults[comp]['sim_comp_posx']['test'] = 'none'
    compresults[comp]['sim_comp_posx']['testval'] = 0.0
    
    compresults[comp]['sim_comp_posy'] = {}
    compresults[comp]['sim_comp_posy']['name'] = 'Simulated comp pos y'
    compresults[comp]['sim_comp_posy']['value'] = complist[comp]['ypix']
    compresults[comp]['sim_comp_posy']['op'] = op
    compresults[comp]['sim_comp_posy']['tol'] = tol
    compresults[comp]['sim_comp_posy']['test'] = 'none'
    compresults[comp]['sim_comp_posy']['testval'] = 0.0

    # Now what we found
    op = 'divf'
    tol = 0.10
    compresults[comp]['dirty_image_max'] = {}
    compresults[comp]['dirty_image_max']['name'] = 'Dirty image max'
    compresults[comp]['dirty_image_max']['value'] = xstat_dirimg['max'][0]
    compresults[comp]['dirty_image_max']['op'] = op
    compresults[comp]['dirty_image_max']['tol'] = tol
    compresults[comp]['dirty_image_max']['test'] = 'sim_comp_flux'
    compresults[comp]['dirty_image_max']['testval'] = compresults[comp]['sim_comp_flux']['value']
    
    op = 'diff'
    tol = 1.0
    compresults[comp]['dirty_image_maxposx'] = {}
    compresults[comp]['dirty_image_maxposx']['name'] = 'Dirty image max pos x'
    compresults[comp]['dirty_image_maxposx']['value'] = xstat_dirimg['maxpos'][0]
    compresults[comp]['dirty_image_maxposx']['op'] = op
    compresults[comp]['dirty_image_maxposx']['tol'] = tol
    compresults[comp]['dirty_image_maxposx']['test'] = 'sim_comp_posx'
    compresults[comp]['dirty_image_maxposx']['testval'] = compresults[comp]['sim_comp_posx']['value']
    
    compresults[comp]['dirty_image_maxposy'] = {}
    compresults[comp]['dirty_image_maxposy']['name'] = 'Dirty image max pos y'
    compresults[comp]['dirty_image_maxposy']['value'] = xstat_dirimg['maxpos'][1]
    compresults[comp]['dirty_image_maxposy']['op'] = op
    compresults[comp]['dirty_image_maxposy']['tol'] = tol
    compresults[comp]['dirty_image_maxposy']['test'] = 'sim_comp_posy'
    compresults[comp]['dirty_image_maxposy']['testval'] = compresults[comp]['sim_comp_posy']['value']
    
    op = 'divf'
    tol = 0.02
    compresults[comp]['clean_image_max'] = {}
    compresults[comp]['clean_image_max']['name'] = 'Restored image max'
    compresults[comp]['clean_image_max']['value'] = xstat_img['max'][0]
    compresults[comp]['clean_image_max']['op'] = op
    compresults[comp]['clean_image_max']['tol'] = tol
    compresults[comp]['clean_image_max']['test'] = 'sim_comp_flux'
    compresults[comp]['clean_image_max']['testval'] = compresults[comp]['sim_comp_flux']['value']

    op = 'diff'
    tol = 1.0
    compresults[comp]['clean_image_maxposx'] = {}
    compresults[comp]['clean_image_maxposx']['name'] = 'Restored image max pos x'
    compresults[comp]['clean_image_maxposx']['value'] = xstat_img['maxpos'][0]
    compresults[comp]['clean_image_maxposx']['op'] = op
    compresults[comp]['clean_image_maxposx']['tol'] = tol
    compresults[comp]['clean_image_maxposx']['test'] = 'sim_comp_posx'
    compresults[comp]['clean_image_maxposx']['testval'] = compresults[comp]['sim_comp_posx']['value']
    
    compresults[comp]['clean_image_maxposy'] = {}
    compresults[comp]['clean_image_maxposy']['name'] = 'Restored image max pos y'
    compresults[comp]['clean_image_maxposy']['value'] = xstat_img['maxpos'][1]
    compresults[comp]['clean_image_maxposy']['op'] = op
    compresults[comp]['clean_image_maxposy']['tol'] = tol
    compresults[comp]['clean_image_maxposy']['test'] = 'sim_comp_posy'
    compresults[comp]['clean_image_maxposy']['testval'] = compresults[comp]['sim_comp_posy']['value']
    
    op = 'logf'
    tol = 0.3
    compresults[comp]['residual_image_max'] = {}
    compresults[comp]['residual_image_max']['name'] = 'Residual image max'
    compresults[comp]['residual_image_max']['value'] = xstat_res['max'][0]
    compresults[comp]['residual_image_max']['op'] = op
    compresults[comp]['residual_image_max']['tol'] = tol
    compresults[comp]['residual_image_max']['test'] = 'value'
    compresults[comp]['residual_image_max']['testval'] = expectedres
    
    compresults[comp]['residual_image_rms'] = {}
    compresults[comp]['residual_image_rms']['name'] = 'Residual image rms'
    compresults[comp]['residual_image_rms']['value'] = xstat_res['sigma'][0]
    compresults[comp]['residual_image_rms']['op'] = op
    compresults[comp]['residual_image_rms']['tol'] = tol
    compresults[comp]['residual_image_rms']['test'] = 'value'
    compresults[comp]['residual_image_rms']['testval'] = expectedres

    op = 'divf'
    tol = expectedres
    compresults[comp]['model_image_sum'] = {}
    compresults[comp]['model_image_sum']['name'] = 'Model image flux'
    compresults[comp]['model_image_sum']['value'] = xstat_mod['sum'][0]
    compresults[comp]['model_image_sum']['op'] = op
    compresults[comp]['model_image_sum']['tol'] = tol
    compresults[comp]['model_image_sum']['test'] = 'sim_comp_flux'
    compresults[comp]['model_image_sum']['testval'] = compresults[comp]['sim_comp_flux']['value']

new_regression['compresults'] = compresults

# Dataset size info
datasize_raw = 1400.0 # MB
datasize_ms =  1400.0 # MB
new_regression['datasize'] = {}
new_regression['datasize']['raw'] = datasize_raw
new_regression['datasize']['ms'] = datasize_ms

# Save timing to regression dictionary
new_regression['timing'] = {}

total = {}
total['wall'] = (endTime - startTime)
total['cpu'] = (endProc - startProc)
total['rate_raw'] = (datasize_raw/(endTime - startTime))
total['rate_ms'] = (datasize_ms/(endTime - startTime))

new_regression['timing']['total'] = total

nstages = 4
new_regression['timing']['nstages'] = nstages

stages = {}
stages[0] = ['simulate',(sim1time-startTime)]
stages[1] = ['invert',(dirty1time-sim1time)]
stages[2] = ['clean',(clean1time-dirty1time)]
stages[3] = ['stats',(stats1time-clean1time)]

new_regression['timing']['stages'] = stages

#
##########################################################################
# Now print results
resultlist = ['clean_dirtyimg_max','clean_imageall_max','clean_imageoff_rms',
              'clean_residall_max','clean_residall_rms','clean_residoff_rms',
              'clean_modelall_flux']

compreslist = ['dirty_image_max','dirty_image_maxposx','dirty_image_maxposy',
              'clean_image_max','clean_image_maxposx','clean_image_maxposy',
              'residual_image_max','residual_image_rms','model_image_sum']

testresults = {}
testcompresults = {}
final_status = 'Passed'

print >>logfile,'---'
print >>logfile,'Results of imaging:'
print >>logfile,'---'
print '---'
print 'Results of imaging:'
print '---'

# Do the component results
for comp in compnames:
    print ""
    print " Component "+comp+" :"
    print >>logfile,""
    print >>logfile," Component "+comp+" :"

    # Set up storage for regression results
    testcompresults[comp] = {}
        
    for keys in compreslist:
        res = compresults[comp][keys]
        new_val = res['value']
        testres = {}
        if res['op']!='none':
           # Test versus another value
            testres['value'] = new_val
            testres['op'] = res['op']
            testres['tol'] = res['tol']
            testres['test'] = res['test']
            prev_val = res['testval']
            if res['op'] == 'divf':
                new_diff = (new_val - prev_val)/prev_val
            elif res['op'] == 'logf':
                new_diff = log10(new_val) - log10(prev_val)
            else:
                new_diff = new_val - prev_val
            
            if abs(new_diff)>res['tol']:
                new_status = 'Failed'
            else:
                new_status = 'Passed'
            
            testres['prev'] = prev_val
            testres['diff'] = new_diff
            testres['status'] = new_status
            
            # Save results
            testcompresults[comp][keys] = testres
                            
            # Print results
            print '--%30s : %12.6f expected %12.6f %4s %12.6f (%6s vs %12.6f)' % ( res['name'], testres['value'], testres['prev'], testres['op'], testres['diff'], testres['status'], testres['tol'] )
            print >>logfile,'--%30s : %12.6f expected %12.6f %4s %12.6f (%6s vs %12.6f)' % ( res['name'], testres['value'], testres['prev'], testres['op'], testres['diff'], testres['status'], testres['tol'] )
                
            if testres['status']=='Failed':
                final_status = 'Failed'
                
        else:
            # Just print results
            print '--%30s : %12.6f ' % ( res['name'], res['value'] )
            print >>logfile,'--%30s : %12.6f ' % ( res['name'], res['value'] )

# Do the overall results
print ""
print " Whole image :"
print >>logfile,""
print >>logfile," Whole image :"
for keys in resultlist:
    res = results[keys]
    new_val = res['value']
    testkey = res['test']
    if res['op']!='none':
        # Test versus another value
        testresults[keys] = {}
        testres = {}
        testres['value'] = new_val
        testres['op'] = res['op']
        testres['tol'] = res['tol']
        testres['test'] = res['test']
        prev_val = res['testval']
        if res['op'] == 'divf':
            new_diff = (new_val - prev_val)/prev_val
        elif res['op'] == 'logf':
            new_diff = log10(new_val) - log10(prev_val)
        else:
            new_diff = new_val - prev_val
        
        if abs(new_diff)>res['tol']:
            new_status = 'Failed'
        else:
            new_status = 'Passed'
            
        testres['prev'] = prev_val
        testres['diff'] = new_diff
        testres['status'] = new_status

        # Save results
        testresults[keys] = testres

        # Print results
        print '--%30s : %12.6f expected %12.6f %4s %12.6f (%6s vs %12.6f)' % ( res['name'], testres['value'], testres['prev'], testres['op'], testres['diff'], testres['status'], testres['tol'] )
        print >>logfile,'--%30s : %12.6f expected %12.6f %4s %12.6f (%6s vs %12.6f)' % ( res['name'], testres['value'], testres['prev'], testres['op'], testres['diff'], testres['status'], testres['tol'] )
        if testres['status']=='Failed':
            final_status = 'Failed'

    else:
        # Just print the current results
        print '--%30s : %12.6f ' % ( res['name'], res['value'] )
        print >>logfile,'--%30s : %12.6f ' % ( res['name'], res['value'] )
    
# Done
                    
# Final regression status
if (final_status == 'Passed'):
    regstate=True
    print >>logfile,'---'
    print >>logfile,'Passed Test for '+mydataset
    print >>logfile,'---'
    print 'Passed Test for '+mydataset
else:
    regstate=False
    print >>logfile,'----FAILED Test for '+mydataset
    print '----FAILED Test for '+mydataset
    
# Write these regression test results to the dictionary
new_regression['testresults'] = testresults
new_regression['testcompresults'] = testcompresults

#
##########################################################################
# Final stats and timing

print ''
print '********* Benchmarking *************************'
print '*                                              *'
print 'Total wall clock time was: %10.3f ' % total['wall']
print 'Total CPU        time was: %10.3f ' % total['cpu']
print 'Raw processing rate MB/s was: %10.3f ' % total['rate_raw']
print 'MS  processing rate MB/s was: %10.3f ' % total['rate_ms']
print ''
print '* Breakdown:                                   *'

print >>logfile,''
print >>logfile,'********* Benchmarking *************************'
print >>logfile,'*                                              *'
print >>logfile,'Total wall clock time was: %10.3f ' % total['wall']
print >>logfile,'Total CPU        time was: %10.3f ' % total['cpu']
print >>logfile,'Raw processing rate MB/s was: %10.3f ' % total['rate_raw']
print >>logfile,'MS  processing rate MB/s was: %10.3f ' % total['rate_ms']
print >>logfile,''
print >>logfile,'* Breakdown:                                   *'

for i in range(nstages):
    print '* %16s * time was: %10.3f ' % tuple(stages[i])
    print >>logfile,'* %16s * time was: %10.3f ' % tuple(stages[i])

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

#
##########################################################################
#
# Save regression results as dictionary using Pickle
#
pickfile = 'out.'+prefix + '.regression.'+datestring+'.pickle'
f = open(pickfile,'w')
p = pickle.Pickler(f)
p.dump(new_regression)
f.close()

print ""
print "Regression result dictionary saved in "+pickfile
print ""
print "Use Pickle to retrieve these"
print ""

# e.g.
# f = open(pickfile)
# u = pickle.Unpickler(f)
# clnmodel = u.load()
# polmodel = u.load()
# f.close()

print >>logfile,""
print >>logfile,"Done with "+mydataset

logfile.close()
print "Results are in "+outfile

print ""
print "Done with "+mydataset

#exit()
#
#=====================================================================