Chi sqrd gradient

montblanc/examples/MS_single_gpu_gradient_example.py

@@ -0,0 +1,134 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Copyright (c) 2016 Marzia Rivi
+#
+# This file is part of montblanc.
+#
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program; if not, see <http://www.gnu.org/licenses/>.
+
+# Compare analytical chi squared gradient computation with respect to 
+# Sersic parameters with the numerical gradent computation. 
+
+# Input: filename - Measurement Set file
+#        ns - number of sersic sources
+
+import sys
+import argparse
+import math
+import numpy as np
+import montblanc
+import montblanc.util as mbu
+
+from  montblanc.config import RimeSolverConfig as Options
+
+ARCS2RAD = np.pi/648000.
+
+#scalelength in arcsec
+minscale = 0.3*ARCS2RAD
+maxscale = 3.5*ARCS2RAD
+
+parser = argparse.ArgumentParser(description='TEST SERSIC GRADIENT')
+parser.add_argument('msfile', help='Input MS filename')
+parser.add_argument('-ns',dest='nssrc', type=int, default=1, help='Number of Sersic Galaxies')
+
+args = parser.parse_args(sys.argv[1:])
+
+# Get the RIME solver
+# Enable gradient computation: sersic_gradient=True
+slvr_cfg = montblanc.rime_solver_cfg(msfile=args.msfile,
+        sources=montblanc.sources(point=0, gaussian=0, sersic=args.nssrc),
+        init_weights=None, weight_vector=False,  
+        sersic_gradient=True, dtype='double', version='v4')
+
+with montblanc.rime_solver(slvr_cfg) as slvr:
+
+    nsrc, nssrc, ntime, nchan = slvr.dim_local_size('nsrc', 'nssrc', 'ntime', 'nchan')
+    np.random.seed(14352)
+
+    # Random source coordinates in the l,m (brightness image) domain
+    l= slvr.ft(np.random.random(nsrc)*100*ARCS2RAD)
+    m= slvr.ft(np.random.random(nsrc)*100*ARCS2RAD) 
+    lm=mbu.shape_list([l,m], shape=slvr.lm.shape, dtype=slvr.lm.dtype)
+    slvr.transfer_lm(lm)
+
+    # Brightness matrix for sources
+    stokes = np.empty(shape=slvr.stokes.shape, dtype=slvr.stokes.dtype)
+    I, Q, U, V = stokes[:,:,0], stokes[:,:,1], stokes[:,:,2], stokes[:,:,3]
+    I[:] = np.ones(shape=I.shape)*70.
+    Q[:] = np.zeros(shape=Q.shape)
+    U[:] = np.zeros(shape=U.shape)
+    V[:] = np.zeros(shape=V.shape)
+    slvr.transfer_stokes(stokes)
+    alpha = slvr.ft(np.ones(nssrc*ntime)*(-0.7)).reshape(nsrc,ntime)
+    slvr.transfer_alpha(alpha)
+
+    # If there are sersic sources, create their
+    # shape matrix and transfer it.
+    mod = slvr.ft(np.random.random(nssrc))*0.3
+    angle = slvr.ft(np.random.random(nssrc))*2*np.pi
+    e1 = mod*np.sin(angle)
+    e2 = mod*np.cos(angle)
+    R = slvr.ft(np.random.random(nssrc)*(maxscale-minscale))+minscale
+    sersic_shape = slvr.ft(np.array([e1,e2,R])).reshape((3,nssrc))
+    slvr.transfer_sersic_shape(sersic_shape)
+    print sersic_shape
+
+    #Set visibility noise variance (muJy)
+    time_acc = 60
+    efficiency = 0.9
+    channel_bandwidth_hz = 20e6
+    SEFD = 400e6
+    sigma = (SEFD*SEFD)/(2.*time_acc*channel_bandwidth_hz*efficiency*efficiency)
+    slvr.set_sigma_sqrd(sigma)
+
+    # Create observed data and upload it to the GPU
+    slvr.solve()
+    with slvr.context:
+        observed = slvr.retrieve_model_vis()
+
+    noiseR = np.random.normal(0,np.sqrt(sigma),size=observed.shape)
+    noiseI = np.random.normal(0,np.sqrt(sigma),size=observed.shape)
+    noise = noiseR +1j*noiseI
+    observed = observed+noise
+    slvr.transfer_observed_vis(observed)
+
+    slvr.solve() 
+
+    print slvr
+    print "analytical: ",slvr.X2_grad
+
+    # Execute the pipeline
+    l1 = slvr.X2
+    dq = np.empty([3,nssrc])
+    for i in xrange(nssrc):
+        e1_inc = sersic_shape.copy()
+        e1_inc[0,i] = e1_inc[0,i]+0.000001
+        slvr.transfer_sersic_shape(e1_inc)
+        slvr.solve()
+        l2 = slvr.X2
+        dq[0,i] = (l2-l1)*1e6
+        e2_inc = sersic_shape.copy()
+        e2_inc[1,i] = e2_inc[1,i]+ 0.000001
+        slvr.transfer_sersic_shape(e2_inc)
+        slvr.solve()
+        l2 = slvr.X2
+        dq[1,i] = (l2-l1)*1e6
+        R_inc = sersic_shape.copy()
+        R_inc[2,i] = R_inc[2,i]+0.00000001
+        slvr.transfer_sersic_shape(R_inc)
+        slvr.solve()
+        l2 = slvr.X2
+        dq[2,i] = (l2-l1)*1e8
+    print "numeric: ",dq

