updated with new main

CC       =  gcc

CXX      =  g++

MPICC    =  mpicc

MPIC++   =  mpiCC

#FFTW_INCL=  -I/homes/gheller/Pleiadi/fftw-3.3.10/local/include/

#FFTW_LIB=  /homes/gheller/Pleiadi/fftw-3.3.10/local/lib/

FFTW_INCL=  -I/homes/gheller/Pleiadi/fftw-3.3.10/local-mpich/include/

FFTW_LIB=  /homes/gheller/Pleiadi/fftw-3.3.10/local-mpich/lib/

FITSIO_INCL= -I/iranet/arcsoft/cfitsio/cfitsio-4.1.0/include/

FITSIO_LIB= /iranet/arcsoft/cfitsio/cfitsio-4.1.0/lib/

NVCC = nvcc

NVFLAGS = -arch=sm_70 -Xcompiler -mno-float128 -std=c++11

NVLIB = -L/cineca/prod/opt/compilers/cuda/10.1/none/lib64/ -lcudart -lcuda

OMP= -fopenmp

CFLAGS +=  -I. $(FFTW_INCL) $(GSL_INCL) $(MPI_INCL) $(FITSIO_INCL)

OPTIMIZE = $(OMP) -O3 -mtune=native

# write the final image

OPT += -DWRITE_IMAGE

# perform w-stacking phase correction

<<<<<<< HEAD

OPT += -DPHASE_ON

# GPU support for FFT using cuFFTMP

OPT += -DCUDA_FFT

# Support CFITSIO

=======

#OPT += -DPHASE_ON

# Support CFITSIO !!! Remember to add the path to the CFITSIO library to LD_LIBRARY_PATH

>>>>>>> main

#OPT += -DFITSIO

# Perform true parallel images writing

#OPT += -DPARALLELIO

# Normalize uvw in case it is not done in the binMS

OPT += -DNORMALIZE_UVW

# Gridding kernel: GAUSS, GAUSS_HI_PRECISION, KAISERBESSEL

OPT += -DGAUSS_HI_PRECISION

ifeq (FITSIO,$(findstring FITSIO,$(OPT)))

	LIBS += -L$(FITSIO_LIB) -lcfitsio

import sys

import pyrap

from pyrap.tables import table, maketabdesc, makecoldesc

from pyrap.tables import tablecolumn

print(sys.argv[1])

#inputMS = '/local/work/gheller/mock/test1hr_t201806301100_SBL150.MS'

#inputMS = 'gauss_t201806301100_SBL180.MS'

#inputMS = 'gauss_t201806301100_SBH300.MS'

#inputMS = 'C_RADIO30MHz500_t201806301100_SBL180.MS'

#inputMS = 'C_sky802_t201806301100_SBL180.MS'

#inputMS = 'onehour_10asec_t201806301100_SBL180.MS'

#oinputMS = 'sixhours_t201806301100_SBL180.MS'

inputMS = sys.argv[1]

t=table(inputMS,readonly=False)

print("     ")

print("Columns in the MS")

print(t.colnames())

print("     ")

t.putcol('DATA', t.getcol('MODEL_DATA'))

t.flush()

#outpath = '/data/gridding/data/shortgauss_t201806301100_SBH255.binMS/'

print(sys.argv[1])

outpath = "/data/gridding/data/Lofarbig/"+sys.argv[1]+".binMS/"

outpath = "/homes/gheller/lofar2/gheller/binMS/"+sys.argv[1]+".binMS/"

os.mkdir(outpath)

# input MS

readtime0 = time.time()

#msfile = "/data/Lofar-data/results/L798046_SB244_uv.uncorr_130B27932t_146MHz.pre-cal.ms"

msfile = "/data/Lofar-Luca/results/"+sys.argv[1]+".ms/"

msfile = "/homes/gheller/lofar2/gheller/MS/"+sys.argv[1]+".ms/"

ms = pt.table(msfile, readonly=True, ack=False)

if rank == 0:

import subprocess

import sys

import time

# input data

datadir = "cgheller@login.m100.cineca.it://m100_scratch/userexternal/cgheller/Lofar/data_model"

cleandir = "cgheller@login.m100.cineca.it://m100_scratch/userexternal/cgheller/Lofar/data_wsclean"

dirtydir = "cgheller@login.m100.cineca.it://m100_scratch/userexternal/cgheller/Lofar/data_wsdirty"

filename_prefix = "RADIO30MHz"

file_ext = ".fits"

# measurement set

workdir = "."

ms_name = "hba-10hrs_t201806301100_SBH255.MS"    #"visibilities_aux.MS"

