diff -r -c -b -B -N cosmomc/params.ini cosmomc_sampler/params.ini
*** cosmomc/params.ini	2011-07-07 15:30:25.000000000 +0200
--- cosmomc_sampler/params.ini	2009-04-30 15:31:00.000000000 +0200
***************
*** 1,97 ****
  #Sample parameters for cosmomc in default parameterization
  
  #Root name for files produced
! file_root = chains/test
  
! #action = 0:  MCMC, action=1: postprocess .data file, action=2: find best fit point only
! action = 0
! 
! #Maximum number of chain steps
! samples = 200000
  
! #Feedback level ( 2=lots,1=chatty,0=none)
! feedback = 1
  
! #Temperature at which to Monte-Carlo
! temperature = 1
  
! #filenames for CMB datasets and SZ templates (added to C_l times parameter(13))
! #Note you may need to change lmax in cmbtypes.f90 to use small scales (e.g. lmax=2100)
! cmb_numdatasets = 1
! cmb_dataset1 = WMAP
! cmb_dataset_SZ1 = data/WMAP_SZ_VBand.dat
! cmb_dataset_SZ_scale1 = 1
! 
! cmb_dataset2 = data/acbar2007_v3_corr.newdat
! cmb_dataset_SZ2 = data/WMAP_SZ_VBand.dat
! cmb_dataset_SZ_scale2 = 0.28
! 
! cmb_dataset3 = data/CBIpol_2.0_final.newdat
! cmb_dataset4 = data/B03_NA_21July05.newdat
! 
! #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
! 
! #filename for supernovae (default SDSS compilation)
! SN_filename = data/supernovae.dataset
! 
! #if true, use HALOFIT for non-linear corrections (astro-ph/0207664).
! #note lyman-alpha (lya) code assumes linear spectrum
! nonlinear_pk = F
! 
! use_CMB = T
! use_HST = F
! use_mpk = F
! use_BAO = F
! use_clusters = F
! use_BBN = F
! use_Age_Tophat_Prior = T
! use_SN = F
! use_lya = F
! use_min_zre = 0
  
! #directory, e.g. window functions in directory windows under data_dir
! data_dir = data/
  
- #Force computation of sigma_8 even if use_mpk = F
- get_sigma8 = F
  
! #1: Simple Metropolis, 2: slice sampling, 3: slice sampling fast parameters, 4: directional gridding
! sampling_method = 1
  
! #if sampling_method =4, iterations per gridded direction
! directional_grid_steps = 20
  
! #use fast-slow parameter distinctions to speed up 
! #(note for basic models WMAP3 code is only ~3x as fast as CAMB)
! use_fast_slow = F
  
! #Can use covariance matrix for proposal density, otherwise use settings below
! #Covariance matrix can be produced using "getdist" program.
! propose_matrix = params_CMB.covmat
  
! #If propose_matrix is blank (first run), can try to use numerical Hessian to 
! #estimate a good propose matrix. As a byproduct you also get an approx best fit point
! estimate_propose_matrix = F
  
! #Tolerance on log likelihood to use when estimating best fit point
! delta_loglike = 2
  
! #Scale of proposal relative to covariance; 2.4 is recommended by astro-ph/0405462 for Gaussians
! #If propose_matrix is much broader than the new distribution, make proportionately smaller
! #Generally make smaller if your acceptance rate is too low
! propose_scale = 2.4
  
! #Increase to oversample fast parameters more, e.g. if space is odd shape
! oversample_fast = 1
  
! #if non-zero number of steps between sample info dumped to file file_root.data
  indep_sample = 0
  
! #number of samples to disgard at start; usually set to zero and remove later
  burn_in = 0
  
  #If zero set automatically
--- 1,60 ----
  #Sample parameters for cosmomc in default parameterization
  
  #Root name for files produced
! file_root = test
  
! samples = 100000
  
! delta_loglike = 2
  
! sampling_method = 1
  
! estimate_propose_matrix = F
  
! propose_scale = 2.4
  
  
! #Temperature at which to Monte-Carlo
! temperature = 1
  
! #action = 0, to MCMC, action=1, postprocess .data file
! action = 0
  
! #Feedback level ( 2=lots,1=chatty,0=less,-1=minimal)
! feedback = 1
  
! #Can re-start from the last line of previous run (.txt file)
! continue_from = 
  
! #Increase to oversample fast parameters,e.g. if space is odd shape
! oversample_fast = 1
  
! #Can use covariance matrix for proposal density, otherwise use settings below
! #Covariance matrix can be produced using "getdist" prorgram.
! propose_matrix = 
  
! #If action = 1
! redo_likelihoods = T
! redo_theory = F
! redo_cls = F
! redo_pk = F
! redo_skip = 0
! redo_outroot = 
! redo_thin = 1
! #If large difference in log likelihoods may need to offset to give sensible weights
! #for exp(difference in likelihoods)
! redo_likeoffset =  0
  
! #Number of distinct points to sample
! #Every accepted point is included
  
! #number of steps between independent samples
! #if non-zero all info is dumped to file file_root.data
  indep_sample = 0
  
! 
! #number of samples to disgard at start
! #May prefer to set to zero and remove later
  burn_in = 0
  
  #If zero set automatically
***************
*** 101,194 ****
  #MPI_Converge_Stop is a (variance of chain means)/(mean of variances) parameter that can be used to stop the chains
  #Set to a negative number not to use this feature. Does not guarantee good accuracy of confidence limits.
  MPI_Converge_Stop = 0.03
! 
! #Do initial period of slice sampling; may be good idea if 
! #cov matrix or widths are likely to be very poor estimates
! MPI_StartSliceSampling  = F
  
  #Can optionally also check for convergence of confidence limits (after MPI_Converge_Stop reached)
- #Can be good idea as small value of MPI_Converge_Stop does not (necessarily) imply good exploration of tails
  MPI_Check_Limit_Converge = F
  
  #if MPI_Check_Limit_Converge = T, give tail fraction to check (checks both tails):
  MPI_Limit_Converge = 0.025
! #permitted quantile chain variance in units of the standard deviation (small values v slow):
! MPI_Limit_Converge_Err = 0.2
  #which parameter's tails to check. If zero, check all parameters:
  MPI_Limit_Param = 0
  
- #if MPI_LearnPropose = T, the proposal density is continally updated from the covariance of samples so far (since burn in)
- MPI_LearnPropose = T
- #can set a value of converge at which to stop updating covariance (so that it becomes rigorously Markovian)
- #e.g. MPI_R_StopProposeUpdate = 0.4 will stop updating when (variance of chain means)/(mean of variances) < 0.4
- MPI_R_StopProposeUpdate = 0
- 
- #If have covmat, R to reach before updating proposal density (increase if covmat likely to be poor)
- #Only used if not varying new parameters that are fixed in covmat
- MPI_Max_R_ProposeUpdate = 2
- #As above, but used if varying new parameters that were fixed in covmat
- MPI_Max_R_ProposeUpdateNew = 30
  
  #if blank this is set from system clock
  rand_seed = 
  
! #If true, generate checkpoint files and terminated runs can be restarted using exactly the same command
! #and chains continued from where they stopped
! #With checkpoint=T note you must delete all chains/file_root.* files if you want new chains with an old file_root
! checkpoint = F
! 
! #file used by CAMB
! highL_unlensed_cl_template = ./camb/HighLExtrapTemplate_lenspotentialCls.dat
! 
! #CAMB parameters
! #If we are including tensors
! compute_tensors = F
! #Initial power spectrum amplitude point (Mpc^{-1})
! pivot_k = 0.05
! #If using tensors, enforce n_T = -A_T/(8A_s)
! inflation_consistency = F
! 
! #Set Y_He from BBN constraint; if false set to fixed value of 0.24 by default.
! bbn_consistency=T 
! 
! #Whether the CMB should be lensed (slows a lot unless also computing matter power)
! CMB_lensing = T
! high_accuracy_default = F
! #increase accuracy_level to run CAMB on higher accuracy
! #(default is about 0.3%, 0.1% if high_accuracy_default =T)
! #accuracy_level can be increased to check stability/higher accuracy, but inefficient
! accuracy_level = 1
! 
! #If action = 1
! redo_likelihoods = T
! redo_theory = F
! redo_cls = F
! redo_pk = F
! redo_skip = 0
! redo_outroot = 
! redo_thin = 1
! redo_add = F
! redo_from_text = F
! #If large difference in log likelihoods may need to offset to give sensible weights
! #for exp(difference in likelihoods)
! redo_likeoffset =  0
  
- #parameter start center, min, max, start width, st. dev. estimate
- param[omegabh2] = 0.0223 0.005 0.1 0.001 0.001
- param[omegadmh2] = 0.105 0.01 0.99 0.01 0.01
- param[theta] = 1.04 0.5 10 0.002 0.002
- param[tau] = 0.09 0.01 0.8 0.03 0.03
- 
- param[omegak] = 0 0 0 0 0
- param[fnu] = 0 0 0 0 0
- param[w] = -1 -1 -1 0 0
- 
- param[ns] = 0.95 0.5 1.5 0.02 0.01
- param[nt] = 0 0 0 0 0
- param[nrun] = 0 0 0 0 0
- 
- #log[10^10 A_s]
- param[logA] = 3 2.7 4 0.01 0.01
- param[r] = 0 0 0 0 0
- #SZ amplitude, as in WMAP analysis
- param[asz]= 1 0 2 0.4 0.4
--- 64,101 ----
  #MPI_Converge_Stop is a (variance of chain means)/(mean of variances) parameter that can be used to stop the chains
  #Set to a negative number not to use this feature. Does not guarantee good accuracy of confidence limits.
  MPI_Converge_Stop = 0.03
! #if MPI_LearnPropose = T, the proposal density is continally updated from the covariance of samples so far (since burn in)
! MPI_LearnPropose = T
  
  #Can optionally also check for convergence of confidence limits (after MPI_Converge_Stop  reached)
  MPI_Check_Limit_Converge = F
  
  #if MPI_Check_Limit_Converge = T, give tail fraction to check (checks both tails):
  MPI_Limit_Converge = 0.025
