Giter Club home page Giter Club logo

commander's Introduction

Documentation Status astropy


cosmoglobe is a python package that interfaces the Cosmoglobe Sky Model with Commander outputs for the purpose of producing astrophysical sky maps.

Features

See the documentation for a more comprehensive guide.

Initialize the Cosmoglobe Sky Model (this downloads and caches a ~800 MB file with the sky model data)

import cosmoglobe

model = cosmoglobe.sky_model(nside=256)

Simulate the sky at 150 GHz in units of MJy/sr, smoothed to 40 arcmin with a gaussian beam:

import astropy.units as u

emission = model(150*u.GHz, fwhm=40*u.arcmin, output_unit="MJy/sr")

Integrate over a bandpass:

import numpy as np
import healpy as hp
import matplotlib.pyplot as plt

# Reading in WMAP K-band bandpass profile.
bandpass_frequencies, bandpass_weights = np.loadtxt(wmap_bandpass.txt, unpack=True)

# The units of the detector must be specified even if the bandpass is pre-normalized.
bandpass_weights *= u.Unit("K_RJ") # Specify K_RJ or K_CMB
bandpass_frequencies *= u.GHz

model.remove_dipole() # Remove the dipole from the CMB component
emission = model(
    freqs=bandpass_frequencies, 
    weights=bandpass_weights, 
    fwhm=0.8*u.deg, 
    output_unit="mK_RJ",
)

hp.mollview(emission[0], hist="norm") # Plotting the intensity
plt.show()

Installation

cosmoglobe can be installed via pip

pip install cosmoglobe

Funding

This work has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreements No 776282 (COMPET-4; BeyondPlanck), 772253 (ERC; bits2cosmology) and 819478 (ERC; Cosmoglobe).

License

GNU GPLv3

commander's People

Contributors

amaurea avatar artembasyrov avatar dncnwtts avatar eirikgje avatar gronx-bot avatar gustavbe avatar haraltho avatar hermda02 avatar hke avatar htihle avatar ingunnkw avatar krisjand avatar lillejohs avatar lmmocanu avatar maksymbrl avatar mariekf avatar metinsa avatar mohanagr avatar raurlien avatar trygvels avatar unfunfunt avatar

Stargazers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

Watchers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

commander's Issues

Ideas for lossless Huffman compression of TOD

It may be possible to perform a lossless Huffman compression of (single precision) TOD, which could save quite significant amounts of space. However, this idea only works if

  • the value range of the data does not include 0
  • the values have a fairly small dynamic range (i.e. maxval/minval is not very far from 1).