montblanc/examples/MS_single_node_gradient_example.py

@@ -0,0 +1,130 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Copyright (c) 2016 Marzia Rivi
+#
+# This file is part of montblanc.
+#
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program; if not, see <http://www.gnu.org/licenses/>.
+
+# Compare analytical chi squared gradient computation with respect to 
+# Sersic parameters with the numerical gradent computation. 
+
+# Input: filename - Measurement Set file
+#        ns - number of sersic sources
+
+import sys
+import argparse
+import math
+import numpy as np
+import montblanc
+import montblanc.util as mbu
+
+from  montblanc.config import RimeSolverConfig as Options
+
+ARCS2RAD = np.pi/648000.
+
+#scalelength in arcsec
+minscale = 0.3*ARCS2RAD
+maxscale = 3.5*ARCS2RAD
+
+parser = argparse.ArgumentParser(description='TEST SERSIC GRADIENT')
+parser.add_argument('msfile', help='Input MS filename')
+parser.add_argument('-ns',dest='nssrc', type=int, default=1, help='Number of Sersic Galaxies')
+
+args = parser.parse_args(sys.argv[1:])
+
+# Get the RIME solver
+# Enable gradient computation: sersic_gradient=True
+slvr_cfg = montblanc.rime_solver_cfg(msfile=args.msfile,
+        sources=montblanc.sources(point=0, gaussian=0, sersic=args.nssrc),
+        init_weights=None, weight_vector=False,  
+        sersic_gradient=True, dtype='double', version='v5')
+
+with montblanc.rime_solver(slvr_cfg) as slvr:
+
+    nsrc, nssrc, ntime, nchan = slvr.dim_local_size('nsrc', 'nssrc', 'ntime', 'nchan')
+    np.random.seed(14352)    
+
+    # Random source coordinates in the l,m (brightness image) domain
+    l= slvr.ft(np.random.random(nsrc)*100*ARCS2RAD)
+    m= slvr.ft(np.random.random(nsrc)*100*ARCS2RAD) 
+    slvr.lm[:] = mbu.shape_list([l,m], shape=slvr.lm.shape, dtype=slvr.lm.dtype)
+
+    # Brightness matrix for sources
+    I, Q, U, V = slvr.stokes[:,:,0], slvr.stokes[:,:,1], slvr.stokes[:,:,2], slvr.stokes[:,:,3]
+    I[:] = np.ones(shape=I.shape)*70.
+    Q[:] = np.zeros(shape=Q.shape)
+    U[:] = np.zeros(shape=U.shape)
+    V[:] = np.zeros(shape=V.shape)
+    slvr.alpha[:] = slvr.ft(np.ones(nssrc*ntime)*(-0.7)).reshape(nsrc,ntime)
+
+    # If there are sersic sources, create their
+    # shape matrix and transfer it.
+    mod = slvr.ft(np.random.random(nssrc))*0.3
+    angle = slvr.ft(np.random.random(nssrc))*2*np.pi
+    e1 = mod*np.sin(angle)
+    e2 = mod*np.cos(angle)
+    R = slvr.ft(np.random.random(nssrc)*(maxscale-minscale))+minscale
+    slvr.sersic_shape[:] = slvr.ft(np.array([e1,e2,R])).reshape((3,nssrc))
+    print slvr.sersic_shape
+
+    #Set visibility noise variance (muJy)
+    time_acc = 60
+    efficiency = 0.9
+    channel_bandwidth_hz = 20e6
+    SEFD = 400e6
+    sigma = (SEFD*SEFD)/(2.*time_acc*channel_bandwidth_hz*efficiency*efficiency)
+    slvr.set_sigma_sqrd(sigma)
+
+    # Create observed data and upload it to the GPU
+    slvr.solve()
+    observed = slvr.model_vis[:].copy() 
+
+    noiseR = np.random.normal(0,np.sqrt(sigma),size=observed.shape)
+    noiseI = np.random.normal(0,np.sqrt(sigma),size=observed.shape)
+    noise = noiseR +1j*noiseI
+    slvr.observed_vis[:] = (observed+noise).copy()
+
+    slvr.solve()
+
+    print slvr 
+    print "analytical: ",slvr.X2_grad
+
+    # Execute the pipeline
+    l1 = slvr.X2
+    sersic_shape = slvr.sersic_shape.copy()
+    dq = np.empty([3,nssrc])
+    for i in xrange(nssrc):
+        e1_inc = sersic_shape.copy()
+        e1_inc[0,i] = e1_inc[0,i]+0.000001
+        slvr.sersic_shape[:] = e1_inc
+        slvr.solve()
+        l2 = slvr.X2
+        dq[0,i] = (l2-l1)*1e6
+        e2_inc = sersic_shape.copy()
+        e2_inc[1,i] = e2_inc[1,i]+ 0.000001
+        slvr.sersic_shape[:] = e2_inc
+        slvr.solve()
+        l2 = slvr.X2
+        dq[1,i] = (l2-l1)*1e6
+        R_inc = sersic_shape.copy()
+        R_inc[2,i] = R_inc[2,i]+0.00000001
+        slvr.sersic_shape[:] = R_inc
+        slvr.solve()
+        l2 = slvr.X2
+        dq[2,i] = (l2-l1)*1e8
+    print "numeric: ",dq
+
+

