Ok now is working. For those interested in what i did here is my Likelihood :
________________________________________________________________________
module boost_estimator
use CosmologyTypes
use CosmoTheory
use Likelihood_Cosmology
implicit none
private
type, extends(TCosmoCalcLikelihood) :: boost_estimatorLikelihood
integer(kind=4) :: lmaxtt = 1024
real(mcp), allocatable :: ell(:)
real(mcp), allocatable :: ClsTT(:)
real(mcp), allocatable :: ClsTTcross(:)
real(mcp), allocatable :: ClsTE(:)
real(mcp), allocatable :: ClsTE_err(:)
real(mcp), allocatable :: ClsEE(:)
real(mcp), allocatable :: ClsEE_err(:)
contains
procedure :: boost_estimator_ReadCls
procedure :: LogLike => boost_estimator_LnLike
end type boost_estimatorLikelihood
public boost_estimatorLikelihood, boost_estimatorLikelihood_Add
contains
!Subroutine area ................... here we want to reade the .ini likelihood file from CosmoMC
subroutine boost_estimatorLikelihood_Add(LikeList, ini)
class(TLikelihoodList) LikeList
class(TSettingIni) Ini
Type (boost_estimatorLikelihood), pointer :: this
character(len=256) :: ClTT_filename = ''
character(LEN=:), allocatable :: string
if (Ini%Read_Logical('boost_use_CMB_TT_only',.false.)) then !EMS there should be use_CMB_only=T inside .ini file
allocate(this)
this%LikelihoodType = 'boost_estimator'
this%name = Ini%Read_String('boost_CMB_TT_only_name')
call Ini%Read('boost_CMB_TT_only_lmax',this%lmaxtt)
string = Ini%Read_String('boost_CMB_TT_only_fname')
if (string /= '') call this%boost_estimator_ReadCls(string)
this%needs_background_functions = .true.
call this%loadParamNames(trim(DataDir)//'boost.paramnames')
call LikeList%Add(this)
end if
end subroutine boost_estimatorLikelihood_Add
!Subroutine area ................... here we want to read the power spectrum
subroutine boost_estimator_ReadCls(this, ClTT_filename)
class(boost_estimatorLikelihood) :: this
character(LEN=:), allocatable :: ClTT_filename
integer j, idx, m
integer :: IO = 0
real(mcp) c1value
real(mcp) c2value
real(mcp) c3value
real(mcp) err
real(mcp) err2
real(mcp) err3
allocate ( this%ClsTT(this%lmaxtt+1) )
allocate ( this%ClsTTcross(this%lmaxtt+1) )
allocate ( this%ClsEE(this%lmaxtt+1) )
allocate ( this%ClsEE_err(this%lmaxtt+1) )
allocate ( this%ClsTE(this%lmaxtt+1) )
allocate ( this%ClsTE_err(this%lmaxtt+1) )
print*,'POWER ',trim(ClTT_filename)
open(unit = 100, file = ClTT_filename, status = 'old', action = 'read')
!do j = 0, this%lmaxtt
!this%ClsTT(j) = -1.d0
!this%ClsTTcross(j) = 0.d0
!this%ClsEE(j) = -1.d0
!this%ClsEE_Err(j) = 0.d0
!this%ClsTE(j) = -1.d0
!this%ClsTE_Err(j) = 0.d0
!end do
do while (IO == 0)
read(100, *, IOSTAT = IO) idx, c1value
if (IO == 0 .and. idx<=this%lmaxtt) then
this%ClsTT(idx) = c1value
!this%ClsTTcross(idx) = err
!print *, 'Power Spectra', idx, this%ClsTT(idx)
end if
end do
!call this%loadParamNames(trim(DataDir)//'boost.paramnames')
end subroutine boost_estimator_ReadCls
!We declare the function that calculates the likelihood a rela valued one
real(mcp) function boost_estimator_LnLike(this, CMB, Theory, DataParams) !the funtion takes parameters from othe modules
class(boost_estimatorLikelihood) :: this
class(CMBParams) CMB
class(TCosmoTheoryPredictions), target :: Theory
real(mcp) :: DataParams(:)
real(mcp) :: boost, theoryval1, theoryval2, noise1, noise2, observedf, coefficientA, coefficientB, theoryf
integer l
integer m
boost_estimator_LnLike = 0.d0
boost=DataParams(1)
associate(Cls=> Theory%Cls(1,1))
do l=0,this%lmaxtt ! loop over multipoles !skip C_0 and C_1
if (l >= 2 .and. this%ClsTT(l) > 0.) then
theoryval1 = CLs%CL(l)/(l*(l+ 1.d0))*2.d0*pi ! already multiplied by l*(l+1)/twopi
theoryval2 = CLs%CL(l + 1)/((l+1.d0)*(l+ 2.d0 ))*2.d0*pi ! already multiplied by l*(l+1)/twopi
noise1 = (14.d-6)**2*exp(l*(l+1.d0)*8.84)
noise2 = (14.d-6)**2*exp((l+1.d0)*(l+2.d0)*8.84)
observedf = this%ClsTTcross(l)
do m=-l,l ! loop over m !skip C_0 and C_1
coefficientA = l*sqrt((l+1.d0)**2-m**2)/sqrt(4*(l + 1.d0)**2-1.d0)
coefficientB = -(l+2)*sqrt((l+1.d0)**2-m**2)/sqrt(4*(l + 1.d0)**2-1.d0)
end do
theoryf = coefficientA*theoryval1 + coefficientB*theoryval2
boost_estimator_LnLike = boost_estimator_LnLike + &
((0.8)**2)*(observedf - (boost)*theoryf)**2 / ((theoryval1 + noise1)*(theoryval1 + noise1))
coefficientA = 0.d0
coefficientB = 0.d0
!print *, 'Power Spectra', boost
end if
end do
end associate
end function boost_estimator_LnLike
!ending the module
end module boost_estimator
________________________________________________________________________
But i still have a the MCMC stops and complains due to some negative Cls, there is something that i am missing? this is the error:
___________________________________________________________________________
Number of MPI processes: 2
file_root:boost_estimator_test
NOTE: num_massive_neutrinos ignored, using specified hierarchy
NOTE: num_massive_neutrinos ignored, using specified hierarchy
Random seeds: 18762, 22676 rand_inst: 2
Random seeds: 18662, 22676 rand_inst: 1
POWER /home/esteban/CosmoMC_test/data/homemade_like/power_cross_test.txt
POWER /home/esteban/CosmoMC_test/data/homemade_like/power_cross_test.txt
Doing non-linear Pk: F
Doing CMB lensing: T
Doing non-linear lensing: T
TT lmax = 2500
EE lmax = 2500
ET lmax = 2500
BB lmax = 2500
PP lmax = 2500
lmax_computed_cl = 2500
Computing tensors: F
max_eta_k = 14000.0000
transfer kmax = 1.00000000
adding parameters for: ECHM
Fast divided into 1 blocks
7 parameters ( 6 slow ( 0 semi-slow), 1 fast ( 0 semi-fast))
starting Monte-Carlo
Chain:0 drag accpt: 1.00000000 fast/slow 2.06666660 slow: 30
Chain:1 drag accpt: 1.00000000 fast/slow 2.20000005 slow: 30
Chain:0 drag accpt: 1.00000000 fast/slow 2.06666660 slow: 60
Chain:1 drag accpt: 1.00000000 fast/slow 2.23333335 slow: 60
Chain:0 drag accpt: 1.00000000 fast/slow 1.93333328 slow: 90
Chain:1 drag accpt: 1.00000000 fast/slow 2.31111121 slow: 90
Chain1, MPI done 'burn', Samples =118, like = 0.00000000
Time: 491.86910210100007 output lines= 116
slow changes 118 power changes 0
MPI_Min_Sample_Update 79 118
Chain2, MPI done 'burn', Samples =119, like = 0.00000000
Time: 500.03799697899990 output lines= 117
slow changes 119 power changes 0
MPI_Min_Sample_Update 79 119
Chain:0 drag accpt: 1.00000000 fast/slow 2.09999990 slow: 120
1 all_burn done
1 DoUpdates
Chain:1 drag accpt: 1.00000000 fast/slow 2.28333330 slow: 120
2 all_burn done
2 DoUpdates
Chain:0 drag accpt: 1.00000000 fast/slow 2.05333328 slow: 150
Chain:1 drag accpt: 1.00000000 fast/slow 2.40000010 slow: 150
Chain:0 drag accpt: 1.00000000 fast/slow 2.07777786 slow: 180
Chain:1 drag accpt: 0.994475126 fast/slow 2.48888898 slow: 180
Chain:1 drag accpt: 0.981308401 fast/slow 2.47619057 slow: 210
Chain:0 drag accpt: 1.00000000 fast/slow 2.00952387 slow: 210
Chain:1 drag accpt: 0.948616624 fast/slow 2.48333335 slow: 240
Chain:0 drag accpt: 1.00000000 fast/slow 1.98333335 slow: 240
Chain:1 drag accpt: 0.950704217 fast/slow 2.45185184 slow: 270
Chain:0 drag accpt: 1.00000000 fast/slow 2.04444456 slow: 270
Chain:1 drag accpt: 0.955413997 fast/slow 2.40000010 slow: 300
Chain:0 drag accpt: 1.00000000 fast/slow 2.09999990 slow: 300
Chain 2 MPI communicating
Chain 1 MPI communicating
Current convergence R-1 = 103.575783 chain steps = 322
slow changes 322 power changes 0
updating proposal density
Chain:1 drag accpt: 0.942857146 fast/slow 2.41212130 slow: 330
Chain:0 drag accpt: 0.985074639 fast/slow 2.10909081 slow: 330
Chain:1 drag accpt: 0.935064912 fast/slow 2.47777772 slow: 360
Chain:0 drag accpt: 0.957446814 fast/slow 2.22222233 slow: 360
Chain:1 drag accpt: 0.926365793 fast/slow 2.51282048 slow: 390
Chain:0 drag accpt: 0.955882370 fast/slow 2.26410246 slow: 390
Chain:1 drag accpt: 0.895522416 fast/slow 2.59761906 slow: 420
Chain:0 drag accpt: 0.952380955 fast/slow 2.31190467 slow: 420
Chain:1 drag accpt: 0.889328063 fast/slow 2.65555549 slow: 450
Chain:0 drag accpt: 0.955413997 fast/slow 2.39333344 slow: 450
Chain:1 drag accpt: 0.893854737 fast/slow 2.69374990 slow: 480
Chain:0 drag accpt: 0.958083808 fast/slow 2.46041656 slow: 480
Chain:1 drag accpt: 0.897887349 fast/slow 2.73333335 slow: 510
Chain:0 drag accpt: 0.956848025 fast/slow 2.48235297 slow: 510
Chain:1 drag accpt: 0.903010011 fast/slow 2.77407408 slow: 540
Chain:0 drag accpt: 0.955752194 fast/slow 2.54814816 slow: 540
Chain:1 drag accpt: 0.907643318 fast/slow 2.82456136 slow: 570
Chain:0 drag accpt: 0.957983196 fast/slow 2.58771920 slow: 570
Chain 2 MPI communicating
Chain 1 MPI communicating
Current convergence R-1 = 13.4338703 chain steps = 602
slow changes 577 power changes 0
updating proposal density
Chain:1 drag accpt: 0.906344414 fast/slow 2.83333325 slow: 600
Chain:0 drag accpt: 0.958466470 fast/slow 2.61333323 slow: 600
Chain:1 drag accpt: 0.907781005 fast/slow 2.85238099 slow: 630
Chain:0 drag accpt: 0.958904088 fast/slow 2.64603186 slow: 630
Chain:1 drag accpt: 0.901639342 fast/slow 2.88484859 slow: 660
Chain:0 drag accpt: 0.953757226 fast/slow 2.67121220 slow: 660
Calculator_CAMB: negative C_l (could edit to silent error here)
MpiStop: 1
--------------------------------------------------------------------------
MPI_ABORT was invoked on rank 1 in communicator MPI_COMM_WORLD
with errorcode 128.
NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.
___________________________________________________________________________
I am attaching the plot that i got using only my likelihood just before the MCMC stopped.