Both conditions should be fulfilled for Planck TOD, correct?
Under these conditions, we can determine the smallest possible distance between floating point values at minval (see, e.g., https://en.wikipedia.org/wiki/Unit_in_the_last_place) and call this d. All numbers in the TOD stream are then, by construction, an exact integer mutiple of d away from minval, and the maximum integer is (maxval-minval)/d. Representing the TOD by these integers is a lossless transformation, and the integers themselves can be differenced and Huffman-compressed.
I'd expect pretty significant space savings when using this method.

input map with (unphysically) high polarisation quadrupole (or dipole) not recovered in runs

Disclaimer: This issue came up during the component separation course AST9240 (so all this is based on that corresponding branch).
Keep in mind that I am completely new to component separation and commander, so I might not be as precise in formulating this as I'd like to be... Nonetheless, I thought it might be good to summarize what I've been doing for our future selves.

The issue

I've stumbled over the issue that for an input map with a high (unphysical) polarisation dipole or similarly for a high polarisation quadrupole we don't recover the input map.

In the next few postings I'll document what modifications I've done and how they affected this issue.

Current situation
(AST9240_2021_sim branch)

The below plots are from a run where the dipole map was added to all CMB maps (i.e. I, Q, and U) instead of just to the CMB I map.

high pol dipole input, not same output

ย 

high pol dipole input, ncorr output

ย 

high pol dipole input, gain output

mpiexec: Error: unknown option "-env"

Lukas recently cloned from AST9240_2022 and got this error when trying to run, having compiled with oneapi.

Loading intel_parallel_studio fixed this issue, but what caused the issue is not clear.

Use NEST ordering for pointing compression?

I'm not absolutely sure, but browsing the code I got the feeling that you are exclusively using Healpix pixel indices in RING scheme. For differencing and Huffman compression it might be quite advantageous to use NESTED scheme instead, since this might lead to a much more concentrated distribution of difference values.
Since the (uncompressed) maps themselves will most likely be stored in RING scheme (otherwise there would be problems with the SHTs), this incurs the overhead of doing quite a few nest2ring conversions whenever uncompressing a scan, but that might be an acceptable price to pay.
I'm sorry if you have already tried this and it didn't produce any appreciable gain. If so, please feel free to close this!

Redundant initialization of bandpass structure

In comm_data_mod

The fields data(n)%tod%label(j) and cpar%ds_tod_dets(i) are both supposed to hold the detector names. Get rid of one.

do j = 1, data(n)%ndet
    data(n)%bp(j)%p => comm_bp(cpar, n, i, detlabel=data(n)%tod%label(j))
end do

if (trim(cpar%ds_tod_type(i)) == 'none') then
    data(n)%bp(0)%p => comm_bp(cpar, n, i, detlabel=data(n)%label)
else
    data(n)%bp(0)%p => comm_bp(cpar, n, i, subdets=cpar%ds_tod_dets(i))
end if

More useful error messages when there is hash table read error

When using an incorrect parameter file, the error message is rather difficult to parse. For example, I'm getting the message

forrtl: severe (408): fort: (7): Attempt to use pointer KEY when it is not associated with a target

Image              PC                Routine            Line        Source
commander3         0000000001D3D690  Unknown               Unknown  Unknown
commander3         0000000000F4CFFA  hashtbl_mp_get_sl         105  hashtbl.f90
commander3         0000000000F505A4  hashtbl_mp_get_ha         177  hashtbl.f90
commander3         0000000000AC94AA  comm_param_mod_mp        2355  comm_param_mod.f90
commander3         0000000000AC8DD4  comm_param_mod_mp        2328  comm_param_mod.f90
commander3         00000000009EF351  comm_param_mod_mp         689  comm_param_mod.f90
commander3         00000000009C6127  comm_param_mod_mp         270  comm_param_mod.f90
commander3         000000000047BAEB  MAIN__                     84  commander.f90
commander3         000000000047B31E  Unknown               Unknown  Unknown
libc-2.17.so       00002B36298BA555  __libc_start_main     Unknown  Unknown
commander3         000000000047B229  Unknown               Unknown  Unknown

It's clear that the issue is that the hash table mod is trying to read a parameter that isn't defined in the parameter file or something along those lines, but it would be nice if there was a more legible output, for instance, the value of KEY that is trying to be read in.

Some rework of baseline subtractions

After discussing with @dncnwtts it seems as if baselines could be isolated and dealt with more 'modularly' than is currently the case. This would lead to cleaner code and as far as I can tell, not an appreciable reduction in speed. So let's see if we can do something about that!

Output files in the input directory

At line 3733 in comm_diffuse_comp_smod we seem to write an output file to the input directory, and not to the chains directory. This is an issue on clusters, where the input directory is sometimes mounted read-only for submitted jobs, and caused a crash on scinet. We should fix this by either removing this code (idk what it does) or fixing it to respect the output directory properly.

BAND_TOD_INIT_FROM_HDF parameter not being read

It appears that it is not possible to initialize sky parameter separately from instrument parameters.

According to the documentation, this should be the case if BAND_TOD_INIT_FROM_HDF&&& = none, but it seems like this is overwritten in the comm_signal_mod.f90 block:

    ! Initialize instrumental parameters
    call update_status(status, "init_chain_inst")
    if (trim(cpar%cs_init_inst_hdf) /= 'none' .or. present(init_from_output)) then
       do i = 1, numband
          if (cpar%ignore_gain_bp) then
             data(i)%gain          = 1.d0
             data(i)%bp(0)%p%delta = 0.d0
          else
             if (trim(cpar%cs_init_inst_hdf) == 'default' .or. present(init_from_output)) then
                call read_hdf(file, trim(adjustl(itext))//'/gain/'//trim(adjustl(data(i)%label)), &
                  & data(i)%gain)

                call get_size_hdf(file, trim(adjustl(itext))//'/bandpass/'//&
                     & trim(adjustl(data(i)%label)), ext)
                if (data(i)%ndet > ext(1)-1) then
                   write(*,*) 'Error -- init HDF file ', trim(chainfile), ' does not contain enough bandpass information'
                   stop
                end if
                allocate(bp_delta(0:ext(1)-1,ext(2)))
                call read_hdf(file, trim(adjustl(itext))//'/bandpass/'//trim(adjustl(data(i)%label)), &
                     & bp_delta)
                do j = 0, data(i)%ndet
                   data(i)%bp(j)%p%delta = bp_delta(j,:)
                end do
                deallocate(bp_delta)
             else
                call get_chainfile_and_samp(cpar%cs_init_inst_hdf, &
                     & chainfile, initsamp2)
                call int2string(initsamp2, itext2)
                call open_hdf_file(chainfile, file2, 'r')
                call read_hdf(file2, trim(adjustl(itext2))//'/gain/'//trim(adjustl(data(i)%label)), &
                     & data(i)%gain)

Potential speedup for Huffman decompression

As far as I understand, the current Huffman decompressor looks at the compressed stream one bit at a time, navigating through the tree depending on the value of that bit. I expect that this could be speeded up quite significantly with the following change:

  • whenever the decompresor starts work on a new value, it looks at the next 8 bits in the stream, not just the next single bit.
  • these 8 bits can be used in several (very small) lookup tables with information
    • whether these 8 bits contain a full key or not
    • if they contain a full key: how many bits the key had and what the corresponding value is
    • if they don't contain a full key: where in the tree to continue searching

After this first modified step, the search would continue exactly as before.
I'm not absolutely sure, but I expect this change could improve performance quite a lot, since most keys are short and since the total number of memory accesses is reduced drastically for the first 8 bits of each key.

Using the first 16 bits instad of 8 might also work, but then the size of the lookup tables may already be uncomfortably large (larger than L2 cache).

Tutorials out of date

The parameter file is pre-defaults in this directory - let's update it to something that works with a more recent version.

Do something sensible about nsides in components and maps

We have recently started smoothing the noise input files, and so we have code that smoothes input maps to a desired nside and fwhm. Should we consider applying these to bands and components as well? Ie. if I have a dust input map at nside 1024 but I want it at nside 16, can we just do that automatically inside commander? You could only do it for lower nsides probably, but it would save a lot of issues with the low-resolutions datasets we are implementing now

Chain statistics files output

In Commander1, there was a file that would output statistics like chi_hiliat, chi_lolat, beta_mean, amp_mean, etc. I believe this was in line ~883 of commander.f90 of Commander1, it would be good to have something like this implemented in Commander3.

Write outputs from missing parameters in a better way

Instead of just killing Commander when a single parameter is missing, loop through the whole parameter file and find the missing parameters, saved to some sort of list. Make these missing parameters written to terminal, only on the root core please.

Commander crush when attempting first run

Hi all,

I am running Commander3 on a cluster. I compiled current master branch with intel compilers. I attempt to run a tutorial parameter file but the job is terminated with the error attached below.

Abort(1090959) on node 0 (rank 0 in comm 0): Fatal error in PMPI_Init: Other MPI error, error stack:
MPIR_Init_thread(176)........: 
MPID_Init(1548)..............: 
MPIDI_OFI_mpi_init_hook(1554): 
(unknown)(): Other MPI error
[unset]: write_line error; fd=-1 buf=:cmd=abort exitcode=1090959
:
system msg for write_line failure : Bad file descriptor
Abort(1090959) on node 0 (rank 0 in comm 0): Fatal error in PMPI_Init: Other MPI error, error stack:
MPIR_Init_thread(176)........: 
MPID_Init(1548)..............: 
MPIDI_OFI_mpi_init_hook(1554): 
(unknown)(): Other MPI error
[unset]: write_line error; fd=-1 buf=:cmd=abort exitcode=1090959
:
system msg for write_line failure : Bad file descriptor
forrtl: severe (174): SIGSEGV, segmentation fault occurred
Image              PC                Routine            Line        Source             
libpthread-2.31.s  0000152B82E37420  Unknown               Unknown  Unknown
libmpi.so.12.0.0   0000152B81233BE1  MPIR_Err_return_c     Unknown  Unknown
libmpi.so.12.0.0   0000152B813D9ED0  MPI_Init              Unknown  Unknown
libmpifort.so.12.  0000152B829D748B  PMPI_INIT             Unknown  Unknown
commander3         000000000049276A  MAIN__                     77  commander.f90
commander3         00000000004923BD  Unknown               Unknown  Unknown
libc-2.31.so       0000152B806DD083  __libc_start_main     Unknown  Unknown
commander3         00000000004922DE  Unknown               Unknown  Unknown
--------------------------------------------------------------------------
Primary job  terminated normally, but 1 process returned
a non-zero exit code. Per user-direction, the job has been aborted.
--------------------------------------------------------------------------
--------------------------------------------------------------------------
mpirun detected that one or more processes exited with non-zero status, thus causing
the job to be terminated. The first process to do so was:

  Process name: [[28425,1],0]
  Exit code:    174
--------------------------------------------------------------------------

The parameter file is from BP10 branch. And the version of MPI I was using is mpirun (Open MPI) 4.1.5 I did not modify much but only change the path of output and data path. Is it an error related to MPI or running out of memory on my cluster? Looking forward to any help.

Point source catalog in chains file

Would be nice to have the chains file containing an array of point source locations in the parameter output; partially to make it easier for the sky model code.

Should be simple, so I'm assigning myself.

Memory corruption leading to segfaults

Currently, the devel branch cannot be run with TOD_NUM_BP_PROPOSALS_PER_ITER >0 - it has to be set to 0, or segfaults will result.

Statement from HKE:
"Yes, I don't think it's worth debugging this quite yet. Been there, done that.. ๐Ÿ™‚ What happens is actually that polres suddenly changes from 3.14e-2 to -7.e-6 or something like that during the tod projection operator -- which has nothing at all to do with the polres variable. But with a very small and negative polres variable, the derived polangle indices goes completely crazy, and the code crashes with a segfault. But that's just a symptom -- the actual problem is a memory corruption bug much further down the line. And the typical culprit for these types of memory corruption errors is passing wrong array sizes to external libraries, which do not perform array boundary checks. And by far the most typical offender is HDF. That's why I would like to have an HDF version with boundary checking turned on. That should hopefully let us see which array is the problematic one. So, my prediction is that this is eventually going to be traced down to some read_hdf call (or similar) with incorrect array size, dimension or type."

Issue with initializing the synchrotron component

Error when attempting to initialize the synchrotron component:

` FITSIO Error Status = 104 : could not open the named file
failed to find or open the following file: (ffopen)

/mn/stornext/u3/hke/xsan/commander3/BP9/data/0.4-Haslam

`

It appears that in the devel version, Commander is attempting to read the band named 0.4-Haslam as a file as opposed to using that band as a prior for the synchrotron monopole.

Merging burned-in chains

There doesn't seem to be an "easy" way (that I'm aware of) to take burned in instrument files from one experiment that was running independently and merge them together. It would be nice to do that since it takes a while for WMAP to burn in by itself.