montblanc/factory.py

@@ -137,12 +137,13 @@ def create_rime_solver_from_ms(slvr_class_type, slvr_cfg):
     autocor = slvr_cfg.get(Options.AUTO_CORRELATIONS)
 
     with MeasurementSetLoader(msfile, auto_correlations=autocor) as loader:
-        ntime, nbl, na, nbands, nchan = loader.get_dims()
+        ntime, nbl, na, nbands, nchan, npol = loader.get_dims()
         slvr_cfg[Options.NTIME] = ntime
         slvr_cfg[Options.NA] = na
         slvr_cfg[Options.NBL] = nbl
         slvr_cfg[Options.NBANDS] = nbands
         slvr_cfg[Options.NCHAN] = nchan
+        slvr_cfg[Options.NPOL] = npol
         slvr = slvr_class_type(slvr_cfg)
         loader.load(slvr, slvr_cfg)
         return slvr
@@ -167,7 +168,8 @@ def rime_solver(slvr_cfg):
     elif version == Options.VERSION_FIVE:
         from montblanc.impl.rime.v5.CompositeRimeSolver \
         import CompositeRimeSolver as RimeSolver
-        slvr_cfg[Options.CONTEXT] = __contexts
+        if slvr_cfg.get(Options.CONTEXT, None) is None: 
+            slvr_cfg[Options.CONTEXT] = __contexts
     else:
         raise ValueError('Invalid version %s' % version)
 