! #And the permitted percentil chain variance in units of the standard deviation (small values v slow):
! MPI_Limit_Converge_Err = 0.3
  #which parameter's tails to check. If zero, check all parameters:
  MPI_Limit_Param = 0
  
  
  #if blank this is set from system clock
  rand_seed = 
  
! #parameter start center, min, max, start width, propose width
! #e.g. for 2D Gaussian
! param1 = 0 -20 20 1 1
! param2 = 0 -20 20 1 1
! param3 = 0 0 0 0 0
! param4 = 0 0 0 0 0
! param5 = 0 0 0 0 0
! param6 = 0 0 0 0 0
! param7 = 0 0 0 0 0
! 
! #fast params
! param8 = 0 -20 20 1 1
! param9 = 0 0 0 0 0
! param10 = 0 0 0 0 0
! param11 = 0 0 0 0 0
! param12 = 0 0 0 0 0
! param13 = 0 0 0 0 0
  
diff -r -c -b -B -N cosmomc/source/calclike.f90 cosmomc_sampler/source/calclike.f90
*** cosmomc/source/calclike.f90	2010-05-10 17:28:00.000000000 +0200
--- cosmomc_sampler/source/calclike.f90	2011-08-03 11:14:37.771243747 +0200
***************
*** 1,25 ****
  module CalcLike
-  use CMB_Cls
-  use cmbtypes
-  use cmbdata
-  use mpk
   use Random
   use settings
   use ParamDef
-  use snovae
-  use WeakLen
-  use Lya
   implicit none
  
-  logical :: Use_Age_Tophat_Prior = .true.
-  logical :: Use_CMB = .true.
-  logical :: Use_BBN = .false.
-  logical :: Use_Clusters = .false.
-   
-  integer :: H0_min = 40, H0_max = 100
-  real :: Omk_min = -0.3, Omk_max = 0.3
-  real :: Use_min_zre = 0
-  integer :: Age_min = 10, Age_max = 20
   real :: Temperature  = 1
  
  contains
--- 1,9 ----
***************
*** 41,72 ****
    end function GenericLikelihoodFunction
  
    
-   function GetLogPrior(CMB, Info) !Get -Ln(Prior)
-     real GetLogPrior
-     real Age
-     Type (CMBParams) CMB
-     Type(ParamSetInfo) Info
- 
-     GetLogPrior = logZero
-  
-     if (.not. generic_mcmc) then
-      if (CMB%H0 < H0_min .or. CMB%H0 > H0_max) return
-      if (CMB%omk < Omk_min .or. CMB%omk > Omk_max .or. CMB%Omv < 0) return
-      if (CMB%zre < Use_min_zre) return
- 
-      Age = GetAge(CMB, Info)
-       !This sets up parameters in CAMB, so do not delete unless you are careful!
-  
-      if (Use_Age_Tophat_Prior .and. (Age < Age_min .or. Age > Age_max) .or. Age < 0) return
-     
-     end if
-     GetLogPrior = 0
-  
-   end function GetLogPrior
- 
    function GetLogLike(Params) !Get -Ln(Likelihood)
      type(ParamSet)  Params 
-     Type (CMBParams) CMB
      real GetLogLike
      real dum(1,1)
   
--- 25,32 ----
***************
*** 81,176 ****
  
      else
  
!      call ParamsToCMBParams(Params%P,CMB)
! 
!      GetLogLike = GetLogLikePost(CMB, Params%Info,dum,.false.)
      end if 
     end function GetLogLike
  
      
-   function GetLogLikePost(CMB, Info, inCls, HasCls) 
-     real GetLogLikePost
-     Type (CMBParams) CMB
-     Type(ParamSetInfo) Info
-     real, intent(in):: inCls(:,:)
-     logical, intent(in) :: HasCls
-     real acl(lmax,num_cls_tot)
-     integer error
- 
-     if (generic_mcmc) stop 'GetLogLikePost: not supported for generic'
- 
-     GetLogLikePost  = GetLogPrior(CMB, Info)
-     if ( GetLogLikePost >= logZero) then
-        GetLogLikePost = logZero
-        
-     else 
-        GetLogLikePost = GetLogLikePost + sum(CMB%nuisance(1:nuisance_params_used)**2)/2
-           !Unit Gaussian prior on all nuisance parameters
-        if (Use_BBN) GetLogLikePost = GetLogLikePost + (CMB%ombh2 - 0.022)**2/(2*0.002**2) 
-           !I'm using increased error bars here
-    
-        if (Use_CMB .or. Use_LSS) then
-           if (HasCls) then
-            acl = inCls
-            error =0
-           else
-            call GetCls(CMB, Info, acl, error)
-           end if
-          if (error /= 0) then
-           GetLogLikePost = logZero 
-          else
-           if (Use_CMB) GetLogLikePost = &
-             CMBLnLike(acl, CMB%norm(norm_freq_ix:norm_freq_ix+num_freq_params-1),CMB%nuisance) + GetLogLikePost
-           if (Use_mpk) GetLogLikePost = GetLogLikePost + LSSLnLike(CMB, Info%theory)
-           if (Use_WeakLen) GetLogLikePost = GetLogLikePost + WeakLenLnLike(CMB, Info%theory)     
-           if (Use_Lya) GetLogLikePost = GetLogLikePost +  LSS_Lyalike(CMB, Info%Theory)
-           if ( GetLogLikePost >= logZero) then
-             GetLogLikePost = logZero
-           end if
-          end if
-          if (Use_SN .and. GetLogLikePost /= logZero ) then
-             if (Info%Theory%SN_loglike /= 0) then
-              GetLogLikePost = GetLogLikePost + Info%Theory%SN_loglike
-             else
-              GetLogLikePost = GetLogLikePost + SN_LnLike(CMB)
-             end if
-                !Assume computed only every time hard parameters change
-          end if
-          if (Use_BAO .and. GetLogLikePost /= logZero ) then
-             if (Info%Theory%BAO_loglike /= 0) then
-              GetLogLikePost = GetLogLikePost + Info%Theory%BAO_loglike
-             else
-              GetLogLikePost = GetLogLikePost + BAO_LnLike(CMB)
-             end if
-                !Assume computed only every time hard parameters change
-          end if
-          if (Use_HST .and. GetLogLikePost /= logZero) then
-             if (Info%Theory%HST_loglike /= 0) then
-              GetLogLikePost = GetLogLikePost + Info%Theory%HST_loglike
-             else
-              GetLogLikePost = GetLogLikePost + HST_LnLike(CMB)
-             end if
-            !!Old: GetLogLikePost = GetLogLikePost + (CMB%H0 - 72)**2/(2*8**2)  !HST 
-          end if  
-      
-          else
-            if (Use_SN) GetLogLikePost = GetLogLikePost + SN_LnLike(CMB)
-            if (Use_BAO) GetLogLikePost = GetLogLikePost + BAO_LnLike(CMB)
-            if (Use_HST) GetLogLikePost = GetLogLikePost + HST_LnLike(CMB)           
-          end if
-          
-      
-       if (Use_Clusters .and. GetLogLikePost /= LogZero) then
-           GetLogLikePost = GetLogLikePost + &
-                  (Info%Theory%Sigma_8-0.9)**2/(2*0.05**2)
-           stop 'Write your cluster prior in calclike.f90 first!'
-       end if
- 
-      if (GetLogLikePost /= LogZero) GetLogLikePost = GetLogLikePost/Temperature
-    
- 
-     end if
  
-   end function GetLogLikePost
  
  end module CalcLike
--- 41,51 ----
  
      else
  
!        stop
      end if 
     end function GetLogLike
  
      
  
  
  end module CalcLike
diff -r -c -b -B -N cosmomc/source/conjgrad_wrapper.f90 cosmomc_sampler/source/conjgrad_wrapper.f90
*** cosmomc/source/conjgrad_wrapper.f90	2011-04-05 11:52:18.000000000 +0200
--- cosmomc_sampler/source/conjgrad_wrapper.f90	2011-08-03 11:14:58.530257636 +0200
***************
*** 4,9 ****
--- 4,10 ----
   use ParamDef
   use CalcLike
   use Random
+  use wrapper
   implicit none
   private 
   real ftol, Bestfit_loglike
diff -r -c -b -B -N cosmomc/source/driver.F90 cosmomc_sampler/source/driver.F90
*** cosmomc/source/driver.F90	2011-07-07 15:26:52.000000000 +0200
--- cosmomc_sampler/source/driver.F90	2011-08-03 11:37:58.804892011 +0200
***************
*** 7,25 ****
          use MonteCarlo
          use ParamDef
          use settings
-         use cmbdata
          use posthoc
-         use WeakLen
          use CalcLike
          use EstCovmatModule
          use ConjGradModule
-         use mpk
          use MatrixUtils
          use IO
          use ParamNames
! #ifdef WMAP_PARAMS
!         use WMAP_OPTIONS
! #endif
          implicit none 
          
          character(LEN=Ini_max_string_len) InputFile, LogFile
--- 7,20 ----
          use MonteCarlo
          use ParamDef
          use settings
          use posthoc
          use CalcLike
          use EstCovmatModule
          use ConjGradModule
          use MatrixUtils
          use IO
          use ParamNames
! 
          implicit none 
          
          character(LEN=Ini_max_string_len) InputFile, LogFile
***************
*** 100,120 ****
  
          propose_scale = Ini_Read_Real('propose_scale',2.4)
  
-         HighAccuracyDefault = Ini_Read_Logical('high_accuracy_default',.false.)
          AccuracyLevel = Ini_Read_Real('accuracy_level',1.)
          
-         if (Ini_HasKey('highL_unlensed_cl_template')) &
-           highL_unlensed_cl_template=  ReadIniFilename(DefIni,'highL_unlensed_cl_template')
- 
          checkpoint = Ini_Read_Logical('checkpoint',.false.)
          if (checkpoint) flush_write = .true.
          