Of course, if this functionality already exists, I retract the issue.

Map noise type not handled uniformly throughout codebase

There are some places where certain operations are applied if the noise type is rms and then does another operation otherwise, or just doesn't do it at all.

Two examples I found so far are in the subroutine compute_chisq (line 76):

          if (trim(data(i)%N%type) == "rms" .and. data(i)%N%nside_chisq_lowres < res%info%nside .and. present(chisq_fullsky) .and. present(lowres_eval)) then
             if (lowres_eval) then
                lowres = .true.
                info_lowres  => comm_mapinfo(data(i)%info%comm, data(i)%N%nside_chisq_lowres, 0, data(i)%info%nmaps, data(i)%info%nmaps==3)

                res_lowres => comm_map(info_lowres)
                res_lowres_temp => comm_map(info_lowres)

                call res%udgrade(res_lowres)
                res_lowres_temp%map = res_lowres%map ! Save temporarily

                call data(i)%N%invN_lowres(res_lowres) ! invN*res
                res_lowres%map = res_lowres_temp%map*res_lowres%map ! res*(invN*res)

                call res_lowres_temp%dealloc(); deallocate(res_lowres_temp)
             end if
          else
             lowres=.false.
             call data(i)%N%sqrtInvN(res)
             res%map = res%map**2 !(sqrtInvN*res)**2 = res*invN*res
          end if

