1!----------------------------------------------------------------------------- 2! 3! Copyright (C) 1997-2013 Krzysztof M. Gorski, Eric Hivon, 4! Benjamin D. Wandelt, Anthony J. Banday, 5! Matthias Bartelmann, Hans K. Eriksen, 6! Frode K. Hansen, Martin Reinecke 7! 8! 9! This file is part of HEALPix. 10! 11! HEALPix is free software; you can redistribute it and/or modify 12! it under the terms of the GNU General Public License as published by 13! the Free Software Foundation; either version 2 of the License, or 14! (at your option) any later version. 15! 16! HEALPix is distributed in the hope that it will be useful, 17! but WITHOUT ANY WARRANTY; without even the implied warranty of 18! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 19! GNU General Public License for more details. 20! 21! You should have received a copy of the GNU General Public License 22! along with HEALPix; if not, write to the Free Software 23! Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA 24! 25! For more information about HEALPix see http://healpix.sourceforge.net 26! 27!----------------------------------------------------------------------------- 28module sub_ngpdf_sho 29 30 31 use healpix_types 32 integer(i4b), parameter :: namax = 20 33 real(DP), dimension(0:namax) :: factorial 34! logical(LGT) :: precomp_done = .false. 35 logical(LGT), save :: precomp_done ! 2013-05-07: edited for G95 36 37 38 private 39 public :: shodev_driver 40 ! subroutine shodev_driver 41 ! subroutine compute_factorial 42 ! function shopdf 43 ! subroutine rand_sho 44 45contains 46 47 48 !-----*----------------------------------------------------------------- 49 50 SUBROUTINE shodev_driver(handle, ns, sigma0, x, nu, bins) 51 52 Use healpix_types 53 Use paramfile_io, ONLY : paramfile_handle, parse_int, parse_double, concatnl, parse_lgt 54 use rngmod, only: rand_init, planck_rng, rand_uni 55 use misc_utils, only : wall_clock_time, assert_alloc 56 Use sky_sub 57 use statistics, only: compute_statistics, print_statistics, tstats 58 59 implicit none 60 61 ! type(paramfile_handle), Intent(In) :: handle 62 integer(I4B), parameter :: KMAP = SP 63 type(paramfile_handle), Intent(InOut) :: handle 64 Integer(I4B), Intent(In) :: ns 65 Real(DP), Intent(Out) :: sigma0 66 Real(KMAP), Intent(Out), dimension(0:ns-1) :: x 67 Real(DP), Intent(Out) :: nu(1:3) 68 Integer, Intent(In), Optional :: bins 69 70 integer :: namax,nsmax 71 integer :: i,j,iseed,na,status 72 real(KMAP) :: xmin,xmax 73 real(DP) :: skewrel, alpha0, mean 74 real(DP), dimension(:), allocatable :: alpha 75 character(len=filenamelen) :: description 76 character(len=20) :: chline 77 logical :: plot, do_bin = .false. 78 character(len = 20) :: form1 = '(a, i1)', form2 = '(a, i1, a)' 79#ifdef PGPLOT 80 Integer, parameter :: npts = 1000 81 real(SP), dimension(1:npts) :: xline,yline 82 real(DP) :: xval,xstep 83 integer :: pgopen, pgdev 84#endif 85 character(len=*), parameter :: code = 'shodev_driver' 86 type(planck_rng) :: rng_handle 87 real(DP), parameter :: one = 1.0_dp, two = 2.0_dp, three = 3.0_dp 88 real(DP), parameter :: four = 4.0_dp, six = 6.0_dp 89 real(DP), allocatable, dimension(:) :: shovar 90 real :: time0, time1, time2 91 integer :: loc(1) 92 type(tstats) :: my_stats 93 94 95 ! initialize variables 96 97 If (present(bins)) do_bin = .True. 98 precomp_done = .false. 99 100 chline = 'iseed' 101 If (do_bin) call add_subscript(chline, bins) 102 description = concatnl( "", & 103 & " Input seed integer iseed:") 104 iseed = parse_int(handle, chline, default=1, descr=description) 105 106 ! --- If necessary use current system time to set seed --- 107 IF (iseed .eq. 0) THEN 108 CALL SYSTEM_CLOCK(COUNT = iseed) 109 !-- Seed should initially be set to a large, odd integer -- 110 IF (MOD(iseed,2) .EQ. 0) iseed = iseed + 1 111 PRINT *," "//code//"> Generating random number seed ", iseed 112 END IF 113 114 call rand_init(rng_handle, iseed) 115 116 chline = 'sigma0' 117 If (do_bin) call add_subscript(chline, bins) 118 description = concatnl( "", & 119 & " Input value of sigma0:") 120! sigma0 = parse_real(handle, chline, default=1.0, descr=description) 121 sigma0 = parse_double(handle, chline, default=1.d0, descr=description) 122 123 chline = 'na' 124 If (do_bin) call add_subscript(chline, bins) 125 description = concatnl( "", & 126 & " Input order of highest non-zero alpha (na):") 127 na = parse_int(handle, chline, default=3, vmin = 0, vmax = 20, descr=description) 128 ! Note maximum value of na is 20, but currently moments can only be calculated analytically 129 ! for na no greater than 3 130 131 132 !---allocate space for arrays--- 133 allocate(shovar(0:ns-1), stat=status) 134 call assert_alloc(status, code, ' cannot allocate space for alpha array') 135 136 allocate(alpha(1:MAX(na,3)), stat=status) 137 call assert_alloc(status, code, ' cannot allocate space for alpha array') 138 alpha = 0.0_dp 139 !---end of allocation for arrays--- 140 141 Do i = 1, na 142 If (i .EQ. 10) Then 143 form1 = "(a, i2)" 144 form2 = "(a, i2, a)" 145 End If 146 Write (chline, form1) "alpha_", i 147 Write (description, form2) ' Input alpha(',i,')' 148 If (do_bin) call add_subscript(chline, bins) 149 description = concatnl( "", description) 150 alpha(i) = parse_double(handle, chline, default = 0.0_dp, vmin = 0.0_dp, & 151 &vmax = sqrt(1.0_dp-SUM(alpha**2)), descr = description) 152 End Do 153 154 ! 155 Write (*,*) 'Drawing pixel values from non-Gaussian pdf...' 156 if (.not. precomp_done) call compute_factorial() 157 call rand_sho(rng_handle, ns, shovar, sigma0, alpha, na) 158 call compute_statistics(shovar, my_stats) 159 call print_statistics(my_stats) 160 xmin = minval(shovar) 161 xmax = maxval(shovar) 162 163 mean = 0.0 164 if (na.le.3) then 165 ! compute moments of this distribution 166 alpha0 = sqrt(1.0 - SUM(alpha**2)) 167 write (*,*) 'alpha0 is ',alpha0 168 nu(1)=sqrt(TWO*sigma0)*(TWO*alpha(1)*alpha(2) + sqrt(TWO)*alpha0*alpha(1) & 169 & + sqrt(SIX)*alpha(2)*alpha(3)) 170 nu(2)=sigma0**2 *(ONE + TWO*alpha(1)**2 + FOUR*alpha(2)**2 + SIX*alpha(3)**2 & 171 & + TWO*sqrt(TWO)*alpha0*alpha(2) + TWO*sqrt(SIX)*alpha(1)*alpha(3)) 172 nu(3)=THREE*alpha0*alpha(1)/sqrt(TWO) + SIX*alpha(1)*alpha(2) + & 173 & 9.0_dp*sqrt(THREE)*alpha(2)*alpha(3)/sqrt(TWO) + sqrt(THREE)*alpha0*alpha(3) 174 nu(3)=(TWO*sigma0**2)**(THREE/TWO)*nu(3) 175 176 skewrel=nu(3)/nu(2)**(THREE/TWO) 177 178 Write(*,*) 'Raw (theoretical) moments of distribution:' 179 write(*,*) 'nu1=', nu(1) 180 write(*,*) 'nu2=', nu(2) 181 write(*,*) 'nu3=', nu(3) 182 write(*,*) 'Relative skewness=', skewrel 183 184 If (nu(1) .NE. 0.0) Then 185 Write(*,*) 'Re-centering distribution on zero' 186 mean = nu(1) 187 shovar = shovar - nu(1) 188 xmin = xmin - nu(1) 189 xmax = xmax - nu(1) 190 nu(3) = nu(3) - THREE*nu(1)*nu(2) + TWO*nu(1)**3 191 nu(2) = nu(2) - nu(1)**2 192 nu(1) = 0.0_dp 193 skewrel=nu(3)/nu(2)**(THREE/TWO) 194 Write (*,*) 'Centered (theoretical) moments are:' 195 Write (*,*) 'nu1=', nu(1) 196 write (*,*) 'nu2=', nu(2) 197 write (*,*) 'nu3=', nu(3) 198 write (*,*) 'Relative skewness=', skewrel 199 End If 200 201 Else 202 nu = 0.0_dp 203 nu(1) = -1.0_dp 204 Write (*,*) 'Cannot calculate moments when highest non-zero alpha > 3' 205 End If 206 207 x = shovar 208 209 ! Find out if plotting is required 210#ifdef PGPLOT 211 description = concatnl( "", & 212 & " Plot histogram of pixel values? (True or False)") 213 plot = parse_lgt(handle, 'plot', default=.false., descr=description) 214 215 If (plot) then 216! call pgbegin(0,'?',1,1) 217111 pgdev = pgopen('?') 218 if (pgdev < 0) goto 111 219 call pghist(ns,real(x,SP),real(xmin,SP),real(xmax,SP),200,0) 220 call pglab('Pixel value','Number of pixels','Non-Gaussian SHO White Noise Map') 221 222 xstep=(xmax-xmin)/real(npts-1, DP) 223 Do i=1,npts 224 xval=xmin+real(i -1, DP)*xstep 225 xline(i)=xval 226 yline(i)=shopdf(xval+mean, sigma0, alpha, na) 227 End Do 228 yline = yline*real(ns)*(xmax-xmin)/200. 229 call pgsci(12) 230 call pgline(npts,xline,yline) 231 call pgend 232 End If !plot 233#endif 234 235 236 END SUBROUTINE shodev_driver 237 238 !================================================================================== 239 subroutine compute_factorial() 240 !------------------------------------------------------------------------------- 241 integer(i4b) :: i 242 243 244 if (precomp_done) return 245 246 factorial(0) = 1.0_dp 247 factorial(1) = 1.0_dp 248 249 do i = 2, namax 250 factorial(i) = factorial(i-1) * real(i, DP) 251 enddo 252 precomp_done = .true. 253 254 return 255 end subroutine compute_factorial 256 257 !-----*----------------------------------------------------------------- 258 259 FUNCTION shopdf(x, sigma0, alpha, na) result(spdf) 260 !------------------------------------------------------------ 261 ! calculates the general simple harmonic oscillator pdf 262 !------------------------------------------------------------ 263 use healpix_types 264 use num_rec, only : othpl 265 implicit none 266 real(DP), intent(in) :: x 267 integer(i4b), intent(in) :: na 268 real(DP), intent(in) :: sigma0 269 real(DP), intent(in) :: alpha(Max(na,3)) 270 real(DP) :: spdf 271 272 integer, parameter :: namax=20 273 integer(i4b) :: i,m,n 274 real(DP), dimension(0:namax) :: pars, hermxp, dhermxp 275 real(DP) :: xp,sumsq,fact 276 !------------------------------------------------------------ 277 278 xp = x/(sqrt(2.d0)*sigma0) 279 280 pars(1:na) = alpha(1:na) 281 sumsq=0d0 282 do i=1,na 283 sumsq=sumsq + pars(i)*pars(i) 284 end do 285 pars(0)=sqrt(1d0-sumsq) 286 287 call othpl(4,na,xp,hermxp,dhermxp) 288 289 if (.not. precomp_done) call compute_factorial 290 291 spdf=0.0_dp 292 do m=0,na 293 do n=0,na 294 fact= sqrt( (2d0**dble(m+n)) * factorial(m) * factorial(n) ) 295 spdf=spdf+ pars(m)*pars(n)*hermxp(m)*hermxp(n) / fact 296 end do 297 end do 298 spdf=spdf* exp(-xp*xp)/sqrt(2d0*pi)/sigma0 299 300 return 301 END FUNCTION shopdf 302 303! ===================================================================== 304 subroutine rand_sho(rng_handle, npix, xvar, sigma0, alpha, na) 305 !------------------------------------------------------------ 306 ! draws npix numbers from a SHO variable 307 !------------------------------------------------------------ 308 use healpix_types 309 use misc_utils, only: fatal_error 310 use rngmod, only: rand_uni, planck_rng 311 implicit none 312 313 type(planck_rng), intent(inout) :: rng_handle 314 integer(i4b), intent(in) :: na, npix 315 real(DP), intent(out) :: xvar(0: npix-1) 316 real(DP), intent(in) :: sigma0 317 real(DP), intent(in) :: alpha(Max(na,3)) 318 319 character(len=*), parameter :: code = 'rand_sho' 320 real(DP), parameter :: threshold = 1.d-12 321 322 real(DP) :: x, xmin, xmax, xstep, xlin, y, yold, yuni 323 real(DP), dimension(:), allocatable :: integral 324 integer(I4B) :: nx, i, p, ilow, ihi, imid 325 !---------------------------------------------------------- 326 327 ! find out x range 328 xmax = 5.0_dp * sigma0 329 xmin = -xmax 330 do 331 if (shopdf(xmax, sigma0, alpha, na) < threshold) exit 332 xmax = xmax + sigma0 333 enddo 334 do 335 if (shopdf(xmin, sigma0, alpha, na) < threshold) exit 336 xmin = xmin - sigma0 337 enddo 338 nx = (xmax - xmin)/sigma0 * 2000 339 340 ! build cumulative table, using trapeze integral 341 allocate(integral(0:nx-1)) 342 xstep = (xmax-xmin)/(nx-1.d0) 343 integral = 0_dp 344 yold = 0.0_dp 345 do i=0, nx - 1 346 x = xmin + i * xstep 347 y = shopdf(x, sigma0, alpha, na) 348 if (i > 0) integral(i) = integral(i-1) + (y+yold)*0.5_dp 349 yold = y 350 enddo 351 integral = integral * xstep ! normalize 352 if (abs(integral(nx-1)-1.d0) > 5.d0*threshold) then 353 call fatal_error(code//': problem with computation of integral') 354 print*,'Error on integral: ',integral(nx-1)-1.d0 355 endif 356 integral(nx-1) = 1.d0 357 358 ! solve y = F(x) 359 do p = 0, npix-1 360 ! uniform variate in [0,1] 361 yuni = rand_uni(rng_handle) 362 ! locate point in tabulated integral by dichotomy 363 ilow = 0 364 ihi = nx-1 365 do 366 imid = (ihi + ilow)/2 367 if ((integral(imid) <= yuni) .eqv. (yuni <= integral(ihi))) then 368 ilow = imid 369 else 370 ihi = imid 371 endif 372 if ((ihi - ilow) == 1) exit 373 enddo 374 ! linear interpolation 375 xlin = (yuni - integral(ilow))/(integral(ihi)-integral(ilow)) * xstep 376 ! return non-gaussian variate 377 xvar(p) = xmin + ilow * xstep + xlin 378 enddo 379 deallocate(integral) 380 381 return 382 end subroutine rand_sho 383 !-----*----------------------------------------------------------------- 384 385end module sub_ngpdf_sho 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408