012_Fisher_calculator_derived.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_derived

    use constants, only: k_b, ev, mass_ratio_he_h
    use precision
    use transfer, only: transfer_sortandindexredshifts, transfer_output_sig8, mt
    use modelparams
    use cosmicfish_types
    use omp_lib
    use cosmicfish_utilities
    use cosmicfish_camb_interface
    use camb

#ifdef COSMICFISH_EFTCAMB
    use eftinitialization
#endif

    implicit none

    private

    public  dimension_derived_parameters, compute_derived_parameters, fisher_derived_param_names, &
        & derived_parameter_derivative, fisher_derived, save_derived_fisher_to_file, save_derived_paramnames_to_file

contains

    ! ---------------------------------------------------------------------------------------------
    subroutine derived_fisher_parameter_compatibility( P, FP )

        implicit none

        Type(cambparams)        :: p
        Type(cosmicfish_params)  :: fp

        integer :: i

        integer :: j
        character(LEN=20) :: str, str2

        ! if we want tomography for the derived parameters fix the transfer functions parameters accorgingly:
        if ( fp%fisher_der%want_sigma8 ) then
            ! write the relevant quantities in Transfer:
            p%PK_WantTransfer = .true.
            p%WantTransfer    = .true.
            p%Transfer%high_precision = .true.
            p%Transfer%PK_num_redshifts = fp%fisher_der%FD_num_redshift
            do i=1, fp%fisher_der%FD_num_redshift
                p%transfer%PK_redshifts(i) = fp%fisher_der%FD_redshift(i)
            end do
            ! ensure ordering in redshift:
            call transfer_sortandindexredshifts(p%Transfer)
            ! copy back the redshift:
            do i=1, fp%fisher_der%FD_num_redshift
                fp%fisher_der%FD_redshift(i) = p%transfer%PK_redshifts(i)
            end do
        end if

        ! we do not want other things to be computed:
        num_redshiftwindows = 0

    end subroutine derived_fisher_parameter_compatibility

    ! ---------------------------------------------------------------------------------------------
    subroutine dimension_derived_parameters( P, FP, num )

        implicit none

        Type(cambparams)      , intent(in)  :: P
        Type(cosmicfish_params), intent(in)  :: FP
        integer               , intent(out) :: num

        ! start with the base derived parameters:
        num = 0
        if ( fp%fisher_der%want_omegab ) then
            num = num + 1
        end if
        if ( fp%fisher_der%want_omegac ) then
            num = num + 1
        end if
        if ( fp%fisher_der%want_omegan ) then
            num = num + 1
        end if
        if ( fp%fisher_der%want_omegav ) then
            num = num + 1
        end if
        if ( fp%fisher_der%want_omegak ) then
            num = num + 1
        end if
        if ( fp%fisher_der%want_omegam ) then
            num = num + 1
        end if
        if ( fp%fisher_der%want_theta  ) then
            num = num + 1
        end if
        if ( fp%fisher_der%want_mnu    ) then
            num = num + 1
        end if
        if ( fp%fisher_der%want_zre    ) then
            num = num + 1
        end if
        if ( fp%fisher_der%want_neff   ) then
            num = num + 1
        end if

        ! add sigma8 if wanted:
        if ( fp%fisher_der%want_sigma8 ) then
            num = num + fp%fisher_der%FD_num_redshift
        end if
        ! add Hubble if wanted:
        if ( fp%fisher_der%want_loghubble ) then
            num = num + fp%fisher_der%FD_num_redshift
        end if
        ! add DA if wanted:
        if ( fp%fisher_der%want_logDA ) then
            num = num + fp%fisher_der%FD_num_redshift
        end if

    end subroutine dimension_derived_parameters

    ! ---------------------------------------------------------------------------------------------
    subroutine compute_derived_parameters( P, FP, derived, num, err )

        implicit none

        Type(cambparams)      , intent(in)  :: P
        Type(cosmicfish_params), intent(in)  :: FP
        real(dl), dimension(:), intent(out) :: derived
        integer               , intent(in)  :: num
        integer, intent(out)    :: err

        ! other definitions:
        integer  :: i, j, nu_i, ind, error
        logical  :: EFTsuccess
        real(dl) :: conv, z
        type(cambdata)  :: outdata
        ! initialization:
        derived = 0._dl
        err     = 0

        ! set the camb parameters:
        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

        ! get the transfer functions for sigma_8
        if ( fp%fisher_der%want_sigma8 ) then
            error = 0
            call camb_gettransfers(p, outdata, error)
            ! test for errors in the computation
            if (error/=0) then
                write (*,*) 'Error result '//trim(global_error_message)
                err = 2
                return
            endif
        end if

        ! write the parameter array:
        ind = 0
        if ( fp%fisher_der%want_omegab ) then
            ind = ind + 1
            derived(ind) = cp%omegab ! Baryons density
        end if
        if ( fp%fisher_der%want_omegac ) then
            ind = ind + 1
            derived(ind) = cp%omegac ! CDM density
        end if
        if ( fp%fisher_der%want_omegan ) then
            ind = ind + 1
            derived(ind) = cp%omegan ! Neutrino density
        end if
        if ( fp%fisher_der%want_omegav ) then
            ind = ind + 1
            derived(ind) = cp%omegav ! DE density
        end if
        if ( fp%fisher_der%want_omegak ) then
            ind = ind + 1
            derived(ind) = cp%omegak ! Curvature density
        end if
        if ( fp%fisher_der%want_omegam ) then
            ind = ind + 1
            derived(ind) = 1-cp%omegak-cp%omegav ! Matter density
        end if
        if ( fp%fisher_der%want_theta  ) then
            ind = ind + 1
            derived(ind) = 100*cosmomctheta()    ! Theta CosmoMC
        end if
        if ( fp%fisher_der%want_mnu    ) then
            ind = ind + 1
            if (cp%Num_Nu_Massive == 1) then
                do nu_i=1, cp%Nu_mass_eigenstates
                    conv = k_b*(8*grhor/grhog/7)**0.25*cp%tcmb/ev * &
                        (cp%nu_mass_degeneracies(nu_i)/cp%nu_mass_numbers(nu_i))**0.25 !approx 1.68e-4
                    derived(ind) = conv*nu_masses(nu_i)
                end do
            else
                derived(ind) = 0._dl ! m_nu in eV
            end if
        end if
        if ( fp%fisher_der%want_zre    ) then
            ind = ind + 1
            derived(ind)  = cp%Reion%redshift ! reionization redshift
        end if
        if ( fp%fisher_der%want_neff   ) then
            ind = ind + 1
            derived(ind) = cp%Num_Nu_massless + cp%Num_Nu_massive
        end if

        ! add sigma8 if wanted:
        if ( fp%fisher_der%want_sigma8 ) then
            do i=1, fp%fisher_der%FD_num_redshift
                ind = ind + 1
                j = p%Transfer%PK_redshifts_index(i)
                derived(ind) = outdata%MTrans%sigma_8(i,1)
            end do
        end if
        ! add Hubble if wanted:
        if ( fp%fisher_der%want_loghubble ) then
            do i=1, fp%fisher_der%FD_num_redshift
                ind = ind + 1
                z = fp%fisher_der%FD_redshift(i)
                derived(ind) = log10(hofz(z))
            end do
        end if
        ! add DA if wanted:
        if ( fp%fisher_der%want_logDA ) then
            do i=1, fp%fisher_der%FD_num_redshift
                ind = ind + 1
                z = fp%fisher_der%FD_redshift(i)
                if (angulardiameterdistance(z)==0._dl) then
                    derived(ind) = -16._dl ! protection against redshift zero...
                else
                    derived(ind) = log10(angulardiameterdistance(z))
                end if
            end do
        end if

        if ( num/=ind ) then
            write(*,*) 'ERROR: Compute_Derived_Parameters num is not: ', num
            write(*,*) 'If you added a parameter change also dimension_derived_parameters'
            stop
        end if

    end subroutine compute_derived_parameters

    ! ---------------------------------------------------------------------------------------------
    subroutine fisher_derived_param_names(P, FP, param_number, param_name, param_name_latex)

        implicit none

        Type(cambparams)      , intent(in)  :: P
        Type(cosmicfish_params), intent(in)  :: FP
        integer               , intent(in)  :: param_number
        character(*)          , intent(out) :: param_name
        character(*)          , intent(out), optional :: param_name_latex

        integer  :: i, j, ind, num_param
        character(len=20) str
        real(dl) :: z

        num_param = 0

        ! standard parameters:

        if ( fp%fisher_der%want_omegab          )  num_param  = num_param + 1 ! omegab
        if ( param_number == num_param          )  then; param_name = 'omegab'; if ( present(param_name_latex) ) param_name_latex = '\Omega_b'; return ;
        end if
        if ( fp%fisher_der%want_omegac          )  num_param  = num_param + 1 ! omegac
        if ( param_number == num_param          )  then; param_name = 'omegac'; if ( present(param_name_latex) ) param_name_latex = '\Omega_c'; return ;
        end if
        if ( fp%fisher_der%want_omegan          )  num_param  = num_param + 1 ! omegan
        if ( param_number == num_param          )  then; param_name = 'omeganu'; if ( present(param_name_latex) ) param_name_latex = '\Omega_{\nu}'; return ;
        end if
        if ( fp%fisher_der%want_omegav          )  num_param  = num_param + 1 ! omegav
        if ( param_number == num_param          )  then; param_name = 'omegav'; if ( present(param_name_latex) ) param_name_latex = '\Omega_{\Lambda}'; return ;
        end if
        if ( fp%fisher_der%want_omegak          )  num_param  = num_param + 1 ! omegak
        if ( param_number == num_param          )  then; param_name = 'omegak'; if ( present(param_name_latex) ) param_name_latex = '\Omega_{K}'; return ;
        end if
        if ( fp%fisher_der%want_omegam          )  num_param  = num_param + 1 ! omegam
        if ( param_number == num_param          )  then; param_name = 'omegam'; if ( present(param_name_latex) ) param_name_latex = '\Omega_{m}'; return ;
        end if
        if ( fp%fisher_der%want_theta           )  num_param  = num_param + 1 ! theta
        if ( param_number == num_param          )  then; param_name = 'theta'; if ( present(param_name_latex) ) param_name_latex = '100\theta_{MC}'; return ;
        end if
        if ( fp%fisher_der%want_mnu             )  num_param  = num_param + 1 ! mnu
        if ( param_number == num_param          )  then; param_name = 'mnu'; if ( present(param_name_latex) ) param_name_latex = '\Sigma m_\nu'; return ;
        end if
        if ( fp%fisher_der%want_zre             )  num_param  = num_param + 1 ! zre
        if ( param_number == num_param          )  then; param_name = 'z_re'; if ( present(param_name_latex) ) param_name_latex = 'z_{\rm re}'; return ;
        end if
        if ( fp%fisher_der%want_neff            )  num_param  = num_param + 1 ! neff
        if ( param_number == num_param          )  then; param_name = 'Neff'; if ( present(param_name_latex) ) param_name_latex = 'N_{\rm eff}'; return ;
        end if

        ind = num_param
        ! add sigma8 if wanted:
        if ( fp%fisher_der%want_sigma8 ) then
            do i=1, fp%fisher_der%FD_num_redshift
                ind = ind + 1
                if ( param_number == ind )  then
                    j = p%Transfer%PK_redshifts_index(i)
                    write (str, *) i
                    str = adjustl(str)
                    param_name = 'sigma8_'//trim(str)
                    write (str,'(F20.2)') p%Transfer%redshifts(j)
                    str = adjustl(str)
                    if ( present(param_name_latex) ) param_name_latex = '\sigma_{8}('//trim(str)//')'
                    return
                end if
            end do
        end if
        ! add Hubble if wanted:
        if ( fp%fisher_der%want_loghubble ) then
            do i=1, fp%fisher_der%FD_num_redshift
                ind = ind + 1
                if ( param_number == ind )  then
                    z = fp%fisher_der%FD_redshift(i)
                    write (str, *) i
                    str = adjustl(str)
                    param_name = 'logHoz_'//trim(str)
                    write (str, '(F20.2)') z
                    str = adjustl(str)
                    if ( present(param_name_latex) ) param_name_latex = '\log_{10}H('//trim(str)//')'
                    return
                end if
            end do
        end if
        ! add DA if wanted:
        if ( fp%fisher_der%want_logDA ) then
            do i=1, fp%fisher_der%FD_num_redshift
                ind = ind + 1
                if ( param_number == ind )  then
                    z = fp%fisher_der%FD_redshift(i)
                    write (str, *) i
                    str = adjustl(str)
                    param_name = 'logDAz_'//trim(str)
                    write (str, '(F20.2)') z
                    str = adjustl(str)
                    if ( present(param_name_latex) ) param_name_latex = '\log_{10}D_{\rm A}('//trim(str)//')'
                    return
                end if
            end do
        end if

    end subroutine fisher_derived_param_names

    ! ---------------------------------------------------------------------------------------------
    subroutine derived_parameter_derivative( P, FP, param_num, &
        & derived_initial, derived_derivative, derived_der_error, &
        & 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)  :: derived_initial
        real(dl), dimension(:), intent(out) :: derived_derivative
        real(dl), dimension(:), intent(out) :: derived_der_error

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

        ! definitions:
        real(dl), dimension(:), allocatable :: der_temp_1, der_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  :: derived_param_size

        ! initialize output variables:
        error_code         = 0
        derived_derivative = 0._dl
        derived_der_error  = 0._dl

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

        call dimension_derived_parameters( p, fp, derived_param_size )

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

        ! 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
        derived_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 compute_derived_parameters( p_temp, fp_temp, der_temp_1, derived_param_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 compute_derived_parameters( p_temp, fp_temp, der_temp_2, derived_param_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
                derived_der_error    = 0.0_dl
                derived_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,:) = (der_temp_1(:) - der_temp_2(:))/(2.0_dl*hh)
        else if (side == 1) then ! left side: increment negative
            a(1,1,:) = (der_temp_2(:) - derived_initial(:))/(-hh)
        else if (side == 2) then ! right side: increment positive
            a(1,1,:) = (der_temp_1(:) - derived_initial(:))/(hh)
        end if

        ! fixed stepsize derivative if adaptivity is not wanted:
        if ( .not. fp%adaptivity ) then
            derived_derivative(:) = a(1,1,:)
            derived_der_error     = 0.0_dl
            deallocate( der_temp_1, der_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(derived_param_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 compute_derived_parameters( p_temp, fp_temp, der_temp_1, derived_param_size, err )
                if ( err /= 0 ) then
                    error_code     = 4
                    derived_der_error  = 0.0_dl
                    derived_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 compute_derived_parameters( p_temp, fp_temp, der_temp_2, derived_param_size, err )
                if ( err /= 0 ) then
                    error_code     = 4
                    derived_der_error  = 0.0_dl
                    derived_derivative = 0.0_dl
                    return
                end if
                ! store the result
                a(1,i,:) = (der_temp_1(:) - der_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 compute_derived_parameters( p_temp, fp_temp, der_temp_1, derived_param_size, err )
                if ( err /= 0 ) then
                    error_code     = 4
                    derived_der_error  = 0.0_dl
                    derived_derivative = 0.0_dl
                    return
                end if
                ! store the result
                a(1,i,:) = (der_temp_1(:) - derived_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 compute_derived_parameters( p_temp, fp_temp, der_temp_1, derived_param_size, err )
                if ( err /= 0 ) then
                    error_code     = 4
                    derived_der_error  = 0.0_dl
                    derived_derivative = 0.0_dl
                    return
                end if
                ! store the result
                a(1,i,:) = (der_temp_1(:) - derived_initial(:))/(hh)
            end if

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

                forall ( k = 1:derived_param_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:derived_param_size, &
                    & errt(k) .le. derived_der_error(k) .and. accept_mask(k) == 0 )

                    derived_der_error(k)  = errt(k)
                    derived_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:derived_param_size, &
                & accept_mask(k) == 0 .and. &
                & (a(i,i,k) - a(i-1,i-1,k)) .ge. safe*derived_der_error(k) )

                accept_mask(k) = 1
            end forall

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

        end do

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

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

        return

    end subroutine derived_parameter_derivative

    ! ---------------------------------------------------------------------------------------------
    subroutine fisher_derived( P_in, FP_in, num_param, num_derived, Derived_Matrix, outroot )

        implicit none

        Type(cambparams)       , intent(in) :: P_in
        Type(cosmicfish_params) , intent(in) :: FP_in

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

        Type(cambparams)       :: P
        Type(cosmicfish_params) :: FP

        real(dl) :: time_1, time_2
        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(:)   :: derived_fiducial, derived_derivative, derived_temp
        real(dl), allocatable, dimension(:,:) :: derived_derivative_array

        ! copy input parameters:
        p  = p_in
        fp = fp_in

        ! enforce parameter compatibility:
        call derived_fisher_parameter_compatibility( p, fp )

        ! warning against 1 parameter derived:
        if ( num_derived==1 ) then
            write(*,*) 'WARNING: for better compatibility with CosmicFish PyLib use at least two derived parameters.'
        end if

        ! 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 )

        allocate( derived_fiducial(num_derived), derived_derivative(num_derived), derived_temp(num_derived), stat = allocatestatus )
        if (allocatestatus /= 0) stop "Allocation failed. Not enough memory to allocate derived fiducial."
        allocate( derived_derivative_array( num_param, num_derived ), stat = allocatestatus)
        if (allocatestatus /= 0) stop "Allocation failed. Not enough memory to allocate derived_derivative_array"

        ! compute the fiducial:
        time_1 = omp_get_wtime()
        call compute_derived_parameters( p, fp, derived_fiducial, num_derived, err )
        time_2 = omp_get_wtime() - time_1

        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 ( fp%cosmicfish_feedback >= 1 ) then
            write(*,'(a)')         '**************************************************************'
            write(*,'(a,i10)')         'Number of derived parameters : ', num_derived
            write(*,'(a,i10)')         'Number of parameters         : ', num_param
            write(*,'(a,f10.3,a)')     'Time for derived 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    : ', time_2*20._dl*num_param, ' (s)'
            else
                write(*,'(a,f10.3,a)') 'Estimated completion time    : ', time_2*2._dl*num_param, ' (s)'
            end if
            write(*,'(a)')         '**************************************************************'
        end if

        ! compute the derivatives of the Cls:
        if ( fp%cosmicfish_feedback >= 1 ) then
            write(*,'(a)') 'Computation of the derived parameters 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 derived_parameter_derivative( p, fp, ind, &
                & derived_fiducial, derived_derivative, derived_temp, &
                & initial_step, &
                & num_param, &
                & err  )

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

            derived_derivative_array(ind,:) = derived_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_derived

                derived_matrix(ind, ind2) = derived_derivative_array(ind, ind2)

            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 derived matrix:', time_2, ' (s)'
        end if

    end subroutine fisher_derived

    ! ---------------------------------------------------------------------------------------------
    subroutine save_derived_fisher_to_file( P_in, FP_in, num_params, num_derived, matrix, filename )

        implicit none

        Type(cambparams)                          :: P_in
        Type(cosmicfish_params)                    :: FP_in
        integer , intent(in)                      :: num_params
        integer , intent(in)                      :: num_derived
        real(dl), intent(in), dimension(num_params, num_derived) :: matrix
        character(len=*), intent(in)              :: filename

        Type(cambparams)        :: P
        Type(cosmicfish_params)  :: FP
        integer :: ind, err
        character(LEN=500) :: param_name
        character(LEN=500) :: param_name_latex

        real(dl), allocatable, dimension(:) :: parameters, parameters_derived

        ! copy input parameters:
        p  = p_in
        fp = fp_in
        ! enforce parameter compatibility:
        call derived_fisher_parameter_compatibility( p, fp )

        allocate( parameters(num_params), parameters_derived(num_derived) )
        call camb_params_to_params_array( p, fp, num_params, parameters )
        call compute_derived_parameters( p, fp, parameters_derived, num_derived, err )
        if ( err /= 0 ) stop 'Error: Compute_Derived_Parameters failed when save_derived_Fisher_to_file'

        open(unit=666, file=filename, action="write", status="replace")

        ! write the header:
        write(666,'(a)') '#'
        write(666,'(a)') '# This file contains a Fisher matrix created with the CosmicFish code.'
        ! write the parameter names:
        write(666,'(a)') '#'
        write(666,'(a)') '# The parameters of this derived parameters Fisher matrix are:'
        write(666,'(a)') '#'

        do ind = 1, num_params
            call fisher_param_names(p, fp, ind, param_name, param_name_latex)
            write(666,'(a,I5,a,a,a,a,a,E25.16)') '# ', ind, '    ', trim(param_name), '    ', trim(param_name_latex), '    ', parameters(ind)
        end do
        ! write the parameter names:
        write(666,'(a)') '#'
        write(666,'(a)') '# The derived parameters of this matrix are:'
        write(666,'(a)') '#'

        do ind = 1, num_derived
            call fisher_derived_param_names(p, fp, ind, param_name, param_name_latex)
            write(666,'(a,I5,a,a,a,a,a,E25.16)') '# ', ind, '    ', trim(param_name), '    ', trim(param_name_latex), '    ', parameters_derived(ind)
        end do
        write(666,'(a)') '#'
        write(666,'(a)')

        ! write the Fisher matrix:
        do ind = 1, num_params
            write(666,'(2x,*(E25.16,2x))') matrix(ind, :)
        end do

        close(666)

    end subroutine save_derived_fisher_to_file

    ! ---------------------------------------------------------------------------------------------
    subroutine save_derived_paramnames_to_file( P_in, FP_in, num_params, num_derived, filename )

        implicit none

        Type(cambparams)                          :: P_in
        Type(cosmicfish_params)                    :: FP_in
        integer , intent(in)                      :: num_params
        integer , intent(in)                      :: num_derived
        character(len=*), intent(in)              :: filename

        Type(cambparams)        :: P
        Type(cosmicfish_params)  :: FP
        integer :: ind, err
        character(LEN=500) :: param_name
        character(LEN=500) :: param_name_latex

        real(dl), allocatable, dimension(:) :: parameters, parameters_derived

        ! copy input parameters:
        p  = p_in
        fp = fp_in
        ! enforce parameter compatibility:
        call derived_fisher_parameter_compatibility( p, fp )
        allocate( parameters(num_params), parameters_derived(num_derived) )
        call camb_params_to_params_array( p, fp, num_params, parameters )

        call compute_derived_parameters( p, fp, parameters_derived, num_derived, err )

        if ( err /= 0 ) stop 'Error: Compute_Derived_Parameters failed when save_derived_Fisher_to_file'

        open(unit=666, file=filename, action="write", status="replace")

        ! write the header:

        write(666,'(a)') '#'
        write(666,'(a)') '# This file contains the parameter names for a derived Fisher matrix.'
        write(666,'(a)') '# Derived parameters are indicated with a * at the end of the name.'
        write(666,'(a)') '#'
        ! write the parameter names:
        do ind = 1, num_params
            call fisher_param_names(p, fp, ind, param_name, param_name_latex)
            write(666,'(a,a,a,a,E25.16)') trim(param_name), '    ', trim(param_name_latex), '    ', parameters(ind)
        end do
        do ind = 1, num_derived
            call fisher_derived_param_names(p, fp, ind, param_name, param_name_latex)
            write(666,'(a,a,a,a,E25.16)') trim(param_name)//'*', '    ', trim(param_name_latex), '    ', parameters_derived(ind)
        end do

        close(666)

    end subroutine save_derived_paramnames_to_file

end module fisher_calculator_derived