It's not clear to me if things like lcut and the in-development rms_qucov should be skipped, or if the else construct is meant to only catch the pixel-pixel qucov noise format.

Another example is in comm_signal_mod.f90, line 447. I've recently added a default class that gives a warning, like this:

          select type (N)
          class is (comm_N_rms)
             if (trim(data(i)%tod%init_from_HDF) == 'default' .or. present(init_from_output)) then
                call data(i)%tod%initHDF(file, initsamp, data(i)%map, rms)
             else
                call get_chainfile_and_samp(data(i)%tod%init_from_HDF, &
                     & chainfile, initsamp2)
                call open_hdf_file(chainfile, file2, 'r')
                call data(i)%tod%initHDF(file2, initsamp2, data(i)%map, rms)
                call close_hdf_file(file2)
             end if
          class default
             write(*,*) 'Noise type is not covered'
          end select

just so that people know that the tod parameters aren't being read in.

In any case, this has had the effect of things that shouldn't have had anything to do with the noise format being affected when the noise format is changed. There's probably some other subtle things elsewhere.

NSIDE check outputs warning even though nsides agree

There is an issue with the if statement in comm_tod_mod that checks whether the nside for the band in the parameter file agrees with that in the hdf file. The output says

Nside= 0 found in parameter file does not match nside= 512 found in data files ,

