#====================================================================
#                                                                    
# Script for imaging 74MHz widefield data
#                                                                    
# Created STM 2009-06-24 from wproj3ddat_regression.py
#
# 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-06-29 STM'
mydataset = '74MHz 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_simtest'
scriptprefix = 'run_widefield_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/demo/3DDAT.fits
fitsdata=pathname+'/data/demo/3DDAT.fits'

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

#
#=====================================================================
#
# Set up user parameters
myfield  = "0"
myspw    = "0"
myimsize = 3000
mycell   = "30arcsec"
myniter  = 300
mythreshold = 0.0
myftmachine = "wproject"
mywprojplanes = 256
myfacets = 3
myweighting = 'natural'
mycyclefactor = 3
mycen    = "J2000 16:00:00.0000 +50.00.00.0000"
mygain = 0.1

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)

# Save the parameters used in clean etc.
params = {}
# Version of script
params['version'] = myscriptvers
# User set params
params['user'] = {}
params['user']['templatems'] = fitsdata
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']['mycyclefactor'] = mycyclefactor
params['user']['mycen'] = mycen
params['user']['mygain'] = mygain

# 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
if benchmarking:
    startTime = time.time()
    startProc = time.clock()

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

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

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

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

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

# Now some flagging
print '--Flagdata--'
flagdata(vis=msfile, mode='manualflag', antenna='VA10', spw='*', field='0')

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

#
#=====================================================================
#
# 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
# Comp A is at pixel 1734,1336 of 3000x3000 image
# This was the max in the original image
# Note that only comp A is centered on a pixel by fiat
complist = {}
complist['A'] = {}
complist['A']['rf'] = 'J2000'
complist['A']['RA'] = '15:48:12.077'
complist['A']['Dec'] = '+48.35.41.038' 
complist['A']['flux'] = 1.0

# Comp B is at 1419,1436
complist['B'] = {}
complist['B']['rf'] = 'J2000'
#complist['B']['RA'] = '16:04:07.904'
#complist['B']['Dec'] = '+49.27.57.066'
# To center on a pixel
complist['B']['RA'] = '16:04:09.263'
complist['B']['Dec'] = '+49.27.43.10'
complist['B']['flux'] = 1.0

# Comp C is near 1575,658
complist['C'] = {}
complist['C']['rf'] = 'J2000'
#complist['C']['RA'] = '15:56:34.204'
#complist['C']['Dec'] = '+42.57.50.469'
# To center on a pixel
complist['C']['RA'] = '15:56:35.020'
complist['C']['Dec'] = '+42.57.43.487'
complist['C']['flux'] = 1.0

# Comp D is near 235,2382 
complist['D'] = {}
complist['D']['rf'] = 'J2000'
#complist['D']['RA'] = '17:16:45.708'
#complist['D']['Dec'] = '+55.59.19.802'
# To center on a pixel
complist['D']['RA'] = '17:16:49.593'
complist['D']['Dec'] = '+55.59.40.561'
complist['D']['flux'] = 1.0

