diff -crBN source_orig_select/driver.F90 source_interp/driver.F90 *** source_orig_select/driver.F90 2010-11-26 14:13:28.000000000 -0500 --- source_interp/driver.F90 2010-11-26 14:15:37.000000000 -0500 *************** *** 200,205 **** --- 200,208 ---- if (Ini_Read_String('use_2dF') /= '') stop 'use_2dF now replaced with use_mpk' Use_Clusters = Ini_Read_Logical('use_clusters',.false.) Use_mpk = Ini_Read_Logical('use_mpk',.false.) ! matter power spectrum, incl 2dF + !Interpolation edit: + mpk_ind = Use_mpk + !end edit Use_HST = Ini_Read_Logical('use_HST',.true.) Use_BBN = Ini_Read_Logical('use_BBN',.false.) Use_Age_Tophat_Prior= Ini_Read_Logical('use_Age_Tophat_Prior',.true.) *************** *** 273,278 **** --- 276,288 ---- numtoget = Ini_Read_Int('samples') + !Interpolation edit: + pts_mult = Ini_Read_Real('factor', 2.0) + interp_order = Ini_Read_Int('interp_order') + cut = Ini_Read_Real('cut',1000.0) + do_interp = Ini_Read_Logical('do_interp', .false.) + frac_cut = Ini_Read_Real('fraction_cut', 1.0) + !End end call Initialize(DefIni,Params) call Ini_Close diff -crBN source_orig_select/interpolation.f90 source_interp/interpolation.f90 *** source_orig_select/interpolation.f90 1969-12-31 19:00:00.000000000 -0500 --- source_interp/interpolation.f90 2010-11-26 14:14:59.000000000 -0500 *************** *** 0 **** --- 1,285 ---- + module interp_mod + use settings + implicit none + + contains + + subroutine interpolate(dataBaseNumberloc) + + use lsq, only: nobs, sserr, startup, includ, sing, tolset, regcf, ss + use CMB_Cls + use cmbtypes + + implicit none + + + !FOR !!ism + #ifdef MPI + integer ierror + integer loc_met, tot_met, tot_accept, tot_databasenumber, entries + integer, allocatable :: num_samples_ar(:), mpirankindex(:), offset_ar(:), recvcounts(:), likeoffset_ar(:) + real par_MaxLike + #endif + integer, intent(in) :: dataBaseNumberloc + real (kind = 8), allocatable :: par_paramdatabase(:,:), par_likedatabase(:), tempDataBase(:,:), par_sigma8database(:), par_r10database(:) + real (kind = 8), allocatable :: xx(:), xx1(:), r10_xx1(:) + logical :: debug_kr = .true. + integer :: i, j, k, cnt, num_samples, ifault1, ifault + real (kind = 8):: yy, s8, wt, one=1.D0 + logical, save :: fst_interp = .true. + + allocate(xx(coeff)) + allocate(xx1(coeff1)) + if (compute_tensors) allocate(r10_xx1(r10_coeff1)) + if (feedback > 1) then + write(*,*) 'Making an interpolation call' + write(logfile_unit,*)'Making an interpolation call' + end if + + ! call startup (coeff, .false.) + wt = one + + !Let's share our samples!! + #ifdef MPI + + CALL MPI_allreduce(dataBaseNumberloc, tot_databasenumber, 1, MPI_integer, MPI_sum, MPI_COMM_WORLD, ierror) + if (ierror /= MPI_SUCCESS) write(*,*) 'ERROR in MPI_reduce' + if (feedback > 1) then + write(*,*) MPIrank, 'db number: loc, tot', dataBaseNumber, tot_databasenumber + write(logfile_unit,*) MPIrank, 'db number: loc, tot', dataBaseNumber, tot_databasenumber + end if + num_samples = tot_databasenumber + + ALLOCATE(par_paramdatabase(num_params, tot_databasenumber)) + ALLOCATE(par_likedatabase(tot_databasenumber)) + ALLOCATE(par_sigma8database(tot_databasenumber)) + if (compute_tensors) ALLOCATE(par_r10database(tot_databasenumber)) + ALLOCATE(num_samples_ar(MPIchains)) + ALLOCATE(mpirankindex(MPIchains)) + ALLOCATE(offset_ar(MPIchains)) + ALLOCATE(recvcounts(MPIchains)) + ALLOCATE(likeoffset_ar(MPIchains)) + !num_samples_ar: Vector containing the number of database entries for each chain. + + !Figure out the max likelihood of the chains + CALL MPI_allreduce(notinterp_MaxLike, par_MaxLike, 1, MPI_real, MPI_MIN, MPI_COMM_WORLD, ierror) + + if (fst_interp) then + fst_interp = .false. + if (feedback > 1) write(logfile_unit,*) 'starting interpolation...' + else + cut = abs(notinterp_MaxLike - par_maxLike) + cut + write (logfile_unit,*) 'cut is:', cut + end if + + notinterp_MaxLike = par_MaxLike + + CALL MPI_allgather(dataBaseNumberloc, 1, MPI_integer, num_samples_ar, 1, MPI_integer, MPI_COMM_WORLD, ierror) + + !Vector containing the offsets: + offset_ar(1) = 0 + do i = 2, MPIchains + offset_ar(i) = offset_ar(i-1) + num_params*num_samples_ar(i-1) + end do + + likeoffset_ar(1) = 0 + do i = 2, MPIchains + likeoffset_ar(i) = likeoffset_ar(i-1) + num_samples_ar(i-1) + end do + + recvcounts(1:MPIchains) = num_params*num_samples_ar(1:MPIchains) + + CALL MPI_allgatherv(paramDataBase, (num_params*dataBaseNumberloc), MPI_double_precision,& + par_paramdatabase, recvcounts, offset_ar, MPI_double_precision, MPI_COMM_WORLD, ierror) + CALL MPI_allgatherv(likeDataBase, dataBaseNumberloc, MPI_double_precision, par_likedatabase,& + num_samples_ar, likeoffset_ar, MPI_double_precision, MPI_COMM_WORLD, ierror) + + if (mpk_ind) CALL MPI_allgatherv(sigma8database, dataBaseNumberloc, MPI_double_precision,& + par_sigma8database, num_samples_ar, likeoffset_ar, MPI_double_precision, MPI_COMM_WORLD, ierror) + + if (compute_tensors) CALL MPI_allgatherv(r10database, dataBaseNumberloc, MPI_double_precision,& + par_r10database, num_samples_ar, likeoffset_ar, MPI_double_precision, MPI_COMM_WORLD, ierror) + #else + + num_samples = dataBaseNumberloc + + allocate(par_paramDataBase(num_params, num_samples)) + allocate(par_likedatabase(num_samples)) + allocate(par_sigma8database(num_samples)) + if (compute_tensors) allocate(par_r10database(num_samples)) + + do i = 1, num_samples + par_paramDataBase(1:num_params, i) = paramDataBase(1:num_params,i) + par_likeDataBase(i) = likeDataBase(i) + end do + + #endif + + allocate(tempDataBase(num_params,num_samples)) + + !Figuring out the parameter means: + mean_ar(1:num_params) = 0 + do i = 1, num_samples + mean_ar(1:num_params) = mean_ar(1:num_params) + par_paramDataBase(1:num_params,i) + end do + mean_ar(1:num_params) = mean_ar(1:num_params)/num_samples + + !Figuring out the standard deviations: + sd_ar(1:num_params) = 0 + do i = 1, num_samples + sd_ar(1:num_params) = sd_ar(1:num_params) + (par_paramDataBase(1:num_params,i) - mean_ar(1:num_params))**2 + end do + sd_ar(1:num_params) = sqrt(sd_ar(1:num_params)/(num_samples)) + + !Normalizing Parameters: + do i = 1, num_samples + do j= 1, num_params + if (sd_ar(j) /= 0) tempDataBase(j,i) =(par_paramDataBase(j,i) - mean_ar(j))/sd_ar(j) + end do + end do + + !Running interpolation + cnt = 0 + + call startup (coeff, .false.) + do i = 1, num_samples + if ( ABS(par_LikeDataBase(I) - notinterp_MaxLike) < cut) then + yy = par_likeDataBase(i) + xx(1:coeff) = 1 + do k = 1, coeff + do j = 1, num_params_used + xx(k)=xx(k)*tempDatabase(params_used(j),i)**emat(k,j) + end do + end do + call includ(wt, xx, yy) + end if + end do + call tolset() + call regcf(PAR,coeff,ifault) + call ss() + if (feedback > 0) write(logfile_unit,*) 'RMSE 1:', sqrt(sserr/(nobs - coeff)) + + call startup (coeff1, .false.) + do i = 1, num_samples + if ( ABS(par_LikeDataBase(I) - notinterp_MaxLike) < cut) then + yy = par_likeDataBase(i) + xx1(1:coeff1) = 1 + do k = 1, coeff1 + do j = 1, num_params_used + xx1(k)=xx1(k)*tempDatabase(params_used(j),i)**emat1(k,j) + end do + end do + call includ(wt, xx1, yy) + end if + end do + call tolset() + call regcf(PAR1,coeff1,ifault) + call ss() + if (feedback > 0) write(logfile_unit,*) 'RMSE 2:', sqrt(sserr/(nobs - coeff1)) + + if (mpk_ind) then + call startup (coeff1, .false.) + do i = 1, num_samples + if ( ABS(par_LikeDataBase(I) - notinterp_MaxLike) < cut) then + yy = par_sigma8database(i) + xx1(1:coeff1) = 1 + do k = 1, coeff1 + do j = 1, num_params_used + xx1(k)=xx1(k)*tempDatabase(params_used(j),i)**emat1(k,j) + end do + end do + call includ(wt, xx1, yy) + end if + end do + call tolset() + call regcf(PAR_s8,coeff1,ifault) + call ss() + if (feedback > 0) write(logfile_unit,*) 'RMSE s8:', sqrt(sserr/(nobs - coeff1)) + endif + + if (compute_tensors) then + call startup (r10_coeff1, .false.) + do i = 1, num_samples + if ( ABS(par_LikeDataBase(I) - notinterp_MaxLike) < cut) then + yy = par_r10database(i) + r10_xx1(1:r10_coeff1) = 1 + do k = 1, r10_coeff1 + do j = 1, num_params_used + r10_xx1(k)=r10_xx1(k)*tempDatabase(params_used(j),i)**r10_emat1(k,j) + end do + end do + call includ(wt, r10_xx1, yy) + end if + end do + call tolset() + call regcf(PAR_r10,r10_coeff1,ifault) + call ss() + if (feedback > 0) write(logfile_unit,*) 'RMSE r10:', sqrt(sserr/(nobs - r10_coeff1)) + endif + + #ifdef MPI + DEALLOCATE(num_samples_ar) + DEALLOCATE(mpirankindex) + DEALLOCATE(offset_ar) + DEALLOCATE(recvcounts) + DEALLOCATE(likeoffset_ar) + #endif + DEALLOCATE(par_paramdatabase, par_likedatabase, par_sigma8database) + if (compute_tensors) deallocate(par_r10database) + + + if (feedback > 0) then + write(logfile_unit, *) 'Peak Likelihood (noninterp)', notinterp_Maxlike + write(logfile_unit,*)'No. observations', nobs + end if + if (nobs < coeff1) write(logfile_unit,*) 'WARNING: not enough samples for interp 1.' + if (nobs == 0 ) then + write(logfile_unit) 'ERROR: DISCONTINUING USE OF INTERPOLATION' + do_interp = .false. + interp_ind = .false. + end if + + DEALLOCATE(xx, xx1) + if (compute_tensors) deallocate(r10_xx1) + call FlushFile(logfile_unit) + + end subroutine interpolate + + !Subroutine sets the arrays for the polynomial's exponents. + !Order is set in settings.f90 or .ini file. + recursive subroutine make_array(level, array, index_ar, order, ar_index, nparams) + + integer, intent(in) :: level, order, nparams + integer, intent(inout) :: array(:,:), index_ar(:) + integer, intent(inout) :: ar_index + + do while (index_ar(level) <= order) + if (level == 1) then + if (sum(index_ar) <= order) then + array(ar_index, 1:nparams)= index_ar(1:nparams) + !array(ar_index, 1:num_params_used)= index_ar(1:num_params_used) + ar_index = ar_index + 1 + index_ar(level) = index_ar(level) + 1 + else + index_ar(level) = index_ar(level) + 1 + endif + else + index_ar(level -1) = 0 + call make_array(level-1, array, index_ar, order, ar_index, nparams) + !call make_array(level-1, array, index_ar, order, ar_index) + index_ar(level) = index_ar(level) + 1 + end if + end do + + end subroutine make_array + + recursive function factorial(n) result(res) + integer*8 :: res, n + if (n==1) then + res = 1 + else + res = n * factorial(n-1) + end if + end function factorial + + end module interp_mod diff -crBN source_orig_select/lsq.f90 source_interp/lsq.f90 *** source_orig_select/lsq.f90 1969-12-31 19:00:00.000000000 -0500 --- source_interp/lsq.f90 2010-11-26 14:14:51.000000000 -0500 *************** *** 0 **** --- 1,994 ---- + MODULE lsq + + ! Module for unconstrained linear least-squares calculations. + ! The algorithm is suitable for updating LS calculations as more + ! data are added. This is sometimes called recursive estimation. + ! Only one dependent variable is allowed. + ! Based upon Applied Statistics algorithm AS 274. + ! Translation from Fortran 77 to Fortran 90 by Alan Miller. + ! A function, VARPRD, has been added for calculating the variances + ! of predicted values, and this uses a subroutine BKSUB2. + + ! Version 1.14, 19 August 2002 - ELF90 compatible version + ! Author: Alan Miller + ! e-mail : amiller @ bigpond.net.au + ! WWW-pages: http://www.ozemail.com.au/~milleraj + ! http://users.bigpond.net.au/amiller/ + + ! Bug fixes: + ! 1. In REGCF a call to TOLSET has been added in case the user had + ! not set tolerances. + ! 2. In SING, each time a singularity is detected, unless it is in the + ! variables in the last position, INCLUD is called. INCLUD assumes + ! that a new observation is being added and increments the number of + ! cases, NOBS. The line: nobs = nobs - 1 has been added. + ! 3. row_ptr was left out of the DEALLOCATE statement in routine startup + ! in version 1.07. + ! 4. In COV, now calls SS if rss_set = .FALSE. 29 August 1997 + ! 5. In TOLSET, correction to accomodate negative values of D. 19 August 2002 + + ! Other changes: + ! 1. Array row_ptr added 18 July 1997. This points to the first element + ! stored in each row thus saving a small amount of time needed to + ! calculate its position. + ! 2. Optional parameter, EPS, added to routine TOLSET, so that the user + ! can specify the accuracy of the input data. + ! 3. Cosmetic change of lsq_kind to dp (`Double precision') + ! 4. Change to routine SING to use row_ptr rather than calculate the position + ! of first elements in each row. + + ! The PUBLIC variables are: + ! dp = a KIND parameter for the floating-point quantities calculated + ! in this module. See the more detailed explanation below. + ! This KIND parameter should be used for all floating-point + ! arguments passed to routines in this module. + + ! nobs = the number of observations processed to date. + ! ncol = the total number of variables, including one for the constant, + ! if a constant is being fitted. + ! r_dim = the dimension of array r = ncol*(ncol-1)/2 + ! vorder = an integer vector storing the current order of the variables + ! in the QR-factorization. The initial order is 0, 1, 2, ... + ! if a constant is being fitted, or 1, 2, ... otherwise. + ! initialized = a logical variable which indicates whether space has + ! been allocated for various arrays. + ! tol_set = a logical variable which is set when subroutine TOLSET has + ! been called to calculate tolerances for use in testing for + ! singularities. + ! rss_set = a logical variable indicating whether residual sums of squares + ! are available and usable. + ! d() = array of row multipliers for the Cholesky factorization. + ! The factorization is X = Q.sqrt(D).R where Q is an ortho- + ! normal matrix which is NOT stored, D is a diagonal matrix + ! whose diagonal elements are stored in array d, and R is an + ! upper-triangular matrix with 1's as its diagonal elements. + ! rhs() = vector of RHS projections (after scaling by sqrt(D)). + ! Thus Q'y = sqrt(D).rhs + ! r() = the upper-triangular matrix R. The upper triangle only, + ! excluding the implicit 1's on the diagonal, are stored by + ! rows. + ! tol() = array of tolerances used in testing for singularities. + ! rss() = array of residual sums of squares. rss(i) is the residual + ! sum of squares with the first i variables in the model. + ! By changing the order of variables, the residual sums of + ! squares can be found for all possible subsets of the variables. + ! The residual sum of squares with NO variables in the model, + ! that is the total sum of squares of the y-values, can be + ! calculated as rss(1) + d(1)*rhs(1)^2. If the first variable + ! is a constant, then rss(1) is the sum of squares of + ! (y - ybar) where ybar is the average value of y. + ! sserr = residual sum of squares with all of the variables included. + ! row_ptr() = array of indices of first elements in each row of R. + ! + !-------------------------------------------------------------------------- + + ! General declarations + + IMPLICIT NONE + + INTEGER, SAVE :: nobs, ncol, r_dim + INTEGER, ALLOCATABLE, SAVE :: vorder(:), row_ptr(:) + LOGICAL, SAVE :: initialized = .false., & + tol_set = .false., rss_set = .false. + + ! Note. dp is being set to give at least 12 decimal digit + ! representation of floating point numbers. This should be adequate + ! for most problems except the fitting of polynomials. dp is + ! being set so that the same code can be run on PCs and Unix systems, + ! which will usually represent floating-point numbers in `double + ! precision', and other systems with larger word lengths which will + ! give similar accuracy in `single precision'. + + INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12,60) + REAL (dp), ALLOCATABLE, SAVE :: d(:), rhs(:), r(:), tol(:), rss(:) + REAL (dp), SAVE :: zero = 0.0_dp, one = 1.0_dp, vsmall + REAL (dp), SAVE :: sserr, toly + + !PUBLIC :: dp, nobs, ncol, r_dim, vorder, row_ptr, & + ! initialized, tol_set, rss_set, & + ! r, d, rhs, tol, rss, sserr + !PUBLIC :: sserr, rss, nobs !KR EDIT + PRIVATE :: zero, one, vsmall + + + CONTAINS + + SUBROUTINE startup(nvar, fit_const) + + ! Allocates dimensions for arrays and initializes to zero + ! The calling program must set nvar = the number of variables, and + ! fit_const = .true. if a constant is to be included in the model, + ! otherwise fit_const = .false. + ! + !-------------------------------------------------------------------------- + + IMPLICIT NONE + INTEGER, INTENT(IN) :: nvar + LOGICAL, INTENT(IN) :: fit_const + + ! Local variable + INTEGER :: i + + vsmall = 10. * TINY(zero) + + + IF (fit_const) THEN + ncol = nvar + 1 + ELSE + ncol = nvar + END IF + + IF (initialized) DEALLOCATE(d, rhs, r, tol, rss, vorder, row_ptr) + r_dim = ncol * (ncol - 1)/2 + ALLOCATE( d(ncol), rhs(ncol), r(r_dim), tol(ncol), rss(ncol), vorder(ncol), & + row_ptr(ncol) ) + + d = zero + rhs = zero + r = zero + sserr = zero + nobs = zero !KR EDIT + + IF (fit_const) THEN + DO i = 1, ncol + vorder(i) = i-1 + END DO + ELSE + DO i = 1, ncol + vorder(i) = i + END DO + END IF ! (fit_const) + + ! row_ptr(i) is the position of element R(i,i+1) in array r(). + + row_ptr(1) = 1 + DO i = 2, ncol-1 + row_ptr(i) = row_ptr(i-1) + ncol - i + 1 + END DO + row_ptr(ncol) = 0 + + initialized = .true. + tol_set = .false. + rss_set = .false. + + RETURN + END SUBROUTINE startup + + + + + SUBROUTINE includ(weight, xrow, yelem) + + ! ALGORITHM AS75.1 APPL. STATIST. (1974) VOL.23, NO. 3 + + ! Calling this routine updates D, R, RHS and SSERR by the + ! inclusion of xrow, yelem with the specified weight. + + ! *** WARNING Array XROW is overwritten. + + ! N.B. As this routine will be called many times in most applications, + ! checks have been eliminated. + ! + !-------------------------------------------------------------------------- + + + IMPLICIT NONE + REAL (dp),INTENT(IN) :: weight, yelem + REAL (dp), DIMENSION(:), INTENT(IN OUT) :: xrow + + ! Local variables + + INTEGER :: i, k, nextr + REAL (dp) :: w, y, xi, di, wxi, dpi, cbar, sbar, xk + + nobs = nobs + 1 + w = weight + y = yelem + rss_set = .false. + nextr = 1 + DO i = 1, ncol + + ! Skip unnecessary transformations. Test on exact zeroes must be + ! used or stability can be destroyed. + + IF (ABS(w) < vsmall) RETURN + xi = xrow(i) + IF (ABS(xi) < vsmall) THEN + nextr = nextr + ncol - i + ELSE + di = d(i) + wxi = w * xi + dpi = di + wxi*xi + cbar = di / dpi + sbar = wxi / dpi + w = cbar * w + d(i) = dpi + DO k = i+1, ncol + xk = xrow(k) + xrow(k) = xk - xi * r(nextr) + r(nextr) = cbar * r(nextr) + sbar * xk + nextr = nextr + 1 + END DO + xk = y + y = xk - xi * rhs(i) + rhs(i) = cbar * rhs(i) + sbar * xk + END IF + END DO ! i = 1, ncol + + ! Y * SQRT(W) is now equal to the Brown, Durbin & Evans recursive + ! residual. + + sserr = sserr + w * y * y + + RETURN + END SUBROUTINE includ + + + + SUBROUTINE regcf(beta, nreq, ifault) + + ! ALGORITHM AS274 APPL. STATIST. (1992) VOL.41, NO. 2 + + ! Modified version of AS75.4 to calculate regression coefficients + ! for the first NREQ variables, given an orthogonal reduction from + ! AS75.1. + ! + !-------------------------------------------------------------------------- + + IMPLICIT NONE + INTEGER, INTENT(IN) :: nreq + INTEGER, INTENT(OUT) :: ifault + REAL (dp), DIMENSION(:), INTENT(OUT) :: beta + + ! Local variables + + INTEGER :: i, j, nextr + + ! Some checks. + + ifault = 0 + IF (nreq < 1 .OR. nreq > ncol) ifault = ifault + 4 + IF (ifault /= 0) RETURN + + IF (.NOT. tol_set) CALL tolset() + + DO i = nreq, 1, -1 + IF (SQRT(d(i)) < tol(i)) THEN + beta(i) = zero + d(i) = zero + ifault = -i + ELSE + beta(i) = rhs(i) + nextr = row_ptr(i) + DO j = i+1, nreq + beta(i) = beta(i) - r(nextr) * beta(j) + nextr = nextr + 1 + END DO ! j = i+1, nreq + END IF + END DO ! i = nreq, 1, -1 + + RETURN + END SUBROUTINE regcf + + + + SUBROUTINE tolset(eps) + + ! ALGORITHM AS274 APPL. STATIST. (1992) VOL.41, NO. 2 + + ! Sets up array TOL for testing for zeroes in an orthogonal + ! reduction formed using AS75.1. + + REAL (dp), INTENT(IN), OPTIONAL :: eps + + ! Unless the argument eps is set, it is assumed that the input data are + ! recorded to full machine accuracy. This is often not the case. + ! If, for instance, the data are recorded to `single precision' of about + ! 6-7 significant decimal digits, then singularities will not be detected. + ! It is suggested that in this case eps should be set equal to + ! 10.0 * EPSILON(1.0) + ! If the data are recorded to say 4 significant decimals, then eps should + ! be set to 1.0E-03 + ! The above comments apply to the predictor variables, not to the + ! dependent variable. + + ! Correction - 19 August 2002 + ! When negative weights are used, it is possible for an alement of D + ! to be negative. + + ! Local variables. + ! + !-------------------------------------------------------------------------- + + ! Local variables + + INTEGER :: col, row, pos + REAL (dp) :: eps1, ten = 10.0, total, work(ncol) + + ! EPS is a machine-dependent constant. + + IF (PRESENT(eps)) THEN + eps1 = MAX(ABS(eps), ten * EPSILON(ten)) + ELSE + eps1 = ten * EPSILON(ten) + END IF + + ! Set tol(i) = sum of absolute values in column I of R after + ! scaling each element by the square root of its row multiplier, + ! multiplied by EPS1. + + work = SQRT(ABS(d)) + DO col = 1, ncol + pos = col - 1 + total = work(col) + DO row = 1, col-1 + total = total + ABS(r(pos)) * work(row) + pos = pos + ncol - row - 1 + END DO + tol(col) = eps1 * total + END DO + + tol_set = .TRUE. + RETURN + END SUBROUTINE tolset + + + + + SUBROUTINE sing(lindep, ifault) + + ! ALGORITHM AS274 APPL. STATIST. (1992) VOL.41, NO. 2 + + ! Checks for singularities, reports, and adjusts orthogonal + ! reductions produced by AS75.1. + + ! Correction - 19 August 2002 + ! When negative weights are used, it is possible for an alement of D + ! to be negative. + + ! Auxiliary routines called: INCLUD, TOLSET + ! + !-------------------------------------------------------------------------- + + INTEGER, INTENT(OUT) :: ifault + LOGICAL, DIMENSION(:), INTENT(OUT) :: lindep + + ! Local variables + + REAL (dp) :: temp, x(ncol), work(ncol), y, weight + INTEGER :: pos, row, pos2 + + ifault = 0 + + work = SQRT(ABS(d)) + IF (.NOT. tol_set) CALL tolset() + + DO row = 1, ncol + temp = tol(row) + pos = row_ptr(row) ! pos = location of first element in row + + ! If diagonal element is near zero, set it to zero, set appropriate + ! element of LINDEP, and use INCLUD to augment the projections in + ! the lower rows of the orthogonalization. + + lindep(row) = .FALSE. + IF (work(row) <= temp) THEN + lindep(row) = .TRUE. + ifault = ifault - 1 + IF (row < ncol) THEN + pos2 = pos + ncol - row - 1 + x = zero + x(row+1:ncol) = r(pos:pos2) + y = rhs(row) + weight = d(row) + r(pos:pos2) = zero + d(row) = zero + rhs(row) = zero + CALL includ(weight, x, y) + ! INCLUD automatically increases the number + ! of cases each time it is called. + nobs = nobs - 1 + ELSE + sserr = sserr + d(row) * rhs(row)**2 + END IF ! (row < ncol) + END IF ! (work(row) <= temp) + END DO ! row = 1, ncol + + RETURN + END SUBROUTINE sing + + + + SUBROUTINE ss() + + ! ALGORITHM AS274 APPL. STATIST. (1992) VOL.41, NO. 2 + + ! Calculates partial residual sums of squares from an orthogonal + ! reduction from AS75.1. + ! + !-------------------------------------------------------------------------- + + ! Local variables + + INTEGER :: i + REAL (dp) :: total + + total = sserr + rss(ncol) = sserr + DO i = ncol, 2, -1 + total = total + d(i) * rhs(i)**2 + rss(i-1) = total + END DO + + rss_set = .TRUE. + RETURN + END SUBROUTINE ss + + + + SUBROUTINE cov(nreq, var, covmat, dimcov, sterr, ifault) + + ! ALGORITHM AS274 APPL. STATIST. (1992) VOL.41, NO. 2 + + ! Calculate covariance matrix for regression coefficients for the + ! first nreq variables, from an orthogonal reduction produced from + ! AS75.1. + + ! Auxiliary routine called: INV + ! + !-------------------------------------------------------------------------- + + INTEGER, INTENT(IN) :: nreq, dimcov + INTEGER, INTENT(OUT) :: ifault + REAL (dp), INTENT(OUT) :: var + REAL (dp), DIMENSION(:), INTENT(OUT) :: covmat, sterr + + ! Local variables. + + INTEGER :: dim_rinv, pos, row, start, pos2, col, pos1, k + REAL (dp) :: total + REAL (dp), ALLOCATABLE :: rinv(:) + + ! Check that dimension of array covmat is adequate. + + IF (dimcov < nreq*(nreq+1)/2) THEN + ifault = 1 + RETURN + END IF + + ! Check for small or zero multipliers on the diagonal. + + ifault = 0 + DO row = 1, nreq + IF (ABS(d(row)) < vsmall) ifault = -row + END DO + IF (ifault /= 0) RETURN + + ! Calculate estimate of the residual variance. + + IF (nobs > nreq) THEN + IF (.NOT. rss_set) CALL ss() + var = rss(nreq) / (nobs - nreq) + ELSE + ifault = 2 + RETURN + END IF + + dim_rinv = nreq*(nreq-1)/2 + ALLOCATE ( rinv(dim_rinv) ) + + CALL INV(nreq, rinv) + pos = 1 + start = 1 + DO row = 1, nreq + pos2 = start + DO col = row, nreq + pos1 = start + col - row + IF (row == col) THEN + total = one / d(col) + ELSE + total = rinv(pos1-1) / d(col) + END IF + DO K = col+1, nreq + total = total + rinv(pos1) * rinv(pos2) / d(k) + pos1 = pos1 + 1 + pos2 = pos2 + 1 + END DO ! K = col+1, nreq + covmat(pos) = total * var + IF (row == col) sterr(row) = SQRT(covmat(pos)) + pos = pos + 1 + END DO ! col = row, nreq + start = start + nreq - row + END DO ! row = 1, nreq + + DEALLOCATE(rinv) + RETURN + END SUBROUTINE cov + + + + SUBROUTINE inv(nreq, rinv) + + ! ALGORITHM AS274 APPL. STATIST. (1992) VOL.41, NO. 2 + + ! Invert first nreq rows and columns of Cholesky factorization + ! produced by AS 75.1. + ! + !-------------------------------------------------------------------------- + + INTEGER, INTENT(IN) :: nreq + REAL (dp), DIMENSION(:), INTENT(OUT) :: rinv + + ! Local variables. + + INTEGER :: pos, row, col, start, k, pos1, pos2 + REAL (dp) :: total + + ! Invert R ignoring row multipliers, from the bottom up. + + pos = nreq * (nreq-1)/2 + DO row = nreq-1, 1, -1 + start = row_ptr(row) + DO col = nreq, row+1, -1 + pos1 = start + pos2 = pos + total = zero + DO k = row+1, col-1 + pos2 = pos2 + nreq - k + total = total - r(pos1) * rinv(pos2) + pos1 = pos1 + 1 + END DO ! k = row+1, col-1 + rinv(pos) = total - r(pos1) + pos = pos - 1 + END DO ! col = nreq, row+1, -1 + END DO ! row = nreq-1, 1, -1 + + RETURN + END SUBROUTINE inv + + + + SUBROUTINE partial_corr(in, cormat, dimc, ycorr, ifault) + + ! Replaces subroutines PCORR and COR of: + ! ALGORITHM AS274 APPL. STATIST. (1992) VOL.41, NO. 2 + + ! Calculate partial correlations after the variables in rows + ! 1, 2, ..., IN have been forced into the regression. + ! If IN = 1, and the first row of R represents a constant in the + ! model, then the usual simple correlations are returned. + + ! If IN = 0, the value returned in array CORMAT for the correlation + ! of variables Xi & Xj is: + ! sum ( Xi.Xj ) / Sqrt ( sum (Xi^2) . sum (Xj^2) ) + + ! On return, array CORMAT contains the upper triangle of the matrix of + ! partial correlations stored by rows, excluding the 1's on the diagonal. + ! e.g. if IN = 2, the consecutive elements returned are: + ! (3,4) (3,5) ... (3,ncol), (4,5) (4,6) ... (4,ncol), etc. + ! Array YCORR stores the partial correlations with the Y-variable + ! starting with YCORR(IN+1) = partial correlation with the variable in + ! position (IN+1). + ! + !-------------------------------------------------------------------------- + + INTEGER, INTENT(IN) :: in, dimc + INTEGER, INTENT(OUT) :: ifault + REAL (dp), DIMENSION(:), INTENT(OUT) :: cormat, ycorr + + ! Local variables. + + INTEGER :: base_pos, pos, row, col, col1, col2, pos1, pos2 + REAL (dp) :: rms(in+1:ncol), sumxx, sumxy, sumyy, work(in+1:ncol) + + ! Some checks. + + ifault = 0 + IF (in < 0 .OR. in > ncol-1) ifault = ifault + 4 + IF (dimc < (ncol-in)*(ncol-in-1)/2) ifault = ifault + 8 + IF (ifault /= 0) RETURN + + ! Base position for calculating positions of elements in row (IN+1) of R. + + base_pos = in*ncol - (in+1)*(in+2)/2 + + ! Calculate 1/RMS of elements in columns from IN to (ncol-1). + + IF (d(in+1) > zero) rms(in+1) = one / SQRT(d(in+1)) + DO col = in+2, ncol + pos = base_pos + col + sumxx = d(col) + DO row = in+1, col-1 + sumxx = sumxx + d(row) * r(pos)**2 + pos = pos + ncol - row - 1 + END DO ! row = in+1, col-1 + IF (sumxx > zero) THEN + rms(col) = one / SQRT(sumxx) + ELSE + rms(col) = zero + ifault = -col + END IF ! (sumxx > zero) + END DO ! col = in+1, ncol-1 + + ! Calculate 1/RMS for the Y-variable + + sumyy = sserr + DO row = in+1, ncol + sumyy = sumyy + d(row) * rhs(row)**2 + END DO ! row = in+1, ncol + IF (sumyy > zero) sumyy = one / SQRT(sumyy) + + ! Calculate sums of cross-products. + ! These are obtained by taking dot products of pairs of columns of R, + ! but with the product for each row multiplied by the row multiplier + ! in array D. + + pos = 1 + DO col1 = in+1, ncol + sumxy = zero + work(col1+1:ncol) = zero + pos1 = base_pos + col1 + DO row = in+1, col1-1 + pos2 = pos1 + 1 + DO col2 = col1+1, ncol + work(col2) = work(col2) + d(row) * r(pos1) * r(pos2) + pos2 = pos2 + 1 + END DO ! col2 = col1+1, ncol + sumxy = sumxy + d(row) * r(pos1) * rhs(row) + pos1 = pos1 + ncol - row - 1 + END DO ! row = in+1, col1-1 + + ! Row COL1 has an implicit 1 as its first element (in column COL1) + + pos2 = pos1 + 1 + DO col2 = col1+1, ncol + work(col2) = work(col2) + d(col1) * r(pos2) + pos2 = pos2 + 1 + cormat(pos) = work(col2) * rms(col1) * rms(col2) + pos = pos + 1 + END DO ! col2 = col1+1, ncol + sumxy = sumxy + d(col1) * rhs(col1) + ycorr(col1) = sumxy * rms(col1) * sumyy + END DO ! col1 = in+1, ncol-1 + + ycorr(1:in) = zero + + RETURN + END SUBROUTINE partial_corr + + + + + SUBROUTINE vmove(from, to, ifault) + + ! ALGORITHM AS274 APPL. STATIST. (1992) VOL.41, NO. 2 + + ! Move variable from position FROM to position TO in an + ! orthogonal reduction produced by AS75.1. + ! + !-------------------------------------------------------------------------- + + INTEGER, INTENT(IN) :: from, to + INTEGER, INTENT(OUT) :: ifault + + ! Local variables + + REAL (dp) :: d1, d2, x, d1new, d2new, cbar, sbar, y + INTEGER :: m, first, last, inc, m1, m2, mp1, col, pos, row + + ! Check input parameters + + ifault = 0 + IF (from < 1 .OR. from > ncol) ifault = ifault + 4 + IF (to < 1 .OR. to > ncol) ifault = ifault + 8 + IF (ifault /= 0) RETURN + + IF (from == to) RETURN + + IF (.NOT. rss_set) CALL ss() + + IF (from < to) THEN + first = from + last = to - 1 + inc = 1 + ELSE + first = from - 1 + last = to + inc = -1 + END IF + + DO m = first, last, inc + + ! Find addresses of first elements of R in rows M and (M+1). + + m1 = row_ptr(m) + m2 = row_ptr(m+1) + mp1 = m + 1 + d1 = d(m) + d2 = d(mp1) + + ! Special cases. + + IF (d1 < vsmall .AND. d2 < vsmall) GO TO 40 + x = r(m1) + IF (ABS(x) * SQRT(d1) < tol(mp1)) THEN + x = zero + END IF + IF (d1 < vsmall .OR. ABS(x) < vsmall) THEN + d(m) = d2 + d(mp1) = d1 + r(m1) = zero + DO col = m+2, ncol + m1 = m1 + 1 + x = r(m1) + r(m1) = r(m2) + r(m2) = x + m2 = m2 + 1 + END DO ! col = m+2, ncol + x = rhs(m) + rhs(m) = rhs(mp1) + rhs(mp1) = x + GO TO 40 + ELSE IF (d2 < vsmall) THEN + d(m) = d1 * x**2 + r(m1) = one / x + r(m1+1:m1+ncol-m-1) = r(m1+1:m1+ncol-m-1) / x + rhs(m) = rhs(m) / x + GO TO 40 + END IF + + ! Planar rotation in regular case. + + d1new = d2 + d1*x**2 + cbar = d2 / d1new + sbar = x * d1 / d1new + d2new = d1 * cbar + d(m) = d1new + d(mp1) = d2new + r(m1) = sbar + DO col = m+2, ncol + m1 = m1 + 1 + y = r(m1) + r(m1) = cbar*r(m2) + sbar*y + r(m2) = y - x*r(m2) + m2 = m2 + 1 + END DO ! col = m+2, ncol + y = rhs(m) + rhs(m) = cbar*rhs(mp1) + sbar*y + rhs(mp1) = y - x*rhs(mp1) + + ! Swap columns M and (M+1) down to row (M-1). + + 40 pos = m + DO row = 1, m-1 + x = r(pos) + r(pos) = r(pos-1) + r(pos-1) = x + pos = pos + ncol - row - 1 + END DO ! row = 1, m-1 + + ! Adjust variable order (VORDER), the tolerances (TOL) and + ! the vector of residual sums of squares (RSS). + + m1 = vorder(m) + vorder(m) = vorder(mp1) + vorder(mp1) = m1 + x = tol(m) + tol(m) = tol(mp1) + tol(mp1) = x + rss(m) = rss(mp1) + d(mp1) * rhs(mp1)**2 + END DO + + RETURN + END SUBROUTINE vmove + + + + SUBROUTINE reordr(list, n, pos1, ifault) + + ! ALGORITHM AS274 APPL. STATIST. (1992) VOL.41, NO. 2 + + ! Re-order the variables in an orthogonal reduction produced by + ! AS75.1 so that the N variables in LIST start at position POS1, + ! though will not necessarily be in the same order as in LIST. + ! Any variables in VORDER before position POS1 are not moved. + + ! Auxiliary routine called: VMOVE + ! + !-------------------------------------------------------------------------- + + INTEGER, INTENT(IN) :: n, pos1 + INTEGER, DIMENSION(:), INTENT(IN) :: list + INTEGER, INTENT(OUT) :: ifault + + ! Local variables. + + INTEGER :: next, i, l, j + + ! Check N. + + ifault = 0 + IF (n < 1 .OR. n > ncol+1-pos1) ifault = ifault + 4 + IF (ifault /= 0) RETURN + + ! Work through VORDER finding variables which are in LIST. + + next = pos1 + i = pos1 + 10 l = vorder(i) + DO j = 1, n + IF (l == list(j)) GO TO 40 + END DO + 30 i = i + 1 + IF (i <= ncol) GO TO 10 + + ! If this point is reached, one or more variables in LIST has not + ! been found. + + ifault = 8 + RETURN + + ! Variable L is in LIST; move it up to position NEXT if it is not + ! already there. + + 40 IF (i > next) CALL vmove(i, next, ifault) + next = next + 1 + IF (next < n+pos1) GO TO 30 + + RETURN + END SUBROUTINE reordr + + + + SUBROUTINE hdiag(xrow, nreq, hii, ifault) + + ! ALGORITHM AS274 APPL. STATIST. (1992) VOL.41, NO. 2 + ! + ! -1 -1 + ! The hat matrix H = x(X'X) x' = x(R'DR) x' = z'Dz + + ! -1 + ! where z = x'R + + ! Here we only calculate the diagonal element hii corresponding to one + ! row (xrow). The variance of the i-th least-squares residual is (1 - hii). + !-------------------------------------------------------------------------- + + INTEGER, INTENT(IN) :: nreq + INTEGER, INTENT(OUT) :: ifault + REAL (dp), DIMENSION(:), INTENT(IN) :: xrow + REAL (dp), INTENT(OUT) :: hii + + ! Local variables + + INTEGER :: col, row, pos + REAL (dp) :: total, wk(ncol) + + ! Some checks + + ifault = 0 + IF (nreq > ncol) ifault = ifault + 4 + IF (ifault /= 0) RETURN + + ! The elements of xrow.inv(R).sqrt(D) are calculated and stored in WK. + + hii = zero + DO col = 1, nreq + IF (SQRT(d(col)) <= tol(col)) THEN + wk(col) = zero + ELSE + pos = col - 1 + total = xrow(col) + DO row = 1, col-1 + total = total - wk(row)*r(pos) + pos = pos + ncol - row - 1 + END DO ! row = 1, col-1 + wk(col) = total + hii = hii + total**2 / d(col) + END IF + END DO ! col = 1, nreq + + RETURN + END SUBROUTINE hdiag + + + + FUNCTION varprd(x, nreq) RESULT(fn_val) + + ! Calculate the variance of x'b where b consists of the first nreq + ! least-squares regression coefficients. + ! + !-------------------------------------------------------------------------- + + INTEGER, INTENT(IN) :: nreq + REAL (dp), DIMENSION(:), INTENT(IN) :: x + REAL (dp) :: fn_val + + ! Local variables + + INTEGER :: ifault, row + REAL (dp) :: var, wk(nreq) + + ! Check input parameter values + + fn_val = zero + ifault = 0 + IF (nreq < 1 .OR. nreq > ncol) ifault = ifault + 4 + IF (nobs <= nreq) ifault = ifault + 8 + IF (ifault /= 0) THEN + WRITE(*, '(1x, a, i4)') 'Error in function VARPRD: ifault =', ifault + RETURN + END IF + + ! Calculate the residual variance estimate. + + var = sserr / (nobs - nreq) + + ! Variance of x'b = var.x'(inv R)(inv D)(inv R')x + ! First call BKSUB2 to calculate (inv R')x by back-substitution. + + CALL BKSUB2(x, wk, nreq) + DO row = 1, nreq + IF(d(row) > tol(row)) fn_val = fn_val + wk(row)**2 / d(row) + END DO + + fn_val = fn_val * var + + RETURN + END FUNCTION varprd + + + + SUBROUTINE bksub2(x, b, nreq) + + ! Solve x = R'b for b given x, using only the first nreq rows and + ! columns of R, and only the first nreq elements of R. + ! + !-------------------------------------------------------------------------- + + INTEGER, INTENT(IN) :: nreq + REAL (dp), DIMENSION(:), INTENT(IN) :: x + REAL (dp), DIMENSION(:), INTENT(OUT) :: b + + ! Local variables + + INTEGER :: pos, row, col + REAL (dp) :: temp + + ! Solve by back-substitution, starting from the top. + + DO row = 1, nreq + pos = row - 1 + temp = x(row) + DO col = 1, row-1 + temp = temp - r(pos)*b(col) + pos = pos + ncol - col - 1 + END DO + b(row) = temp + END DO + + RETURN + END SUBROUTINE bksub2 + + + END MODULE lsq diff -crBN source_orig_select/Makefile source_interp/Makefile *** source_orig_select/Makefile 2010-11-26 14:17:22.000000000 -0500 --- source_interp/Makefile 2010-11-26 14:21:48.000000000 -0500 *************** *** 128,134 **** OBJFILES= utils.o ParamNames.o Matrix_utils.o settings.o IO.o cmbtypes.o Planck_like.o \ cmbdata.o WeakLen.o bbn.o bao.o lrggettheory.o mpk.o supernovae.o HST.o SDSSLy-a-v3.o \ $(CLSFILE) paramdef.o $(PROPOSE) $(PARAMETERIZATION) calclike.o \ ! conjgrad_wrapper.o EstCovmat.o postprocess.o MCMC.o driver.o ifeq ($(EXTDATA),LRG) --- 128,135 ---- OBJFILES= utils.o ParamNames.o Matrix_utils.o settings.o IO.o cmbtypes.o Planck_like.o \ cmbdata.o WeakLen.o bbn.o bao.o lrggettheory.o mpk.o supernovae.o HST.o SDSSLy-a-v3.o \ $(CLSFILE) paramdef.o $(PROPOSE) $(PARAMETERIZATION) calclike.o \ ! conjgrad_wrapper.o EstCovmat.o postprocess.o MCMC.o driver.o \ ! lsq.o interpolation.o ifeq ($(EXTDATA),LRG) *************** *** 183,188 **** --- 184,191 ---- postprocess.o: calclike.o MCMC.o: calclike.o driver.o: MCMC.o + paramdef.o: interpolation.o + interpolation.o: lsq.o export FFLAGS export F90C diff -crBN source_orig_select/MCMC.f90 source_interp/MCMC.f90 *** source_orig_select/MCMC.f90 2010-11-26 14:14:04.000000000 -0500 --- source_interp/MCMC.f90 2010-11-26 14:15:06.000000000 -0500 *************** *** 462,467 **** --- 462,483 ---- subroutine MCMCsample(Params, samples_to_get) + + !Interpolation Edit + use CMB_Cls + use cmbtypes + use interp_mod + + !Variables for interpolation: + integer i,j,k + real (kind=8) :: interp, temptrial(1:num_params) + integer :: dum, thresh_val, r10params, r10_interp_order + integer*8 :: tcoeff, tr10_coeff, tnum_params_used, tinterp_order, tr10_interp_order + real (kind = 8):: CambLike, iLike(2), i_sigma8, temp + Type (CMBParams) CMB + logical :: i_ind + integer, allocatable :: index_ar(:), r10_index_ar(:) + integer samples_to_get Type(ParamSet) Trial, Params, CurParams real Like, CurLike *************** *** 474,479 **** --- 490,496 ---- character(LEN=128) logLine MaxLike = LogZero + notinterp_maxlike = LogZero !Interpolation edit CurLike = StartLike testCurLike = StartLike CurParams = Params *************** *** 485,490 **** --- 502,565 ---- mc_update_burn_num = mc_update_burn mc_updates = 0 + !Interpolation Edit: WRITING TO LOGFILE AND ALLOCATING PAR ARRAYS + + tnum_params_used = num_params_used + tinterp_order = interp_order + tr10_interp_order = 2 + + i_r10a = 0 + + tcoeff = factorial(tnum_params_used + tinterp_order)/(factorial(tinterp_order)*factorial(tnum_params_used)) + coeff1 = tcoeff + tcoeff = factorial(tnum_params_used + tinterp_order-1)/(factorial(tinterp_order-1)*factorial(tnum_params_used)) + coeff = tcoeff + + if (compute_tensors) then + tr10_coeff = factorial(tnum_params_used + tr10_interp_order)/(factorial(tr10_interp_order)*factorial(tnum_params_used)) + r10_coeff1 = tr10_coeff + endif + + + ALLOCATE(paramdatabase(num_params, samples_to_get), likedatabase(samples_to_get), sigma8database(samples_to_get)) + ALLOCATE(PAR(coeff), PAR1(coeff1), par_s8(coeff1)) + ALLOCATE(emat(coeff,num_params_used), emat1(coeff1, num_params_used)) + ALLOCATE(index_ar(num_params_used)) + if (compute_tensors) ALLOCATE(r10database(samples_to_get), r10_emat1(r10_coeff1,num_params_used), par_r10(r10_coeff1), r10_index_ar(num_params_used)) + + thresh_val = Int(pts_mult)*coeff1 + + write(*,*) 'Includes Interpolation Patch' + write(logfile_unit,*) 'Includes Interpolation Patch' + write(*,*) 'sampling method', sampling_method + if (Feedback > 1) then + write(logfile_unit,*) 'interp_order:', interp_order + write(logfile_unit,*) 'threshold to start interpolation :', thresh_val + write(logfile_unit,*) 'first cut is:', cut + end if + !Initialize Exponent Arrays: + if (feedback > 1) write(*,*) 'Initializing exponent array...' + + dum = 1 + index_ar(1:num_params_used) = 0 + call make_array(num_params_used, emat1, index_ar, interp_order, dum, num_params_used) + + dum = 1 + index_ar(1:num_params_used) = 0 + call make_array(num_params_used, emat, index_ar, interp_order-1, dum, num_params_used) + + if (compute_tensors) then + dum = 1 + r10_index_ar(1:num_params_used) = 0 + r10_interp_order = tr10_interp_order + call make_array(num_params_used, r10_emat1, r10_index_ar, r10_interp_order, dum, num_params_used) !2nd order + endif + + i_ind = .false. + interp_ind = .false. + dataBaseNumber = 0 + !End Interpolation Edit + do while (num <= samples_to_get) num = num + 1 *************** *** 534,541 **** call GetProposal(CurParams, Trial) end if ! Like = GetLogLike(Trial) if (Feedback > 1) write (*,*) 'Likelihood: ', Like, 'Current Like:', CurLike !!! accpt = (Like /= logZero) .and. (CurLike > Like .or. randexp1() > Like - CurLike) --- 609,705 ---- call GetProposal(CurParams, Trial) end if ! !Like = GetLogLike(Trial) + !Interpolation Edit + if (interp_ind) then + + if (feedback > 1 ) write(*,*) 'Calculating Interpolated Likelihood...' + iLike(1) = 0 + iLike(2) = 0 + temptrial(1:num_params) = dble(Trial%P(1:num_params)) + + !Normalizing: + + do i=1,num_params + if (sd_ar(i) /= 0) temptrial(i) = (temptrial(i) - mean_ar(i))/sd_ar(i) + end do + + !Compute interpolated likelihoods: + + do j=1,coeff + temp = 1 + do k=1,num_params_used + temp = temp*temptrial(params_used(k))**emat(j,k) + end do + iLike(1) = ilike(1) + PAR(j)*temp + end do + do j=1, coeff1 + temp = 1 + do k=1, num_params_used + temp = temp*temptrial(params_used(k))**emat1(j,k) + end do + iLike(2) = iLike(2) + PAR1(j)*temp + end do + + !!! accpt = (Like /= logZero) .and. (CurLike > Like .or. randexp1() > Like - CurLike) + !Include the min() so that compilers not doing optimal compilation don't complain + if (mpk_ind) then + i_sigma8 = 0 + do j=1,coeff1 + temp = 1 + do k=1, num_params_used + temp = temp*temptrial(params_used(k))**emat1(j,k) + end do + i_sigma8 = i_sigma8 + PAR_S8(j)*temp + end do + if (feedback > 1) write(*,*) 'interpolated sigma8:', i_sigma8 + end if + + if (compute_tensors) then + i_r10 = 0 + do j=1,r10_coeff1 + temp = 1 + do k=1, num_params_used + temp = temp*temptrial(params_used(k))**r10_emat1(j,k) + end do + i_r10 = i_r10 + PAR_r10(j)*temp + end do + i_r10 = i_r10*Trial%P(12) + if (feedback > 1) write(*,*) 'interpolated r10:', i_r10 + end if + + !NOTE: Have not checked boundaries yet. This is done with the LogPrior call. + if (feedback > 1) write(*,*) 'Likelihoods: ', iLike(1), ', ', iLike(2), 'Current Like:', CurLike + + if ((abs((iLike(1)-iLike(2))/(ilike(2) - notinterp_maxlike)) > frac_cut) .or.& + abs(iLike(2) - notinterp_maxlike) > cut) then + if (Feedback > 1) write(*,*) 'NOT Using Interpolation.' + Like = GetLogLike(Trial) + i_ind = .false. + elseif(any(Trial%P > Scales%PMax) .or. any(Trial%P < Scales%PMin)) then + Like = logZero + else + call ParamsToCMBParams(Trial%P,CMB) + Camblike = GetLogPrior(CMB, Trial%Info) + if (Camblike >= logZero) then + Like = logZero + else + Like = iLike(2) + if (mpk_ind) Trial%Info%Theory%sigma_8 = i_sigma8 + i_ind = .true. + if (Feedback > 1) write(*,*) mpirank, 'Using Interpolation.' + end if + end if + + else + + !If the interpolation has not been calculated yet: + Like = GetLogLike(Trial) !Original + i_ind = .false. + end if + + !End Interpolation Edit if (Feedback > 1) write (*,*) 'Likelihood: ', Like, 'Current Like:', CurLike !!! accpt = (Like /= logZero) .and. (CurLike > Like .or. randexp1() > Like - CurLike) *************** *** 582,587 **** --- 746,760 ---- if (rmult /= 0.) then output_lines = output_lines +1 call WriteParams(CurParams,rmult,CurLike) + + !if trial has been accepted and is an interpolated value then store i_r10 in i_r10a. otherwise calculate proper r10 value. + !Interpolation Edit + if (i_ind .and. compute_tensors) then + i_r10a = i_r10 + else + i_r10a =Trial%Info%Theory%cl_tensor(10,1)/Trial%Info%Theory%cl(10,1) + end if + !End Interpolation Edit end if end if *************** *** 590,596 **** CurParams = Trial num_accept = num_accept + 1 num_metropolis_accept = num_metropolis_accept + 1 ! if (Like < MaxLike) MaxLike = Like mult=1 if (Feedback > 1) write (*,*) num_metropolis, ' accepting. ratio:', real(num_metropolis_accept)/num_metropolis --- 763,785 ---- CurParams = Trial num_accept = num_accept + 1 num_metropolis_accept = num_metropolis_accept + 1 ! ! if (Like < MaxLike) MaxLike = Like ! !Interpolation Edit ! ! if (.not.(i_ind)) then ! dataBaseNumber = dataBaseNumber + 1 ! paramDataBase(1:num_params, dataBaseNumber) = Trial%P(1:num_params) ! likeDataBase(dataBaseNumber) = Like ! if (mpk_ind) sigma8database(databaseNumber) = Trial%Info%Theory%Sigma_8 ! if (compute_tensors) r10database(databaseNumber) = Trial%Info%Theory%cl_tensor(10,1)/Trial%Info%Theory%cl(10,1)/Trial%P(12) ! end if ! ! if (Like < MaxLike) then ! MaxLike = Like ! if (.not. (i_ind)) notinterp_MaxLike = Like ! endif ! !End Interpolation Edit ! mult=1 if (Feedback > 1) write (*,*) num_metropolis, ' accepting. ratio:', real(num_metropolis_accept)/num_metropolis *************** *** 609,622 **** end if else mult=mult+1 if (Feedback > 1) write (*,*) num_metropolis,' rejecting. ratio:', real(num_metropolis_accept)/num_metropolis end if end if !not slicing if (CurLike /= logZero) then ! if (sampling_method /= sampling_slowgrid) call AddMPIParams(CurParams%P,CurLike) if (mod(num,100)==0) call CheckParamChange if (sampling_method == sampling_multicanonical .and. & mod(num,mc_update_steps + mc_updates*mc_steps_inc)==0) then --- 798,846 ---- end if else mult=mult+1 + + !Interpolation Edit + !FILTER FOR VALUES OF 1E30 + if(Like < 1.0E29 .and. .not.(i_ind)) then + dataBaseNumber = dataBaseNumber + 1 + paramDataBase(1:num_params, dataBaseNumber) = Trial%P(1:num_params) + likeDataBase(dataBaseNumber) = Like + if (mpk_ind) sigma8database(databaseNumber) = Trial%Info%Theory%Sigma_8 + if (compute_tensors) r10database(databaseNumber) = Trial%Info%Theory%cl_tensor(10,1)/Trial%Info%Theory%cl(10,1)/Trial%P(12) + elseif (.not.(i_ind) .and. Trial%P(12) < 0) then + endif + !End Interpolation Edit + if (Feedback > 1) write (*,*) num_metropolis,' rejecting. ratio:', real(num_metropolis_accept)/num_metropolis end if end if !not slicing if (CurLike /= logZero) then ! ! if (sampling_method /= sampling_slowgrid) call AddMPIParams(CurParams%P,CurLike) if (mod(num,100)==0) call CheckParamChange + if (sampling_method /= sampling_slowgrid) then + !Interpolation edit + loc_accept = num_accept + call AddMPIParams(CurParams%P,CurLike) !INTERPOLATION CALL IS IN PARAMDEF.F90 + #ifndef MPI + if (do_interp) then + if ((mod(num,100) == 1) .and. (interp_ind .eqv. false)) then + k = 0 + do i = 1, num_samples + if (ABS(par_LikeDataBase(i) - notinterp_maxlike) < cut) k = k+1 + end do + if (k >= thresh_val) then + interp_ind = .true. + call interpolate(databasenumber) + end if + elseif ( (interp_ind .eqv. .true.) .and. (mod(num,500) == 1) ) then + call interpolate(databasenumber) + endif + endif + #endif + end if + !end edit if (sampling_method == sampling_multicanonical .and. & mod(num,mc_update_steps + mc_updates*mc_steps_inc)==0) then *************** *** 647,652 **** --- 871,883 ---- end do + !Interpolation Edit + DEALLOCATE(PAR, PAR1, par_s8) + DEALLOCATE(emat, emat1) + Deallocate(paramdatabase, likedatabase, sigma8database, index_ar) + if (compute_tensors) DEALLOCATE(r10database, r10_emat1, par_r10, r10_index_ar) + !End Interpolation Edit + if (Feedback > 0) then write(*,*) 'Stopping as have ',samples_to_get ,' samples. ' write(*,*) MPIRank, ' Slow proposals: ',slow_proposals diff -crBN source_orig_select/paramdef.F90 source_interp/paramdef.F90 *** source_orig_select/paramdef.F90 2010-11-26 14:13:19.000000000 -0500 --- source_interp/paramdef.F90 2010-11-26 14:15:27.000000000 -0500 *************** *** 9,14 **** --- 9,15 ---- use cmbdata use Random use settings + use interp_mod !Interpolation edit implicit none Type ParamSet *************** *** 405,411 **** integer index, STATUS(MPI_STATUS_SIZE),STATs(MPI_STATUS_SIZE*(MPIChains-1)) logical flag ! real, allocatable, dimension(:), save ::MPIMean real, allocatable, dimension(:,:,:) ::MPICovmats real, allocatable, dimension(:,:) ::MPIMeans --- 406,413 ---- integer index, STATUS(MPI_STATUS_SIZE),STATs(MPI_STATUS_SIZE*(MPIChains-1)) logical flag ! integer maxpoints, thresh_val !Interpolation edit ! real, allocatable, dimension(:), save ::MPIMean real, allocatable, dimension(:,:,:) ::MPICovmats real, allocatable, dimension(:,:) ::MPIMeans *************** *** 426,431 **** --- 428,434 ---- logical, save :: all_burn = .false., done_check = .false., DoUpdates = .false. character(LEN=128) logLine + thresh_val = int(pts_mult)*coeff1 !Interpolation edit !Read in checkpoing stuff at restart if (checkpoint .and. present(checkpoint_start)) then *************** *** 614,619 **** --- 617,639 ---- call MPI_BCAST(MPIMeans(:,i),Size(MPIMean),MPI_REAL,j,MPI_COMM_WORLD,ierror) end do + !Interpolation Edit + if (do_interp) then + if (interp_ind .eqv. .false.) then + j=0 + do i = 1, databasenumber + if (ABS(LikeDataBase(i) - notinterp_maxlike) < cut) j = j+1 + end do + call MPI_allreduce(j, maxpoints, 1,MPI_integer, MPI_sum, MPI_COMM_WORLD, ierror) !TOTAL COLLECTED > THRESH + if (maxpoints >= thresh_val) then + call interpolate(dataBaseNumber) + interp_ind =.true. + end if + elseif (interp_ind .eqv. .true.) then + call interpolate(dataBaseNumber) + end if + endif + !End edit ! These should be better but don't work *************** *** 794,797 **** ! \ No newline at end of file --- 814,817 ---- ! diff -crBN source_orig_select/params_CMB.f90 source_interp/params_CMB.f90 *** source_orig_select/params_CMB.f90 2010-11-26 14:13:49.000000000 -0500 --- source_interp/params_CMB.f90 2010-11-26 14:15:52.000000000 -0500 *************** *** 253,259 **** call ParamsToCMBParams(P%P,CMB) if (lmax_tensor /= 0 .and. compute_tensors) then ! r10 = P%Info%Theory%cl_tensor(10,1)/P%Info%Theory%cl(10,1) else r10 = 0 end if --- 253,266 ---- call ParamsToCMBParams(P%P,CMB) if (lmax_tensor /= 0 .and. compute_tensors) then ! !Interpolation edit ! ! r10 = P%Info%Theory%cl_tensor(10,1)/P%Info%Theory%cl(10,1) ! if (fst) then ! i_r10a = P%Info%Theory%cl_tensor(10,1)/P%Info%Theory%cl(10,1) ! fst = .false. ! endif ! r10 = i_r10a ! !End interpolation edit else r10 = 0 end if diff -crBN source_orig_select/settings.f90 source_interp/settings.f90 *** source_orig_select/settings.f90 2010-11-26 14:13:36.000000000 -0500 --- source_interp/settings.f90 2010-11-26 14:15:45.000000000 -0500 *************** *** 67,72 **** --- 67,84 ---- logical :: Use_LSS = .true. + !Interpolation Edit + real (kind=8), allocatable:: PAR(:), PAR1(:), PAR_S8(:), PAR_r10(:) + integer :: interp_order, dataBaseNumber + integer, allocatable :: emat(:,:), emat1(:,:), r10_emat1(:,:) + real (kind=8), allocatable :: paramDataBase(:,:), likedatabase(:), sigma8database(:), r10database(:) + real (kind=8) :: sd_ar(1:num_params), mean_ar(1:num_params) + logical :: interp_ind = .false., mpk_ind, do_interp, fst = .true. + integer :: coeff, coeff1, r10_coeff1, loc_accept=0 + real :: notinterp_MaxLike, i_r10, i_r10a, cut, frac_cut, pts_mult=2 + integer :: kr_temp, kr_temp1 + !End Edit + integer :: logfile_unit = 0 integer :: outfile_handle = 0 integer :: indepfile_handle = 0