Giter VIP home page Giter VIP logo

commander's Introduction

Commander is an Optimal Monte-carlo Markov chAiN Driven EstimatoR which implements fast and efficient end-to-end CMB posterior exploration through Gibbs sampling.

| Main features | Installation | Usage | License | Projects | Citation |

Main Features

The latest version - Commander3 - brings together critical features such as:

  • Modern Linear Solver
  • Map Making
  • Parallelism
  • Sky and instrumental modelling
  • CMB Component Separation

Commander3 is written using modern Fortran standards such as modules, sub modules, and object oriented derived types. The code is highly tuned and optimized to run on High Performance Computing (HPC) facilities, but it can also be run on your local machine.

The previous incarnation of Commander, - Commander2 - is now an internal part of Commander3, while the first version of the code, - Commander1 - is used mainly for debugging and/or legacy purposes. However, Commander1 has not been officially released; thus, it doesn't support CMake installation, as described in official documentation.

Installation

For the complete installation guide please refer to the official documentation, where you can find how to compile and run Commander on different platforms, including HPCs such as NERSC, UNINETT Sigma2, OWLs etc. Below you can find the short summary of how to compile it from source.

Prerequisites

To successfully run Commander, you need the following libraries:

  • MPI - required regardless of installation type;
  • OpenMP - required regardless of installation type;
  • BLAS - required regardless of installation type;
  • LAPACK - required regardless of installation type;
  • HDF5 - required only if compiled via Makefile;
  • FFTW - required only if compiled via Makefile;
  • Sharp2 - required only if compiled via Makefile;
  • Healpix - required only if compiled via Makefile;
  • CFitsio - required only if compiled via Makefile;

In addition you may want to install/update the following packages:

  • Automake version 1.16 or higher - required regardless of installation type;
  • Autoconf version 2.69 or higher - required regardless of installation type;
  • Libtool version 2.4.6 or higher - required regardless of installation type;

Compile using CMake

CMake is a tool which allows you to compile your code on various platform, via generation of build files (e.g. on Linux are Makefiles). It is configured to scan your system and identify present/missing libraries to download and install the missing ones. So, please install CMake before proceeding by this installation type.

Once CMake is installed, the Commander3 installation procedure is quite simple and consists of the following steps:

$ git clone https://github.com/hke/Commander.git
$ cd Commander
$ mkdir build
$ cd build

then to configure Commander3 compillation with, e.g. Intel Fortran compilers, use:

$ cmake -DCMAKE_C_COMPILER=icc -DCMAKE_CXX_COMPILER=icpc -DCMAKE_Fortran_COMPILER=ifort ..

wait while configuration is finished and then run:

$ cmake --build . -j n

where n is the amount of processors you wish to use to speed up the installation.

Because Commander3 is usually run on HPCs, where users do not have the sudo/root previleges, the default installation path is configured to be inside /path/to/cloned/repo/Commander/build/install/, where Commander3 binary can be found inside bin folder, under the name commander3.

Compile from source

After you cloned the latest version of the repo, you need to do the following:

  1. Determine the locations of your MPI compilers (mpif90, mpif77, mpicc, etc), and ensure that they function correctly.
  2. Determine the locations of the CFITSIO and LAPACK libraries, and how to link to these libraries.
  3. Look in the config/ directory and see if a configuration already exists which is similar to your machine. Copy the config (or the config.example) to a new file in the same directory. Call this new file config.<machine> where <machine> is the name of the system you will be building the software on.
  4. Edit the config file and specify all the options for your system.
  5. cd into the top level commander directory and set the COMMANDER environment variable to the string you used above as the <machine>. For example, if your config is named config.mylaptop, then you would set the environment variable in .bashrc as:
    • For Bourne-like shells
    $ export COMMANDER=mylaptop
    • For csh-like shells
    % setenv COMMANDER mylaptop
  6. To view some additional help about the available make targets:
$ make help
  1. To build and install the software:
$ make
$ make install

Usage

In short, to run Commander3, you need to:

$ export OMP_NUM_THREADS=1

then create the chain directory as specified in parameter file:

$ mkdir chains_dir

copy parameter file into it:

$ cp param_file.txt chains_dir/

and, finally, run Commander3 via following command:

$ mpirun -np num_proc ~/Commander/src/commander/commander param_file.txt 2>&1 | tee chains_dir/slurm.txt

Here, num_proc is the number of processors to use, slurm.txt is the file to store output logs.

As stated previously, Commander1 has not been officially released and is used primarily for debugging. If you wish to run it, however, you can compile it with Makefile using:

$ cd commander1
$ make

and then run the follwoing command:

$ mpirun -n num_proc ~/Commander/commander1/src/commander/commander param_file.txt 2>&1 | tee chains_dir/slurm.txt

Projects

Commander framework is part of the following projects:

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

Citation

If used for published results, please cite these papers:

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

Forkers

tskisner akarakci dncnwtts shazalvi mardom artembasyrov maksymbrl lukashergt simrannerval knut0815 fazlurahman moriguchir gracospa windshirely

commander's Issues

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.

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!

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

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.

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

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.

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.

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.

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)

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.

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.

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..

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.

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).

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.

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.

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.

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."

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.

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.

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.

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

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.

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).

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.