# Which of these to use
compnames = ['A','B','C','D']
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 = 20
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 = "J2000 "+complist[comp]['RA']+" "+complist[comp]['Dec']
    myflux = complist[comp]['flux']
    cl.addcomponent(dir=mydirection, flux=myflux, freq='230.0GHz', 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

# Use listobs to summarize the MS
print '--Listobs--'
listobs(msfile)
#
#=====================================================================
# MeasurementSet Name:
# /home/sandrock3/smyers/Testing3/Patch4/Widefield/widefield_regression.ms
# MS Version 2
#=====================================================================
#    Observer: APLOW     Project:   
# Observation: VLA
# Data records: 29700       Total integration time = 5790 seconds
#    Observed from   12-Oct-1998/20:59:10.0   to   12-Oct-1998/22:35:40.0 (TAI)
# 
#    ObservationID = 0         ArrayID = 0
#   Date        Timerange (TAI)          Scan  FldId FieldName    nVis   Int(s)   SpwIds
#   12-Oct-1998/20:59:10.0 - 21:04:30.0     1      0 POSA         9900   10       [0]
#               21:45:00.0 - 21:50:20.0     2      0 POSA         9900   10       [0]
#               22:30:20.0 - 22:35:40.0     3      0 POSA         9900   10       [0]
#            (nVis = Total number of time/baseline visibilities per scan) 
# Fields: 1
#   ID   Code Name         RA            Decl           Epoch   nVis   
#   0         POSA         16:00:00.0000 +50.00.00.0000 J2000   29700  
#    (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          12 TOPO  73.1347168  122.070312  1464.84375  73.8        RR  LL  
# The SOURCE table is empty: see the FIELD table
# Antennas: 25:
#   ID   Name  Station   Diam.    Long.         Lat.         
#   0    VA01  VLA:W16   25.0 m   -107.37.57.4  +33.53.33.0  
#   1    VA02  VLA:N16   25.0 m   -107.37.10.9  +33.54.48.0  
#   2    VA03  VLA:N12   25.0 m   -107.37.09.0  +33.54.30.0  
#   3    VA04  VLA:W8    25.0 m   -107.37.21.6  +33.53.53.0  
#   5    VA06  VLA:E28   25.0 m   -107.34.39.3  +33.53.04.9  
#   6    VA07  VLA:E24   25.0 m   -107.35.13.4  +33.53.18.1  
#   7    VA08  VLA:N8    25.0 m   -107.37.07.5  +33.54.15.8  
#   8    VA09  VLA:E8    25.0 m   -107.36.48.9  +33.53.55.1  
#   9    VA10  VLA:W4    25.0 m   -107.37.10.8  +33.53.59.1  
#   10   VA11  VLA:N24   25.0 m   -107.37.16.1  +33.55.37.7  
#   11   VA12  VLA:E16   25.0 m   -107.36.09.8  +33.53.40.0  
#   12   VA13  VLA:W36   25.0 m   -107.40.32.6  +33.52.06.0  
#   14   VA15  VLA:N28   25.0 m   -107.37.18.7  +33.56.02.5  
#   15   VA16  VLA:W32   25.0 m   -107.39.54.8  +33.52.27.2  
#   16   VA17  VLA:W28   25.0 m   -107.39.20.2  +33.52.46.6  
#   17   VA18  VLA:E20   25.0 m   -107.35.43.6  +33.53.29.9  
#   18   VA19  VLA:N4    25.0 m   -107.37.06.5  +33.54.06.1  
#   19   VA20  VLA:E36   25.0 m   -107.33.20.2  +33.52.34.3  
#   20   VA21  VLA:E12   25.0 m   -107.36.31.7  +33.53.48.5  
#   21   VA22  VLA:W12   25.0 m   -107.37.37.4  +33.53.44.2  
#   22   VA23  VLA:W24   25.0 m   -107.38.49.0  +33.53.04.0  
#   23   VA24  VLA:N32   25.0 m   -107.37.22.0  +33.56.33.6  
#   24   VA25  VLA:W20   25.0 m   -107.38.21.4  +33.53.19.5  
#   25   VA26  VLA:E32   25.0 m   -107.34.01.5  +33.52.50.3  
#   27   VA28  VLA:E4    25.0 m   -107.37.00.8  +33.53.59.7  
#=====================================================================

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

#
#=====================================================================
#
# Now image the data
#
print '--Widefield--'
default('widefield')

# Pick up our split source data
vis = msfile

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

# This is a single-source MS with one spw
field = myfield
spw = myspw

# Set center explicitly
phasecenter =  mycen

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

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

# Do a Clark clean during minor cycles
psfmode = 'clark'

cell = [mycell,mycell]

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

# Set up the weighting
weighting = myweighting

# Set up cycles, ftmachine, wprojection, facets, etc.
ftmachine=myftmachine
wprojplanes=mywprojplanes
facets=myfacets
cyclefactor=mycyclefactor

# No clean mask for now
mask = ''

print "First make dirty image"
dirtimag = prefix + ".dirty"
imagename          =  dirtimag
niter              =  0
widefield()

print "Output dirty images are called "+imagename

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

# Set maximum number of iterations
niter = myniter

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

# Set clean mask using boxes
mask = clnbox
print "Will be cleaning using box(es) "+str(mask)

#Clean image
contimag = prefix + ".clean"
imagename          =  contimag
niter              =  myniter
print "Now actually clean the image using niter = "+str(niter)

widefield()

print "Output clean images are called "+imagename

dirimage = dirtimag+'.image'
clnimage = contimag+'.image'
#modimage = contimag+'.model'
# 2.4.0 still has bug that model image has no postfix
modimage = contimag
resimage = contimag+'.residual'

# 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 = 50
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 = 1.0
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 = 1.0
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 = 1.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 = 1.0
    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 =  9.7 # MB
datasize_ms =  25.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 = 6
new_regression['timing']['nstages'] = nstages

stages = {}
stages[0] = ['import',(import1time-startTime)]
stages[1] = ['flag',(flag1time-import1time)]
stages[2] = ['simulate',(sim1time-flag1time)]
stages[3] = ['invert',(dirty1time-sim1time)]
stages[4] = ['clean',(clean1time-dirty1time)]
stages[5] = ['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()
#
#=====================================================================