! #ifdef WMAP_PARAMS
!         use_TT_beam_ptsrc = Ini_read_Logical('use_TT_beam_ptsrc')
!         use_TE = Ini_read_Logical('use_TE')
!         use_TT = Ini_Read_Logical('use_TT')
!         print *, 'WMAP beam TE TT', use_TT_beam_ptsrc, use_TE, use_TT
! #endif
  
  #ifdef MPI 
          
--- 95,106 ----
  
          propose_scale = Ini_Read_Real('propose_scale',2.4)
  
          AccuracyLevel = Ini_Read_Real('accuracy_level',1.)
          
          checkpoint = Ini_Read_Logical('checkpoint',.false.)
          if (checkpoint) flush_write = .true.
          
! 
  
  #ifdef MPI 
          
***************
*** 206,246 ****
            call InitRandom()
          end if
  
-         use_nonlinear = Ini_Read_Logical('nonlinear_pk',.false.)
-         pivot_k = Ini_Read_Real('pivot_k',0.05)
-         inflation_consistency = Ini_read_Logical('inflation_consistency',.false.)
-         bbn_consistency = Ini_Read_Logical('bbn_consistency',.true.)
- 
-         w_is_w = Ini_Read_Logical ('w_is_w',.true.)     
          oversample_fast = Ini_Read_Int('oversample_fast',1)
          use_fast_slow = Ini_read_Logical('use_fast_slow',.true.)
   
          if (Ini_Read_Logical('cmb_hyperparameters', .false.)) &
              call DoAbort( 'Hyperparameters not supported any more')
  
-         if (Ini_Read_String('use_2dF') /= '') stop 'use_2dF now replaced with use_mpk'
-         Use_Clusters = Ini_Read_Logical('use_clusters',.false.)
-         Use_mpk = Ini_Read_Logical('use_mpk',.false.) ! matter power spectrum, incl 2dF
-         Use_HST = Ini_Read_Logical('use_HST',.true.)
-         Use_BBN = Ini_Read_Logical('use_BBN',.false.)
-         Use_Age_Tophat_Prior= Ini_Read_Logical('use_Age_Tophat_Prior',.true.)
-         Use_SN = Ini_Read_Logical('use_SN',.false.)
-         if (Use_SN) SN_filename = ReadIniFileName(DefIni,'SN_filename')
-         Use_BAO = Ini_Read_Logical('use_BAO',.false.)
-         Use_CMB = Ini_Read_Logical('use_CMB',.true.)
-         Use_WeakLen = Ini_Read_Logical('use_WeakLen',.false.)
-         Use_min_zre = Ini_Read_Double('use_min_zre',0.d0) 
-         Use_Lya = Ini_Read_logical('use_lya',.false.)
-        
- 
-         if (Use_Lya .and. use_nonlinear) &
-              call DoAbort('Lya.f90 assumes LINEAR power spectrum input')
- 
-         !flag to force getting sigma8 even if not using LSS data 
-         use_LSS = Ini_Read_Logical('get_sigma8',.false.)
-         ! use_LSS = Use_2dF .or. Use_Clusters .or. Use_WeakLen
-         use_LSS = Use_LSS .or. Use_mpk .or. Use_Clusters .or. Use_WeakLen .or. Use_Lya
- 
          Temperature = Ini_Read_Real('temperature',1.)
          
          num_threads = Ini_Read_Int('num_threads',0)
--- 192,203 ----
***************
*** 258,297 ****
  
          Ini_fail_on_not_found = .true.
  
-         numsets = Ini_Read_Int('cmb_numdatasets')
          num_points = 0
          nuisance_params_used = 0
-         if (Use_CMB) then
-          do i= 1, numsets
-           filename(i) = ReadIniFileName(DefIni,numcat('cmb_dataset',i)) 
-           call ReadDataset(filename(i))
-           num_points = num_points + datasets(i)%num_points
-           keyname=numcat('cmb_dataset_SZ',i)
-           SZTemplate(i) = ''
-           if (Ini_HasKey(KeyName)) SZTemplate(i) = Ini_Read_String(keyname, .false.) 
-           if (SZTemplate(i)/='') then
-            SZScale(i) = Ini_read_Real(numcat('cmb_dataset_SZ_scale',i),1.0)
-            call ReadSZTemplate(datasets(i), SZTemplate(i),SZScale(i))
-           end if
-           nuisance_params_used = nuisance_params_used + datasets(i)%nuisance_parameters
-          end do
-          if (Feedback > 1) write (*,*) 'read datasets'
-         end if
  
          Ini_fail_on_not_found = .true.
          
-         nummpksets = Ini_Read_Int('mpk_numdatasets',0)
-         if (Use_mpk) then
-          do i= 1, nummpksets
-           mpk_filename(i) = ReadIniFileName(DefIni,numcat('mpk_dataset',i)) 
-           call ReadMpkDataset(mpk_filename(i))
-          end do
-          if (Feedback>1) write(*,*) 'read mpk datasets'
-         end if
- 
-         if(Use_BAO .and. use_dr7lrg) &
-          call MpiStop('DR7 LRG and BAO are based on the same dataset. You cannot use both.')
- 
          numtoget = Ini_Read_Int('samples')
  
          call Initialize(DefIni,Params)
--- 215,225 ----
diff -r -c -b -B -N cosmomc/source/EstCovmat.f90 cosmomc_sampler/source/EstCovmat.f90
*** cosmomc/source/EstCovmat.f90	2011-04-05 11:53:35.000000000 +0200
--- cosmomc_sampler/source/EstCovmat.f90	2011-08-03 11:22:03.070110984 +0200
***************
*** 12,17 ****
--- 12,18 ----
   use CalcLike
   use settings
   use Matrixutils
+  use wrapper
   implicit none
  
   real, dimension(:,:,:,:), allocatable :: Lgrid
