diff --git a/scripts/getStructure_from_phases.py b/scripts/getStructure_from_phases.py
deleted file mode 100755
index 73f2bb60602a6656d74866ef576ddb322919941b..0000000000000000000000000000000000000000
--- a/scripts/getStructure_from_phases.py
+++ /dev/null
@@ -1,253 +0,0 @@
-#!/usr/bin/python3
-
-##########################################################################
-#
-# Written my Maaijke Mevius -- see Radio Science publication DOI: 10.1002/2016RS006028
-#
-##########################################################################
-# History
-# 20/07/2016 Outlier rejection and argparse added -- Tim Shimwell
-# 06/06/2019 Changed for the use with h5parms -- Alexander Drabent
-
-import matplotlib as mpl
-mpl.use("Agg")
-from losoto.h5parm import h5parm
-from losoto.lib_operations import *
-from pylab import *
-import scipy
-from scipy.optimize import leastsq
-import argparse
-import logging
-import numpy
-
-def mad(arr):
- """ Median Absolute Deviation: a "Robust" version of standard deviation.
- Indices variabililty of the sample.
- https://en.wikipedia.org/wiki/Median_absolute_deviation
- """
- arr = np.ma.array(arr).compressed() # should be faster to not use masked arrays.
- med = np.median(arr)
- return np.median(np.abs(arr - med))
-
-def model(t, coeffs):
- coeffs[0] = abs(coeffs[0])
- return coeffs[0] + (t**coeffs[2])/coeffs[1]
-def residuals(coeffs, y, t):
- coeffs[0] = abs(coeffs[0])
- return y - model(t, coeffs)
-
-def main(h5parmfile,solset='sol000',soltab='phase000',nr_grid=1,doplot=True,outbasename='ionosphere',output_dir='./',dofilter=True):
- data = h5parm(h5parmfile)
- solset = data.getSolset(solset)
- soltab = solset.getSoltab(soltab)
- soltype = soltab.getType()
- if soltype != 'phase':
- logging.warning("Soltab is of type " + soltype + ", but should be phase. Skipping.")
- data.close()
- return(0)
- outfile = open(str(output_dir) + '/' + outbasename +'_structure.txt','w')
- vals = soltab.val
- station_names = soltab.ant[:]
- stations = solset.getAnt()
- phx = []
- phy = []
- allposx = []
- allposy = []
- nrtimes = len(soltab.time)
- nrfreqs = len(soltab.freq)
- timestep=int(nrtimes/nr_grid)
- for vals, coord, selection in soltab.getValuesIter(returnAxes=soltab.getAxesNames(), weight=False):
- try:
- vals = reorderAxes( vals, soltab.getAxesNames(), ['dir', 'pol', 'ant', 'time', 'freq'])
- vals = vals[0]
- except:
- vals = reorderAxes( vals, soltab.getAxesNames(), ['pol', 'ant', 'time', 'freq'])
-
- valsx = vals[0]
- valsy = vals[1]
- for i,station_name in enumerate(station_names):
- if "CS" not in station_name:
- continue
- valx = valsx[i,:,0].flatten()
- valy = valsy[i,:,0].flatten()
- phx.append(np.nan_to_num(valx))
- phy.append(np.nan_to_num(valy))
- if type(station_name) != str and type(station_name) != numpy.str_:
- allposx.append(stations[station_name.decode()])
- allposy.append(stations[station_name.decode()])
- else:
- allposx.append(stations[station_name])
- allposy.append(stations[station_name])
-
-
- phx = np.array(phx)
- phy = np.array(phy)
- allposx = np.array(allposx)
- allposy = np.array(allposy)
- D = allposx[:,np.newaxis,:]-allposx[np.newaxis,:]
- D2 = np.sqrt(np.sum(D**2,axis=-1))
- Dy = allposy[:,np.newaxis,:]-allposy[np.newaxis,:]
- D2y = np.sqrt(np.sum(Dy**2,axis=-1))
- S0s = []
- betas = []
- for itime in range(nr_grid+1):
- tm=[0,1000000000]
- if itime<nr_grid:
- tm[0]=int(itime*timestep)
- tm[1]=int(tm[0]+timestep)
- dphx=phx[:,np.newaxis,tm[0]:tm[1]]-phx[np.newaxis,:,tm[0]:tm[1]]
- dphy=phy[:,np.newaxis,tm[0]:tm[1]]-phy[np.newaxis,:,tm[0]:tm[1]]
- dphx=np.unwrap(np.remainder(dphx-dphx[:,:,0][:,:,np.newaxis]+np.pi,2*np.pi))
- dvarx=np.var(dphx,axis=-1)
- dphy=np.unwrap(np.remainder(dphy-dphy[:,:,0][:,:,np.newaxis]+np.pi,2*np.pi))
- dvary=np.var(dphy,axis=-1)
- myselect = [False]
- var_tolerance = 0
- flagselect = [[]]
- x = []
- while (len(np.unique(myselect)) == 1 and not np.unique(myselect)[0]) or len(flagselect[0]) == len(x):
- var_tolerance += 0.1
- print("INFO: Trying to use variance tolerance level of: ", var_tolerance)
- myselect=np.logical_and(D2>0,np.logical_and(np.any(np.logical_and(dvarx>1e-7,dvarx<var_tolerance),axis=0)[np.newaxis],np.any(np.logical_and(dvarx>1e-7,dvarx<var_tolerance),axis=1)[:,np.newaxis]))
- x = D2y[myselect]
- y = dvary[myselect]
- flagselect = np.where(y > 1.0)
- if dofilter:
- x=D2[myselect]
- y=dvarx[myselect]
- # Seems like real values dont occur above 1.0
- flagselect = np.where(y > 1.0)
- xplot = np.delete(x,flagselect)
- yplot = np.delete(y,flagselect)
- bins = np.logspace(np.log10(np.min(xplot)),np.log10(np.max(xplot)),10)
- binys = []
- binxs = []
- for i in range(0,len(bins)-1):
- binmin = bins[i]
- binmax = bins[i+1]
- inbin = np.intersect1d(np.where(xplot > binmin),np.where(xplot < binmax))
- binxs.append((binmin+binmax)/2.0)
- binys.append(np.median(yplot[inbin])+10*mad(yplot[inbin]))
- x0 = [0.1,1.0,1.0]
- xfit, flag = scipy.optimize.leastsq(residuals, x0, args=(binys,binxs))
- flagselect = np.where(y > model(x,xfit))
- print(xfit,'fitting')
- if doplot:
- x=D2[myselect]
- y=dvarx[myselect]
- subplot(3,1,1)
- scatter(x,y,color='b')
- if dofilter:
- y2 = y[flagselect]
- x2 = x[flagselect]
- scatter(x2,y2,color='r')
- scatter(x,model(x,xfit),color='gray',linestyle='-')#,'gray',linestyle='-',linewidth=3)
- x=D2[np.logical_and(D2>1e3,myselect)]
- y=dvarx[np.logical_and(D2>1e3,myselect)]
- if dofilter:
- flagselect = np.where(y > model(x,xfit))
- x = np.delete(x,flagselect)
- y = np.delete(y,flagselect)
- A=np.ones((2,x.shape[0]),dtype=float)
- A[1,:]=np.log10(x)
- par=np.dot(np.linalg.inv(np.dot(A,A.T)),np.dot(A,np.log10(y)))
- S0=10**(-1*np.array(par)[0]/np.array(par)[1])
- if doplot:
- plot(x,10**(par[0]+par[1]*np.log10(x)),'r-')
- xscale("log")
- yscale("log")
- xlim(30,4000)
- xlabel('baseline length [m]')
- ylabel('XX phase variance [rad$^2$]')
- title('XX diffractive scale: %3.1f km'%(float(S0)/1000.))
- myselect = [False]
- var_tolerance = 0
- flagselect = [[]]
- x = []
- while (len(np.unique(myselect)) == 1 and not np.unique(myselect)[0]) or len(flagselect[0]) == len(x):
- var_tolerance += 0.1
- print("INFO: Trying to use variance tolerance level of: ", var_tolerance)
- myselect=np.logical_and(D2>0,np.logical_and(np.any(np.logical_and(dvarx>1e-7,dvarx<var_tolerance),axis=0)[np.newaxis],np.any(np.logical_and(dvarx>1e-7,dvarx<var_tolerance),axis=1)[:,np.newaxis]))
- x = D2y[myselect]
- y = dvary[myselect]
- flagselect = np.where(y > 1.0)
- if dofilter:
- x=D2y[myselect]
- y=dvary[myselect]
- # seems like real values dont occur above 1.0
- flagselect = np.where(y > 1.0)
- xplot = np.delete(x,flagselect)
- yplot = np.delete(y,flagselect)
- bins = np.logspace(np.log10(np.min(xplot)),np.log10(np.max(xplot)),20)
- binys = []
- binxs = []
- for i in range(0,len(bins)-1):
- binmin = bins[i]
- binmax = bins[i+1]
- inbin = np.intersect1d(np.where(xplot > binmin),np.where(xplot < binmax))
- binxs.append((binmin+binmax)/2.0)
- binys.append(np.median(yplot[inbin])+10*mad(yplot[inbin]))
- x0 = [0.1,1.0,1.0]
- xfit, flag = scipy.optimize.leastsq(residuals, x0, args=(binys,binxs))
- flagselect = np.where(y > model(x,xfit))
- print(xfit,'fitting')
- if doplot:
- x=D2y[myselect]
- y=dvary[myselect]
- subplot(3,1,3)
- scatter(x,y,color='b')
- if dofilter:
- y2 = y[flagselect]
- x2 = x[flagselect]
- scatter(x2,y2,color='r')
- scatter(x,model(x,xfit),color='gray',linestyle='-')#,'gray',linestyle='-',linewidth=3)
- x=D2y[np.logical_and(D2y>1e3,myselect)]
- y=dvary[np.logical_and(D2y>1e3,myselect)]
- if dofilter:
- flagselect = np.where(y > model(x,xfit))
- x = np.delete(x,flagselect)
- y = np.delete(y,flagselect)
- A=np.ones((2,x.shape[0]),dtype=float)
- A[1,:]=np.log10(x)
- pary=np.dot(np.linalg.inv(np.dot(A,A.T)),np.dot(A,np.log10(y)))
- S0y=10**(-1*np.array(pary)[0]/np.array(pary)[1])
- if doplot:
- plot(x,10**(pary[0]+pary[1]*np.log10(x)),'r-')
- xscale("log")
- yscale("log")
- xlim(30,4000)
- xlabel('baseline length [m]')
- ylabel('YY phase variance [rad$^2$]')
- title('YY diffractive scale: %3.1f km'%(float(S0y)/1000.))
- savefig(output_dir + '/' + outbasename + '_structure.png')
- close()
- cla()
- S0s.append([S0,S0y])
- betas.append([par[1],pary[1]])
- outfile.write('S0s ****%s**** %s beta %s %s\n'%(S0,S0y,par[1],pary[1]))
- data.close()
- return(0)
-
-
-if __name__ == "__main__":
- parser = argparse.ArgumentParser()
- parser.add_argument('h5parmfile',help='Name of input h5parm file')
- parser.add_argument('--solset',help='Name of solution set to use')
- parser.add_argument('--soltab',help='Name of solution table to use')
- parser.add_argument('--outbasename',help='Outfile base name')
- parser.add_argument('--output_dir', help='Output directory', default='./')
- parser.add_argument('-p',dest='doplot',default=True,help='Do plotting')
- parser.add_argument('-f',dest='dofilter',default=True,help='Attempt to filter out odd solutions')
- parser.add_argument('-nr_grid',dest='nr_grid',default=1,help='Time step')
-
- args = parser.parse_args()
- h5parmfile = args.h5parmfile
- solset = args.solset
- soltab = args.soltab
- dofilter = args.dofilter
- doplot = args.doplot
- outbasename = args.outbasename
- output_dir = args.output_dir
- nr_grid = args.nr_grid
-
- main(h5parmfile,solset,soltab,nr_grid,doplot,outbasename,output_dir,dofilter)
diff --git a/steps/LoSoTo.Structure.cwl b/steps/LoSoTo.Structure.cwl
new file mode 100644
index 0000000000000000000000000000000000000000..e441532fb1821434c0faac24fa127b606f35a764
--- /dev/null
+++ b/steps/LoSoTo.Structure.cwl
@@ -0,0 +1,69 @@
+#!/usr/bin/env cwl-runner
+
+class: CommandLineTool
+cwlVersion: v1.2
+id: losoto_structure
+
+$namespaces:
+ lofar: https://git.astron.nl/eosc/ontologies/raw/master/schema/lofar.owl
+doc: |
+ Find the structure function from phase solutions of core stations.
+ WEIGHT: compliant
+
+
+requirements:
+ InlineJavascriptRequirement:
+ expressionLib:
+ - { $include: utils.js}
+ InitialWorkDirRequirement:
+ listing:
+ - entryname: 'parset.config'
+ entry: $(get_losoto_config('STRUCTURE').join('\n'))
+
+ - entryname: $(inputs.input_h5parm.basename)
+ entry: $(inputs.input_h5parm)
+ writable: true
+
+baseCommand: "losoto"
+
+arguments:
+ - '--verbose'
+ - $(inputs.input_h5parm.basename)
+ - parset.config
+
+hints:
+ DockerRequirement:
+ dockerPull: astronrd/linc
+
+inputs:
+ - id: input_h5parm
+ type: File
+ format: lofar:#H5Parm
+ - id: soltab
+ type: string
+ doc: "Solution table"
+ - id: doUnwrap
+ type: boolean?
+ - id: refAnt
+ type: string?
+ doc: Reference antenna, by default the first.
+ - id: plotName
+ type: string?
+ doc: Plot file name, by default no plot.
+ - id: ndiv
+ type: int?
+
+outputs:
+ - id: output_plot
+ type: File
+ outputBinding:
+ glob: '$(inputs.plotName)'
+ - id: logfile
+ type: File[]
+ outputBinding:
+ glob: '$(inputs.input_h5parm.basename)-losoto*.log'
+
+stdout: $(inputs.input_h5parm.basename)-losoto.log
+stderr: $(inputs.input_h5parm.basename)-losoto_err.log
+$schema:
+ - https://git.astron.nl/eosc/ontologies/raw/master/schema/lofar.owl
\ No newline at end of file
diff --git a/steps/structure_function.cwl b/steps/structure_function.cwl
deleted file mode 100644
index 389d271efd49a21054e6ab28b99d0bc0d6c68e9a..0000000000000000000000000000000000000000
--- a/steps/structure_function.cwl
+++ /dev/null
@@ -1,60 +0,0 @@
-class: CommandLineTool
-cwlVersion: v1.2
-$namespaces:
- lofar: 'https://git.astron.nl/eosc/ontologies/raw/master/schema/lofar.owl'
-id: structure_function
-baseCommand:
- - getStructure_from_phases.py
-inputs:
- - format: 'lofar:#H5Parm'
- id: h5parmFile
- type: File
- inputBinding:
- position: 0
- doc: List of h5parm files
- - default: 'target'
- id: solset
- type: string?
- inputBinding:
- position: 0
- prefix: '--solset'
- doc: Input solset name
- - default: 'phase000'
- id: soltab
- type: string?
- inputBinding:
- position: 0
- prefix: '--soltab'
- doc: Name of the soltab
- - default: 'POINTING'
- id: outbasename
- type: string?
- inputBinding:
- position: 0
- prefix: '--outbasename'
- doc: Namebase of the output files
-outputs:
- - id: structure_plot
- doc: Output plot
- type: File
- outputBinding:
- glob: $(inputs.outbasename+'_structure.png')
- - id: structure_txt
- doc: Output text
- type: File
- outputBinding:
- glob: $(inputs.outbasename+'_structure.txt')
- - id: log
- type: File[]
- outputBinding:
- glob: 'structure_function*.log'
-label: structure_function
-requirements:
- - class: InlineJavascriptRequirement
-hints:
- - class: DockerRequirement
- dockerPull: astronrd/linc
-stdout: structure_function.log
-stderr: structure_function_err.log
-$schema:
- - 'https://git.astron.nl/eosc/ontologies/raw/master/schema/lofar.owl'
diff --git a/workflows/HBA_target.cwl b/workflows/HBA_target.cwl
index 94ce324b5cf29120560ccd9e5af88b563a9647bc..ba9d64d7a6540992ef556b64ec020a84a6fbd4c5 100644
--- a/workflows/HBA_target.cwl
+++ b/workflows/HBA_target.cwl
@@ -142,9 +142,6 @@ inputs:
- id: clipAteam
type: boolean?
default: true
- - id: ncores
- type: int?
- default: 16
- id: lbfgs_historysize
type: int?
default: 10
@@ -268,8 +265,6 @@ steps:
source: clip_sources
- id: clipAteam
source: clipAteam
- - id: ncores
- source: ncores
- id: gsmcal_step
source: gsmcal_step
- id: lbfgs_historysize
diff --git a/workflows/linc_target.cwl b/workflows/linc_target.cwl
index cd38254ef0dc8ea023310e5e9e522133d57346a0..c1f5d5b8b7bbc5afb90a5626f39acae21f51b733 100644
--- a/workflows/linc_target.cwl
+++ b/workflows/linc_target.cwl
@@ -142,9 +142,6 @@ inputs:
- id: clipAteam
type: boolean?
default: true
- - id: ncores
- type: int?
- default: 16
- id: lbfgs_historysize
type: int?
default: 10
@@ -367,8 +364,6 @@ steps:
source: demix
- id: demix_sources
source: demix_sources
- - id: ncores
- source: ncores
- id: removed_bands
source: gsmcal/removed_bands
- id: min_unflagged_fraction
diff --git a/workflows/linc_target/finalize.cwl b/workflows/linc_target/finalize.cwl
index cb2b945aa44a25b818e8c1e8705104c5f5ce3715..c46be139e2a68b9528d53732c6796fbcf5018634 100644
--- a/workflows/linc_target/finalize.cwl
+++ b/workflows/linc_target/finalize.cwl
@@ -57,9 +57,6 @@ inputs:
default: '[CR]S*&'
- id: flags
type: File[]
- - id: ncores
- type: int?
- default: 16
- id: targetname
type: string?
default: 'pointing'
@@ -84,7 +81,7 @@ outputs:
linkMerge: merge_flattened
- id: inspection
outputSource:
- - structure_function/structure_plot
+ - structure_function/output_plot
- wsclean/image
- uvplot/output_image
type: File[]
@@ -216,8 +213,6 @@ steps:
- bad_antennas
- compare_stations_filter
valueFrom: $(self.join(''))
- - id: structure_file
- source: structure_function/structure_txt
- id: Ateam_separation_file
source: check_Ateam_separation.json
- id: solutions
@@ -361,22 +356,20 @@ steps:
label: h5parm_pointingname
- id: structure_function
in:
- - id: h5parmFile
+ - id: input_h5parm
source: write_solutions/outh5parm
- - id: solset
- default: target
- id: soltab
source:
- skymodel_source
- gsmcal_step
- valueFrom: $(self.join(''))
- - id: outbasename
+ valueFrom: $('target/'+self.join(''))
+ - id: plotName
source: targetname
+ valueFrom: $(self+'_structure.png')
out:
- - id: structure_plot
- - id: structure_txt
- - id: log
- run: ../../steps/structure_function.cwl
+ - id: output_plot
+ - id: logfile
+ run: ../../steps/LoSoTo.Structure.cwl
label: structure_function
- id: concat_logfiles_applygsm
in:
@@ -409,8 +402,7 @@ steps:
- id: file_list
linkMerge: merge_flattened
source:
- - structure_function/log
- - structure_function/structure_txt
+ - structure_function/logfile
- id: file_prefix
source: targetname
valueFrom: $(self+'_structure')