even though the nside in the parameter file is 512.

Write parameter files to output chain file

@MetinSa wrote a commit to master that would output a parameters field to the output hdf5 file. Part of the issue is that this requires changing the comm_hdf_mod.f90.in file, which requires some fixes for the character output fields.

As it stands, hdf5 elements with character arrays as the output fields are not output correctly.

Unnecessary FFT work in comm_tod_noise_mod

In several places in comm_tod_noise_mod.f90 there is code analogous to

          dt(1:ntod)           = d_prime(:)
          dt(2*ntod:ntod+1:-1) = dt(1:ntod)
          call sfftw_execute_dft_r2c(plan_fwd, dt, dv)

The array dt represents an even function, i.e. the FFT transforms a function symmetric around the origin. This is in fact a discrete cosine transform of type 2, which can in principle be computed twice as quickly than the real-valued FFT used by the code.

If a significant fraction of overall run time is spent on these FFTs, it's probably worth switching to DCTs. They are supported by FFTW (see http://www.fftw.org/fftw3_doc/Real-even_002fodd-DFTs-_0028cosine_002fsine-transforms_0029.html).

Retired system referenced in documentation

The documentation references a retired system, Cori at NERSC. Can you consider retiring the build instructions for that system please?

If you want to continue supporting users of Commander at NERSC you may want to update your build procedure to target Perlmutter instead. Thanks.