diff -r -c -b -B -N cosmomc/source/inifile.f90 cosmomc_sampler/source/inifile.f90
*** cosmomc/source/inifile.f90	1970-01-01 01:00:00.000000000 +0100
--- cosmomc_sampler/source/inifile.f90	2011-08-03 11:22:18.669373253 +0200
***************
*** 0 ****
--- 1,671 ----
+ !Module to read in name/value pairs from a file, with each line of the form line 'name = value'
+ !Should correctly interpret FITS headers
+ !Antony Lewis (http://cosmologist.info/). Released to the public domain.
+ !This version Oct 2009.
+ 
+ module IniFile
+  implicit none
+  public
+ 
+   integer, parameter :: Ini_max_name_len = 128
+ 
+   integer, parameter :: Ini_max_string_len = 1024
+   logical :: Ini_fail_on_not_found = .false.
+ 
+   logical :: Ini_Echo_Read = .false.
+ 
+   type TNameValue
+    !no known way to make character string pointers..
+     character(Ini_max_name_len)  :: Name
+     character(Ini_max_string_len):: Value
+   end type TNameValue
+ 
+   type TNameValue_pointer
+      Type(TNameValue), pointer :: P
+   end type TNameValue_pointer
+ 
+   Type TNameValueList
+     integer Count
+     integer Delta
+     integer Capacity
+     type(TNameValue_pointer), dimension(:), pointer :: Items
+   end Type TNameValueList
+ 
+   Type TIniFile
+    logical SlashComments
+    Type (TNameValueList) :: L, ReadValues
+   end Type TIniFile
+  
+   Type(TIniFile) :: DefIni
+ 
+ contains
+ 
+    subroutine TNameValueList_Init(L)
+     Type (TNameValueList) :: L
+     
+      L%Count = 0
+      L%Capacity = 0
+      L%Delta = 128
+      nullify(L%Items)
+ 
+    end subroutine TNameValueList_Init
+ 
+    subroutine TNameValueList_Clear(L)
+     Type (TNameValueList) :: L
+     integer i, status
+      
+     do i=L%count,1,-1
+      deallocate (L%Items(i)%P, stat = status)
+     end do
+     deallocate (L%Items, stat = status)
+     call TNameValueList_Init(L)
+ 
+    end subroutine TNameValueList_Clear
+ 
+    subroutine TNameValueList_ValueOf(L, AName, AValue)
+      Type (TNameValueList), intent(in) :: L
+      character(LEN=*), intent(in) :: AName
+      CHARACTER(LEN=*), intent(out) :: AValue
+      integer i
+ 
+      do i=1, L%Count
+        if (L%Items(i)%P%Name == AName) then
+           AValue = L%Items(i)%P%Value 
+           return
+        end if
+      end do
+      AValue = ''
+ 
+    end subroutine TNameValueList_ValueOf
+ 
+    function TNameValueList_HasKey(L, AName) result (AValue)
+      Type (TNameValueList), intent(in) :: L
+      character(LEN=*), intent(in) :: AName
+      logical :: AValue
+      integer i
+ 
+      do i=1, L%Count
+        if (L%Items(i)%P%Name == AName) then
+           AValue = .true.
+           return
+        end if
+      end do
+      AValue = .false.
+      
+    end function TNameValueList_HasKey
+     
+    subroutine TNameValueList_Add(L, AName, AValue)
+     Type (TNameValueList) :: L
+     character(LEN=*), intent(in) :: AName, AValue
+ 
+     if (L%Count == L%Capacity) call TNameValueList_SetCapacity(L, L%Capacity + L%Delta)
+     L%Count = L%Count + 1
+     allocate(L%Items(L%Count)%P)
+     L%Items(L%Count)%P%Name = AName
+     L%Items(L%Count)%P%Value = AValue
+ 
+    end subroutine TNameValueList_Add
+ 
+    subroutine TNameValueList_SetCapacity(L, C)
+     Type (TNameValueList) :: L
+     integer C
+     type(TNameValue_pointer), dimension(:), pointer :: TmpItems
+     
+     if (L%Count > 0) then
+       if (C < L%Count) stop 'TNameValueList_SetCapacity: smaller than Count'
+       allocate(TmpItems(L%Count))
+       TmpItems = L%Items(1:L%Count)
+       deallocate(L%Items)
+       allocate(L%Items(C))
+       L%Items(1:L%Count) = TmpItems
+       deallocate(TmpItems)
+     else
+      allocate(L%Items(C))
+     end if  
+     L%Capacity = C
+   
+    end subroutine TNameValueList_SetCapacity
+ 
+    subroutine TNameValueList_Delete(L, i)
+     Type (TNameValueList) :: L
+     integer, intent(in) :: i
+      
+      deallocate(L%Items(i)%P)
+      if (L%Count > 1) L%Items(i:L%Count-1) = L%Items(i+1:L%Count)
+      L%Count = L%Count -1
+      
+    end subroutine TNameValueList_Delete
+ 
+   subroutine Ini_NameValue_Add(Ini,AInLine)
+     Type(TIniFile) :: Ini
+     character (LEN=*), intent(IN) :: AInLine
+     integer EqPos, slashpos, lastpos
+     character (LEN=len(AInLine)) :: AName, S, InLine
+ 
+       InLine=trim(adjustl(AInLine))
+       EqPos = scan(InLine,'=')
+       if (EqPos/=0 .and. InLine(1:1)/='#' .and. InLine(1:7) /= 'COMMENT' ) then
+    
+          AName = trim(InLine(1:EqPos-1))
+          
+          S = adjustl(InLine(EqPos+1:)) 
+            if (Ini%SlashComments) then
+            slashpos=scan(S,'/')
+            if (slashpos /= 0) then
+               S  = S(1:slashpos-1)
+            end if
+          end if
+          lastpos=len_trim(S)
+          if (lastpos>1) then
+           if (S(1:1)=='''' .and. S(lastpos:lastpos)=='''') then
+            S = S(2:lastpos-1)
+           end if
+          end if
+          call TNameValueList_Add(Ini%L, AName, S)
+ 
+       end if
+ 
+   end subroutine Ini_NameValue_Add
+ 
+   subroutine Ini_Open(filename, unit_id,  error, slash_comments)
+      character (LEN=*), intent(IN) :: filename
+      integer, intent(IN) :: unit_id
+      logical, optional, intent(OUT) :: error
+      logical, optional, intent(IN) :: slash_comments
+      logical aerror
+ 
+      call TNameValueList_Init(DefIni%L)
+      call TNameValueList_Init(DefIni%ReadValues)
+           
+      if (present(slash_comments)) then
+       call Ini_Open_File(DefIni,filename,unit_id,aerror,slash_comments)
+      else
+       call Ini_Open_File(DefIni,filename,unit_id,aerror)
+      end if
+ 
+      if (present(error)) then
+        error = aerror
+      else
+       if (aerror) then
+         write (*,*) 'Ini_Open: Error opening file ' // trim(filename)
+         stop
+       end if
+      end if
+ 
+   end subroutine Ini_Open
+ 
+ 
+ 
+   subroutine Ini_Open_File(Ini, filename, unit_id,  error, slash_comments)
+      Type(TIniFile) :: Ini
+ 
+      character (LEN=*), intent(IN) :: filename
+      integer, intent(IN) :: unit_id
+      logical, intent(OUT) :: error
+      logical, optional, intent(IN) :: slash_comments
+      character (LEN=120) :: InLine
+     
+ 
+     call TNameValueList_Init(Ini%L)
+     call TNameValueList_Init(Ini%ReadValues)
+ 
+     if (present(slash_comments)) then
+      Ini%SlashComments = slash_comments
+     else
+      Ini%SlashComments = .false.
+     end if
+          
+     open(unit=unit_id,file=filename,form='formatted',status='old', err=500)
+    
+     do 
+       read (unit_id,'(a)',end=400) InLine
+       if (InLine == 'END') exit;
+       if (InLine /= '') call Ini_NameValue_Add(Ini,InLine) 
+     end do
+ 
+ 400 close(unit_id)
+     error=.false.
+     return
+ 
+ 500 error=.true.
+ 
+   end subroutine Ini_Open_File
+ 
+   subroutine Ini_Open_Fromlines(Ini, Lines, NumLines, slash_comments)
+     Type(TIniFile) :: Ini
+ 
+     integer, intent(IN) :: NumLines
+     character (LEN=*), dimension(NumLines), intent(IN) :: Lines
+     logical, intent(IN) :: slash_comments
+     integer i
+ 
+     call TNameValueList_Init(Ini%L)
+     call TNameValueList_Init(Ini%ReadValues)
+ 
+     Ini%SlashComments = slash_comments
+ 
+     do i=1,NumLines
+        call Ini_NameValue_Add(Ini,Lines(i))
+     end do
+ 
+   end  subroutine Ini_Open_Fromlines
+ 
+   subroutine Ini_Close
+ 
+     call Ini_close_File(DefIni)
+ 
+   end subroutine Ini_Close
+ 
+ 
+   subroutine Ini_Close_File(Ini)
+     Type(TIniFile) :: Ini
+    
+     call TNameValueList_Clear(Ini%L)
+     call TNameValueList_Clear(Ini%ReadValues)
+ 
+   end  subroutine Ini_Close_File
+   
+ 
+ 
+   function Ini_Read_String(Key, NotFoundFail) result(AValue)
+    character (LEN=*), intent(IN) :: Key
+    logical, optional, intent(IN) :: NotFoundFail
+    character(LEN=Ini_max_string_len) :: AValue
+ 
+      if (present(NotFoundFail)) then
+       AValue = Ini_Read_String_File(DefIni, Key, NotFoundFail)
+      else
+       AValue = Ini_Read_String_File(DefIni, Key)
+      end if
+ 
+   end function Ini_Read_String
+ 
+ 
+   function Ini_Read_String_File(Ini, Key, NotFoundFail) result(AValue)
+    Type(TIniFile) :: Ini
+    character (LEN=*), intent(IN) :: Key
+    logical, optional, intent(IN) :: NotFoundFail
+    character(LEN=Ini_max_string_len) :: AValue
+ 
+    call TNameValueList_ValueOf(Ini%L, Key, AValue)
+ 
+    if (AValue/='') then
+ 
+     call  TNameValueList_Add(Ini%ReadValues, Key, AValue)
+     if (Ini_Echo_Read) write (*,*) trim(Key)//' = ',trim(AValue)
+     return
+ 
+    end if
+    if (present(NotFoundFail)) then
+       if (NotFoundFail) then
+          write(*,*) 'key not found : '//trim(Key)
+          stop
+       end if
+    else if (Ini_fail_on_not_found) then
+       write(*,*) 'key not found : '//trim(Key)
+       stop
+    end if
+ 
+   end function Ini_Read_String_File
+   function Ini_HasKey(Key) result(AValue)
+    character (LEN=*), intent(IN) :: Key
+    logical AValue
+ 
+    AValue = Ini_HasKey_File(DefIni, Key)
+ 
+   end function Ini_HasKey
+ 
+   function Ini_HasKey_File(Ini, Key) result(AValue)
+    type(TIniFile), intent(in) :: Ini
+    character (LEN=*), intent(IN) :: Key
+    logical AValue
+ 
+       Avalue = TNameValueList_HasKey(Ini%L, Key)
+       
+   end function Ini_HasKey_File
+ 
+  function Ini_Key_To_Arraykey(Key, index)  result(AValue)
+     character (LEN=*), intent(IN) :: Key
+     integer, intent(in) :: index
+     character(LEN=Ini_max_string_len) :: AValue
+     
+     character(LEN=32) :: numstr
+     write (numstr,*) index 
+     numstr=adjustl(numstr)
+     AValue = trim(Key) // '(' // trim(numStr) // ')' 
+  
+  end function Ini_Key_To_Arraykey
+ 
+   function Ini_Read_String_Array(Key, index, NotFoundFail) result(AValue)
+    character (LEN=*), intent(IN) :: Key
+    integer, intent(in) :: index
+    logical, optional, intent(IN) :: NotFoundFail
+    character(LEN=Ini_max_string_len) :: AValue
+ 
+      if (present(NotFoundFail)) then
+       AValue = Ini_Read_String_Array_File(DefIni, Key, index, NotFoundFail)
+      else
+       AValue = Ini_Read_String_Array_File(DefIni, Key, index)
+      end if
+ 
+   end function Ini_Read_String_Array
+ 
+   function Ini_Read_String_Array_File(Ini, Key, index, NotFoundFail) result(AValue)
+    Type(TIniFile) :: Ini
+    integer, intent(in) :: index
+    character (LEN=*), intent(IN) :: Key
+    logical, optional, intent(IN) :: NotFoundFail
+    character(LEN=Ini_max_string_len) :: AValue
+    character(LEN=Ini_max_string_len) :: ArrayKey
+    
+      ArrayKey = Ini_Key_To_Arraykey(Key,index)
+      if (present(NotFoundFail)) then
+       AValue = Ini_Read_String_File(Ini, ArrayKey, NotFoundFail)
+      else
+       AValue = Ini_Read_String_File(Ini, ArrayKey)
+      end if
+    
+   end function Ini_Read_String_Array_File
+ 
+   function Ini_Read_Int_Array(Key, index, Default)
+      integer, optional, intent(IN) :: Default
+      integer, intent(in) :: index
+      character (LEN=*), intent(IN) :: Key
+      integer Ini_Read_Int_Array
+ 
+      if (present(Default)) then
+       Ini_Read_Int_Array = Ini_Read_Int_Array_File(DefIni, Key, index, Default)
+      else
+       Ini_Read_Int_Array = Ini_Read_Int_Array_File(DefIni, Key, index)
+      end if
+      
+    end function Ini_Read_Int_Array
+ 
+   function Ini_Read_Int_Array_File(Ini,Key, index, Default)
+   !Reads Key(1), Key(2), etc.
+    Type(TIniFile) :: Ini
+    integer Ini_Read_Int_Array_File 
+    integer, optional, intent(IN) :: Default
+    integer, intent(in) :: index
+    character (LEN=*), intent(IN) :: Key
+    character(LEN=Ini_max_string_len) :: ArrrayKey
+      ArrrayKey = Ini_Key_To_Arraykey(Key,index)
+      if (present(Default)) then
+       Ini_Read_Int_Array_File = Ini_Read_Int_File(Ini, ArrrayKey, Default)
+      else
+       Ini_Read_Int_Array_File = Ini_Read_Int_File(Ini, ArrrayKey)
+      end if
+   end function Ini_Read_Int_Array_File
+ 
+ 
+   function Ini_Read_Int(Key, Default)
+      integer, optional, intent(IN) :: Default
+      character (LEN=*), intent(IN) :: Key
+      integer Ini_Read_Int
+ 
+      if (present(Default)) then
+       Ini_Read_Int = Ini_Read_Int_File(DefIni, Key, Default)
+      else
+       Ini_Read_Int = Ini_Read_Int_File(DefIni, Key)
+      end if
+   end function Ini_Read_Int
+ 
+   function Ini_Read_Int_File(Ini, Key, Default)
+    Type(TIniFile) :: Ini
+    integer Ini_Read_Int_File
+    integer, optional, intent(IN) :: Default
+    character  (LEN=*), intent(IN) :: Key
+    character(LEN=Ini_max_string_len) :: S
+    
+    S = Ini_Read_String_File(Ini, Key,.not. present(Default))
+    if (S == '') then
+       if (.not. present(Default)) then
+         write(*,*) 'no value for key: '//Key
+         stop
+       end if
+       Ini_Read_Int_File = Default
+       write (S,*) Default
+       call  TNameValueList_Add(Ini%ReadValues, Key, S)
+    else
+     if (verify(trim(S),'-+0123456789') /= 0) goto 10
+     read (S,*, err = 10) Ini_Read_Int_File
+    end if
+   return
+ 10 write (*,*) 'error reading integer for key: '//Key
+    stop
+   
+   end function Ini_Read_Int_File
+ 
+   function Ini_Read_Double(Key, Default)
+      double precision, optional, intent(IN) :: Default
+      character (LEN=*), intent(IN) :: Key
+      double precision Ini_Read_Double
+ 
+      if (present(Default)) then
+       Ini_Read_Double = Ini_Read_Double_File(DefIni, Key, Default)
+      else
+       Ini_Read_Double = Ini_Read_Double_File(DefIni, Key)
+      end if
+   
+   end function Ini_Read_Double
+ 
+   function Ini_Read_Double_File(Ini,Key, Default)
+    Type(TIniFile) :: Ini
+    double precision Ini_Read_Double_File 
+    double precision, optional, intent(IN) :: Default
+    character (LEN=*), intent(IN) :: Key
+    character(LEN=Ini_max_string_len) :: S
+    
+    S = Ini_Read_String_File(Ini,Key,.not. present(Default))
+    if (S == '') then
+       if (.not. present(Default)) then
+         write(*,*) 'no value for key: '//Key
+         stop
+       end if
+       Ini_Read_Double_File = Default
+       write (S,*) Default
+ 
+       call  TNameValueList_Add(Ini%ReadValues, Key, S)
+ 
+    else
+     read (S,*, err=10) Ini_Read_Double_File
+    end if
+ 
+   return
+ 
+ 10 write (*,*) 'error reading double for key: '//Key
+    stop
+ 
+   end function Ini_Read_Double_File
+ 
+ 
+ 
+     function Ini_Read_Double_Array(Key, index, Default)
+      double precision, optional, intent(IN) :: Default
+      integer, intent(in) :: index
+      character (LEN=*), intent(IN) :: Key
+      double precision Ini_Read_Double_Array
+ 
+      if (present(Default)) then
+       Ini_Read_Double_Array = Ini_Read_Double_Array_File(DefIni, Key, index, Default)
+      else
+       Ini_Read_Double_Array = Ini_Read_Double_Array_File(DefIni, Key, index)
+      end if
+      
+     end function Ini_Read_Double_Array
+ 
+ 
+   function Ini_Read_Double_Array_File(Ini,Key, index, Default)
+ 
+   !Reads Key(1), Key(2), etc.
+ 
+    Type(TIniFile) :: Ini
+ 
+    double precision Ini_Read_Double_Array_File 
+    double precision, optional, intent(IN) :: Default
+    integer, intent(in) :: index
+    character (LEN=*), intent(IN) :: Key
+    character(LEN=Ini_max_string_len) ::  ArrrayKey
+ 
+      ArrrayKey = Ini_Key_To_Arraykey(Key,index)
+      if (present(Default)) then
+ 
+       Ini_Read_Double_Array_File = Ini_Read_Double_File(Ini, ArrrayKey, Default)
+      else
+       Ini_Read_Double_Array_File = Ini_Read_Double_File(Ini, ArrrayKey)
+      end if
+   end function Ini_Read_Double_Array_File
+ 
+     function Ini_Read_Real(Key, Default)
+      real, optional, intent(IN) :: Default
+      character (LEN=*), intent(IN) :: Key
+      real Ini_Read_Real
+ 
+      if (present(Default)) then
+       Ini_Read_Real = Ini_Read_Real_File(DefIni, Key, Default)
+      else
+       Ini_Read_Real = Ini_Read_Real_File(DefIni, Key)
+      end if
+ 
+     end function Ini_Read_Real
+ 
+     function Ini_Read_Real_File(Ini,Key, Default)
+     Type(TIniFile) :: Ini
+     real Ini_Read_Real_File 
+     real, optional, intent(IN) :: Default
+     character (LEN=*), intent(IN) :: Key
+     character(LEN=Ini_max_string_len) :: S
+    
+     S = Ini_Read_String_File(Ini,Key,.not. present(Default))
+     if (S == '') then
+       if (.not. present(Default)) then
+         write(*,*) 'no value for key: '//Key
+         stop
+       end if
+       Ini_Read_Real_File = Default
+       write (S,*) Default
+       call  TNameValueList_Add(Ini%ReadValues, Key, S)
+ 
+    else
+     read (S,*, err=10) Ini_Read_Real_File
+    end if
+ 
+   return
+ 
+ 10 write (*,*) 'error reading double for key: '//Key
+    stop
+ 
+   end function Ini_Read_Real_File
+ 
+ 
+ 
+    function Ini_Read_Real_Array(Key, index, Default)
+      real, optional, intent(IN) :: Default
+      integer, intent(in) :: index
+      character (LEN=*), intent(IN) :: Key
+      real Ini_Read_Real_Array
+ 
+      if (present(Default)) then
+       Ini_Read_Real_Array = Ini_Read_Real_Array_File(DefIni, Key, index, Default)
+      else
+       Ini_Read_Real_Array = Ini_Read_Real_Array_File(DefIni, Key, index)
+      end if
+    end function Ini_Read_Real_Array
+ 
+   function Ini_Read_Real_Array_File(Ini,Key, index, Default)
+   !Reads Key(1), Key(2), etc.
+    Type(TIniFile) :: Ini
+    real Ini_Read_Real_Array_File 
+    real, optional, intent(IN) :: Default
+    integer, intent(in) :: index
+    character (LEN=*), intent(IN) :: Key
+    character(LEN=Ini_max_string_len) :: ArrrayKey
+    
+      ArrrayKey = Ini_Key_To_Arraykey(Key,index)
+      if (present(Default)) then
+       Ini_Read_Real_Array_File = Ini_Read_Real_File(Ini, ArrrayKey, Default)
+      else
+       Ini_Read_Real_Array_File = Ini_Read_Real_File(Ini, ArrrayKey)
+      end if
+   end function Ini_Read_Real_Array_File
+ 
+   function Ini_Read_Logical(Key, Default)
+      Logical, optional, intent(IN) :: Default
+      character (LEN=*), intent(IN) :: Key
+     logical Ini_Read_Logical
+ 
+      if (present(Default)) then
+       Ini_Read_Logical = Ini_Read_Logical_File(DefIni, Key, Default)
+      else
+       Ini_Read_Logical = Ini_Read_Logical_File(DefIni, Key)
+      end if
+   end function Ini_Read_Logical
+ 
+   function Ini_Read_Logical_File(Ini, Key, Default)
+    Type(TIniFile) :: Ini
+ 
+    logical Ini_Read_Logical_File
+    logical, optional, intent(IN) :: Default
+    character  (LEN=*), intent(IN) :: Key
+   
+    character(LEN=Ini_max_string_len) :: S
+    
+    S = Ini_Read_String_File(Ini,Key,.not. present(Default))
+    if (S == '') then
+       if (.not. present(Default)) then
+         write(*,*) 'no value for key: '//Key
+         stop
+       end if
+       Ini_Read_Logical_File = Default
+       write (S,*) Default
+ 
+       call  TNameValueList_Add(Ini%ReadValues, Key, S)
+ 
+    else
+ 
+     if (verify(trim(S),'10TF') /= 0) goto 10  
+     read (S,*, err = 10) Ini_Read_Logical_File
+    end if
+ 
+   return
+ 
+ 10 write (*,*) 'error reading logical for key: '//Key
+    stop
+   end function Ini_Read_Logical_File
+ 
+ 
+ 
+   subroutine Ini_SaveReadValues(afile,unit_id)
+    character(LEN=*)  :: afile
+    integer, intent(in) :: unit_id
+ 
+    call Ini_SaveReadValues_File(DefIni, afile, unit_id)
+ 
+   end subroutine Ini_SaveReadValues
+ 
+ 
+ 
+   subroutine Ini_SaveReadValues_File(Ini, afile, unit_id)
+    Type(TIniFile) :: Ini
+    character(LEN=*), intent(in) :: afile
+    integer, intent(in) :: unit_id
+    integer i
+ 
+    open(unit=unit_id,file=afile,form='formatted',status='replace', err=500)
+ 
+    do i=1, Ini%ReadValues%Count
+ 
+     write (unit_id,'(a)') trim(Ini%ReadValues%Items(i)%P%Name) // ' = ' &
+                         //trim(Ini%ReadValues%Items(i)%P%Value)
+ 
+    end do
+ 
+    close(unit_id)
+    return
+ 
+ 500 write(*,*) 'Ini_SaveReadValues_File: Error creating '//trim(afile)
+ 
+   end subroutine Ini_SaveReadValues_File
+ 
+ end module IniFile
+ 
diff -r -c -b -B -N cosmomc/source/Makefile cosmomc_sampler/source/Makefile
*** cosmomc/source/Makefile	2011-07-11 18:00:59.000000000 +0200
--- cosmomc_sampler/source/Makefile	2011-08-03 11:41:35.164647841 +0200
***************
*** 1,218 ****
! #You may need to edit the library paths for MKL for Intel
! #Beware of using optmizations that lose accuracy - may give errors when running
  
! ##Uncomment the next line to include dr7 LRG
! EXTDATA =
! #EXTDATA = LRG
! 
! #set WMAP empty not to compile with WMAP
! WMAP = /home/aml1005/WMAP7/likelihood_v4
! 
! #Only needed for WMAP
! cfitsio = /usr/local/cfitsio/intel10/64/3.040
! 
! #GSL only needed for DR7 LRG
! GSLPATH = /home/aml1005/libs/gsl
! 
! IFLAG = -I
! INCLUDE=
! 
! #Intel MPI (assuming mpif77 set to point to ifort)
! #these settings for ifort 11.1 and higher; may need to add explicit link directory otherwise
! F90C     = mpif90
! FFLAGS = -O2 -ip -W0 -WB -openmp -fpp -DMPI -vec_report0 -mkl=parallel
! LAPACKL = -lmkl_lapack
! 
! #HPCF settings. Use Inteal 9 or 10.1+, not 10.0
! #F90C     = mpif90
! #FFLAGS = -O2 -Vaxlib -W0 -WB -openmp -fpp -DMPI -vec_report0
! #LAPACKL = -L/usr/local/Cluster-Apps/intel/mkl/10.2.2.025/lib/em64t -lmkl_lapack -lmkl -lguide -lpthread
! #GSLPATH = /usr/local/Cluster-Apps/gsl/1.9
! #cfitsio = /usr/local/Cluster-Users/cpac/cmb/2.1.0/cfitsio
! #INCLUDE=
! 
! #COSMOS: use "module load cosmolib latest"
! #use "runCosmomc" (globally installed) to run, defining required memory usage
! ifeq ($(COSMOHOST),cosmos)
! F90C = ifort
! FFLAGS = -openmp -O3 -w -fpp2 -DMPI
! LAPACKL = -mkl=sequential -lmkl_lapack -lmpi
! cfitsio = $(CFITSIO)
! WMAP = $(COSMOLIB)/WMAP7
! GSLPATH = $(GSL_ROOT)
  endif
  
! #Intel fortran 8, check you have the latest update from the Intel web pages
! #See Makefile_intel for ifc 7.1 or lower (some versions have problems)
! #F90C     = ifort
! #FFLAGS = -O2 -Vaxlib -ip -W0 -WB -openmp -fpp
! #LAPACKL = -L/opt/intel/mkl61/lib/32 -lmkl_lapack -lmkl_ia32 -lguide -lpthread
! 
! #G95; make sure LAPACK and MPI libs also compiled with g95
! #F90C     = mpif90
! #FFLAGS = -O2 -cpp -DMPI
! #LAPACKL = /LAPACK/lapack_LINUX.a /LAPACK/blas_LINUX.a
! 
! #GFortran: if pre v4.3 add -D__GFORTRAN__
! #F90C     = gfortran
! #FFLAGS = -O2 -ffree-form -x f95-cpp-input
! #LAPACKL = -Wl,-framework -Wl,accelerate
! #may need to delete -Wl,accelerate, and/or add -shared
! 
! #SGI, -mp toggles multi-processor. Use -O2 if -Ofast gives problems.
! #Not various versions of the compiler are buggy giving erroneous seg faults with -mp.
! #Version 7.3 is OK, currently version 7.4 is bugged, as are some earlier versions.
! #F90C     = f90
! #LAPACKL = -lcomplib.sgimath
! #FFLAGS  = -Ofast -mp
! 
! #Digital/Compaq fortran, -omp toggles multi-processor
! #F90C    = f90
! #FFLAGS  = -omp -O4 -arch host -math_library fast -tune host -fpe1
! #LAPACKL  = -lcxml
! 
! #Absoft ProFortran, single processor, set -cpu:[type] for your local system
! #F90C     = f95
! #FFLAGS = -O2 -s -cpu:athlon -lU77 -w -YEXT_NAMES="LCS" -YEXT_SFX="_"
! #LAPACKL =  -llapack -lblas -lg2c
! #IFLAG = -p
! 
! #NAGF95, single processor:
! #F90C     = f95
! #FFLAGS = -DNAGF95 -O3
! #LAPACKL = -llapack -lblas -lg2c
! 
! #PGF90
! #F90C = pgf90
! #FFLAGS = -O2 -DESCAPEBACKSLASH
! #LAPACKL = -llapack -lblas
! 
! 
! #Sun, single processor:
! #F90C     = f90
! #FFLAGS = -fast -ftrap=%none
! #LAPACKL = -lsunperf -lfsu
! #LAPACKL = -lsunperf -lfsu -lsocket -lm
! #IFLAG = -M
! 
! #Sun MPI
! #F90C = mpf90
! #FFLAGS =  -O4 -openmp -ftrap=%none -dalign -DMPI
! #LAPACKL = -lsunperf -lfsu -lmpi_mt
! #IFLAG = -M
! 
! #Sun parallel enterprise:
! #F90C     = f95
! #FFLAGS =  -O4 -xarch=native64 -openmp -ftrap=%none
! #LAPACKL = -lsunperf -lfsu
! #IFLAG = -M
! 
! 
! #IBM XL Fortran, multi-processor (run "module load lapack" then run "gmake")
! # See also http://cosmocoffee.info/viewtopic.php?t=326
! #F90C     = xlf90_r $(LAPACK)
! #FFLAGS  = -WF,-DIBMXL -qsmp=omp -qsuffix=f=f90:cpp=F90 -O3 -qstrict -qarch=pwr3 -qtune=pwr3
! #INCLUDE = -lessl
! #LAPACKL =
  
! PROPOSE = propose.o
! CLSFILE = CMB_Cls_simple.o
  
! #Can use params_H if you prefer more generic parameters
! PARAMETERIZATION = params_CMB.o
  
! F90FLAGS = -DMATRIX_SINGLE $(FFLAGS) $(IFLAG)../camb $(INCLUDE)
! LINKFLAGS = -L../camb -lcamb $(LAPACKL)
  
- DISTFILES = ParamNames.o Matrix_utils.o settings.o IO.o GetDist.o
  
! OBJFILES=  ParamNames.o Matrix_utils.o settings.o IO.o cmbtypes.o Planck_like.o  \
! 	cmbdata.o WeakLen.o bbn.o bao.o lrggettheory.o mpk.o supernovae.o HST.o SDSSLy-a-v3.o \
! 	$(CLSFILE) paramdef.o $(PROPOSE) $(PARAMETERIZATION) calclike.o \
! 	conjgrad_wrapper.o EstCovmat.o postprocess.o MCMC.o driver.o
  
  
! ifeq ($(EXTDATA),LRG)
! F90FLAGS += -DDR71RG
! OBJFILES += bsplinepk.o
! GSLINC = -I$(GSLPATH)/include
! LINKFLAGS += -L$(GSLPATH)/lib -lgsl -lgslcblas
! endif
  
! ifneq ($(WMAP),)
! F90FLAGS += $(IFLAG)$(cfitsio)/include $(IFLAG)$(WMAP)
! LINKFLAGS +=  -L$(cfitsio)/lib -L$(WMAP) -lcfitsio
! 
! OBJFILES += $(WMAP)/read_archive_map.o \
! 	$(WMAP)/read_fits.o \
! 	$(WMAP)/healpix_types.o \
! 	$(WMAP)/WMAP_7yr_options.o \
! 	$(WMAP)/WMAP_7yr_util.o \
! 	$(WMAP)/WMAP_7yr_tt_pixlike.o \
! 	$(WMAP)/WMAP_7yr_teeebb_pixlike.o \
! 	$(WMAP)/WMAP_7yr_likelihood.o \
! 	$(WMAP)/WMAP_7yr_gibbs.o \
! 	$(WMAP)/WMAP_7yr_tt_beam_ptsrc_chisq.o \
! 	$(WMAP)/br_mod_dist.o
! else
! F90FLAGS += -DNOWMAP
! endif
  
! default: cosmomc
  
! all : cosmomc getdist
  
! settings.o: ../camb/libcamb.a
! cmbtypes.o: settings.o
! Planck_like.o: cmbtypes.o
! cmbdata.o: Planck_like.o  $(WMAPOBJS)
! WeakLen.o: cmbtypes.o
! bbn.o: settings.o
! bao.o: cmbtypes.o
! mpk.o: cmbtypes.o lrggettheory.o
! HST.o: cmbtypes.o
! supernovae.o: cmbtypes.o
! SDSSLy-a-v3.o: cmbtypes.o
! $(CLSFILE): cmbtypes.o mpk.o HST.o bao.o IO.o
! paramdef.o: $(CLSFILE)
! $(PROPOSE): paramdef.o
! $(PARAMETERIZATION): paramdef.o
! calclike.o: $(PARAMETERIZATION)
! conjgrad_wrapper.o: calclike.o
! EstCovmat.o: conjgrad_wrapper.o
! postprocess.o: calclike.o
! MCMC.o: calclike.o
! driver.o: MCMC.o
  
! export FFLAGS
! export F90C
  
  .f.o:
  	f77 $(F90FLAGS) -c $<
  
- %.o: %.c
- 	$(CC) $(GSLINC) -c $*.c
- 
  %.o: %.f90
  	$(F90C) $(F90FLAGS) -c $*.f90
  
  %.o: %.F90
  	$(F90C) $(F90FLAGS) -c $*.F90
  
  
! cosmomc: camb $(OBJFILES)
! 	$(F90C) -o ../cosmomc $(OBJFILES) $(LINKFLAGS) $(F90FLAGS)
  
! 
! clean: cleancosmomc
! 	rm -f ../camb/*.o ../camb/*.obj ../camb/*.mod
! 
! cleancosmomc:
  	rm -f *.o *.mod *.d *.pc *.obj ../core
  
! 
! getdist: camb $(DISTFILES)
! 	$(F90C) -o ../getdist $(DISTFILES) $(LINKFLAGS) $(F90FLAGS)
! 
! camb:
! 	cd ../camb && $(MAKE) --file=Makefile_main libcamb.a
--- 1,85 ----
! # >>> DESIGNED FOR GMAKE <<<
  
! # Unified Systems makefile for GENEMC
! # Add FLAGS -DMPI for using MPI
! 
! ext=$(shell uname | cut -c1-3)
! 
! ifeq ($(ext),IRI)
! F90C= f90
! FFLAGS= -Ofast -mp -n32 -LANG:recursive=ON -lmpi -DMPI
! LAPACKL = -lcomplib.sgimath_mp
! INCLUDE = 
  endif
  
! ifeq ($(ext),Lin)
! F90C=gfortran
! FFLAGS= -O -fopenmp 
! LAPACKL = -llapack
! INCLUDE = -I/usr/include 
! endif 
  
! ifeq ($(ext),OSF)
! F90C=f90
! FFLAGS= -omp -O -arch host -math_library fast -tune host -fpe1
! WMAPFLAGS= $(FFLAGS)
! LAPACKL  = -lcxml
! INCLUDE = 
! endif
  
! ifeq ($(ext),Sun)
! F90C=f90
! FFLAGS= -O4 -xarch=native64 -openmp -ftrap=%none
! WMAPFLAGS= $(FFLAGS)
! LAPACKL = -lsunperf -lfsu
! INCLUDE =
! endif
  
! ifeq ($(ext),AIX)
! F90C    = mpxlf90_r
! FFLAGS  = -O4 -WF,-DIBMXL,-DMPI -qstrict -qsmp=omp -qmaxmem=-1 -qsuffix=f=f90:cpp=F90
! LAPACKL = -lessl
! INCLUDE = 
! endif
  
  
! PROPOSE = propose.o
  
+ PARAMETERIZATION = params_EXT.o
  
! MCMCLIB = $(LAPACKL)
  
! DISTLIB = $(LAPACKL)
  
! F90FLAGS = -DMATRIX_SINGLE $(FFLAGS) $(INCLUDE)
  
! DISTFILES = utils.o inifile.o ParamNames.o Matrix_utils.o settings.o GetDist.o
  
! OBJFILES= inifile.o utils.o ParamNames.o Matrix_utils.o  settings.o IO.o \
! 	paramdef.o wrapper.o $(PROPOSE) $(PARAMETERIZATION) calclike.o \
! 	conjgrad_wrapper.o EstCovmat.o postprocess.o MCMC.o driver.o
  
! default: genemc.$(ext)
! 
! all : genemc.$(ext) getdist.$(ext)
  
  .f.o:
  	f77 $(F90FLAGS) -c $<
  
  %.o: %.f90
  	$(F90C) $(F90FLAGS) -c $*.f90
  
  %.o: %.F90
  	$(F90C) $(F90FLAGS) -c $*.F90
  
+ utils.o:
+ 	$(F90C) $(F90FLAGS) -c ../camb/utils.F90
  
! genemc.$(ext): $(OBJFILES)
! 	$(F90C) -o ../$@ $(OBJFILES) $(MCMCLIB) $(F90FLAGS)
  
! clean:
  	rm -f *.o *.mod *.d *.pc *.obj ../core
  
! getdist.$(ext): $(DISTFILES)
! 	$(F90C) -o ../$@ $(DISTFILES) $(DISTLIB) $(F90FLAGS) 
diff -r -c -b -B -N cosmomc/source/MCMC.f90 cosmomc_sampler/source/MCMC.f90
*** cosmomc/source/MCMC.f90	2010-05-10 17:28:00.000000000 +0200
--- cosmomc_sampler/source/MCMC.f90	2011-08-03 11:22:30.778800558 +0200
***************
*** 8,13 ****
--- 8,14 ----
   use CalcLike
   use Random
   use propose
+  use wrapper
   use IO
   implicit none
  
diff -r -c -b -B -N cosmomc/source/paramdef.F90 cosmomc_sampler/source/paramdef.F90
*** cosmomc/source/paramdef.F90	2011-04-08 13:21:21.000000000 +0200
--- cosmomc_sampler/source/paramdef.F90	2011-08-03 11:39:27.200706843 +0200
***************
*** 1,16 ****
  !Defines a parameterization, computes likelihoods and defined proposal density
  !Change this file to change these things, and MPI updating
- !The Cls are computed using routines defined in the CMB_Cls module
  
  module ParamDef
    !module defines a parameterization, and works out power spectra
!  use CMB_Cls
!  use cmbtypes
!  use cmbdata
   use Random
   use settings
   implicit none
  
   Type ParamSet
      real P(num_params)
      Type(ParamSetInfo) Info
--- 1,20 ----
  !Defines a parameterization, computes likelihoods and defined proposal density
  !Change this file to change these things, and MPI updating
  
  module ParamDef
    !module defines a parameterization, and works out power spectra
! ! use cmbtypes
! ! use cmbdata
   use Random
   use settings
+  use IO
   implicit none
  
+  Type ParamSetInfo  
+     integer :: infoInt
+     real :: infoReal
+  end Type ParamSetInfo
+ 
   Type ParamSet
      real P(num_params)
      Type(ParamSetInfo) Info
***************
*** 152,158 ****
          
          Ini_fail_on_not_found = .false.
  
!         call CMB_Initialize(Params%Info)
  
          prop_mat = trim(Ini_Read_String_File(Ini,'propose_matrix'))
          has_propose_matrix = prop_mat /= ''
--- 156,162 ----
          
          Ini_fail_on_not_found = .false.
  
! !        call CMB_Initialize(Params%Info)
  
          prop_mat = trim(Ini_Read_String_File(Ini,'propose_matrix'))
          has_propose_matrix = prop_mat /= ''
***************
*** 369,410 ****
  end subroutine SetProposeMatrix
  
  
- 
-  
-   subroutine WriteCMBParams(CMB,Theory,mult,like, with_data)     
-      Type (CosmoTheory) Theory 
-      real, intent(in) :: mult, like
-      logical, intent(in), optional :: with_data
-      Type(ParamSet) P
-      Type(CMBParams) CMB
- 
- 
-      call CMBParamsToParams(CMB,P%P)
- 
-      P%Info%Theory = Theory
-      if (present(with_data)) then
-       if (with_data) then
-       call WriteParamsAndDat(P, mult,like)
-       else
-       call WriteParams(P, mult,like)
-       end if
-      else
-       call WriteParams(P, mult,like)
-      end if
- 
-   end subroutine WriteCMBParams
-    
    subroutine WriteIndepSample(P, like)
      Type(ParamSet) P
      real like
!     Type(CMBParams) C
      if (indepfile_handle ==0) return
!     call ParamsToCMBParams(P%P,C)
!     call WriteModel(indepfile_handle, C,P%Info%Theory,like)
     end subroutine WriteIndepSample
  
  
     subroutine AddMPIParams(P,like, checkpoint_start)     
       real, intent(in) ::like
       real P(:)
       logical, intent(in), optional :: checkpoint_start
--- 373,390 ----
  end subroutine SetProposeMatrix
  
  
    subroutine WriteIndepSample(P, like)
      Type(ParamSet) P
      real like
! 
      if (indepfile_handle ==0) return
! 
     end subroutine WriteIndepSample
  
  
     subroutine AddMPIParams(P,like, checkpoint_start)     
+      use IO
+      use MAtrixUtils
       real, intent(in) ::like
       real P(:)
       logical, intent(in), optional :: checkpoint_start
***************
*** 739,744 ****
--- 719,725 ----
  #ifdef MPI
  
   subroutine CheckLimitsConverge(L)
+    use IO
    !Check limits from last half chains agree well enough across chains to be confident of result
    !Slowly explored tails will cause problems (long time till stops)
    Type(TList_RealArr), intent(in) :: L
diff -r -c -b -B -N cosmomc/source/params_EXT.f90 cosmomc_sampler/source/params_EXT.f90
*** cosmomc/source/params_EXT.f90	1970-01-01 01:00:00.000000000 +0100
--- cosmomc_sampler/source/params_EXT.f90	2011-08-03 11:27:47.503819235 +0200
***************
*** 0 ****
--- 1,63 ----
+    subroutine SetParamNames(Names)
+     use settings
+     use ParamNames
+     Type(TParamNames) :: Names
+  
+     if (generic_mcmc) then
+      Names%nnames=0
+      if (Feedback>0) write (*,*) 'edit SetParamNames in params_EXT.f90 if you want to use named params'
+ !     call ParamNames_init(Names, './params_EXT.paramnames')
+   end if
+ end subroutine SetParamNames
+ 
+ 
+   subroutine WriteParams(P, mult, like)
+      use settings 
+      use ParamDef
+      use IO
+      implicit none
+      Type(ParamSet) P
+      real, intent(in) :: mult, like
+      real r10
+      real, allocatable :: output_array(:)
+   
+     if (outfile_handle ==0) return
+   
+     if (generic_mcmc) then
+ 
+       call IO_OutputChainRow(outfile_handle, mult, like, P%P)
+      
+     else
+     
+        stop
+      
+     end if
+ 
+   end  subroutine WriteParams
+ 
+ 
+ 
+ 
+   subroutine WriteParamsAndDat(P, mult, like)
+      use settings
+      use ParamDef
+      use IO
+      implicit none
+     Type(ParamSet) P
+     real, intent(in) :: mult, like
+     character(LEN =30) fmt
+     real r10
+     real,allocatable :: output_array(:)
+        
+     if (outfile_handle ==0) return
+ 
+     allocate(output_array(num_real_params))
+     output_array(1:num_real_params) =  P%P(1:num_real_params)
+      
+     call IO_OutputChainRow(outfile_handle, mult, like, output_array)
+     deallocate(output_array)
+     
+   end  subroutine WriteParamsAndDat
+ 
+ 
+  
diff -r -c -b -B -N cosmomc/source/postprocess.f90 cosmomc_sampler/source/postprocess.f90
*** cosmomc/source/postprocess.f90	2010-05-10 17:28:00.000000000 +0200
--- cosmomc_sampler/source/postprocess.f90	2011-08-03 11:22:57.187551578 +0200
***************
*** 7,15 ****
  
  module posthoc
    use settings
-   use cmbtypes
-   use CMB_Cls
-   use cmbdata
    use CalcLike
    implicit none
  
--- 7,12 ----
***************
*** 71,287 ****
     subroutine postprocess(InputFile)
          use IO
          character(LEN=*), intent(INOUT):: InputFile
-         Type(CMBParams) LastCMB,CMB, newCMB
-         Type(CosmoTheory) Theory, CorrectTheory
-         Type(ParamSetInfo) Info
-         real Cls(lmax,num_cls_tot), truelike,mult,like
-         real weight_min, weight_max, mult_sum, mult_ratio, mult_max,weight
  
!         real max_like, max_truelike
!         integer error,num, debug
!         character (LEN=120) :: post_root
!         real Params(num_params)
!         integer infile_handle, outdata_handle
! 
!         flush_write = .false.
!         weight_min= 1e30
!         weight_max = -1e30
!         mult_sum = 0
!         mult_ratio = 0
!         mult_max = -1e30
!         max_like = 1e30
! 
!         max_truelike =1e30
! 
!         debug = 0
!         Info%Theory%Sn_LogLike = 0
!         Info%Theory%HST_LogLike = 0
!         Info%Theory%BAO_LogLike = 0
! 
!         infile_handle = 0
!         Temperature = PostParams%redo_temperature 
! 
!         if (Feedback>0 .and. PostParams%redo_change_like_only) &
!               write (*,*) 'Warning: only changing likelihoods not weights'
! 
!         if (PostParams%redo_datafile /= '') InputFile = PostParams%redo_datafile
! 
!         if (PostParams%redo_from_text) then
! 
!          infile_handle = IO_OpenChainForRead(trim(InputFile)//'.txt')
!          
!          if (.not. PostParams%redo_theory) write (*,*) '**You probably want to set redo_theory**'
!          if (.not. PostParams%redo_cls .and. Use_CMB) write (*,*) '**You probably want to set redo_cls**'
!          if (.not. PostParams%redo_pk .and. Use_mpk) write (*,*) '**You probably want to set redo_pk**'
! 
!          if (PostParams%redo_thin>1) write (*,*) 'redo_thin only OK with redo_from_text if input weights are 1'
! 
!         else
!          infile_handle = IO_OpenDataForRead(trim(InputFile)//'.data')
!         end if
! 
!          
!         if (PostParams%redo_outroot == '') then
!           post_root = trim(ExtractFilePath(InputFile))//'post_'// trim(ExtractFileName(InputFile))
!         else
!           post_root = PostParams%redo_outroot
!           if (instance /= 0) post_root = numcat(trim(post_root)//'_',instance)
!         end if
! 
!         if (Feedback > 0) then
!           if (PostParams%redo_from_text) then
!            write (*,*) 'reading from: ' //  trim(InputFile)//'.txt'
!           else
!            write (*,*) 'reading from: ' //  trim(InputFile)//'.data'
!           end if
!            write (*,*) 'writing to: ' // trim(post_root)//'.*'
!         end if
! 
!         write (*,*) 'Using temperature: ', Temperature
! 
!         outfile_handle = IO_OutputOpenForWrite(trim(post_root)//'.txt')
!         
!         outdata_handle = IO_DataOpenForWrite(trim(post_root)//'.data')
!      
!         num = 0
!         do
! 
!         if (PostParams%redo_from_text) then
!           error = 0
!           if (.not. IO_ReadChainRow(infile_handle, mult, like, Params, num_params)) exit
!           call ParamsToCMBParams(Params, CMB)
! 
!         else
!           call ReadModel(infile_handle,CMB,Theory,mult,like, error)
!         end if
! 
!            
!         if (error ==1) then
!           if (num==0) call MpiStop('Error reading data file.')
! 
!           exit
!         end if
!         num=num+1
!         if (num>PostParams%redo_skip .and. mod(num,PostParams%redo_thin) == 0) then
!      
!         LastCMB = CMB
! 
!         newCMB = CMB  
!         if (PostParams%redo_theory) then
!     
!          call GetClsInfo(newCMB, CorrectTheory, error, PostParams%redo_cls, PostParams%redo_pk)
!     
!          if (PostParams%redo_cls) then
!             Theory%cl = CorrectTheory%cl
!             Theory%cl_tensor = CorrectTheory%cl_tensor
!             !!In last version redo_cls just for going to higher l on temperature
! 
!          end if
!   
!          if (PostParams%redo_pk) then
!             Theory%sigma_8 = CorrectTheory%sigma_8
!             Theory%Matter_Power = CorrectTheory%Matter_Power
!          end if
! 
!          Theory%Age = CorrectTheory%Age
!      
!         else
!           error = 0
!         end if
! 
! 
!           CorrectTheory = Theory
! 
! 
!          if (error ==0) then
!  
!            if (PostParams%redo_like) then
!                
! 
!               if (Use_LSS .and. CorrectTheory%sigma_8==0) &
!                  call MpiStop('Matter power/sigma_8 have not been computed. Use redo_theory and redo_pk.')
! 
! 
!               call ClsFromTheoryData(CorrectTheory, newCMB, Cls)
!               Info%Theory = CorrectTheory
!             
!               if (any(Cls(2:lmax,1) < 0)) then
!                  write (*,*) 'WARNING: bad model with C_l < 0 being rejected'
!                  write (*,*) 'in '//trim(InputFile)
!                  !This shouldn't happen. 
!                  !But good sanity check, esp when playing around with extended models
!                  !Think have fixed problems with Om_k \sim -7e-4, w\sim -0.7.
!                  write (*,*) newCMB
!                  write (*,*)
!                  truelike = logZero
!               else
!                truelike = GetLogLikePost(newCMB, Info, Cls,.true.)
!               end if
!               if (truelike == logZero) then
!                weight = 0
!               else
!                if (PostParams%redo_add) truelike = like + truelike
! 
!                weight = exp(like-truelike+PostParams%redo_likeoffset)
!               end if
!            
!              if (.not. PostParams%redo_change_like_only)  mult = mult*weight
!            else
!               truelike = like
!               weight = 1
!            end if
! 
! 
!            max_like = min(max_like,like)
! 
!            max_truelike = min(max_truelike,truelike)
! 
!            mult_ratio = mult_ratio + weight
!            mult_sum = mult_sum + mult
!           
!            if (mult /= 0) then        
!             call WriteCMBParams(newCMB, CorrectTheory, mult, truelike,txt_theory)
!             call WriteModel(outdata_handle, newCMB, CorrectTheory,truelike,mult)
! 
!            else 
! 
!             if (Feedback >1 ) write (*,*) 'Zero weight: new like = ', truelike
!            end if
!           
!            if (Feedback > 1) write (*,*) 'done ', num, 'mult= ', mult,' weight = ', weight
!            weight_max = max(weight,weight_max)       
!            weight_min = min(weight,weight_min)
!            mult_max = max(mult_max,mult)
! 
!         end if
!         end if
!         
!         end do
! 
!         call IO_Close(infile_handle)
!         call IO_Close(outfile_handle)
!         call IO_DataCloseWrite(outdata_handle)
!       
!         num = (num - PostParams%redo_skip) / PostParams%redo_thin
!         if (Feedback>0) then 
!            write(*,*) 'finished. Processed ',num,' models'
!            write (*,*) 'max weight= ',weight_max, ' min weight = ',weight_min
!            write (*,*) 'mean mult  = ', mult_sum/num
!            write (*,*) 'mean importance weight (approx evidence ratio) = ',mult_ratio/num
!            write (*,*) 'effective number of samples =',mult_sum/mult_max
!            write (*,*) 'Best redo_likeoffset = ',max_truelike - max_like
!          end if
!         
!         if ((mult_ratio < 1e-6 .or. mult_ratio > 1e8) & 
! 
!             .and. .not.PostParams%redo_change_like_only) then
! 
!          write (*,*) 'WARNING: use redo_likeoffset to rescale likelihoods'
!         end if
! 
! 
! 
!         return
  
          end subroutine postprocess
  
--- 68,75 ----
     subroutine postprocess(InputFile)
          use IO
          character(LEN=*), intent(INOUT):: InputFile
  
!      stop 'postprocess not implemented!'
  
          end subroutine postprocess
  
diff -r -c -b -B -N cosmomc/source/settings.f90 cosmomc_sampler/source/settings.f90
*** cosmomc/source/settings.f90	2011-03-30 19:50:48.000000000 +0200
--- cosmomc_sampler/source/settings.f90	2011-08-03 11:23:04.257217218 +0200
***************
*** 43,49 ****
  
    real, parameter :: cl_norm = 1e-10 !units for As
  
!   logical, parameter ::  generic_mcmc= .false. 
      !set to true to not call CAMB, etc.
      !write GenericLikelihoodFunction in calclike.f90   
  
--- 43,49 ----
  
    real, parameter :: cl_norm = 1e-10 !units for As
  
!   logical, parameter ::  generic_mcmc= .true. 
      !set to true to not call CAMB, etc.
      !write GenericLikelihoodFunction in calclike.f90   
  
diff -r -c -b -B -N cosmomc/source/wrapper.f90 cosmomc_sampler/source/wrapper.f90
*** cosmomc/source/wrapper.f90	1970-01-01 01:00:00.000000000 +0100
--- cosmomc_sampler/source/wrapper.f90	2009-04-30 16:51:19.000000000 +0200
***************
*** 0 ****
--- 1,23 ----
+ module wrapper
+   use paramdef, only : paramset, paramsetinfo
+   implicit none
+   
+ 
+ contains
+ 
+ 
+   subroutine AcceptReject(accpt, CurParams, Trial)
+     logical, intent(in) :: accpt
+     Type(ParamSetInfo) CurParams, Trial
+     
+     
+     if (accpt) then
+        print *,'accepted'
+     else
+        print *,'rejected'
+     end if
+     
+   end subroutine AcceptReject
+ 
+ end module wrapper
+ 
