009_Fisher_calculator_SN.f90

Go to the documentation of this file.

!----------------------------------------------------------------------------------------
!
! This file is part of CosmicFish.
!
! Copyright (C) 2015-2016 by the CosmicFish authors
!
! The CosmicFish code is free software;
! You can use it, redistribute it, and/or modify it under the terms
! of the GNU General Public License as published by the Free Software Foundation;
! either version 3 of the License, or (at your option) any later version.
! The full text of the license can be found in the file LICENSE at
! the top level of the CosmicFish distribution.
!
!----------------------------------------------------------------------------------------


!----------------------------------------------------------------------------------------


module fisher_calculator_sn

    use precision
    use modelparams
    use cosmicfish_types
    use omp_lib
    use cosmicfish_utilities
    use cosmicfish_camb_interface

#ifdef COSMICFISH_EFTCAMB
    use eftinitialization
#endif

    implicit none

    private

    public observed_sn_magnitude, generate_sn_mock_data, generate_sn_mock_covariance, save_sn_mock_to_file, &
        & observed_sn_magnitude_derivative, fisher_sn

contains

    ! ---------------------------------------------------------------------------------------------
    subroutine observed_sn_magnitude( P, FP, sn_magnitude, redshift, color, stretch, num, err )

        implicit none

        Type(cambparams)      , intent(in)  :: P
        Type(cosmicfish_params), intent(in)  :: FP
        real(dl), dimension(:), intent(out) :: sn_magnitude
        real(dl), dimension(:), intent(in)  :: redshift
        real(dl), dimension(:), intent(in)  :: color
        real(dl), dimension(:), intent(in)  :: stretch
        integer               , intent(in)  :: num
        integer, intent(out)    :: err

        ! other definitions:
        integer  :: i
        real(dl) :: alpha, beta, zero_mag
        logical  :: EFTsuccess
        ! initialization:
        sn_magnitude = 0._dl
        err          = 0

        call cambparams_set(p)

#ifdef COSMICFISH_EFTCAMB
        ! check stability if EFTCAMB is used
        if (p%EFTflag /= 0) then
            call eftcamb_initialization(eftsuccess)
            if (.not. eftsuccess) then
                write(*,*) 'EFTCAMB: unstable model'
                err = 4
                return
            end if
        end if