Writing huffman compressed data from Commander simulations runs

Currently, what I've implemented in the write_huffman branch just writes uncompressed data, and I have written a python script that will compress the uncompressed data. It would be nice if this were done in Commander, but several things need to be done:

  • Create the huffman tree structure in commander
  • Encode the symbols as a binary string
  • Write the tree and the symbols into the output hdf5 file.

The first two just require fixing my buggy implementation, but the last one requires deleting a dataset in the existing hdf5 file, which as far as I can tell, is not actually implemented in the Fortran implementation of the hdf library.

Component 08 not found in default parameter files?

I've uploaded a new BP9 parameter file now. However, in this parameter setup the CMB relativistic quadrupole (component 08) still has explicit numbering in the included file. When trying to set this to ?? instead of 08, it fails..

Commander timing module

Here's a proposal for a Commander timing module. Any comments or suggestions? The idea is essentially to introduce a lot of timers in a new class, and one can start and stop each of them individually (like this: call timer%start(TOT_RUNTIME)), and ask for an ASCII report to be dumped to a given file. There are two types of timers, namely global and per-channel.

module comm_timing_mod
  use comm_utils
  implicit none
  
  ! Global parameters                                                                                                                                                                         
  integer(i4b), parameter :: NUM_GLOBAL    =  9
  integer(i4b), parameter :: TOT_RUNTIME   =  1
  integer(i4b), parameter :: TOT_AMPSAMP   =  2
  integer(i4b), parameter :: TOT_TODPROC   =  3
  integer(i4b), parameter :: TOT_SPECIND   =  4
  integer(i4b), parameter :: TOT_INIT      =  5
  integer(i4b), parameter :: TOT_FFT       =  6
  integer(i4b), parameter :: TOT_SHT       =  7
  integer(i4b), parameter :: TOT_OUTPUT    =  8
  integer(i4b), parameter :: TOT_INIT      =  9
  
  ! Channel specific parameters                                                                                                                                                               
  integer(i4b), parameter :: NUM_TOD       = 13
  integer(i4b), parameter :: TOD_INIT      =  1
  integer(i4b), parameter :: TOD_SL_PRE    =  2
  integer(i4b), parameter :: TOD_SL_INT    =  3
  integer(i4b), parameter :: TOD_PROJECT   =  4
  integer(i4b), parameter :: TOD_ORBITAL   =  5
  integer(i4b), parameter :: TOD_DECOMP    =  6
  integer(i4b), parameter :: TOD_ABSCAL    =  7
  integer(i4b), parameter :: TOD_RELCAL    =  8
  integer(i4b), parameter :: TOD_DELTAG    =  9
  integer(i4b), parameter :: TOD_NCORR     = 10
  integer(i4b), parameter :: TOD_XI_N      = 11
  integer(i4b), parameter :: TOD_MAPBIN    = 12
  integer(i4b), parameter :: TOD_MAPSOLVE  = 13

  private
  public comm_timing

  type comm_timing
     integer(i4b) :: numband, numsamp, comm, myid, n_tot
     real(dp),     allocatable, dimension(:)   :: t     ! Accumulated time for global timers (NUM_GLOBAL+NUM_TOD*numband)                                                                     
     real(dp),     allocatable, dimension(:)   :: t1    ! Start time for currently active timers for global timers (NUM_GLOBAL+NUM_TOD*numband)                                               
   contains
     procedure :: start        => comm_timer_start
     procedure :: start        => comm_timer_start
     procedure :: stop         => comm_timer_stop
     procedure :: incr_numsamp => comm_timer_incr_numsamp
     procedure :: dumpASCII    => comm_timer_dumpASCII
  end type comm_timing

  interface comm_timing
     procedure constructor
  end interface comm_timing

