trand.py - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
 (HTM) git clone git://src.adamsgaard.dk/pism
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) LICENSE
       ---
       trand.py (2898B)
       ---
            1 ############################################################################
            2 #
            3 #  This file is a part of siple.
            4 #
            5 #  Copyright 2010 David Maxwell
            6 #
            7 #  siple is free software: you can redistribute it and/or modify
            8 #  it under the terms of the GNU General Public License as published by
            9 #  the Free Software Foundation, either version 2 of the License, or
           10 #  (at your option) any later version.
           11 # 
           12 ############################################################################
           13 
           14 import numpy as np
           15 
           16 # This module contains utility functions for generating vectors
           17 # with randomly generated coefficients.
           18 
           19 def random_vector(u,scale=1,seed=None):
           20   """Returns a vector in the same function space as u with normally
           21   distributed coefficients with mean zero and standard deviation 'scale'
           22   """
           23   if not (seed is None):
           24     np.random.seed(seed)
           25 
           26   h = u.vector_like()
           27   h.set(np.random.normal(scale=scale, size=u.size()))
           28   return h
           29   
           30 
           31 def random_direction(u,scale=1,seed=None):
           32   """
           33   Returns a random vector with normally distributed entries with standard 
           34   deviation ||u||_\infty * scale.
           35   
           36   Useful for obtaining a random direction vector with size on the order of u.
           37   """
           38   scale *= u.norm('linf')
           39   
           40   return random_vector(u,scale=scale,seed=seed)
           41 
           42 
           43 def addnoise(u,relative_percent_error,seed=None):
           44   """
           45   Add gaussian noise to a vector.  The standard deviation of the noise
           46   is (relative_percent_error)\% of the maximum value of the vector u.
           47   """
           48   if not (seed is None):
           49     np.random.seed(seed)
           50 
           51   noisy_u = u.copy()    
           52   max_u = u.norm('linf')
           53 
           54   scale = (relative_percent_error/100.*max_u)
           55   noisy_u += random_vector(u,scale=scale,seed=seed)
           56 
           57   return (nu,scale)
           58 
           59 def noiseForGrid(N,M,dx,dy,minWaveLength):
           60   """
           61   Returns noise on an N by M grid with spacing dx and dy with a frequency spectrum as follows:
           62 
           63     1) frequencies corresponding to wavelengths small than minWaveLength are all zero
           64     2) amplitude of each nonzero component in normally distributed
           65     3) phase of each nonzero component is uniformly distributed
           66 
           67     The noise has zero mean and is scaled to have RMS=1.
           68   """
           69   x=np.arange(0,N,1)
           70   y=np.arange(0,M,1)
           71 
           72   (Y,X)=np.meshgrid(y,x)
           73 
           74   r = np.random.randn(N,M)
           75 
           76   theta = np.random.rand(N,M)*2*np.pi*(1J)
           77   freqs = r*np.exp(theta)
           78 
           79   a=N*dx/minWaveLength
           80   b=M*dy/minWaveLength
           81 
           82   mask1 = np.sqrt(X**2/a**2+Y**2/b**2)>1
           83   mask2 = np.sqrt((X-N)**2/a**2+Y**2/b**2)>1
           84 
           85   freqs[mask1&mask2]=0
           86 
           87   # This could be done faster with fancy slicing
           88   for i in range(1,N/2):
           89     for j in range(1,M/2):
           90       freqs[N-i,M-j] = np.conj(freqs[i,j])
           91       freqs[i,M-j] = np.conj(freqs[N-i,j])
           92   for i in range(1,N/2):
           93     freqs[N-i,0]=np.conj(freqs[i,0])
           94   for j in range(1,M/2):
           95     freqs[0,M-j]=np.conj(freqs[0,j])
           96   # Ensure we have mean zero
           97   freqs[0,0]=0
           98 
           99   noise = np.real(np.fft.ifft2(freqs))
          100 
          101   rms = np.sqrt(sum(sum(noise*noise))/(N*M))
          102   noise/=rms
          103 
          104   return noise