CosmicFish  Reference documentation for version 1.0
Looking into future Cosmology
008_Fisher_calculator_Cls.f90
Go to the documentation of this file.
1 !----------------------------------------------------------------------------------------
2 !
3 ! This file is part of CosmicFish.
4 !
5 ! Copyright (C) 2015-2016 by the CosmicFish authors
6 !
7 ! The CosmicFish code is free software;
8 ! You can use it, redistribute it, and/or modify it under the terms
9 ! of the GNU General Public License as published by the Free Software Foundation;
10 ! either version 3 of the License, or (at your option) any later version.
11 ! The full text of the license can be found in the file LICENSE at
12 ! the top level of the CosmicFish distribution.
13 !
14 !----------------------------------------------------------------------------------------
15 
18 
19 !----------------------------------------------------------------------------------------
21 
23 
24 module fisher_calculator_cls
25 
26  use precision
27  use modelparams
28  use cosmicfish_types
29  use cambmain
30  use camb
31  use constants, only: const_pi
32  use cosmicfish_camb_interface
33  use lensnoise
34  use fisher_calculator_cls_windows
35  use omp_lib
36  use cosmicfish_utilities
37  use matrix_utilities
38 
39 #ifdef COSMICFISH_EFTCAMB
40  use eftinitialization
41 #endif
42 
43  implicit none
44 
45  private
46 
47  public dimension_cl_covariance, cl_covariance_out, save_cl_covariance_to_file, &
48  & cls_derivative, add_noise_cls_to_fiducial, fisher_cls
49 
50 contains
51 
52  ! ---------------------------------------------------------------------------------------------
55  subroutine dimension_cl_covariance( P, FP, num )
56 
57  implicit none
58 
59  Type(cambparams) , intent(in) :: P
60  Type(cosmicfish_params), intent(in) :: FP
61  integer , intent(out) :: num
62 
63  integer :: i
64 
65  num = 0
66 
67  if ( fp%fisher_cls%Fisher_want_CMB_T ) num = num + 1
68  if ( fp%fisher_cls%Fisher_want_CMB_E ) num = num + 1
69  if ( fp%fisher_cls%Fisher_want_CMB_lensing ) num = num + 1
70 
71  do i = 1, num_redshiftwindows
72  if ( redshift_w(i)%kind == window_counts .and. fp%fisher_cls%Fisher_want_LSS_counts ) then
73  num = num + 1
74  else if ( redshift_w(i)%kind == window_lensing .and. fp%fisher_cls%Fisher_want_LSS_lensing ) then
75  num = num + 1
76  end if
77  end do
78 
79  if ( fp%fisher_cls%Fisher_want_CMB_B ) num = num + 1
80 
81  end subroutine dimension_cl_covariance
82 
83  ! ---------------------------------------------------------------------------------------------
86  subroutine cl_covariance_out( P, FP, cl_dim, l_min, l_max, cl_out, err )
87 
88  implicit none
89 
90  Type(cambparams) :: P
91  Type(cosmicfish_params) :: FP
92  integer, intent(in) :: cl_dim
93  integer, intent(in) :: l_min
94  integer, intent(in) :: l_max
95 
96  real(dl) :: cl_out(cl_dim, cl_dim, l_max-l_min+1)
97  integer, intent(out) :: err
103 
104  ! definitions:
105  logical :: EFTsuccess
106  real(dl) :: ell
107  real(dl), dimension(:,:), allocatable :: outarr
108  integer , dimension(:) , allocatable :: array_select
109  integer , dimension(cl_dim) :: l_maxes
110 
111  integer :: i, j, k, l, m, temp, stat, ind
112  real(dl) :: temp_factor
113 
114  ! initialization:
115  cl_out = 0.0_dl
116  err = 0
117 
118  ! set the window function:
119  if ( fp%fisher_cls%LSS_number_windows > 0 ) call init_camb_sources_windows( fp )
120 
121  ! set the camb parameters:
122  call cambparams_set(p)
123 
124 #ifdef COSMICFISH_EFTCAMB
125  ! check stability if EFTCAMB is used:
126  if (p%EFTflag /= 0) then
127  call eftcamb_initialization(eftsuccess)
128  if (.not. eftsuccess) then
129  write(*,*) 'EFTCAMB: unstable model'
130  err = 4
131  return
132  end if
133  end if
134 #endif
135 
136  ! call the CAMB_GetResults to compute the Cls.
137  call camb_getresults(p)
138 
139  ! test for errors in the computation
140  if (global_error_flag/=0) then
141  write (*,*) 'Error result '//trim(global_error_message)
142  err = 2
143  return
144  endif
145 
146  ! use the scal Cls array as a starting point for the Cls covariance:
147 
148  allocate(outarr(1:3+num_redshiftwindows,1:3+num_redshiftwindows))
149  allocate(array_select(3+num_redshiftwindows))
150 
151  ! decide which columns and rows of the scalar matrix to mantain:
152  array_select = 0
153 
154  do i=1, 3+num_redshiftwindows
155 
156  if ( fp%fisher_cls%Fisher_want_CMB_T ) array_select(1) = 1
157  if ( fp%fisher_cls%Fisher_want_CMB_E ) array_select(2) = 1
158  if ( fp%fisher_cls%Fisher_want_CMB_lensing ) array_select(3) = 1
159 
160  do j = 1, num_redshiftwindows
161  if ( redshift_w(j)%kind == window_counts .and. fp%fisher_cls%Fisher_want_LSS_counts ) then
162  array_select(3+j) = 1
163  else if ( redshift_w(j)%kind == window_lensing .and. fp%fisher_cls%Fisher_want_LSS_lensing ) then
164  array_select(3+j) = 1
165  end if
166  end do
167 
168  end do
169 
170  ! cycle over the scalar Cls matrix.
171  do i=l_min, l_max
172 
173  ell = REAL(i)
174  temp_factor = 2._dl*pi/(ell*(ell+1._dl))
175 
176  outarr = temp_factor*cl_scalar_array(i,1,1:3+num_redshiftwindows,1:3+num_redshiftwindows)
177  outarr(1:2,:)=sqrt(fp%output_factor)*outarr(1:2,:)
178  outarr(:,1:2)=sqrt(fp%output_factor)*outarr(:,1:2)
179 
180  l = 1
181  m = 1
182 
183  do j=1, 3+num_redshiftwindows
184  if ( array_select(j)==1 ) then
185  m = 1
186  do k=1, 3+num_redshiftwindows
187  if ( array_select(k)==1 ) then
188  cl_out(l, m, i-l_min+1) = outarr( j, k )
189  m = m + 1
190  end if
191  end do
192  l = l + 1
193  end if
194  end do
195 
196  end do
197 
198  ! now replace specific ones with lensed or add B mode power spectrum
199  if ( p%WantTensors .and. .not. p%DoLensing) then
200  ! if want tensor not lensing out = total
201  do i = l_min, min( cp%Max_l_tensor, l_max )
202  ell = REAL(i)
203  temp_factor = 2*const_pi/(ell*(ell+1))
204  j = i-l_min+1
205 
206  ! TT
207  if ( fp%fisher_cls%Fisher_want_CMB_T ) cl_out(1,1,j) = temp_factor*fp%output_factor*(cl_scalar(i, 1, c_temp)+ cl_tensor(i, 1, c_temp))
208  ! EE and TE
209  if ( fp%fisher_cls%Fisher_want_CMB_T .and. fp%fisher_cls%Fisher_want_CMB_E ) then
210  cl_out(2,2,j) = temp_factor*fp%output_factor*(cl_scalar(i, 1, c_e)+ cl_tensor(i, 1, c_e))
211  cl_out(1,2,j) = temp_factor*fp%output_factor*(cl_scalar(i, 1, c_cross) + cl_tensor(i, 1, ct_cross))
212  end if
213  ! if only EE
214  if ( .not. fp%fisher_cls%Fisher_want_CMB_T .and. fp%fisher_cls%Fisher_want_CMB_E ) then
215  cl_out(1,1,j) = temp_factor*fp%output_factor*(cl_scalar(i, 1, c_e)+ cl_tensor(i, 1, c_e))
216  end if
217 
218  ! BB
219  if ( fp%fisher_cls%Fisher_want_CMB_B ) cl_out( cl_dim, cl_dim, j ) = temp_factor*fp%output_factor*cl_tensor(i, 1, ct_b)
220 
221  end do
222  do i = cp%Max_l_tensor+1,l_max
223  ell = REAL(i)
224  temp_factor = 2*const_pi/(ell*(ell+1))
225  j = i-l_min+1
226 
227  ! TT
228  if ( fp%fisher_cls%Fisher_want_CMB_T ) cl_out(1,1,j) = temp_factor*fp%output_factor*cl_scalar(i, 1, c_temp)
229  ! EE and TE
230  if ( fp%fisher_cls%Fisher_want_CMB_T .and. fp%fisher_cls%Fisher_want_CMB_E ) then
231  cl_out(2,2,j) = temp_factor*fp%output_factor*cl_scalar(i, 1, c_e)
232  cl_out(1,2,j) = temp_factor*fp%output_factor*cl_scalar(i, 1, c_cross)
233  end if
234  ! if only EE
235  if ( .not. fp%fisher_cls%Fisher_want_CMB_T .and. fp%fisher_cls%Fisher_want_CMB_E ) then
236  cl_out(1,1,j) = temp_factor*fp%output_factor*cl_scalar(i, 1, c_e)
237  end if
238 
239  ! BB
240  if ( fp%fisher_cls%Fisher_want_CMB_B ) cl_out( cl_dim, cl_dim, j ) = 0._dl
241 
242  end do
243  else if (.not. p%WantTensors .and. p%DoLensing) then
244  ! if want lensing and not tensor out = lensed
245  do i = l_min, min( lmax_lensed, l_max )
246  ell = REAL(i)
247  temp_factor = 2*const_pi/(ell*(ell+1))
248  j = i-l_min+1
249 
250  ! TT
251  if ( fp%fisher_cls%Fisher_want_CMB_T ) cl_out(1,1,j) = temp_factor*fp%output_factor*cl_lensed(i, 1, ct_temp)
252  ! EE and TE
253  if ( fp%fisher_cls%Fisher_want_CMB_T .and. fp%fisher_cls%Fisher_want_CMB_E ) then
254  cl_out(2,2,j) = temp_factor*fp%output_factor*cl_lensed(i, 1, ct_e)
255  cl_out(1,2,j) = temp_factor*fp%output_factor*cl_lensed(i, 1, ct_cross)
256  end if
257  ! if only EE
258  if ( .not. fp%fisher_cls%Fisher_want_CMB_T .and. fp%fisher_cls%Fisher_want_CMB_E ) then
259  cl_out(1,1,j) = temp_factor*fp%output_factor*cl_lensed(i, 1, ct_e)
260  end if
261 
262  ! BB
263  if ( fp%fisher_cls%Fisher_want_CMB_B ) cl_out( cl_dim, cl_dim, j ) = temp_factor*fp%output_factor*cl_lensed(i, 1, ct_b)
264  end do
265  else if (p%WantTensors .and. p%DoLensing) then
266  ! if want lensing and want tensor out = total lensed
267  do i = l_min, min( cp%Max_l_tensor, lmax_lensed, l_max )
268  ell = REAL(i)
269  temp_factor = 2*const_pi/(ell*(ell+1))
270  j = i-l_min+1
271 
272  ! TT
273  if ( fp%fisher_cls%Fisher_want_CMB_T ) cl_out(1,1,j) = temp_factor*fp%output_factor*(cl_lensed(i, 1, ct_temp) +cl_tensor(i, 1, ct_temp))
274  ! EE and TE
275  if ( fp%fisher_cls%Fisher_want_CMB_T .and. fp%fisher_cls%Fisher_want_CMB_E ) then
276  cl_out(2,2,j) = temp_factor*fp%output_factor*(cl_lensed(i, 1, ct_e) +cl_tensor(i, 1, ct_e))
277  cl_out(1,2,j) = temp_factor*fp%output_factor*(cl_lensed(i, 1, ct_cross) +cl_tensor(i, 1, ct_cross))
278  end if
279  ! if only EE
280  if ( .not. fp%fisher_cls%Fisher_want_CMB_T .and. fp%fisher_cls%Fisher_want_CMB_E ) then
281  cl_out(1,1,j) = temp_factor*fp%output_factor*(cl_lensed(i, 1, ct_e) +cl_tensor(i, 1, ct_e))
282  end if
283 
284  ! BB
285  if ( fp%fisher_cls%Fisher_want_CMB_B ) cl_out( cl_dim, cl_dim, j ) = temp_factor*fp%output_factor*(cl_lensed(i, 1, ct_b) +cl_tensor(i, 1, ct_b))
286  end do
287  do i = min( cp%Max_l_tensor, lmax_lensed, l_max )+1, min( lmax_lensed, l_max )
288  ell = REAL(i)
289  temp_factor = 2*const_pi/(ell*(ell+1))
290  j = i-l_min+1
291 
292  ! TT
293  if ( fp%fisher_cls%Fisher_want_CMB_T ) cl_out(1,1,j) = temp_factor*fp%output_factor*cl_lensed(i, 1, ct_temp)
294  ! EE and TE
295  if ( fp%fisher_cls%Fisher_want_CMB_T .and. fp%fisher_cls%Fisher_want_CMB_E ) then
296  cl_out(2,2,j) = temp_factor*fp%output_factor*cl_lensed(i, 1, ct_e)
297  cl_out(1,2,j) = temp_factor*fp%output_factor*cl_lensed(i, 1, ct_cross)
298  end if
299  ! if only EE
300  if ( .not. fp%fisher_cls%Fisher_want_CMB_T .and. fp%fisher_cls%Fisher_want_CMB_E ) then
301  cl_out(1,1,j) = temp_factor*fp%output_factor*cl_lensed(i, 1, ct_e)
302  end if
303 
304  ! BB
305  if ( fp%fisher_cls%Fisher_want_CMB_B ) cl_out( cl_dim, cl_dim, j ) = temp_factor*fp%output_factor*cl_lensed(i, 1, ct_b)
306  end do
307  end if
308 
309  ! apply the l_max cut:
310 
311  l_maxes = 1
312  ind = 1
313 
314  if ( fp%fisher_cls%Fisher_want_CMB_T ) then
315  l_maxes(ind) = fp%fisher_cls%l_max_TT
316  ind = ind + 1
317  end if
318 
319  if ( fp%fisher_cls%Fisher_want_CMB_E ) then
320  l_maxes(ind) = fp%fisher_cls%l_max_EE
321  ind = ind + 1
322  end if
323 
324  if ( fp%fisher_cls%Fisher_want_CMB_lensing ) then
325  l_maxes(ind) = min( fp%fisher_cls%l_max_TT, &
326  & fp%fisher_cls%l_max_EE, &
327  & fp%fisher_cls%l_max_BB )
328  ind = ind + 1
329  end if
330 
331  if ( fp%fisher_cls%Fisher_want_CMB_B ) then
332  l_maxes(cl_dim) = fp%fisher_cls%l_max_BB
333  end if
334 
335  do k = 1, fp%fisher_cls%LSS_number_windows
336  l_maxes(ind) = fp%fisher_cls%LSS_lmax(k)
337  ind = ind + 1
338  end do
339 
340  ! apply the l cut:
341  do i=l_min, l_max
342  do j=1, cl_dim
343  do k=1, j
344 
345  if ( i <= min(l_maxes(j), l_maxes(k)) ) then
346  cl_out( k, j, i-l_min+1) = cl_out( k, j, i-l_min+1)
347  else
348  cl_out( k, j, i-l_min+1) = 0._dl
349  end if
350 
351  end do
352  end do
353  end do
354 
355  ! ensure the fact that the covariance is symmetric:
356  do i=l_min, l_max
357  do j=1, cl_dim
358  do k=1, j
359 
360  cl_out( j, k, i-l_min+1) = cl_out( k, j, i-l_min+1)
361 
362  end do
363  end do
364  end do
365 
366  ! remove cross correlation if it is not wanted:
367  if ( .not. fp%fisher_cls%Fisher_want_XC ) then
368  do i=l_min, l_max
369  do j=1, cl_dim
370  do k=1, cl_dim
371 
372  if ( j/=k ) then
373  cl_out( j, k, i-l_min+1) = 0._dl
374  end if
375 
376  end do
377  end do
378  end do
379  end if
380 
381  end subroutine cl_covariance_out
382 
383  ! ---------------------------------------------------------------------------------------------
385  subroutine save_cl_covariance_to_file( cl_cov_in, cl_dim, l_min, l_max, filename )
386 
387  implicit none
388 
389  integer , intent(in) :: cl_dim
390  integer , intent(in) :: l_min
391  integer , intent(in) :: l_max
392  real(dl), dimension(cl_dim, cl_dim, l_max-l_min+1), intent(in) :: cl_cov_in
393  character(len=*), intent(in) :: filename
394 
395  integer :: i, j, k
396  real(dl) :: factor
397  real(dl), dimension(cl_dim, cl_dim) :: dummy
398 
399  open(unit=666, file=filename, action="write", status="replace")
400  ! cycle over l
401  do i = 1, l_max-l_min+1
402  ! change this factor if want to compensate for factors l()...
403  factor = 1._dl
404  dummy(:,:) = cl_cov_in(:,:,i)
405  write(666,'(E15.5)',advance='no') REAL(i+l_min-1)
406  ! cycle over the Cl cov matrix at a fixed l
407  do j = 1, cl_dim
408  do k = 1, cl_dim
409  write(666,'(ES20.5E3)',advance='no') factor*dummy(j,k)
410  end do
411  end do
412  ! go to newline
413  write(666,'(a)') ' '
414  end do
415  close(666)
416 
417  end subroutine
418 
419  ! ---------------------------------------------------------------------------------------------
422  subroutine cls_derivative( P, FP, param_num, &
423  & cl_initial, cl_derivative, cl_error, &
424  & initial_step, &
425  & out_dim, l_min, l_max, num_param, &
426  & error_code )
427 
428  implicit none
429 
430  Type(cambparams) , intent(in) :: P
431  Type(cosmicfish_params) , intent(in) :: FP
432 
433  integer , intent(in) :: param_num
434  integer , intent(in) :: out_dim
435  integer , intent(in) :: l_min
436  integer , intent(in) :: l_max
437  integer , intent(in) :: num_param
438 
439  real(dl), dimension(out_dim, out_dim, l_max-l_min+1), intent(in) :: cl_initial
440  real(dl), dimension(out_dim, out_dim, l_max-l_min+1), intent(out) :: cl_derivative
441  real(dl), dimension(out_dim, out_dim, l_max-l_min+1), intent(out) :: cl_error
442 
443  real(dl), intent(in) :: initial_step
444  integer , intent(out) :: error_code
449 
450  ! definitions:
451  real(dl), dimension(:,:,:), allocatable :: cl_temp_1, cl_temp_2, errt
452  integer , dimension(:,:,:), allocatable :: accept_mask
453  real(dl), dimension(:,:,:,:,:), allocatable :: a
454  real(dl), dimension(num_param) :: param_vector, fiducial_param_vector, temp_param_vector
455 
456  ! parameters of the algorithm:
457  real(dl), parameter :: con = 1.4_dl ! parameter that decides the decrease in the stepsize h
458  real(dl), parameter :: big = 1.d+30 ! a big number used to fill the error array initially
459  real(dl), parameter :: safe = 2._dl ! parameter used to define the stopping cryterium
460  integer , parameter :: ntab = 30 ! maximum number of Cls computation
461  real(dl) :: con2 = con*con
462  integer :: i, j, k, l, m, side, err, AllocateStatus
463  real(dl) :: fac, hh
464  real(dl) :: temp1, temp2, temp3
465  Type(cambparams) :: P_temp
466  Type(cosmicfish_params) :: FP_temp
467  integer :: accepted
468  real(dl) :: acceptance
469 
470  ! initialize output variables:
471  error_code = 0
472  cl_derivative = 0._dl
473  cl_error = 0._dl
474 
475  ! check input values are correct
476  if ( initial_step == 0._dl ) then
477  stop 'cls_derivative initial stepsize is zero'
478  end if
479 
480  ! allocate the arrays:
481 
482  if ( fp%adaptivity ) then
483  allocate( cl_temp_1(out_dim, out_dim, l_max-l_min+1), &
484  & cl_temp_2(out_dim, out_dim, l_max-l_min+1), &
485  & errt(out_dim, out_dim, l_max-l_min+1), &
486  & accept_mask(out_dim, out_dim, l_max-l_min+1),&
487  & stat = allocatestatus )
488  if (allocatestatus /= 0) stop "cls_derivative: Allocation failed. Not enough memory"
489  allocate( a(ntab, ntab, out_dim, out_dim, l_max-l_min+1) , stat = allocatestatus )
490  if (allocatestatus /= 0) stop "cls_derivative: Allocation failed. Not enough memory for the derivative array"
491  else
492  allocate( cl_temp_1(out_dim, out_dim, l_max-l_min+1), &
493  & cl_temp_2(out_dim, out_dim, l_max-l_min+1), &
494  & stat = allocatestatus )
495  if (allocatestatus /= 0) stop "cls_derivative: Allocation failed. Not enough memory"
496  allocate( a(1,1, out_dim, out_dim, l_max-l_min+1) , stat = allocatestatus )
497  if (allocatestatus /= 0) stop "cls_derivative: Allocation failed. Not enough memory for the derivative array"
498  end if
499 
500  ! initialize the parameter array:
501  p_temp = p
502  fp_temp = fp
503 
504  call camb_params_to_params_array( p_temp, fp, num_param, param_vector )
505 
506  fiducial_param_vector = param_vector
507  temp_param_vector = param_vector
508 
509  ! initialize the computation:
510  hh = initial_step
511  cl_error = big
512  side = 0 ! both sides
513  if ( fp%adaptivity ) accept_mask = 0
514 
515  ! do the initial step on the right:
516  param_vector(param_num) = fiducial_param_vector(param_num) + hh
517  call params_array_to_camb_param( p_temp, fp_temp, num_param, param_vector )
518  call cl_covariance_out( p_temp, fp_temp, out_dim, l_min, l_max, cl_temp_1, err )
519  ! if on the right we cannot go we go to the left:
520  if ( err /= 0 ) side = 1
521  ! do the initial step on the left:
522  temp_param_vector(param_num) = fiducial_param_vector(param_num) - hh
523  call params_array_to_camb_param( p_temp, fp_temp, num_param, temp_param_vector )
524  call cl_covariance_out( p_temp, fp_temp, out_dim, l_min, l_max, cl_temp_2, err )
525 
526  ! if on the left we cannot go decide what to do:
527  if ( err /= 0 ) then
528  if ( side == 1 ) then
529  ! we cannot go to the left nor to the right. Return derivative zero with error code.
530  error_code = 4
531  cl_error = 0.0_dl
532  cl_derivative = 0.0_dl
533  return
534  else if ( side == 0 ) then
535  ! we cannot go to the left but we can go to the right:
536  side = 2
537  end if
538  end if
539 
540  ! store the results based on the side we decided to go:
541  if (side == 0) then ! both sides
542  a(1,1,:,:,:) = (cl_temp_1(:,:,:) - cl_temp_2(:,:,:))/(2.0_dl*hh)
543  else if (side == 1) then ! left side: increment negative
544  a(1,1,:,:,:) = (cl_temp_2(:,:,:) - cl_initial(:,:,:))/(-hh)
545  else if (side == 2) then ! right side: increment positive
546  a(1,1,:,:,:) = (cl_temp_1(:,:,:) - cl_initial(:,:,:))/(hh)
547  end if
548 
549  ! fixed stepsize derivative if adaptivity is not wanted:
550  if ( .not. fp%adaptivity ) then
551  cl_derivative(:,:,:) = a(1,1,:,:,:)
552  cl_error = 0.0_dl
553  deallocate( cl_temp_1, cl_temp_2, a )
554  return
555  end if
556 
557  ! start the derivative computation:
558  do i = 2, ntab
559 
560  ! feedback for the adaptive derivative:
561  if ( fp%cosmicfish_feedback >= 2 ) then
562  accepted = sum(accept_mask)
563  acceptance = REAL(accepted)/REAL(out_dim*out_dim*(l_max-l_min+1))
564  print*, 'iteration:', i
565  print*, ' stepsize:', hh
566  print*, ' accepted elements', accepted, 'accepted ratio', acceptance
567  end if
568 
569  ! decrease the stepsize:
570  hh = hh/con
571 
572  ! compute the finite difference:
573  if (side == 0) then
574  param_vector(param_num) = fiducial_param_vector(param_num) + hh
575  temp_param_vector(param_num) = fiducial_param_vector(param_num) - hh
576  ! compute for + hh
577  call params_array_to_camb_param( p_temp, fp_temp, num_param, param_vector )
578  call cl_covariance_out( p_temp, fp_temp, out_dim, l_min, l_max, cl_temp_1, err )
579  if ( err /= 0 ) then
580  error_code = 4
581  cl_error = 0.0_dl
582  cl_derivative = 0.0_dl
583  return
584  end if
585  ! compute for - hh
586  call params_array_to_camb_param( p_temp, fp_temp, num_param, temp_param_vector )
587  call cl_covariance_out( p_temp, fp_temp, out_dim, l_min, l_max, cl_temp_2, err )
588  if ( err /= 0 ) then
589  error_code = 4
590  cl_error = 0.0_dl
591  cl_derivative = 0.0_dl
592  return
593  end if
594  ! store the result
595  a(1,i,:,:,:) = (cl_temp_1(:,:,:) - cl_temp_2(:,:,:))/(2.0_dl*hh)
596  else if (side==1) then
597  param_vector(param_num) = fiducial_param_vector(param_num) - hh
598  call params_array_to_camb_param( p_temp, fp_temp, num_param, param_vector )
599  call cl_covariance_out( p_temp, fp_temp, out_dim, l_min, l_max, cl_temp_1, err )
600  if ( err /= 0 ) then
601  error_code = 4
602  cl_error = 0.0_dl
603  cl_derivative = 0.0_dl
604  return
605  end if
606  ! store the result
607  a(1,i,:,:,:) = (cl_temp_1(:,:,:) - cl_initial(:,:,:))/(-hh)
608  else if (side==2) then
609  param_vector(param_num) = fiducial_param_vector(param_num) + hh
610  call params_array_to_camb_param( p_temp, fp_temp, num_param, param_vector )
611  call cl_covariance_out( p_temp, fp_temp, out_dim, l_min, l_max, cl_temp_1, err )
612  if ( err /= 0 ) then
613  error_code = 4
614  cl_error = 0.0_dl
615  cl_derivative = 0.0_dl
616  return
617  end if
618  ! store the result
619  a(1,i,:,:,:) = (cl_temp_1(:,:,:) - cl_initial(:,:,:))/(hh)
620  end if
621 
622  ! now do the extrapolation table:
623  fac=con2
624  do j = 2, i
625 
626  forall ( k = 1:out_dim, l = 1:out_dim, m = 1:l_max-l_min+1, &
627  & accept_mask(k,l,m) == 0 )
628  ! compute the interpolation table:
629  a(j,i,k,l,m) = (a(j-1,i,k,l,m)*fac-a(j-1,i-1,k,l,m))/(fac-1._dl)
630  ! estimate the error:
631  errt(k,l,m) = max( abs(a(j,i,k,l,m)-a(j-1,i,k,l,m)), &
632  & abs(a(j,i,k,l,m)-a(j-1,i-1,k,l,m)))
633  end forall
634 
635  fac=con2*fac
636 
637  ! if there is an improvement save the result:
638 
639  forall ( k = 1:out_dim, l = 1:out_dim, m = 1:l_max-l_min+1, &
640  & errt(k,l,m) .le. cl_error(k,l,m) .and. accept_mask(k,l,m) == 0 )
641 
642  cl_error(k,l,m) = errt(k,l,m)
643  cl_derivative(k,l,m) = a(j,i,k,l,m)
644  end forall
645 
646  end do
647 
648  ! accept a value of the derivative if the value of the extrapolation is stable
649  ! Store this information into the mask matrix
650  ! so that the value contained in cl_error and cl_derivative is not overwritten.
651  forall ( k = 1:out_dim, l = 1:out_dim, m = 1:l_max-l_min+1, &
652  & accept_mask(k,l,m) == 0 .and. &
653  & (a(i,i,k,l,m) - a(i-1,i-1,k,l,m)) .ge. safe*cl_error(k,l,m) )
654 
655  accept_mask(k,l,m) = 1
656  end forall
657 
658  ! if all the entries match the required accuracy exit the loop
659  if ( sum(accept_mask) == out_dim*out_dim*(l_max-l_min+1) ) exit
660 
661  end do
662 
663  ! quality check on the results:
664  if ( sum(accept_mask) /= out_dim*out_dim*(l_max-l_min+1) ) then
665  error_code = 3
666  end if
667 
668  ! deallocate the arrays:
669  if ( fp%adaptivity ) then
670  deallocate( cl_temp_1, cl_temp_2, errt, accept_mask, a )
671  end if
672 
673  return
674 
675  end subroutine cls_derivative
676 
677  ! ---------------------------------------------------------------------------------------------
680  subroutine add_noise_cls_to_fiducial( P, FP, cl_fiducial, out_dim, l_min, l_max, error_code)
681 
682  implicit none
683 
684  Type(cambparams) , intent(in) :: P
685  Type(cosmicfish_params) , intent(in) :: FP
686 
687  integer , intent(in) :: out_dim
688  integer , intent(in) :: l_min
689  integer , intent(in) :: l_max
690 
691  real(dl), dimension(out_dim, out_dim, l_max-l_min+1), intent(inout) :: cl_fiducial
692 
693  integer , intent(out) :: error_code
698 
699  ! definitions:
700 
701  integer :: jjmaxTmax = 200
702 
703  integer :: ind, i, j, k
704  real(dl) :: temp1, temp2, ell, factor
705  real(dl) :: Compute_CMB_TT_Noise, Compute_CMB_EE_Noise, Compute_CMB_BB_Noise
706  real(dl), allocatable :: NoiseVar(:), NoiseVarPol(:), sigma2(:)
707  real(dl), allocatable :: incls(:,:), inlcls(:,:), clnoise(:,:), nldd(:)
708  real(dl), allocatable :: fskyes(:)
709  logical :: EFTsuccess
710 
711  ! initialization:
712  error_code = 0
713  ind = 1
714 
715  ! initialize CMB stuff:
716  if ( fp%fisher_cls%Fisher_want_CMB_T .or. fp%fisher_cls%Fisher_want_CMB_E .or. &
717  & fp%fisher_cls%Fisher_want_CMB_B .or. fp%fisher_cls%Fisher_want_CMB_lensing ) then
718 
719  allocate( noisevar(fp%fisher_cls%CMB_n_channels), &
720  & noisevarpol(fp%fisher_cls%CMB_n_channels), &
721  & sigma2(fp%fisher_cls%CMB_n_channels) )
722 
723  temp1 = p%Tcmb*const_pi/180._dl/60._dl
724  temp2 = 180._dl*60._dl*sqrt(8._dl*log(2._dl))/const_pi
725 
726  forall( i = 1:fp%fisher_cls%CMB_n_channels )
727  noisevar(i) = ( fp%fisher_cls%CMB_temp_sens(i)*fp%fisher_cls%CMB_fwhm(i)*temp1 )**2
728  noisevarpol(i) = ( fp%fisher_cls%CMB_pol_sens(i)*fp%fisher_cls%CMB_fwhm(i)*temp1 )**2
729  sigma2(i) = -1._dl*( fp%fisher_cls%CMB_fwhm(i)/temp2 )**2
730  end forall
731 
732  end if
733 
734  ! add TT noise:
735  if ( fp%fisher_cls%Fisher_want_CMB_T ) then
736 
737  do j = l_min, l_max
738  ! compute noise:
739  ell = REAL(j)
740  compute_cmb_tt_noise = 0._dl
741  do i=1, fp%fisher_cls%CMB_n_channels
742  compute_cmb_tt_noise = compute_cmb_tt_noise &
743  & + 1._dl/noisevar(i)*exp(ell*(ell+1)*sigma2(i) )
744  end do
745  compute_cmb_tt_noise = 1._dl/compute_cmb_tt_noise
746  ! add noise to fiducial:
747  k = j - l_min + 1
748  cl_fiducial( ind, ind, k ) = cl_fiducial( ind, ind, k ) + compute_cmb_tt_noise
749  end do
750  ind = ind + 1
751 
752  end if
753  ! add EE noise:
754  if ( fp%fisher_cls%Fisher_want_CMB_E ) then
755 
756  do j = l_min, l_max
757  ! compute noise:
758  ell = REAL(j)
759  compute_cmb_ee_noise = 0._dl
760  do i=1, fp%fisher_cls%CMB_n_channels
761  compute_cmb_ee_noise = compute_cmb_ee_noise &
762  & + 1._dl/noisevarpol(i)*exp(ell*(ell+1)*sigma2(i) )
763  end do
764  compute_cmb_ee_noise = 1._dl/compute_cmb_ee_noise
765  ! add noise to fiducial:
766  k = j - l_min + 1
767  cl_fiducial( ind, ind, k ) = cl_fiducial( ind, ind, k ) + compute_cmb_ee_noise
768  end do
769  ind = ind + 1
770 
771  end if
772  ! add lensing noise:
773  if ( fp%fisher_cls%Fisher_want_CMB_lensing ) then
774 
775  ! we need scalar and lensed scalar Cls. We have to call camb and we want the model to be stable
776 #ifdef COSMICFISH_EFTCAMB
777  ! check stability if EFTCAMB is used:
778  if (p%EFTflag /= 0) then
779  call cambparams_set(p)
780  call eftcamb_initialization(eftsuccess)
781  if (.not. eftsuccess) then
782  write(*,*) 'EFTCAMB: unstable fiducial model. Calculation cannot proceed.'
783  stop
784  end if
785  end if
786 #endif
787  call camb_getresults(p)
788 
789  ! test for errors in the computation:
790  if (global_error_flag/=0) then
791  write (*,*) 'Error result '//trim(global_error_message)
792  write (*,*) 'add_noise_cls_to_fiducial: calculation of fiducial model failed. Calculation cannot proceed.'
793  stop
794  endif
795 
796  allocate(incls(5, 1:l_max))
797  allocate(inlcls(4, 1:l_max))
798  allocate(clnoise(2, 1:l_max))
799  allocate(nldd(1:l_max))
800  ! initialization of the Cls matrices:
801  incls = 0._dl
802  inlcls = 0._dl
803  clnoise = 0._dl
804  nldd = 0._dl
805 
806  ! write the Cls in the proper places:
807  do i = l_min, l_max
808 
809  ell = REAL(i)
810  factor = 2*const_pi/(ell*(ell+1))
811  ! write the cls in the format that FUTURECMB likes:
812  incls(1,i) = fp%output_factor*factor*cl_scalar_array(i,1,1,1) ! TT
813  incls(2,i) = fp%output_factor*factor*cl_scalar_array(i,1,2,2) ! EE
814  incls(3,i) = fp%output_factor*factor*cl_scalar_array(i,1,1,2) ! TE
815  incls(4,i) = 2*const_pi*cl_scalar_array(i,1,3,3) ! dd
816  incls(5,i) = sqrt(fp%output_factor)*sqrt(ell*(ell+1))*factor*cl_scalar_array(i,1,1,3) ! Td
817 
818  inlcls(1,i) = fp%output_factor*factor*cl_lensed(i, 1, ct_temp) ! TT ! micro K^2
819  inlcls(2,i) = fp%output_factor*factor*cl_lensed(i, 1, ct_e) ! EE ! micro K^2
820  inlcls(3,i) = fp%output_factor*factor*cl_lensed(i, 1, ct_cross) ! TE ! micro K^2
821  inlcls(4,i) = fp%output_factor*factor*cl_lensed(i, 1, ct_b) ! BB ! micro K^2
822 
823  ell = REAL(i)
824  compute_cmb_tt_noise = 0._dl
825  do j=1, fp%fisher_cls%CMB_n_channels
826  compute_cmb_tt_noise = compute_cmb_tt_noise &
827  & + 1._dl/noisevar(j)*exp(ell*(ell+1)*sigma2(j) )
828  end do
829  compute_cmb_tt_noise = 1._dl/compute_cmb_tt_noise
830  clnoise(1,i) = compute_cmb_tt_noise
831 
832  compute_cmb_ee_noise = 0._dl
833  do j=1, fp%fisher_cls%CMB_n_channels
834  compute_cmb_ee_noise = compute_cmb_ee_noise &
835  & + 1._dl/noisevarpol(j)*exp(ell*(ell+1)*sigma2(j) )
836  end do
837  compute_cmb_ee_noise = 1._dl/compute_cmb_ee_noise
838  clnoise(2,i) = compute_cmb_ee_noise
839 
840  end do
841 
842  ! compute the lensing noise with FUTURCMB:
843  call calc_nldd(incls, inlcls, clnoise, nldd, jjmaxtmax, l_max)
844 
845  ! add noise to the fiducial:
846  do i = l_min, l_max
847  ell = REAL(i)
848  factor = 1._dl/(ell*(ell+1))
849  k = i - l_min + 1
850  cl_fiducial( ind, ind, k ) = cl_fiducial( ind, ind, k ) + factor*nldd(i)
851  end do
852 
853  ind = ind + 1
854 
855  end if
856  ! add BB noise:
857  if ( fp%fisher_cls%Fisher_want_CMB_b ) then
858 
859  do j = l_min, l_max
860  ! compute noise:
861  ell = REAL(j)
862  compute_cmb_bb_noise = 0._dl
863  do i=1, fp%fisher_cls%CMB_n_channels
864  compute_cmb_bb_noise = compute_cmb_bb_noise &
865  & + 1._dl/noisevarpol(i)*exp(ell*(ell+1)*sigma2(i) )
866  end do
867  compute_cmb_bb_noise = 1._dl/compute_cmb_bb_noise
868  ! add noise to fiducial:
869  k = j - l_min + 1
870  cl_fiducial( out_dim, out_dim, k ) = cl_fiducial( out_dim, out_dim, k ) + compute_cmb_bb_noise
871  end do
872 
873  end if
874 
875  ! add counts and lensing windows noises:
876  j = 1
877  do k = 1, num_redshiftwindows
878  ! counts noise
879  if ( redshift_w(k)%kind == window_counts .and. fp%fisher_cls%Fisher_want_LSS_counts ) then
880  do i = 1, l_max-l_min+1
881  cl_fiducial( ind, ind, i ) = &
882  & cl_fiducial( ind, ind, i ) + 1._dl/fp%fisher_cls%LSS_num_galaxies(j)
883  end do
884  ind = ind + 1
885  j = j + 1
886  ! lensing noise
887  else if ( redshift_w(k)%kind == window_lensing .and. fp%fisher_cls%Fisher_want_LSS_lensing ) then
888  do i = 1, l_max-l_min+1
889  cl_fiducial( ind, ind, i ) = &
890  & cl_fiducial( ind, ind, i ) &
891  & +fp%fisher_cls%LSS_intrinsic_ellipticity(j)**2/fp%fisher_cls%LSS_num_galaxies(j)
892  end do
893  ind = ind + 1
894  j = j + 1
895  end if
896  end do
897 
898  ! multiply by f_sky:
899 
900  ! prepare an f_sky vector:
901  allocate( fskyes(out_dim) )
902  fskyes = 1._dl
903 
904  ind = 1
905 
906  if ( fp%fisher_cls%Fisher_want_CMB_T ) then
907  fskyes(ind) = fp%fisher_cls%CMB_TT_fsky
908  ind = ind + 1
909  end if
910 
911  if ( fp%fisher_cls%Fisher_want_CMB_E ) then
912  fskyes(ind) = fp%fisher_cls%CMB_EE_fsky
913  ind = ind + 1
914  end if
915 
916  if ( fp%fisher_cls%Fisher_want_CMB_lensing ) then
917  fskyes(ind) = min( fp%fisher_cls%CMB_TT_fsky, &
918  & fp%fisher_cls%CMB_EE_fsky, &
919  & fp%fisher_cls%CMB_BB_fsky )
920  ind = ind + 1
921  end if
922 
923  if ( fp%fisher_cls%Fisher_want_CMB_B ) then
924  fskyes(out_dim) = fp%fisher_cls%CMB_BB_fsky
925  end if
926 
927  do k = 1, fp%fisher_cls%LSS_number_windows
928  fskyes(ind) = fp%fisher_cls%LSS_fsky(k)
929  ind = ind + 1
930  end do
931 
932  ! multiply in the appropriate manner the fsky vector to the Cls matrix:
933  do j=1, out_dim
934  do k=1, out_dim
935 
936  cl_fiducial(j, k, :) = cl_fiducial(j, k, :)/sqrt(sqrt(fskyes(j)*fskyes(k)))
937 
938  end do
939  end do
940 
941  end subroutine add_noise_cls_to_fiducial
942 
943  ! ---------------------------------------------------------------------------------------------
945  subroutine fisher_cls( P, FP, num_param, Fisher_Matrix, outroot )
946 
947  implicit none
948 
949  Type(cambparams) , intent(in) :: P
950  Type(cosmicfish_params) , intent(in) :: FP
951 
952  integer , intent(in) :: num_param
953  real(dl), dimension(num_param,num_param), intent(out) :: Fisher_Matrix
954  character(len=*), intent(in), optional :: outroot
955 
956  real(dl) :: time_1, time_2
957  real(dl) :: initial_step, temp, temp1, temp2
958  integer :: cl_dim, l_min, l_max, err, AllocateStatus, ind, ind2
959  integer :: i, j, k, l
960  real(dl), allocatable :: param_array(:)
961  real(dl), dimension(:,:,:), allocatable :: cl_fiducial, cl_derivative, cl_temp
962  real(dl), dimension(:,:,:,:), allocatable :: cl_derivative_array
963  real(dl), dimension(:,:) , allocatable :: dummy_matrix_1, dummy_matrix_2
964  real(dl), dimension(:,:) , allocatable :: dummy_matrix_3, dummy_matrix_4, dummy_matrix_5
965 
966  ! allocate a parameter array:
967  allocate( param_array(num_param) )
968  ! write the parameters in the parameter array:
969  call camb_params_to_params_array( p, fp, num_param, param_array )
970  ! get the size of the cls covariance matrix:
971  call dimension_cl_covariance( p, fp, cl_dim )
972  ! protect against errors:
973  if ( cl_dim == 0 ) then
974  write(*,*) 'WARNING: the Cls matrix is of zero dimension'
975  fisher_matrix = 0._dl
976  return
977  end if
978 
979  ! decide l_min and l_max:
980  l_min = 2
981  l_max = maxval( fp%fisher_cls%LSS_lmax, fp%fisher_cls%LSS_number_windows )
982  l_max = max( l_max, fp%fisher_cls%l_max_TT, &
983  & fp%fisher_cls%l_max_EE, fp%fisher_cls%l_max_BB )
984 
985  ! allocation:
986  allocate(cl_fiducial(cl_dim, cl_dim, l_max-l_min+1), stat = allocatestatus)
987  if (allocatestatus /= 0) stop "Allocation failed. Not enough memory to allocate cl_fiducial"
988  allocate(cl_temp(cl_dim, cl_dim, l_max-l_min+1), stat = allocatestatus)
989  if (allocatestatus /= 0) stop "Allocation failed. Not enough memory to allocate cl_temp"
990  allocate(cl_derivative(cl_dim, cl_dim, l_max-l_min+1), stat = allocatestatus)
991  if (allocatestatus /= 0) stop "Allocation failed. Not enough memory to allocate cl_derivative"
992 
993  ! computation of the fiducial model:
994  time_1 = omp_get_wtime()
995  call cl_covariance_out( p, fp, cl_dim, l_min, l_max, cl_fiducial, err )
996  time_2 = omp_get_wtime() - time_1
997 
998  if ( present(outroot) ) call save_cl_covariance_to_file( cl_fiducial, cl_dim, l_min, l_max, filename=trim(outroot)//'fiducial.dat' )
999 
1000  if ( err == 0 ) then
1001  else if ( err == 4 ) then
1002  write(*,*) 'Fiducial model unstable. Calculation cannot proceed.'
1003  stop
1004  else
1005  write(*,*) 'Error in computing the fiducial. Calculation cannot proceed.'
1006  stop
1007  end if
1008 
1009  ! print some feedback
1010  if ( fp%cosmicfish_feedback >= 1 ) then
1011  write(*,'(a)') '**************************************************************'
1012  write(*,'(a,i10)') 'Dimension of the Cls matrix : ', cl_dim
1013  if ( fp%fisher_cls%Fisher_want_XC ) then
1014  write(*,'(a)') 'Using cross correlation.'
1015  else
1016  write(*,'(a)') 'Not using cross correlation.'
1017  end if
1018  write(*,'(a,i10)') 'Minimum l : ', l_min
1019  write(*,'(a,i10)') 'Maximum l : ', l_max
1020  write(*,'(a,i10)') 'Number of parameters : ', num_param
1021  write(*,'(a,f10.3,a)') 'Time taken for Cl calculation: ', time_2, ' (s)'
1022  write(*,'(a,i10)') 'Number of omp processes : ', omp_get_max_threads()
1023  if ( fp%adaptivity ) then
1024  write(*,'(a,f10.3,a)') 'Estimated completion time : ', time_2*20._dl*num_param, ' (s)'
1025  else
1026  write(*,'(a,f10.3,a)') 'Estimated completion time : ', time_2*2._dl*num_param, ' (s)'
1027  end if
1028  write(*,'(a)') '**************************************************************'
1029  end if
1030 
1031  ! compute the derivatives of the Cls:
1032  if ( fp%cosmicfish_feedback >= 1 ) then
1033  write(*,'(a)') '**************************************************************'
1034  write(*,'(a)') 'Computation of the Cls derivatives'
1035  write(*,'(a)') '**************************************************************'
1036  end if
1037 
1038  allocate(cl_derivative_array(num_param,cl_dim, cl_dim, l_max-l_min+1), stat = allocatestatus)
1039  if (allocatestatus /= 0) stop "Allocation failed. Not enough memory to allocate cl_derivative_array"
1040 
1041  time_1 = omp_get_wtime()
1042  do ind = 1, num_param
1043 
1044  if ( fp%cosmicfish_feedback >= 1 ) then
1045  write(*,'(a,i5,a,E15.5)') 'Doing parameter number: ', ind, ' value: ', param_array(ind)
1046  end if
1047 
1048  if ( param_array(ind) /= 0._dl ) then
1049  initial_step = param_array(ind)*3.0_dl/100._dl
1050  else
1051  initial_step = 3.0_dl/100._dl
1052  end if
1053 
1054  call cls_derivative( p, fp, ind, &
1055  & cl_fiducial, cl_derivative, cl_temp, &
1056  & initial_step, &
1057  & cl_dim, l_min, l_max, num_param, &
1058  & err )
1059 
1060  if (err /= 0) print*, 'WARNING: something went wrong with the Cls derivative of parameter:', ind
1061 
1062  if ( present(outroot) .and. fp%cosmicfish_feedback >= 2 ) then
1063  call save_cl_covariance_to_file( cl_derivative, cl_dim, l_min, l_max, filename=trim(outroot)//'derivative_'//trim(adjustl(integer_to_string(ind)))//'.dat' )
1064  end if
1065 
1066  cl_derivative_array(ind,:,:,:) = cl_derivative(:,:,:)
1067 
1068  end do
1069  time_2 = omp_get_wtime() - time_1
1070 
1071  if ( fp%cosmicfish_feedback >= 1 ) then
1072  write(*,*)
1073  write(*,'(a,f8.3,a)') 'Time to compute all the derivatives:', time_2, ' (s)'
1074  end if
1075 
1076  ! add the noise to the fiducial model:
1077  if ( fp%cosmicfish_feedback >= 1 ) then
1078  write(*,'(a)') '**************************************************************'
1079  write(*,'(a)') 'Adding experimental noise and computing Fisher:'
1080  write(*,'(a)') '**************************************************************'
1081  end if
1082 
1083  time_1 = omp_get_wtime()
1084  ! save the fiducial
1085  cl_temp = cl_fiducial
1086  ! add noise to the Cls
1087  call add_noise_cls_to_fiducial( p, fp, cl_fiducial, cl_dim, l_min, l_max, err)
1088  if ( err /= 0 ) stop 'Something went wrong with the computation of the cls noise'
1089 
1090  time_2 = omp_get_wtime() - time_1
1091 
1092  if ( present(outroot) .and. fp%cosmicfish_feedback >= 2 ) then
1093  call save_cl_covariance_to_file( cl_fiducial, cl_dim, l_min, l_max, filename=trim(outroot)//'fid_noise.dat' )
1094  end if
1095 
1096  if ( fp%cosmicfish_feedback >= 1 ) then
1097  write(*,'(a,f8.3,a)') 'Time to compute the noise Cls:', time_2, ' (s)'
1098  end if
1099 
1100  ! compute the Fisher matrix:
1101  allocate( dummy_matrix_1(cl_dim, cl_dim), dummy_matrix_2(cl_dim, cl_dim) )
1102  allocate( dummy_matrix_3(cl_dim, cl_dim), dummy_matrix_4(cl_dim, cl_dim) )
1103  allocate( dummy_matrix_5(cl_dim, cl_dim) )
1104 
1105  ! invert the fiducial + noise:
1106  do ind = 1, l_max-l_min+1
1107 
1108  forall ( k = 1:cl_dim, l = 1:cl_dim )
1109  dummy_matrix_1(k,l) = cl_fiducial(k,l, ind)
1110  end forall
1111 
1112  call sym_matrix_inversion( dummy_matrix_1, err )
1113  if ( err /= 0 ) then
1114  write(*,*) 'Matrix inversion failed. Calculation cannot proceed. Error at l= ', ind
1115  stop
1116  end if
1117 
1118  forall ( k = 1:cl_dim, l = 1:cl_dim )
1119  cl_fiducial(k,l, ind) = dummy_matrix_1(k,l)
1120  end forall
1121 
1122  end do
1123 
1124  ! compute the Fisher:
1125  do ind = 1, num_param
1126  do ind2 = 1, num_param
1127 
1128  temp2 = 0._dl
1129 
1130  do i = 1, l_max-l_min+1
1131  ! Cosmic variance term:
1132  temp = REAL(i+l_min) -0.5_dl
1133  ! Initialize all the matrix to use:
1134  forall ( j = 1:cl_dim, k = 1:cl_dim )
1135  dummy_matrix_1(j, k) = cl_derivative_array(ind, j, k, i)
1136  dummy_matrix_2(j, k) = cl_derivative_array(ind2, j, k, i)
1137  dummy_matrix_3(j, k) = cl_fiducial(j, k, i)
1138  end forall
1139  ! do the products:
1140  dummy_matrix_1 = matmul(dummy_matrix_1, dummy_matrix_3)
1141  dummy_matrix_1 = matmul(dummy_matrix_3, dummy_matrix_1)
1142  dummy_matrix_4 = matmul(dummy_matrix_2, dummy_matrix_1)
1143  ! take the trace of the result:
1144  temp1 = 0._dl
1145  do j = 1, cl_dim
1146  temp1 = temp1 + dummy_matrix_4(j,j)
1147  end do
1148  ! sum to build up the Fisher matrix entry:
1149  temp2 = temp2 + temp*temp1
1150 
1151  end do
1152 
1153  fisher_matrix(ind, ind2) = temp2
1154 
1155  end do
1156  end do
1157  time_2 = omp_get_wtime() - time_1
1158 
1159  if ( fp%cosmicfish_feedback >= 1 ) then
1160  write(*,'(a,f8.3,a)') 'Time to compute the Fisher matrix:', time_2, ' (s)'
1161  end if
1162 
1163  end subroutine fisher_cls
1164 
1165 end module fisher_calculator_cls
matrix_utilities::sym_matrix_inversion
subroutine, public sym_matrix_inversion(A, error_code)
This subroutine computes the inverse of a square symmetric matrix by means of LAPACK LU decomposition...
Definition: 001_Matrix_utilities.f90:42
cosmicfish_camb_interface::params_array_to_camb_param
subroutine, public params_array_to_camb_param(P, FP, num_param, params_array)
This subroutine updates the values of the parameters P and FP based on the values contained in the pa...
Definition: 006_cosmicfish_camb_interface.f90:1187
fisher_calculator_cls::save_cl_covariance_to_file
subroutine, public save_cl_covariance_to_file(cl_cov_in, cl_dim, l_min, l_max, filename)
This subroutine saves to file the input Cls covariance matrix.
Definition: 008_Fisher_calculator_Cls.f90:386
fisher_calculator_cls::cls_derivative
subroutine, public cls_derivative(P, FP, param_num, cl_initial, cl_derivative, cl_error, initial_step, out_dim, l_min, l_max, num_param, error_code)
This subroutine computes the derivative of the Cls covariance matrix with a fixed stepsize finite dif...
Definition: 008_Fisher_calculator_Cls.f90:427
cosmicfish_camb_interface::camb_params_to_params_array
subroutine, public camb_params_to_params_array(P, FP, num_param, params_array)
This subroutine updates the values of the parameters array based on the values of P and FP...
Definition: 006_cosmicfish_camb_interface.f90:787
cosmicfish_types
This module contains the definitions of derived data types used to store the informations about the f...
Definition: 004_cosmicfish_types.f90:26
fisher_calculator_cls_windows
This module contains the definition of the window functions that are used by CAMB sources...
Definition: 008_Fisher_calculator_Cls_windows.f90:26
matrix_utilities
This module contains the subroutine and functions to manipulate and modify matrices.
Definition: 001_Matrix_utilities.f90:26
cosmicfish_camb_interface
This module contains the interface layer between CosmicFish and CAMB.
Definition: 006_cosmicfish_camb_interface.f90:25
cosmicfish_utilities
This module contains several general purpose utilities.
Definition: 000_Utilities.f90:24
fisher_calculator_cls::add_noise_cls_to_fiducial
subroutine, public add_noise_cls_to_fiducial(P, FP, cl_fiducial, out_dim, l_min, l_max, error_code)
This subroutine computes and adds the noise cls to the fiducial cls. On output the fiducial cls will ...
Definition: 008_Fisher_calculator_Cls.f90:681
lensnoise
This module contains the code to compute the noise for CMB lensing.
Definition: 007_FUTURCMB_lensnoise.f90:25
fisher_calculator_cls::fisher_cls
subroutine, public fisher_cls(P, FP, num_param, Fisher_Matrix, outroot)
This subroutine computes the Fisher matrix for Cls.
Definition: 008_Fisher_calculator_Cls.f90:946
cosmicfish_utilities::integer_to_string
character(10) function, public integer_to_string(number)
This function converts an integer to a string. Usefull for numbered files output. ...
Definition: 000_Utilities.f90:39
lensnoise::calc_nldd
subroutine, public calc_nldd(incls, inlcls, clnoise, nldd, jjmaxTmax, lmaxM)
This subroutine computes the noise for the Okamoto & Hu's estimator of the CMB lensing.
Definition: 007_FUTURCMB_lensnoise.f90:53
fisher_calculator_cls
This module contains the code that computes the Fisher matrix for Cls.
Definition: 008_Fisher_calculator_Cls.f90:24
fisher_calculator_cls::cl_covariance_out
subroutine, public cl_covariance_out(P, FP, cl_dim, l_min, l_max, cl_out, err)
This subroutine that gives the covariance matrix of the Cls for the desired model. Calls internally CAMB_get results.
Definition: 008_Fisher_calculator_Cls.f90:87
fisher_calculator_cls_windows::init_camb_sources_windows
subroutine, public init_camb_sources_windows(FP)
This function initialises the camb sources number counts and lensing window function. The default is the gaussian window that is included in camb sources. If another window is chosen it will be initialised by the CosmicFish value.
Definition: 008_Fisher_calculator_Cls_windows.f90:48
fisher_calculator_cls::dimension_cl_covariance
subroutine, public dimension_cl_covariance(P, FP, num)
This subroutine computes the number of entrance of the cls covariance matrix The matrix of the Cls is...
Definition: 008_Fisher_calculator_Cls.f90:56
cosmicfish_types::cosmicfish_params
This derived data type contains all the informations that the cosmicfish code needs to run...
Definition: 004_cosmicfish_types.f90:215
Generated on Tue Jun 21 2016 03:18:43 for CosmicFish by   doxygen 1.8.9.1