COSMOMC (v2011 Oct): mpk.f90+lrg error

[Savvas Nesseris]

Hi,

I am trying to test the latest COSMOMC with the matter power spectrum datasets and the BAO data but I am getting the following error when I compile out of the box with EXTDATA = LRG (and of course adding the appropriate paths to GSL/LAPACK/CFITSIO etc):

$ ./cosmomc.exe ./params.ini
Random seeds: 20197, 23610 rand_inst: 0
reading: lrg_2009
reading: data/lrgDR7_kbands_kmax02kmin02_maxLv2.ALL_MAGCOVv3.txt
reading: data/lrgDR7_windows_kmax02kmin02_maxLv2.ALL_MAGCOVv3.txt
trying to read this many points 250
reading: data/lrgDR7_zerowinin_kmax02kmin02_maxLv2.ALL_MAGCOVv3.txt
reading: data/lrgDR7_zerowindat_kmax02kmin02_maxLv2.ALL_MAGCOVv3.txt
reading: data/lrgDR7_cov_kmax02kmin02_maxLv2.ALL_MAGCOVv3.txt
Failed a1, a2: -1.1482000000000001 8.2409999999999997
pos 8.2409999999999997 8.2397702500000012

This error seems to be coming from the subroutine LSS_LRG_mpklike_init() found in mpk.f90.

The relevant parts of the params.ini file are (everything not shown, please consider it to be the default setting):

Code: Select all

#filenames for matter power spectrum datasets, incl twodf
mpk_numdatasets = 1
mpk_dataset1 = data/lrgDR7kmax02kmin02newmaxLv2ALL_MAGCOVv3.dataset
#mpk_dataset1 = data/sdss_lrgDR4.dataset
#mpk_dataset1 = data/2df_2005.dataset

#filenames for BAO data sets
bao_numdatasets = 1
bao_dataset1 = data/wigglez_bao.dataset
#bao_dataset1 = data/sdss_bao.dataset

use_CMB = T
use_HST = F
use_mpk = T
use_BAO = T
use_clusters = F
use_BBN = F
use_Age_Tophat_Prior = F
use_SN = F
use_lya = F
use_min_zre = 0

If I compile without the LRG option in the Makefile and then include BAOs in the param.ini then I get an error

mpk: edit makefile to have "EXTDATA = LRG" to inlude LRGs

Just some more details, I am compiling on Windows 7 under Cygwin, with gcc 4.3.4 but I get the same error on a linux machine with gcc 4.1.2. When I use just the CMB/SnIa/HST data Cosmomc works as expected.

Since I am compiling and running Cosmomc out of the box with all the default params etc I am not sure what is going on. Am I missing something?

Many thanks in advance for any help!

Cheers,
Savvas

[Jason Dossett]

I cannot address the issue you are getting out of the LSS_LRG_mpklike_init() routine. However it appears that when compiling without the LRG option you are leaving use_mpk=T.

When compiling without the LRG option you need to either set use_mpk=F or change the mpk_dataset1 to a data set that is not lrgDR7. That should resolve the issue you were having when trying to run use the BAOs.
Hope that helps.


Best,
Jason

[Savvas Nesseris]

Hi Jason,

thanks for your reply, but compiling (out of the box) without the LRG option and setting "use_mpk=F" in the params.ini, I get the error (!)

Random seeds: 17375, 21835 rand_inst: 0
reading BAO data set: wigglez_2011
STOP Error reading mpk file

The relevant parts of the params.ini are:

Code: Select all

#filenames for matter power spectrum datasets, incl twodf
mpk_numdatasets = 1
#mpk_dataset1 = data/lrgDR7kmax02kmin02newmaxLv2ALL_MAGCOVv3.dataset
#mpk_dataset1 = data/sdss_lrgDR4.dataset
#mpk_dataset1 = data/2df_2005.dataset

#filenames for BAO data sets
bao_numdatasets = 1
bao_dataset1 = data/wigglez_bao.dataset
#bao_dataset2 = data/sdss_bao.dataset