montblanc/impl/common/loaders/loaders.py

@@ -31,6 +31,7 @@
 ANTENNA_TABLE = 'ANTENNA'
 SPECTRAL_WINDOW = 'SPECTRAL_WINDOW'
 FIELD_TABLE = 'FIELD'
+POL_TABLE = 'POLARIZATION'
 
 class MeasurementSetLoader(BaseLoader):
     LOG_PREFIX = 'LOADER:'
@@ -43,6 +44,7 @@ def __init__(self, msfile, auto_correlations=False):
         self.antfile = '::'.join((self.msfile, ANTENNA_TABLE))
         self.freqfile = '::'.join((self.msfile, SPECTRAL_WINDOW))
         self.fieldfile = '::'.join((self.msfile, FIELD_TABLE))
+        self.polfile = '::'.join((self.msfile, POL_TABLE))
 
         montblanc.log.info("{lp} Opening Measurement Set {ms}.".format(
             lp=self.LOG_PREFIX, ms=self.msfile))
@@ -58,6 +60,7 @@ def __init__(self, msfile, auto_correlations=False):
         # (2) baseline (ANTENNA1, ANTENNA2)
         # (3) band (SPECTRAL_WINDOW_ID via DATA_DESC_ID)
         ordering_query = ' '.join(["SELECT FROM $ms",
+            #"WHERE ANTENNA1 != ANTENNA2",
             "WHERE FIELD_ID={fid}".format(fid=field_id),
             "" if auto_correlations else "AND ANTENNA1 != ANTENNA2",
             "ORDERBY TIME, ANTENNA1, ANTENNA2, "
@@ -70,6 +73,7 @@ def __init__(self, msfile, auto_correlations=False):
         self.tables['ant']  = at = pt.table(self.antfile, ack=False, readonly=True)
         self.tables['freq'] = ft = pt.table(self.freqfile, ack=False, readonly=True)
         self.tables['field'] = fit = pt.table(self.fieldfile, ack=False, readonly=True)
+        self.tables['pol'] = polt = pt.table(self.polfile, ack=False, readonly=True)
 
         self.nrows = ordered_ms.nrows()
         self.na = at.nrows()
@@ -80,6 +84,9 @@ def __init__(self, msfile, auto_correlations=False):
         bl_query = "SELECT FROM $ordered_ms ORDERBY UNIQUE ANTENNA1, ANTENNA2"
         self.nbl = pt.taql(bl_query).nrows()
 
+        # Number of polarizations (assuming to be the same for all spectral windows)
+        self.npol = polt.getcol('NUM_CORR')[0]
+
         # Number of channels per band
         chan_per_band = ft.getcol('NUM_CHAN')
 
@@ -105,6 +112,8 @@ def __init__(self, msfile, auto_correlations=False):
             aeq=autocor_eq, ebl=expected_nbl, astr=autocor_str, nbl=self.nbl))
         self.log("Found {nb} band(s), containing {cpb} channels.".format(
             nb=self.nbands, nc=chan_per_band[0], cpb=chan_per_band))
+        self.log("Found {npol} polarization(s) in the POLARIZATION table.".format(
+            npol=self.npol))
 
         # Sanity check computed rows vs actual rows
         computed_rows = self.ntime*self.nbl*self.nbands
@@ -116,13 +125,14 @@ def __init__(self, msfile, auto_correlations=False):
                         nbl=self.nbl, nb=self.nbands,
                         cr=computed_rows, msr=self.nrows))
 
+
     def get_dims(self):
         """
         Returns a tuple with the number of
         timesteps, baselines, antenna, channel bands and channels
         """
         # Determine the problem dimensions
-        return self.ntime, self.nbl, self.na, self.nbands, self.nchan
+        return self.ntime, self.nbl, self.na, self.nbands, self.nchan, self.npol
 
     def log(self, msg, *args, **kwargs):
         montblanc.log.info('{lp} {m}'.format(lp=self.LOG_PREFIX, m=msg),