Commit f83d1e47 authored by Vigeesh Gangadharan's avatar Vigeesh Gangadharan
Browse files

Merge branch '50-apply-user-supplied-weights' into 'develop'

Resolve "Add an option to apply user-supplied weights"

See merge request !35
parents d79b2851 8f4c86bc
Pipeline #3638 passed with stage
in 2 minutes and 52 seconds
......@@ -91,13 +91,21 @@ export LD_LIBRARY_PATH=$CONDA_PREFIX/lib:$LD_LIBRARY_PATH
```sh
Usage: vfisv [OPTIONS]
SDC:GRIS Inversion Pipeline code
Options:
-p, --path TEXT Path to the fits files
-d, --id INTEGER Observation ID
-o, --out TEXT Output fits file
-n, --numproc INTEGER Number of processors to run on
-l, --line FLOAT Wavelength of the line (Angstrom)
-w, --width FLOAT Wavelength range +/- (Angstrom)
-l, --line FLOAT Wavelength of the line Angstrom
-w, --width FLOAT Wavelength range in Angstrom
--weights INTEGER... Weights (e.g --weights 1 2 2 2
--preview TEXT Filename to save the plot
--errors TEXT Filename to save the uncertainties
--diagnose TEXT Filename to write the inversion fits
--log TEXT Filename to write the log file
--version Show the version and exit.
-h, --help Show this message and exit.
```
......@@ -159,6 +167,21 @@ mpiexec -n 1 vfisv \
--log='test_log.txt'
```
Another example with user supplied weights,
```shell
mpiexec -n 1 vfisv \
--path='/dat/sdc/gris/20150919/level1_split/' \
--id=3 --line=15648.514 \
--numproc=20 \
--width=1.8 \
--weights 1 2 2 2\
--out='test_inversion.fits' \
--preview='test_preview.png' \
--errors='test_errors.fits' \
--diagnose='test_diagnose.fits' \
--log='test_log.txt'
```
### Python interface
The pipeline can be also be run within a python script. To do so, you need to call the python using `mpiexec` with one processor.
......
......@@ -166,6 +166,7 @@ def vfisv(
line: Optional[float] = 15662.017,
width: Optional[float] = 1.8,
numproc: Optional[int] = 2,
weights = None
):
"""
Calls VFISV fortran
......@@ -212,16 +213,29 @@ def vfisv(
vfisv_data = VFISVpackage(path, id, line, width)
# If weights not given, set to 0 for default vfisv weights
if not weights:
weights = [0, 0, 0, 0]
fweights = np.array(
weights,
dtype="int32",
)
# Spawn jobs to the grisinv:master+slaves
pycomm = MPI.COMM_SELF.Spawn(get_executable(), maxprocs=numproc)
# Merge intercommunicator
common_comm = pycomm.Merge(False)
for irank in range(numproc):
common_comm.Send([vfisv_data.line_props, MPI.FLOAT], dest=irank + 1, tag=5550)
common_comm.Send([vfisv_data.map_props, MPI.INT], dest=irank + 1, tag=5551)
common_comm.Send([vfisv_data.WAVE, MPI.FLOAT], dest=irank + 1, tag=5552)
common_comm.Send([fweights, MPI.INT], dest=irank + 1, tag=5553)
# May cause issue: large buffers may create problem
for irank in range(numproc):
......@@ -1507,7 +1521,7 @@ def create_preview(data, header, preview):
return None
def write_log_file(log, header, data, stokes, total_time, numproc):
def write_log_file(log, header, data, total_time, numproc,weights=None):
"""
Create the log file of the run
......@@ -1519,8 +1533,6 @@ def write_log_file(log, header, data, stokes, total_time, numproc):
Astropy header object
data : dict of numpy arrays
Array containing the inversion results and uncertainties
stokes : dict of numpy arrays
Array containing the stokes vector and fits
total_time : float
time spent on the inversion part
numproc : int
......@@ -1555,11 +1567,14 @@ def write_log_file(log, header, data, stokes, total_time, numproc):
f.write(f"Number of processors : {numproc}\n")
f.write(f"Total profiles inverted : {data['vlos'].size}\n\n")
if weights:
f.write(f"Weights used : {weights} \n\n")
# Some info about the uncertainties
f.write(f"Mean chi_sqr : {data['chsq'].value.mean():5.2f}\n")
f.write(f"Best chi_sqr : {data['chsq'].value.min():5.2f}\n")
f.write(f"Worst chi_sqr : {data['chsq'].value.max():5.2f}\n")
f.close()
......@@ -2144,13 +2159,14 @@ class VFISVpackage:
@click.option("-n", "--numproc", type=int, help="Number of processors to run on")
@click.option("-l", "--line", type=float, help="Wavelength of the line Angstrom")
@click.option("-w", "--width", type=float, help="Wavelength range in Angstrom")
@click.option("--weights", nargs=4, type=int, help="Weights (e.g --weights 1 2 2 2")
@click.option("--preview", type=str, help="Filename to save the plot")
@click.option("--errors", type=str, help="Filename to save the uncertainties")
@click.option("--diagnose", type=str, help="Filename to write the inversion fits")
@click.option("--log", type=str, help="Filename to write the log file")
@click.version_option(get_version())
def main(path, id, out, numproc, line, width, preview, errors, diagnose, log):
def main(path, id, out, numproc, line, width, weights, preview, errors, diagnose, log):
"""SDC:GRIS Inversion Pipeline code """
# Get the MPI Communicator
......@@ -2179,7 +2195,7 @@ def main(path, id, out, numproc, line, width, preview, errors, diagnose, log):
sys.exit()
# Run VFISV
data, stokes, header = vfisv(path, id, line, width, numproc)
data, stokes, header = vfisv(path, id, line, width,numproc,weights=weights)
# Get the MPI timer
total_time = MPI.Wtime() - start_time
......@@ -2202,7 +2218,7 @@ def main(path, id, out, numproc, line, width, preview, errors, diagnose, log):
# write log file
if log:
write_log_file(log, header, data, stokes, total_time, numproc)
write_log_file(log, header, data, total_time, numproc,weights=weights if weights else None)
if __name__ == "__main__":
......
SUBROUTINE INVERT (K,RES,SFIT)
SUBROUTINE INVERT (K,RES,SFIT,USERWEIGHTS)
!
! J M Borrero
! Jan 20, 2020
......@@ -20,6 +20,7 @@ SUBROUTINE INVERT (K,RES,SFIT)
REAL(SP), DIMENSION(5) :: ERR
REAL(SP), DIMENSION(PY,16) :: RES
REAL(SP), DIMENSION(4,PY,NUMW) :: SFIT
integer :: USERWEIGHTS(4)
!---------------------------------------------------------------------
REAL(SP), DIMENSION(4,NUMW) :: OBS, SCAT
REAL(SP), DIMENSION(4) :: WEIGHTS
......@@ -56,7 +57,13 @@ SUBROUTINE INVERT (K,RES,SFIT)
!
ICVAL=ICONT(K,J)
CALL GET_TOTPOL(OBS,TOTPOL)
CALL GET_WEIGHT(OBS,ICVAL,ICMEAN,RATUMQS,WEIGHTS)
IF (USERWEIGHTS(1) == 0) THEN
CALL GET_WEIGHT(OBS,ICVAL,ICMEAN,RATUMQS,WEIGHTS)
ELSE
WEIGHTS = USERWEIGHTS*NOISE
END IF
!-------------------------
! Initializie guess model
!-------------------------
......
......@@ -45,6 +45,7 @@ PROGRAM VFISV_SPEC
integer :: MYTASK1,NTASKS1
REAL(SP) :: LINE_PROP(6)
integer :: MAP_PROP(8)
integer :: USERWEIGHTS(4)
INTEGER :: MPI_IS_INITIALIZED
......@@ -116,9 +117,12 @@ PROGRAM VFISV_SPEC
CALL READ_SPEC_LINE(NTASKS)
! Initialize Wavelength grid
CALL MPI_RECV(WAVE,NUMW,MPI_REAL,0,5552,intracomm,STATUS,IERR)
! Get weights
CALL MPI_RECV(USERWEIGHTS,4,MPI_INTEGER,0,5553,intracomm,STATUS,IERR)
! Keep for allocation
CALL READ_SPEC_MAP(NTASKS)
!Receive the Stokes data
CALL MPI_RECV(SI,LX*LY*NUMW,MPI_REAL,0,7770,intracomm,STATUS,IERR)
CALL MPI_RECV(SQ,LX*LY*NUMW,MPI_REAL,0,7771,intracomm,STATUS,IERR)
......@@ -163,7 +167,7 @@ PROGRAM VFISV_SPEC
IF (NTASKS /= 1) THEN
! Slave recieves from Master
CALL MPI_RECV(ROLL1,1,MPI_INTEGER,1,1000,intracomm,STATUS,IERR)
CALL INVERT(ROLL1,RES,SFIT)
CALL INVERT(ROLL1,RES,SFIT,USERWEIGHTS)
! Slave sends to Master
CALL MPI_SEND(RES,PY*16,MPI_REAL,1,2000,intracomm,IERR)
CALL MPI_SEND(SFIT,PY*NUMW*4,MPI_REAL,1,3000,intracomm,IERR)
......@@ -196,7 +200,7 @@ PROGRAM VFISV_SPEC
IF ((NTASKS /= 1).AND.(NTASKS <= LX-ITERR*(MYTASK-2))) THEN
! Slave recieves from Master
CALL MPI_RECV(ROLL1,1,MPI_INTEGER,1,1000,intracomm,STATUS,IERR)
CALL INVERT(ROLL1,RES,SFIT)
CALL INVERT(ROLL1,RES,SFIT,USERWEIGHTS)
! Slave sends to Master
CALL MPI_SEND(RES,PY*16,MPI_REAL,1,2000,intracomm,IERR)
CALL MPI_SEND(SFIT,PY*NUMW*4,MPI_REAL,1,3000,intracomm,IERR)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment