[Dart-dev] [5610] DART/branches/development/random_nr: Extended the tests to include the standard deviation ( for runs longer than 1)
nancy at ucar.edu
nancy at ucar.edu
Mon Mar 19 09:11:09 MDT 2012
Revision: 5610
Author: thoar
Date: 2012-03-19 09:11:09 -0600 (Mon, 19 Mar 2012)
Log Message:
-----------
Extended the tests to include the standard deviation (for runs longer than 1)
and writes out the first 100 random numbers generated from the new seed.
postprocessing and plotting routines as well ...
Modified Paths:
--------------
DART/branches/development/random_nr/test/input.nml
DART/branches/development/random_nr/test_reseed.f90
Added Paths:
-----------
DART/branches/development/random_nr/test/postprocess_test_reseed.csh
DART/branches/development/random_nr/test/random_proof.m
-------------- next part --------------
Modified: DART/branches/development/random_nr/test/input.nml
===================================================================
--- DART/branches/development/random_nr/test/input.nml 2012-03-19 14:50:48 UTC (rev 5609)
+++ DART/branches/development/random_nr/test/input.nml 2012-03-19 15:11:09 UTC (rev 5610)
@@ -1,5 +1,12 @@
+&test_reseed_nml
+ time_days = 148000,
+ time_seconds = 0,
+ time_step_seconds = 21600
+ /
+
&utilities_nml
- TERMLEVEL = 1,
- module_details = .false.
- logfilename = 'dart_log.out' /
+ TERMLEVEL = 1,
+ module_details = .false.,
+ logfilename = 'dart_log.out'
+ /
Added: DART/branches/development/random_nr/test/postprocess_test_reseed.csh
===================================================================
--- DART/branches/development/random_nr/test/postprocess_test_reseed.csh (rev 0)
+++ DART/branches/development/random_nr/test/postprocess_test_reseed.csh 2012-03-19 15:11:09 UTC (rev 5610)
@@ -0,0 +1,41 @@
+#!/bin/csh
+#
+# DART software - Copyright 2004 - 2011 UCAR. This open source software is
+# provided by UCAR, "as is", without charge, subject to all terms of use at
+# http://www.image.ucar.edu/DAReS/DART/DART_download
+#
+# $Id$
+
+foreach SEEDSEP ( 2 60 3600 21600 43200 86400 )
+
+ cd `printf seeds_%05d $SEEDSEP`
+ set fname = `printf ran_test_r8_%05d_full.out $SEEDSEP`
+
+ # focus on just the following trials
+
+ grep 'seed' $fname >! seeds
+ grep 'uniform 1000' $fname >! u1000
+ grep 'gaussian 1000' $fname >! g1000
+ grep 'gaussian 1000000' $fname >! g1000000
+
+ # strip out the string and leave the numbers to ingest into matlab
+
+ sed s/"trial "// seeds >! bob
+ sed s/"new seed: "// bob >! seeds
+ sed s/"uniform 1000"// u1000 >! unif_1000
+ sed s/"gaussian 1000"// g1000 >! gauss_1000
+ sed s/"gaussian 1000000"// g1000000 >! gauss_1000000
+
+ rm u1000 g1000 g1000000 bob
+
+ cd ..
+
+end
+
+exit $status
+
+# <next few lines under version control, do not edit>
+# $URL$
+# $Revision$
+# $Date$
+
Property changes on: DART/branches/development/random_nr/test/postprocess_test_reseed.csh
___________________________________________________________________
Added: svn:executable
+ *
Added: svn:mime-type
+ text/plain
Added: svn:keywords
+ Date Rev Author HeadURL Id
Added: svn:eol-style
+ native
Added: DART/branches/development/random_nr/test/random_proof.m
===================================================================
--- DART/branches/development/random_nr/test/random_proof.m (rev 0)
+++ DART/branches/development/random_nr/test/random_proof.m 2012-03-19 15:11:09 UTC (rev 5610)
@@ -0,0 +1,167 @@
+%% DART:random_proof checks the integrity of the random number generation
+% under usage patterns common to DART.
+%
+% Run test_reseed.f90; postprocess_test_reseed.csh
+
+%% DART software - Copyright 2004 - 2011 UCAR. This open source software is
+% provided by UCAR, "as is", without charge, subject to all terms of use at
+% http://www.image.ucar.edu/DAReS/DART/DART_download
+%
+% <next few lines under version control, do not edit>
+% $URL$
+% $Id$
+% $Revision$
+% $Date$
+
+seedsep = [2, 60, 3600, 21600, 43200, 86400];
+nseeds = length(seedsep);
+nreseeds = 1000;
+xax = 1:nreseeds;
+
+u1.mean = zeros(nreseeds,nseeds);
+u1.std = zeros(nreseeds,nseeds);
+g1.mean = zeros(nreseeds,nseeds);
+g1.std = zeros(nreseeds,nseeds);
+gM.mean = zeros(nreseeds,nseeds);
+gM.std = zeros(nreseeds,nseeds);
+
+for i = 1:nseeds,
+
+ u = load(sprintf('seeds_%05d/unif_1000' ,seedsep(i)));
+ g = load(sprintf('seeds_%05d/gauss_1000' ,seedsep(i)));
+ G = load(sprintf('seeds_%05d/gauss_1000000',seedsep(i)));
+
+ u1.mean(:,i) = u(:,1);
+ u1.std( :,i) = u(:,2);
+ g1.mean(:,i) = g(:,1);
+ g1.std( :,i) = g(:,2);
+ gM.mean(:,i) = G(:,1);
+ gM.std( :,i) = G(:,2);
+
+end
+
+%% uniform(0,1) length 1000
+
+figure(1); clf
+
+nbins = 50;
+hm = zeros(nbins,nseeds); xm = zeros(nbins,nseeds);
+hs = zeros(nbins,nseeds); xs = zeros(nbins,nseeds);
+legstr = cell(1,nseeds);
+
+for i = 1:nseeds,
+ [hm(:,i), xm(:,i)] = hist(u1.mean(:,i),nbins);
+ [hs(:,i), xs(:,i)] = hist(u1.std( :,i),nbins);
+
+ legstr{i} = sprintf('date separation %d',seedsep(i));
+end
+
+subplot(2,1,1)
+ plot(xm,hm,'-*')
+ h = legend(legstr);
+ legend(h,'boxoff')
+ title({'distribution of uniform(0,1) for a length = 1000 series',...
+ 'state vector times separated by 2 seconds'})
+ ylabel('1000 seeds')
+ xlabel('mean')
+
+subplot(2,1,2)
+ plot(xs,hs,'-*')
+ ylabel('1000 seeds')
+ xlabel('standard deviation')
+
+orient tall
+print -dpdf uniform_1000.pdf
+
+%% gaussian(0,1) length 1000
+
+figure(2); clf
+
+for i = 1:nseeds,
+ [hm(:,i), xm(:,i)] = hist(g1.mean(:,i),nbins);
+ [hs(:,i), xs(:,i)] = hist(g1.std( :,i),nbins);
+end
+
+subplot(2,1,1)
+ plot(xm,hm,'-*')
+ h = legend(legstr);
+ legend(h,'boxoff')
+ title({'distribution of gaussian(0,1) for a length = 1000 series', ...
+ 'state vector times separated by 2 seconds'})
+ ylabel('1000 seeds')
+ xlabel('mean')
+
+subplot(2,1,2)
+ plot(xs,hs,'-*')
+ ylabel('1000 seeds')
+ xlabel('standard deviation')
+
+orient tall
+print -dpdf gaussian_1000.pdf
+
+
+%% gaussian(0,1) length 1000000
+figure(3); clf
+
+for i = 1:nseeds,
+ [hm(:,i), xm(:,i)] = hist(g1.mean(:,i),nbins);
+ [hs(:,i), xs(:,i)] = hist(g1.std( :,i),nbins);
+end
+
+subplot(2,1,1)
+ plot(xm,hm,'-*')
+ h = legend(legstr);
+ legend(h,'boxoff')
+ title({'distribution of gaussian(0,1) for a length = 1,000,000 series', ...
+ 'state vector times separated by 2 seconds'})
+ ylabel('1000 seeds')
+ xlabel('mean')
+
+subplot(2,1,2)
+ plot(xs,hs,'-*')
+ ylabel('1000 seeds')
+ xlabel('standard deviation')
+
+orient tall
+print -dpdf gaussian_1000000.pdf
+
+%% plot of the first 100 values - gaussian(0,1) distrubution
+
+datmat = zeros(100,1000);
+n = 2;
+
+for i = 1:nseeds,
+
+ figure(i); clf
+
+ fname = sprintf('seeds_%05d/ran_test_r8_%05d_numbers.out',seedsep(i),seedsep(i));
+ datmat(:) = load(fname);
+ first = datmat(1,:);
+ dupes = find(first == first(1));
+ fprintf('There are %d first values identical to the first value.\n',length(dupes))
+
+ sname = sprintf('seeds_%05d/seeds',seedsep(i));
+ bob = load(sname);
+ seeds = bob(:,2); clear bob;
+
+ disp('These are the seeds that generate the identical first values:')
+ for j = 1:length(dupes),
+ fprintf(' %d\n',seeds(dupes(j)))
+ end
+
+ subplot('position',[0.1 0.1 0.76 0.8])
+ plot(datmat(1:n,:)','.')
+ h = title({sprintf('First %d values for each reseeding - gaussian(0,1) distribution',n),fname});
+ set(h,'Interpreter','none')
+ xlabel('reseeding number')
+ ylabel('gaussian(0,1)')
+ axlims = axis;
+
+ subplot('position',[0.87 0.1 0.1 0.8])
+ [N,Y] = hist(datmat(:),nbins);
+ plot(N,Y)
+ axis([-Inf Inf axlims(3) axlims(4)])
+ set(gca,'YAxisLocation','right')
+
+end
+
Property changes on: DART/branches/development/random_nr/test/random_proof.m
___________________________________________________________________
Added: svn:mime-type
+ text/plain
Added: svn:keywords
+ Date Rev Author HeadURL Id
Added: svn:eol-style
+ native
Modified: DART/branches/development/random_nr/test_reseed.f90
===================================================================
--- DART/branches/development/random_nr/test_reseed.f90 2012-03-19 14:50:48 UTC (rev 5609)
+++ DART/branches/development/random_nr/test_reseed.f90 2012-03-19 15:11:09 UTC (rev 5610)
@@ -10,13 +10,15 @@
! $Revision$
! $Date$
-use types_mod, only : r4, r8, digits12, i8
-use utilities_mod, only : register_module, error_handler, E_ERR, &
- initialize_utilities, finalize_utilities
+use types_mod, only : r4, r8, i8
+use utilities_mod, only : register_module, error_handler, E_ERR, &
+ initialize_utilities, finalize_utilities, &
+ open_file, close_file, find_namelist_in_file, &
+ check_namelist_read
use time_manager_mod, only : time_type, operator(+), set_time, get_time, &
set_calendar_type, print_time, print_date
-use random_seq_mod, only : random_seq_type, init_random_seq, &
- random_uniform, random_gaussian
+use random_seq_mod, only : random_seq_type, init_random_seq, &
+ random_uniform, random_gaussian
implicit none
@@ -28,44 +30,75 @@
type (random_seq_type) :: seq
type (time_type) :: state_time, delta_time
-integer :: i, l, rep
+
+integer :: i, l, rep, nextseed
integer :: maxreps = 1000
integer, parameter :: ntests = 5
integer :: nloops(ntests) = (/ 1, 10, 100, 1000, 1000000 /)
-real(r8) :: next_val, mean
-integer :: nextseed
+real(r8), allocatable, dimension(:) :: ranarray
+real(r8) :: mean, var
+integer :: iunit1,iunit2,iunit3, io
+character(len=32) :: filename1,filename2,filename3
+! Namelist with default values
+
+integer :: time_days = 148000
+integer :: time_seconds = 0
+integer :: time_step_seconds = 3600
+
+namelist /test_reseed_nml/ time_days, time_seconds, time_step_seconds
+
+!---------------------------------------------------------------
+
call initialize_utilities('test_reseed')
call register_module(source,revision,revdate)
+call set_calendar_type('GREGORIAN')
+! Read the namelist entry
+call find_namelist_in_file("input.nml", "test_reseed_nml", iunit1)
+read(iunit1, nml = test_reseed_nml, iostat = io)
+call check_namelist_read(iunit1, io, "test_reseed_nml")
-call set_calendar_type('GREGORIAN')
+state_time = set_time(time_seconds, time_days)
+delta_time = set_time(time_step_seconds)
-state_time = set_time(0, 148000)
-delta_time = set_time(43200)
-
print *, ' '
call print_date(state_time, 'setting start time')
call print_time(delta_time, 'delta time')
print *, ' '
+print *, 'defining r4 to be ',r4
+print *, 'defining r8 to be ',r8
print *, 'doing ', maxreps, ' repeats of ', ntests, ' counts:'
do i=1, ntests
- print *, 'mean for ', nloops(i), ' uniform distribution'
+ print *, 'mean[,std] for ', nloops(i), ' uniform distribution'
enddo
do i=1, ntests
- print *, 'mean for ', nloops(i), ' gaussian distribution'
+ print *, 'mean[,std] for ', nloops(i), ' gaussian distribution'
enddo
print *, ' '
+write(filename1,'(''ran_test_r'',I1,''_'',I5.5,''_compact.out'')') r8, time_step_seconds
+write(filename2,'(''ran_test_r'',I1,''_'',I5.5,''_full.out'' )') r8, time_step_seconds
+write(filename3,'(''ran_test_r'',I1,''_'',I5.5,''_numbers.out'')') r8, time_step_seconds
+
+iunit1 = open_file(filename1, form='formatted', action='write')
+iunit2 = open_file(filename2, form='formatted', action='write')
+iunit3 = open_file(filename3, form='formatted', action='write')
+
+write(*,*)'creating compact output file ',trim(filename1)
+write(*,*)'creating output file ',trim(filename2)
+write(*,*)'creating output file ',trim(filename3)
+
do rep=1, maxreps
nextseed = generate_seed(state_time)
- print *, ' '
- print *, 'new seed: ', nextseed
+ write( * ,'('' trial '',i5,'' new seed: '', i14)')rep, nextseed
+ write(iunit1,'('' trial '',i5,'' new seed: '', i14)')rep, nextseed
+ write(iunit2,'('' trial '',i5,'' new seed: '', i14)')rep, nextseed
! ---------
@@ -73,21 +106,38 @@
call init_random_seq(seq, nextseed)
+ allocate(ranarray(nloops(l)))
+
mean = 0.0_r8
+ var = 0.0_r8
do i=1, nloops(l)
-
- next_val = random_uniform(seq)
- mean = mean + next_val
-
+ ranarray(i) = random_uniform(seq)
+ mean = mean + ranarray(i)
enddo
- print *, mean/nloops(l), nloops(l), ' uniform'
+ mean = mean/real(nloops(l),r8)
+ do i=1, nloops(l)
+ var = var + (mean - ranarray(i))**2
+ enddo
+
+ if (l /= 1) then
+ write(iunit1,'(''uniform '',i7,1x,f9.6, 1x,f9.6 )') nloops(l), mean, sqrt(var/(nloops(l)-1))
+ write(iunit2,'(''uniform '',i7,1x,f19.16,1x,f19.16)') nloops(l), mean, sqrt(var/(nloops(l)-1))
+ else
+ write(iunit1,'(''uniform '',i7,1x,f9.6 )') nloops(l), mean
+ write(iunit2,'(''uniform '',i7,1x,f19.16 )') nloops(l), mean
+ endif
+
+ deallocate(ranarray)
+
enddo
! ---------
- print *, ' '
+ write( * ,*)
+ write(iunit1,*)
+ write(iunit2,*)
! ---------
@@ -95,21 +145,44 @@
call init_random_seq(seq, nextseed)
+ allocate(ranarray(nloops(l)))
+
mean = 0.0_r8
+ var = 0.0_r8
do i=1, nloops(l)
-
- next_val = random_gaussian(seq, 0.0_r8, 1.0_r8)
- mean = mean + next_val
-
+ ranarray(i) = random_gaussian(seq, 0.0_r8, 1.0_r8)
+ mean = mean + ranarray(i)
enddo
- print *, mean/nloops(l), nloops(l), ' gaussian'
+ mean = mean/real(nloops(l),r8)
+ do i=1, nloops(l)
+ var = var + (mean - ranarray(i))**2
+ enddo
+
+ if (l /= 1) then
+ write(iunit1,'(''gaussian '',i7,1x,f9.6, 1x,f9.6 )') nloops(l), mean, sqrt(var/(nloops(l)-1))
+ write(iunit2,'(''gaussian '',i7,1x,f19.16,1x,f19.16)') nloops(l), mean, sqrt(var/(nloops(l)-1))
+ else
+ write(iunit1,'(''gaussian '',i7,1x,f9.6 )') nloops(l), mean
+ write(iunit2,'(''gaussian '',i7,1x,f19.16 )') nloops(l), mean
+ endif
+
+ if ( nloops(l) == 100 ) then
+ do i=1, nloops(l)
+ write(iunit3,*) ranarray(i)
+ enddo
+ endif
+
+ deallocate(ranarray)
+
enddo
! ---------
- print *, ' '
+ write( * ,*)
+ write(iunit1,*)
+ write(iunit2,*)
! ---------
@@ -117,6 +190,10 @@
enddo
+call close_file(iunit1)
+call close_file(iunit2)
+call close_file(iunit3)
+
call finalize_utilities()
contains
@@ -124,7 +201,7 @@
!-------------------------------------------------------------------------
function generate_seed(state_time)
-! use the state time to set the seed for the (repeatable) random sequence
+! use the state time to set the seed for the (repeatable) random sequence
type(time_type), intent(in) :: state_time
integer :: generate_seed
More information about the Dart-dev
mailing list