use_CMB = T
use_HST = F
use_mpk = F
use_BAO = T
use_clusters = F
use_BBN = F
use_Age_Tophat_Prior = F
use_SN = F
use_lya = F
use_min_zre = 0

Does it work OK in your machine?

EDIT: I have just compiled the previous version (Aug 2011) without the LRGs and it works just fine with the following settings (and everything else default)

Code: Select all

use_mpk = F
use_BAO = T

Of course the BAOs are just SDSS, so probably something got broken when the WiggleZ were added???

EDIT 2: Again in the Aug 2011 version, using

Code: Select all

use_mpk = T
use_BAO = F

and

Code: Select all

mpk_dataset1 = data/lrgDR7kmax02kmin02newmaxLv2ALL_MAGCOVv3.dataset

gives the original error (Failed a1 a2....) while with

Code: Select all

mpk_dataset1 = data/sdss_lrgDR4.dataset

it works OK.

EDIT 3: Just confirmed all this with g95 as well, so I guess it's not a compiler issue...

[Savvas Nesseris]

I have (partly) solved the problem with the BAO data:

Random seeds: 17375, 21835 rand_inst: 0
reading BAO data set: wigglez_2011
STOP Error reading mpk file

The error message "Error reading mpk file" comes from the files "bao.f90" and "mpk.f90". In both cases the problem comes from the lines:

Code: Select all

if (iopb.ne.0) then
   stop 'Error reading mpk file'
endif

which checks if the data have been read in successfully:

Code: Select all

read (tmp_file_unit,*, iostat=iopb) bset%bao_invcov(i,:)

In my case (both the linux and windows/cygwin machines) iostat doesn't work properly (first case) or doesn't exist (second case). So, commenting out these lines that contain "if (iopb.ne.0)" makes the whole thing work like a charm. Of course, a better solution should be found (@Antony). Also, platform specific apps like iostat should be avoided as much as possible, that's programming 101, right?

As soon as I solve the other problem:

Failed a1, a2: −1.1482000000000001 8.2409999999999997
pos 8.2409999999999997 8.2397702500000012

I will again update the thread.


References:
1) iostat missing in Cygwin:
http://cygwin.com/ml/cygwin/2002-04/msg01661.html

[Jason Dossett]

Thanks for catching that I had forgot to change mpk to bao in that error line by the way.

[Ben Gold]

In my case (both the linux and windows/cygwin machines) iostat doesn't work properly (first case) or doesn't exist (second case). So, commenting out these lines that contain "if (iopb.ne.0)" makes the whole thing work like a charm. Of course, a better solution should be found (@Antony). Also, platform specific apps like iostat should be avoided as much as possible, that's programming 101, right?

References:
1) iostat missing in Cygwin:
http://cygwin.com/ml/cygwin/2002-04/msg01661.html

You are confusing the standalone executable "iostat" with the FORTRAN "iostat" specifier to a read statement. They have nothing to do with each other. The former is indeed missing in Cygwin, but this should have absolutely nothing to do with COSMOMC which is not using it.

The "iostat" specifier to a FORTRAN read statement is part of the language specification, should always exist, and should always return zero if there was no error. Any behavior otherwise is a compiler bug, and probably isn't something Antony can help you with.

It's curious that commenting that line out fixes something, but it really means you're ignoring a file error and could quite possibly be setting yourself up for a lot of trouble later.

[Savvas Nesseris]

You are right that I mixed up the two iostats, sorry :(

It's curious that commenting that line out fixes something, but it really means you're ignoring a file error and could quite possibly be setting yourself up for a lot of trouble later.

Actually as a test, I have just changed the part of the code to this:

Code: Select all

if (iopb.ne.0) then
	write(*,*)'Error type is',iopb
	stop
   !stop 'Error reading mpk file'
        endif

So, instead of the "hard stop", I can now see what the error is.