#endif

        alpha    = fp%fisher_SN%alpha_SN
        beta     = fp%fisher_SN%alpha_SN
        zero_mag = fp%fisher_SN%M0_SN

        ! compute the supernova magnitude array:
        do i=1, num
            sn_magnitude(i) = 5._dl*log10( luminositydistance(redshift(i)) ) &
                & -alpha*stretch(i) +beta*color(i) + zero_mag -25._dl
        end do

    end subroutine observed_sn_magnitude

    ! ---------------------------------------------------------------------------------------------
    subroutine generate_sn_mock_data( P, FP, sn_mock )

        use random_utilities

        implicit none

        Type(cambparams)        , intent(in)  :: P
        Type(cosmicfish_params)  , intent(in)  :: FP
        real(dl), dimension(:,:), intent(out) :: sn_mock

        ! other definitions:
        integer  :: i, j, k

        ! initialize the random number generator:
        rand_feedback = 0
        call initrandom()

        ! generate the mock data:
        i = 1
        do j=1, fp%fisher_SN%number_SN_windows
            do k=1, fp%fisher_SN%SN_number(j)
                sn_mock(1,i) = random_uniform( fp%fisher_SN%SN_redshift_start(j), fp%fisher_SN%SN_redshift_end(j)) ! mock redshift
                sn_mock(2,i) = random_gaussian( 0._dl, fp%fisher_SN%color_dispersion )                             ! mock color
                sn_mock(3,i) = random_gaussian( 0._dl, fp%fisher_SN%stretch_dispersion )                           ! mock stretch
                i = i + 1
            end do
        end do

    end subroutine generate_sn_mock_data

    ! ---------------------------------------------------------------------------------------------
    subroutine generate_sn_mock_covariance( P, FP, sn_mock, sn_mock_covariance )

        use constants, only : c

        implicit none

        Type(cambparams)        , intent(in)  :: P
        Type(cosmicfish_params)  , intent(in)  :: FP
        real(dl), dimension(:,:), intent(in)  :: sn_mock
        real(dl), dimension(:)  , intent(out) :: sn_mock_covariance

        ! other definitions:
        integer  :: i, j, k

        ! generate the mock covariance first terms:
        do i=1, fp%fisher_SN%total_SN_number
            sn_mock_covariance(i) = (fp%fisher_SN%magnitude_sigma)**2 &
                & +(fp%fisher_SN%sigma_lens_0*sn_mock(1,i))**2 & ! lensing term
                & +(5._dl/sn_mock(1,i)/log(10._dl)*fp%fisher_SN%c_sigmaz/c*1000._dl)**2 ! redshift term
        end do

        ! add the systematic effects to the covariance:
        do i=1, fp%fisher_SN%total_SN_number
            sn_mock_covariance(i) = sn_mock_covariance(i) &
                & + fp%fisher_SN%beta_SN**2*(fp%fisher_SN%dcolor_offset+fp%fisher_SN%dcolor_zcorr*sn_mock(1,i)**2)**2  &
                & + fp%fisher_SN%alpha_SN**2*(fp%fisher_SN%dshape_offset+fp%fisher_SN%dshape_zcorr*sn_mock(1,i)**2)**2  &
                & + 2._dl*fp%fisher_SN%alpha_SN*(fp%fisher_SN%cov_ms_offset+fp%fisher_SN%cov_ms_zcorr*sn_mock(1,i)**2)     &
                & - 2._dl*fp%fisher_SN%beta_SN*(fp%fisher_SN%cov_mc_offset+fp%fisher_SN%cov_mc_zcorr*sn_mock(1,i)**2)     &
                & - 2._dl*fp%fisher_SN%alpha_SN*fp%fisher_SN%beta_SN*(fp%fisher_SN%cov_sc_offset+fp%fisher_SN%cov_sc_zcorr*sn_mock(1,i)**2)
        end do

    end subroutine generate_sn_mock_covariance

    ! ---------------------------------------------------------------------------------------------
    subroutine save_sn_mock_to_file( mock_size, sn_mock_data, sn_mock_covariance, filename, magnitude )

        implicit none

        integer                 , intent(in) :: mock_size
        real(dl), dimension(:,:), intent(in) :: sn_mock_data
        real(dl), dimension(:)  , intent(in) :: sn_mock_covariance
        character(len=*)        , intent(in) :: filename

        real(dl), dimension(:)  , intent(in), optional :: magnitude

        ! other definitions:
        integer  :: i

        ! write the mock:
        open(unit=666, file=filename, action="write", status="replace")

        if ( present(magnitude) ) then

            write(666,'(a)') "# number of the SN, redshift, magnitude, color, stretch, covariance "
            do i = 1, mock_size
                write(666,'(I10,5E15.5)') i, sn_mock_data(1,i), magnitude(i), sn_mock_data(2,i), sn_mock_data(3,i), sn_mock_covariance(i)
            end do

        else

            write(666,'(a)') "# number of the SN, redshift, color, stretch, covariance "
            do i = 1, mock_size
                write(666,'(I10,4E15.5)') i, sn_mock_data(1,i), sn_mock_data(2,i), sn_mock_data(3,i), sn_mock_covariance(i)
            end do

        end if

        close(666)

    end subroutine save_sn_mock_to_file

    ! ---------------------------------------------------------------------------------------------
    subroutine observed_sn_magnitude_derivative( P, FP, param_num, &
        & mag_initial, mag_derivative, mag_der_error, &
        & sn_mock, &
        & initial_step, &
        & num_param, &
        & error_code  )

        implicit none

        Type(cambparams)       , intent(in) :: P
        Type(cosmicfish_params) , intent(in) :: FP

        integer , intent(in)  :: param_num
        integer , intent(in)  :: num_param

        real(dl), dimension(:), intent(in)  :: mag_initial
        real(dl), dimension(:), intent(out) :: mag_derivative
        real(dl), dimension(:), intent(out) :: mag_der_error

        real(dl), dimension(:,:), intent(in) :: sn_mock

        real(dl), intent(in)  :: initial_step
        integer , intent(out) :: error_code

        ! definitions:
        real(dl), dimension(:), allocatable :: sn_temp_1, sn_temp_2, errt
        integer , dimension(:), allocatable :: accept_mask
        real(dl), dimension(:,:,:), allocatable :: a
        real(dl), dimension(num_param) :: param_vector, fiducial_param_vector, temp_param_vector
        ! parameters of the algorithm:
        real(dl), parameter :: con  = 1.4_dl  ! parameter that decides the decrease in the stepsize h
        real(dl), parameter :: big  = 1.d+30  ! a big number used to fill the error array initially
        real(dl), parameter :: safe = 2._dl   ! parameter used to define the stopping cryterium
        integer , parameter :: ntab = 30      ! maximum number of Cls computation
        real(dl) :: con2 = con*con
        integer  :: i, j, k, l, m, side, err, AllocateStatus
        real(dl) :: fac, hh
        real(dl) :: temp1, temp2, temp3
        Type(cambparams) :: P_temp
        Type(cosmicfish_params) :: FP_temp
        integer  :: accepted
        real(dl) :: acceptance
        integer  :: mock_sn_size
        real(dl), dimension(:), allocatable :: redshift, color, stretch

        ! initialize output variables:
        error_code     = 0
        mag_derivative = 0._dl
        mag_der_error  = 0._dl

        ! check input values are correct
        if ( initial_step == 0._dl ) then
            stop 'Observed_SN_Magnitude_derivative initial stepsize is zero'
        end if

        mock_sn_size = fp%fisher_SN%total_SN_number

        ! allocate the arrays:
        if ( fp%adaptivity ) then
            allocate( sn_temp_1(mock_sn_size),  &
                & sn_temp_2(mock_sn_size),  &
                & errt(mock_sn_size),       &
                & accept_mask(mock_sn_size),&
                & stat = allocatestatus )
            if (allocatestatus /= 0) stop "Observed_SN_Magnitude_derivative: Allocation failed. Not enough memory"
            allocate( a(ntab, ntab, mock_sn_size) , stat = allocatestatus )
            if (allocatestatus /= 0) stop "Observed_SN_Magnitude_derivative: Allocation failed. Not enough memory for the derivative array"
        else
            allocate( sn_temp_1(mock_sn_size),  &
                & sn_temp_2(mock_sn_size),  &
                & stat = allocatestatus )
            if (allocatestatus /= 0) stop "Observed_SN_Magnitude_derivative: Allocation failed. Not enough memory"
            allocate( a(1,1, mock_sn_size) , stat = allocatestatus )
            if (allocatestatus /= 0) stop "Observed_SN_Magnitude_derivative: Allocation failed. Not enough memory for the derivative array"
        end if

        allocate( redshift(mock_sn_size),       &
            & color(mock_sn_size),       &
            & stretch(mock_sn_size),     &
            & stat = allocatestatus )

        ! copy the data from the sn mock:
        redshift(:) = sn_mock(1,:)
        color(:)    = sn_mock(2,:)
        stretch(:)  = sn_mock(3,:)

        ! initialize the parameter array:
        p_temp  = p
        fp_temp = fp

        call camb_params_to_params_array( p_temp, fp, num_param, param_vector )

        fiducial_param_vector = param_vector
        temp_param_vector     = param_vector

        ! initialize the computation:
        hh            = initial_step
        mag_der_error = big
        side          = 0 ! both sides
        if ( fp%adaptivity ) accept_mask = 0

        ! do the initial step on the right:
        param_vector(param_num)      = fiducial_param_vector(param_num) + hh
        call params_array_to_camb_param( p_temp, fp_temp, num_param, param_vector )
        call observed_sn_magnitude( p_temp, fp_temp, sn_temp_1, redshift, color, stretch, mock_sn_size, err )
        ! if on the right we cannot go we go to the left:
        if ( err /= 0 ) side = 1
        ! do the initial step on the left:
        temp_param_vector(param_num) = fiducial_param_vector(param_num) - hh
        call params_array_to_camb_param( p_temp, fp_temp, num_param, temp_param_vector )
        call observed_sn_magnitude( p_temp, fp_temp, sn_temp_2, redshift, color, stretch, mock_sn_size, err )

        ! if on the left we cannot go decide what to do:
        if ( err /= 0 ) then
            if ( side == 1 ) then ! we cannot go to the left nor to the right. Return derivative zero with error code.
                error_code      = 4
                mag_der_error   = 0.0_dl
                mag_derivative   = 0.0_dl
                return
            else if ( side == 0 ) then ! we cannot go to the left but we can go to the right:
                side = 2
            end if
        end if

        ! store the results based on the side we decided to go:
        if (side == 0) then ! both sides
            a(1,1,:) = (sn_temp_1(:) - sn_temp_2(:))/(2.0_dl*hh)
        else if (side == 1) then ! left side: increment negative
            a(1,1,:) = (sn_temp_2(:) - mag_initial(:))/(-hh)
        else if (side == 2) then ! right side: increment positive
            a(1,1,:) = (sn_temp_1(:) - mag_initial(:))/(hh)
        end if

        ! fixed stepsize derivative if adaptivity is not wanted:
        if ( .not. fp%adaptivity ) then
            mag_derivative(:) = a(1,1,:)
            mag_der_error     = 0.0_dl
            deallocate( sn_temp_1, sn_temp_2, a )
            return
        end if

        ! start the derivative computation:
        do i = 2, ntab

            ! feedback for the adaptive derivative:
            if ( fp%cosmicfish_feedback >= 2 ) then
                accepted   = sum(accept_mask)
                acceptance = REAL(accepted)/REAL(mock_sn_size)
                print*, 'iteration:', i
                print*, '  stepsize:', hh
                print*, '  accepted elements', accepted, 'accepted ratio', acceptance
            end if

            ! decrease the stepsize:
            hh = hh/con

            ! compute the finite difference:
            if (side == 0) then
                param_vector(param_num)      = fiducial_param_vector(param_num) + hh
                temp_param_vector(param_num) = fiducial_param_vector(param_num) - hh
                ! compute for + hh
                call params_array_to_camb_param( p_temp, fp_temp, num_param, param_vector )
                call observed_sn_magnitude( p_temp, fp_temp, sn_temp_1, redshift, color, stretch, mock_sn_size, err )
                if ( err /= 0 ) then
                    error_code     = 4
                    mag_der_error  = 0.0_dl
                    mag_derivative = 0.0_dl
                    return
                end if
                ! compute for - hh
                call params_array_to_camb_param( p_temp, fp_temp, num_param, temp_param_vector )
                call observed_sn_magnitude( p_temp, fp_temp, sn_temp_2, redshift, color, stretch, mock_sn_size, err )
                if ( err /= 0 ) then
                    error_code     = 4
                    mag_der_error  = 0.0_dl
                    mag_derivative = 0.0_dl
                    return
                end if
                ! store the result
                a(1,i,:) = (sn_temp_1(:) - sn_temp_2(:))/(2.0_dl*hh)
            else if (side==1) then
                param_vector(param_num)      = fiducial_param_vector(param_num) - hh
                call params_array_to_camb_param( p_temp, fp_temp, num_param, param_vector )
                call observed_sn_magnitude( p_temp, fp_temp, sn_temp_1, redshift, color, stretch, mock_sn_size, err )
                if ( err /= 0 ) then
                    error_code     = 4
                    mag_der_error  = 0.0_dl
                    mag_derivative = 0.0_dl
                    return
                end if
                ! store the result
                a(1,i,:) = (sn_temp_1(:) - mag_initial(:))/(-hh)
            else if (side==2) then
                param_vector(param_num)      = fiducial_param_vector(param_num) + hh
                call params_array_to_camb_param( p_temp, fp_temp, num_param, param_vector )
                call observed_sn_magnitude( p_temp, fp_temp, sn_temp_1, redshift, color, stretch, mock_sn_size, err )
                if ( err /= 0 ) then
                    error_code     = 4
                    mag_der_error  = 0.0_dl
                    mag_derivative = 0.0_dl
                    return
                end if
                ! store the result
                a(1,i,:) = (sn_temp_1(:) - mag_initial(:))/(hh)
            end if

            ! now do the extrapolation table:
            fac=con2
            do j = 2, i

                forall ( k = 1:mock_sn_size, accept_mask(k) == 0 )
                    ! compute the interpolation table:
                    a(j,i,k) = (a(j-1,i,k)*fac-a(j-1,i-1,k))/(fac-1._dl)
                    ! estimate the error:
                    errt(k)  = max( abs(a(j,i,k)-a(j-1,i,k)), &
                        & abs(a(j,i,k)-a(j-1,i-1,k)))
                end forall

                fac=con2*fac

                ! if there is an improvement save the result:

                forall ( k = 1:mock_sn_size, &
                    & errt(k) .le. mag_der_error(k) .and. accept_mask(k) == 0 )

                    mag_der_error(k) = errt(k)
                    mag_derivative(k) = a(j,i,k)
                end forall

            end do
            ! accept a value of the derivative if the value of the extrapolation is stable
            ! Store this information into the mask matrix
            ! so that the value contained in cl_error and cl_derivative is not overwritten.
            forall ( k = 1:mock_sn_size, &
                & accept_mask(k) == 0 .and. &
                & (a(i,i,k) - a(i-1,i-1,k)) .ge. safe*mag_der_error(k) )

                accept_mask(k) = 1
            end forall

            ! if all the entries match the required accuracy exit the loop
            if ( sum(accept_mask) == mock_sn_size ) exit

        end do

        ! quality check on the results:
        if ( sum(accept_mask) /= mock_sn_size ) then
            error_code = 3
        end if

        ! deallocate the arrays:
        if ( fp%adaptivity ) then
            deallocate( sn_temp_1, sn_temp_2, errt, accept_mask, a )
        end if

        return

    end subroutine observed_sn_magnitude_derivative

    ! ---------------------------------------------------------------------------------------------
    subroutine fisher_sn( P, FP, num_param, Fisher_Matrix, outroot )

        implicit none

        Type(cambparams)       , intent(in) :: P
        Type(cosmicfish_params) , intent(in) :: FP

        integer , intent(in)  :: num_param
        real(dl), dimension(num_param,num_param), intent(out) :: Fisher_Matrix
        character(len=*), intent(in), optional :: outroot

        real(dl) :: time_1, time_2, time_1_MC, time_2_MC
        real(dl) :: initial_step, temp, temp1, temp2
        integer  :: AllocateStatus, err, ind, ind2
        integer  :: i, j, k, l, number
        real(dl), allocatable :: param_array(:)
        real(dl), allocatable, dimension(:,:) :: sn_mock
        real(dl), allocatable, dimension(:)   :: sn_mock_covariance
        real(dl), allocatable, dimension(:)   :: sn_fiducial, sn_derivative, sn_temp
        real(dl), allocatable, dimension(:,:) :: sn_derivative_array

        integer  :: MC_i
        real(dl), dimension(num_param,num_param) :: MC_Fisher_Matrix

        ! allocate a parameter array:
        allocate( param_array(num_param) )
        ! write the parameters in the parameter array:
        call camb_params_to_params_array( p, fp, num_param, param_array )

        number = fp%fisher_SN%total_SN_number

        ! allocation:
        allocate( sn_mock(3,number), stat = allocatestatus )
        if (allocatestatus /= 0) stop "Allocation failed. Not enough memory to allocate SN mock."
        allocate( sn_mock_covariance(number), stat = allocatestatus )
        if (allocatestatus /= 0) stop "Allocation failed. Not enough memory to allocate SN mock covariance."
        allocate( sn_fiducial(number), sn_derivative(number), sn_temp(number), stat = allocatestatus )
        if (allocatestatus /= 0) stop "Allocation failed. Not enough memory to allocate SN fiducial."
        allocate( sn_derivative_array( num_param, number ), stat = allocatestatus)
        if (allocatestatus /= 0) stop "Allocation failed. Not enough memory to allocate sn_derivative_array"

        ! do the Monte Carlo:

        mc_fisher_matrix = 0._dl

        time_1_mc = omp_get_wtime()
        do mc_i=1, fp%fisher_SN%SN_Fisher_MC_samples

            ! generate the SN mock:
            call generate_sn_mock_data( p, fp, sn_mock )

            ! generate the mock covariance:
            call generate_sn_mock_covariance( p, fp, sn_mock, sn_mock_covariance )

            ! compute the fiducial:
            time_1 = omp_get_wtime()
            call observed_sn_magnitude( p, fp, sn_fiducial, sn_mock(1,:), sn_mock(2,:), sn_mock(3,:), &
                & number, err )
            time_2 = omp_get_wtime() - time_1

            if ( present(outroot) ) call save_sn_mock_to_file( number, sn_mock, sn_mock_covariance, &
                & filename=trim(outroot)//'sn_fiducial_'//trim(adjustl(integer_to_string(mc_i)))//'.dat', magnitude = sn_fiducial )

            if ( err == 0 ) then
            else if ( err == 4 ) then
                write(*,*) 'Fiducial model unstable. Calculation cannot proceed.'
                stop
            else
                write(*,*) 'Error in computing the fiducial. Calculation cannot proceed.'
                stop
            end if

            ! print some feedback
            if ( mc_i==1 .and. fp%cosmicfish_feedback >= 1 ) then
                write(*,'(a)')         '**************************************************************'
                write(*,'(a,i10)')         'Number of supernovae         : ', number
                write(*,'(a,i10)')         'Number of parameters         : ', num_param
                write(*,'(a,i10)')         'Number of MC repetitions     : ', fp%fisher_SN%SN_Fisher_MC_samples
                write(*,'(a,f10.3,a)')     'Time taken for SN calculation: ', time_2, ' (s)'
                write(*,'(a,i10)')         'Number of omp processes      : ', omp_get_max_threads()
                if ( fp%adaptivity ) then
                    write(*,'(a,f10.3,a)') 'Estimated completion time    : ', fp%fisher_SN%SN_Fisher_MC_samples*time_2*20._dl*num_param, ' (s)'
                else
                    write(*,'(a,f10.3,a)') 'Estimated completion time    : ', fp%fisher_SN%SN_Fisher_MC_samples*time_2*2._dl*num_param, ' (s)'
                end if
                write(*,'(a)')         '**************************************************************'
            end if

            if ( mc_i>1 .and. fp%cosmicfish_feedback >= 1 ) then
                write(*,'(a)') '**************************************************************'
                write(*,'(a,i10)') 'Doing MC computation of SN Fisher sample: ', mc_i
                write(*,'(a)') '**************************************************************'
            end if

            ! compute the derivatives of the Cls:
            if ( fp%cosmicfish_feedback >= 1 ) then
                write(*,'(a)') 'Computation of the SN magnitude derivatives'
            end if

            time_1 = omp_get_wtime()
            do ind = 1, num_param

                if ( fp%cosmicfish_feedback >= 1 ) then
                    write(*,'(a,i5,a,E15.5)') 'Doing parameter number: ', ind, ' value: ', param_array(ind)
                end if

                if ( param_array(ind) /= 0._dl ) then
                    initial_step = param_array(ind)*3.0_dl/100._dl
                else
                    initial_step = 3.0_dl/100._dl
                end if

                call observed_sn_magnitude_derivative( p, fp, ind, &
                    & sn_fiducial, sn_derivative, sn_temp, &
                    & sn_mock, &
                    & initial_step, &
                    & num_param, &
                    & err  )

                if (err /= 0) print*, 'WARNING: something went wrong with the SN derivative of parameter:', ind, 'Error code:', err

                if ( mc_i==1 .and. present(outroot) .and. fp%cosmicfish_feedback >= 2 ) then
                    call save_sn_mock_to_file( number, sn_mock, sn_mock_covariance, &
                        & filename=trim(outroot)//'sn_derivative_'//trim(adjustl(integer_to_string(ind)))//'.dat', magnitude = sn_derivative )
                end if

                sn_derivative_array(ind,:) = sn_derivative(:)

            end do
            time_2 = omp_get_wtime() - time_1

            if ( fp%cosmicfish_feedback >= 1 ) then
                write(*,*)
                write(*,'(a,f8.3,a)') 'Time to compute all the derivatives:', time_2, ' (s)'
            end if

            ! compute the Fisher:
            do ind = 1, num_param
                do ind2 = 1, num_param

                    temp2 = 0._dl

                    do k = 1, number
                        temp2 = temp2 + sn_derivative_array(ind,k)*sn_derivative_array(ind2,k)/sn_mock_covariance(k)
                    end do

                    mc_fisher_matrix(ind, ind2) = mc_fisher_matrix(ind, ind2) + temp2

                end do
            end do
            time_2 = omp_get_wtime() - time_1

            if ( fp%cosmicfish_feedback >= 1 ) then
                write(*,'(a,f8.3,a)') 'Time to compute the Fisher matrix:', time_2, ' (s)'
            end if

        end do
        time_2_mc = omp_get_wtime() - time_1

        fisher_matrix = mc_fisher_matrix/fp%fisher_SN%SN_Fisher_MC_samples

        if ( fp%cosmicfish_feedback >= 1 ) then
            write(*,'(a)') '**************************************************************'
            write(*,'(a,f8.3,a)') 'Time to compute the MC Fisher matrix:', time_2_mc, ' (s)'
            write(*,'(a)') '**************************************************************'
        end if

    end subroutine fisher_sn

end module fisher_calculator_sn