contains
  
  function constructor(numband, comm) result(res)
    !                                                                                                                                                                                         
    ! Constructor routine for timer object                                                                                                                                                    
    !                                                                                                                                                                                         
    ! Input variables:                                                                                                                                                                        
    !    numband  = number of frequency channels in current run                                                                                                                               
    !    comm     = MPI communicator                                                                                                                                                          
    !                                                                                                                                                                                         
    ! Output variable:                                                                                                                                                                        
    !    res      = pointer to comm_timing object                                                                                                                                             
    !                                                                                                                                                                                         
    implicit none
    integer(i4b),               intent(in) :: numband, comm
    type(comm_timing), pointer             :: res

    integer(i4b) :: ierr

    allocate(res)
    res%numband = numband
    res%n_tot   = NUM_GLOBAL + NUM_TOD * numband
    res%comm    = comm
    call mpi_comm_rank(comm, res%myid, ierr)

    allocate(res%t(n_tot), res%t1(n_tot))
    res%numsamp = 0
    res%t       = 0.d0
    res%t1      = 0.d0

  end function constructor

  subroutine comm_timer_start(self, timer_id, band)
    !                                                                                                                                                                                         
    ! Routine for starting timers                                                                                                                                                             
    !                                                                                                                                                                                         
    ! Input variables:                                                                                                                                                                        
    !    self     = comm_timing object                                                                                                                                                        
    !    timer_id = timer ID                                                                                                                                                                  
    !    bands    = band ID (optional)                                                                                                                                                        
    !                                                                                                                                                                                         
    implicit none
    type(comm_timing), intent(inout)          :: self
    integer(i4b),      intent(in)             :: timer_id
    integer(i4b),      intent(in),   optional :: band

    integer(i4b) :: i, timer
    real(dp)     :: t1

    call wall_time(t1)
    timer = timer_id; if (present(band)) timer = NUM_GLOBAL + NUM_TOD*(band-1) + timer_id
    self%t1(timer) = t1

  end subroutine comm_timer_start

  subroutine comm_timer_stop(self, timer_id, band)
    !                                                                                                                                                                                         
    ! Routine for stopping timers; difference between t2 and t1 will be                                                                                                                       
    ! accumulated into self%t_tot. Only currently active timers are accumulated                                                                                                               
    !                                                                                                                                                                                         
    ! Input variables:                                                                                                                                                                        
    !    self     = comm_timing object                                                                                                                                                        
    !    timer_id = timer ID                                                                                                                                                                  
    !    bands    = band ID (optional)                                                                                                                                                        
    !                                                                                                                                                                                         
    implicit none
    type(comm_timing), intent(inout)           :: self
    integer(i4b),      intent(in)              :: timer_id
    integer(i4b),      intent(in),   optional  :: band

    integer(i4b) :: i, timer
    real(dp)     :: t2

    timer = timer_id; if (present(band)) timer = NUM_GLOBAL + NUM_TOD*(band-1) + timer_id
    if (self%t1(timer) > 0) then
       call wall_time(t2)
       self%t(timer) = t2 - self%t1(timer)
    end if

  end subroutine comm_timer_stop

  subroutine comm_timer_incr_numsamp(self)
    !                                                                                                                                                                                         
    ! Routine for incrementing sample counter; used to output time per sample                                                                                                                 
    !                                                                                                                                                                                         
    ! Input variables:                                                                                                                                                                        
    !    self     = comm_timing object                                                                                                                                                        
    !                                                                                                                                                                                         
    implicit none
    type(comm_timing),               intent(inout) :: self

    self%numsamp = self%numsamp + 1

  end subroutine comm_timer_incr_numsamp

  subroutine comm_timer_dumpASCII(self, filename)
    !                                                                                                                                                                                         
    ! Routine for outputting timing information                                                                                                                                               
    !                                                                                                                                                                                         
    ! Input variables:                                                                                                                                                                        
    !    self     = comm_timing object                                                                                                                                                        
    !    filename = output filename                                                                                                                                                           
    !                                                                                                                                                                                         
    implicit none
    type(comm_timing),               intent(in) :: self
    character(len=*),                intent(in) :: filename

    integer(i4b) :: unit, ierr, band, b
    real(dp), dimension(NUM_TIMER) :: t

    if (self%numsamp == 0) return

    allocate(t(self%n_tot))
    call mpi_reduce(self%t, t, self%ntot, MPI_DOUBLE_PRECISION, &
         & MPI_SUM, 0, self%comm, ierr)

    if (self%myid == 0) then
       call getlun(unit)
       open(unit,file=trim(filename), recl=1024)
       write(unit,*) 'Timing summary'
       write(unit,*) ''
       write(unit,*) '   Global total timers:'
       write(unit,*) '       Number of samples         = ', self%numsamp
       write(unit,*) '       Total runtime             = ', t(TOT_RUNTIME)
       write(unit,*) '       Initialization            = ', t(TOT_INIT)
       write(unit,*) ''
       write(unit,*) '   Global per-sample timers:'
       write(unit,*) '       Chain output              = ', t(TOT_OUTPUT)  / self%numsamp
       write(unit,*) '       Amplitude sampling        = ', t(TOT_AMPSAMP) / self%numsamp
       write(unit,*) '       Spectral index sampling   = ', t(TOT_SPECIND) / self%numsamp
       write(unit,*) '       TOD processing            = ', t(TOT_TODPROC) / self%numsamp
       write(unit,*) '       Total FFT                 = ', t(TOT_FFT)     / self%numsamp
       write(unit,*) '       Total SHT                 = ', t(TOT_SHT)     / self%numsamp
       write(unit,*) ''
       write(unit,*) '   Channel-specific global timers:'

       do band = 1, self%numband
          b = NUM_GLOBAL + (band-1)*NUM_TOD
          if (all(t(b+1:b+NUM_TOD) == 0.d0)) cycle
          write(unit,*)
          write(unit,*) '     Channel ID                = ', band
          write(unit,*) '     TOD initialization        = ', t(b+TOD_INIT)
          write(unit,*) '     TOD sidelobe precomputation  = ', t(b+TOD_SL_PRE)   / self%numsamp
          write(unit,*) '     TOD sidelobe interpolation   = ', t(b+TOD_SL_INT)   / self%numsamp
          write(unit,*) '     TOD sky-to-tod projection    = ', t(b+TOD_PROJECT)  / self%numsamp
          write(unit,*) '     TOD orbital dipole           = ', t(b+TOD_ORBITAL)  / self%numsamp
          write(unit,*) '     TOD decompression            = ', t(b+TOD_DECOMP)   / self%numsamp
          write(unit,*) '     TOD absolute calibration     = ', t(b+TOD_ABSCAL)   / self%numsamp
          write(unit,*) '     TOD relative calibration     = ', t(b+TOD_RELCAL)   / self%numsamp
          write(unit,*) '     TOD delta G calibration      = ', t(b+TOD_DELTAG)   / self%numsamp
          write(unit,*) '     TOD correlated noise         = ', t(b+TOD_NCORR)    / self%numsamp
          write(unit,*) '     TOD corr noise PSD           = ', t(b+TOD_XI_N)     / self%numsamp
          write(unit,*) '     TOD binning                  = ', t(b+TOD_MAPBIN)   / self%numsamp
          write(unit,*) '     TOD map solution             = ', t(b+TOD_MAPSOLVE) / self%numsamp
       end do
       close(unit)
    end if

  end subroutine comm_timer_dumpASCII

end module comm_timing_mod

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    ๐Ÿ–– Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. ๐Ÿ“Š๐Ÿ“ˆ๐ŸŽ‰

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google โค๏ธ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.