ms_tar = workdir+"/"+ms_name+"-template.tar"

# general parameters

resolution = "2.0asec"

size = "2000"

# command for creating visibilities

vis_command = "wsclean"

# parameters for the noise

noise_command = "../useful_scripts/noise.py"

noise_factor = "0.0001"

# command for imaging

img_command = "wsclean"

img_parallel = "4"

img_iter = "10000"

#img_options = "-weight briggs 0 -taper-gaussian 60 -apply-primary-beam -reorder -niter 10000 -mgain 0.8 -auto-threshold 5"

#img_options = "-apply-primary-beam -reorder -niter 10000 -mgain 0.8 -auto-threshold 5"

## create empty MS

## LBA 30 MHz

#/opt/losito/bin/synthms --name ms_name --tobs 1 --ra 1.0 --dec 1.57 --lofarversion 1 --station LBA -minfreq 30e6 --maxfreq 30e6

## HBA 150 MHz

#/opt/losito/bin/synthms --name dataset --tobs 1 --ra 1.0 --dec 1.57 --lofarversion 1 --station HBA --minfreq 150e6 --maxfreq 150e6

#

#tar cvf ms_tar ms_name

tinverse = 0.0

tcopy = 0.0

tnoise = 0.0

timaging = 0.0

tget = 0.0

tput = 0.0

for id in range(1,2):

    imgname = filename_prefix+str(id)

    #fitsfilename = datadir+'/'+imgname+file_ext

    fitsfilename = datadir+'/'+filename_prefix+file_ext

    commandline = "scp "+datadir+"/"+imgname+file_ext+" ."

    t0 = time.time()

    error = subprocess.call(commandline,shell=True)

    error = subprocess.run(["cp",imgname+file_ext,"test-model.fits"])

    tget = tget+(time.time()-t0)

    # untar the template MS

    error = subprocess.run(["tar","xvf",ms_tar])

    # create the visibilities

    t0 = time.time()

    error = subprocess.run([vis_command,"-j",img_parallel,"-scale",resolution,"--predict","-use-idg","-grid-with-beam","-name", "test", ms_name])

    tinverse = tinverse+(time.time()-t0)

    # Copy the MODEL_DATA column int the DATA column

    t0 = time.time()

    error = subprocess.run(["python","../useful_scripts/copycolumn.py",ms_name])

    tcopy = tcopy+(time.time()-t0)

    # Add noise

    t0 = time.time()

    error = subprocess.run([noise_command,"--msin",ms_name,"--factor",noise_factor])

    tnoise = tnoise+(time.time()-t0)

    # Imaging

    t0 = time.time()

    error = subprocess.run([img_command,"-j",img_parallel,"-apply-primary-beam","-reorder","-niter",img_iter,\

                            "-mgain", "0.8", "-auto-threshold", "5", "-size",size,size,"-scale",resolution,\

                            "-name",imgname,ms_name])

    timaging = timaging+(time.time()-t0)

    # remove at the end the measurement set

    error = subprocess.run(["rm","-fr",workdir+"/"+ms_name])

    error = subprocess.run(["rm","test-model.fits"])

    error = subprocess.run(["mkdir",imgname])

    #commandline = "rm *beam*.fits"

    #error = subprocess.Popen(commandline, shell=True,stdout=subprocess.PIPE,stderr=subprocess.PIPE)

    t0 = time.time()

    commandline = "scp "+imgname+"-dirty.fits"+" "+dirtydir

    error = subprocess.call(commandline, shell=True)

    commandline = "scp "+imgname+"-image-pb.fits"+" "+cleandir

    error = subprocess.call(commandline, shell=True)

    tput = tput+(time.time()-t0)

    commandline = "mv *.fits "+imgname

    error = subprocess.call(commandline, shell=True)

##    error = subprocess.run(["tar","cvf",imgname+".tar",imgname])

    error = subprocess.run(["rm","-fr",imgname])

    #commandline = "tar cvf "+imgname+".tar *.fits"

    #error = subprocess.Popen(commandline, shell=True,stdout=subprocess.PIPE,stderr=subprocess.PIPE)

    #commandline = "rm *.fits"

    #error = subprocess.Popen(commandline, shell=True,stdout=subprocess.PIPE,stderr=subprocess.PIPE)

print("Summary of timings: ")

print("Get data (sec)     = ",tget)

print("Put data (sec)     = ",tput)

print("Inverse time (sec) = ",tinverse)

print("Copy time (sec)    = ",tcopy)

print("Noise time (sec)   = ",tnoise)

print("Imaging time (sec) = ",timaging)