danh
/
pywafo
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Avoid using "non-Pythonic" variable names

Browse Source
master
Per A Brodtkorb 8 years ago
parent e26325dfc3
commit 9930a88572
1 changed files with 84 additions and 82 deletions
Show Stats Download Patch File Download Diff File

166
wafo/covariance/core.py
Unescape Escape View File

@ -176,7 +176,7 @@ class CovData1D(PlotData):
nugget effect to ensure that round off errors do not result in nugget effect to ensure that round off errors do not result in
negative spectral estimates. Good choice might be 10^-12. negative spectral estimates. Good choice might be 10^-12.
trunc : scalar, real trunc : scalar, real
truncates all spectral values where S/max(S) < trunc truncates all spectral values where spec/max(spec) < trunc
0 <= trunc <1 This is to ensure that high frequency 0 <= trunc <1 This is to ensure that high frequency
noise is not added to the spectrum. (default 1e-5) noise is not added to the spectrum. (default 1e-5)
fast : bool fast : bool
@ -185,7 +185,7 @@ class CovData1D(PlotData):
Returns Returns
-------- --------
S : SpecData1D object spec : SpecData1D object
spectral density spectral density
NB! This routine requires that the covariance is evenly spaced NB! This routine requires that the covariance is evenly spaced
@ -205,14 +205,14 @@ class CovData1D(PlotData):
>>> S0 = R0.tospecdata() >>> S0 = R0.tospecdata()
>>> Sj = sm.Jonswap() >>> Sj = sm.Jonswap()
>>> S = Sj.tospecdata() >>> spec = Sj.tospecdata()
>>> R2 = S.tocovdata() >>> R2 = spec.tocovdata()
>>> S1 = R2.tospecdata() >>> S1 = R2.tospecdata()
>>> abs(S1.data-S.data).max() < 1e-4 >>> abs(S1.data-spec.data).max() < 1e-4
True True
S1.plot('r-') S1.plot('r-')
S.plot('b:') spec.plot('b:')
pylab.show() pylab.show()
See also See also
@ -252,26 +252,26 @@ class CovData1D(PlotData):
nf = int(nfft // 2) # number of frequencies nf = int(nfft // 2) # number of frequencies
acf = r_[acf, zeros(nfft - 2 * n + 2), acf[n - 2:0:-1]] acf = r_[acf, zeros(nfft - 2 * n + 2), acf[n - 2:0:-1]]
Rper = (fft(acf, nfft).real).clip(0) # periodogram r_per = (fft(acf, nfft).real).clip(0) # periodogram
RperMax = Rper.max() r_per_max = r_per.max()
Rper = where(Rper < trunc * RperMax, 0, Rper) r_per = where(r_per < trunc * r_per_max, 0, r_per)
S = abs(Rper[0:(nf + 1)]) * dt / pi spec = abs(r_per[0:(nf + 1)]) * dt / pi
w = linspace(0, pi / dt, nf + 1) w = linspace(0, pi / dt, nf + 1)
So = _wafospec.SpecData1D(S, w, type=spectype, freqtype=ftype) spec_out = _wafospec.SpecData1D(spec, w, type=spectype, freqtype=ftype)
So.tr = self.tr spec_out.tr = self.tr
So.h = self.h spec_out.h = self.h
So.norm = self.norm spec_out.norm = self.norm
if method != 'fft' and rate > 1: if method != 'fft' and rate > 1:
So.args = linspace(0, pi / dt, nf * rate) spec_out.args = linspace(0, pi / dt, nf * rate)
if method == 'stineman': if method == 'stineman':
So.data = stineman_interp(So.args, w, S) spec_out.data = stineman_interp(spec_out.args, w, spec)
else: else:
intfun = interpolate.interp1d(w, S, kind=method) intfun = interpolate.interp1d(w, spec, kind=method)
So.data = intfun(So.args) spec_out.data = intfun(spec_out.args)
So.data = So.data.clip(0) # clip negative values to 0 spec_out.data = spec_out.data.clip(0) # clip negative values to 0
return So return spec_out
def sampling_period(self): def sampling_period(self):
''' '''
@ -303,7 +303,7 @@ class CovData1D(PlotData):
Parameters Parameters
---------- ----------
ns : scalar ns : scalar
number of simulated points. (default length(S)-1=n-1). number of simulated points. (default length(spec)-1=n-1).
If ns>n-1 it is assummed that R(k)=0 for all k>n-1 If ns>n-1 it is assummed that R(k)=0 for all k>n-1
cases : scalar cases : scalar
number of replicates (default=1) number of replicates (default=1)
@ -336,8 +336,8 @@ class CovData1D(PlotData):
Example: Example:
>>> import wafo.spectrum.models as sm >>> import wafo.spectrum.models as sm
>>> Sj = sm.Jonswap() >>> Sj = sm.Jonswap()
>>> S = Sj.tospecdata() #Make spec >>> spec = Sj.tospecdata() #Make spec
>>> R = S.tocovdata() >>> R = spec.tocovdata()
>>> x = R.sim(ns=1000,dt=0.2) >>> x = R.sim(ns=1000,dt=0.2)
See also See also
@ -365,7 +365,7 @@ class CovData1D(PlotData):
n = acf.size n = acf.size
acf.shape = (n, 1) acf.shape = (n, 1)
dT = self.sampling_period() dt = self.sampling_period()
x = zeros((ns, cases + 1)) x = zeros((ns, cases + 1))
@ -391,10 +391,10 @@ class CovData1D(PlotData):
acf = r_[acf, zeros((nfft - m2, 1)), acf[n - 1:1:-1, :]] acf = r_[acf, zeros((nfft - m2, 1)), acf[n - 1:1:-1, :]]
# m2=2*n-2 # m2=2*n-2
S = fft(acf, nfft, axis=0).real # periodogram spec = fft(acf, nfft, axis=0).real # periodogram
I = S.argmax() I = spec.argmax()
k = flatnonzero(S < 0) k = flatnonzero(spec < 0)
if k.size > 0: if k.size > 0:
_msg = ''' _msg = '''
Not able to construct a nonnegative circulant vector from ACF. Not able to construct a nonnegative circulant vector from ACF.
@ -404,7 +404,7 @@ class CovData1D(PlotData):
# truncating negative values to zero to ensure that # truncating negative values to zero to ensure that
# that this noise is not added to the simulated timeseries # that this noise is not added to the simulated timeseries
S[k] = 0. spec[k] = 0.
ix = flatnonzero(k > 2 * I) ix = flatnonzero(k > 2 * I)
if ix.size > 0: if ix.size > 0:
@ -413,13 +413,13 @@ class CovData1D(PlotData):
# that high frequency noise is not added to # that high frequency noise is not added to
# the simulated timeseries. # the simulated timeseries.
ix0 = k[ix[0]] ix0 = k[ix[0]]
S[ix0:-ix0] = 0.0 spec[ix0:-ix0] = 0.0
trunc = 1e-5 trunc = 1e-5
maxS = S[I] max_spec = spec[I]
k = flatnonzero(S[I:-I] < maxS * trunc) k = flatnonzero(spec[I:-I] < max_spec * trunc)
if k.size > 0: if k.size > 0:
S[k + I] = 0. spec[k + I] = 0.
# truncating small values to zero to ensure that # truncating small values to zero to ensure that
# that high frequency noise is not added to # that high frequency noise is not added to
# the simulated timeseries # the simulated timeseries
@ -430,18 +430,18 @@ class CovData1D(PlotData):
# randn = np.random.randn # randn = np.random.randn
epsi = randn(nfft, cases2) + 1j * randn(nfft, cases2) epsi = randn(nfft, cases2) + 1j * randn(nfft, cases2)
Ssqr = sqrt(S / (nfft)) # sqrt(S(wn)*dw ) sqrt_spec = sqrt(spec / (nfft)) # sqrt(spec(wn)*dw )
ephat = epsi * Ssqr # [:,np.newaxis] ephat = epsi * sqrt_spec # [:,np.newaxis]
y = fft(ephat, nfft, axis=0) y = fft(ephat, nfft, axis=0)
x[:, 1:cases + 1] = hstack((y[2:ns + 2, 0:cases2].real, x[:, 1:cases + 1] = hstack((y[2:ns + 2, 0:cases2].real,
y[2:ns + 2, 0:cases1].imag)) y[2:ns + 2, 0:cases1].imag))
x[:, 0] = linspace(0, (ns - 1) * dT, ns) # (0:dT:(dT*(np-1)))' x[:, 0] = linspace(0, (ns - 1) * dt, ns) # (0:dt:(dt*(np-1)))'
if derivative: if derivative:
Ssqr = Ssqr * \ sqrt_spec = sqrt_spec * \
r_[0:(nfft / 2 + 1), -(nfft / 2 - 1):0] * 2 * pi / nfft / dT r_[0:(nfft / 2 + 1), -(nfft / 2 - 1):0] * 2 * pi / nfft / dt
ephat = epsi * Ssqr # [:,newaxis] ephat = epsi * sqrt_spec # [:,newaxis]
y = fft(ephat, nfft, axis=0) y = fft(ephat, nfft, axis=0)
xder[:, 1:(cases + 1)] = hstack((y[2:ns + 2, 0:cases2].imag - xder[:, 1:(cases + 1)] = hstack((y[2:ns + 2, 0:cases2].imag -
y[2:ns + 2, 0:cases1].real)) y[2:ns + 2, 0:cases1].real))
@ -486,24 +486,26 @@ class CovData1D(PlotData):
else: else:
return acf[:n] return acf[:n]
def _split_cov(self, sigma, i_known, i_unknown): @staticmethod
def _split_cov(sigma, i_known, i_unknown):
''' '''
Split covariance matrix between known/unknown observations Split covariance matrix between known/unknown observations
Returns Returns
------- -------
Soo covariance between known observations soo covariance between known observations
S11 = covariance between unknown observations s1o = covariance between known and unknown obs
S1o = covariance between known and unknown obs s11 = covariance between unknown observations
''' '''
Soo, So1 = sigma[i_known][:, i_known], sigma[i_known][:, i_unknown] soo, so1 = sigma[i_known][:, i_known], sigma[i_known][:, i_unknown]
S11 = sigma[i_unknown][:, i_unknown] s11 = sigma[i_unknown][:, i_unknown]
return Soo, So1, S11 return soo, so1, s11
def _update_window(self, idx, i_unknown, num_x, num_acf, @staticmethod
def _update_window(idx, i_unknown, num_x, num_acf,
overlap, nw, num_restored): overlap, nw, num_restored):
Nsig = len(idx) num_sig = len(idx)
start_max = num_x - Nsig start_max = num_x - num_sig
if (nw == 0) and (num_restored < len(i_unknown)): if (nw == 0) and (num_restored < len(i_unknown)):
# move to the next missing data # move to the next missing data
start_ix = min(i_unknown[num_restored + 1] - overlap, start_max) start_ix = min(i_unknown[num_restored + 1] - overlap, start_max)
@ -594,38 +596,38 @@ class CovData1D(PlotData):
if method.startswith('exac'): if method.startswith('exac'):
# exact but slow. It also may not return any result # exact but slow. It also may not return any result
if num_acf > 0.3 * num_x: if num_acf > 0.3 * num_x:
Sigma = toeplitz(hstack((acf, zeros(num_x - num_acf)))) sigma = toeplitz(hstack((acf, zeros(num_x - num_acf))))
else: else:
acf[0] = acf[0] * 1.00001 acf[0] = acf[0] * 1.00001
Sigma = sptoeplitz(hstack((acf, zeros(num_x - num_acf)))) sigma = sptoeplitz(hstack((acf, zeros(num_x - num_acf))))
Soo, So1, S11 = self._split_cov(Sigma, i_known, i_unknown) soo, so1, s11 = self._split_cov(sigma, i_known, i_unknown)
if issparse(Sigma): if issparse(sigma):
So1 = So1.todense() so1 = so1.todense()
S11 = S11.todense() s11 = s11.todense()
S1o_Sooinv = spsolve(Soo + Soo.T, 2 * So1).T s1o_sooinv = spsolve(soo + soo.T, 2 * so1).T
else: else:
Sooinv_So1, _res, _rank, _s = lstsq(Soo + Soo.T, 2 * So1, sooinv_so1, _res, _rank, _s = lstsq(soo + soo.T, 2 * so1,
cond=1e-4) cond=1e-4)
S1o_Sooinv = Sooinv_So1.T s1o_sooinv = sooinv_so1.T
mu1o = S1o_Sooinv.dot(x[i_known]) mu1o = s1o_sooinv.dot(x[i_known])
Sigma1o = S11 - S1o_Sooinv.dot(So1) sigma1o = s11 - s1o_sooinv.dot(so1)
if (diag(Sigma1o) < 0).any(): if (diag(sigma1o) < 0).any():
raise ValueError('Failed to converge to a solution') raise ValueError('Failed to converge to a solution')
mu1o_std = sqrt(diag(Sigma1o)) mu1o_std = sqrt(diag(sigma1o))
sample[:] = rndnormnd(mu1o, Sigma1o, cases=1).ravel() sample[:] = rndnormnd(mu1o, sigma1o, cases=1).ravel()
elif method.startswith('appr'): elif method.startswith('appr'):
# approximating by only condition on the closest points # approximating by only condition on the closest points
Nsig = min(2 * num_acf, num_x) num_sig = min(2 * num_acf, num_x)
Sigma = toeplitz(hstack((acf, zeros(Nsig - num_acf)))) sigma = toeplitz(hstack((acf, zeros(num_sig - num_acf))))
overlap = int(Nsig / 4) overlap = int(num_sig / 4)
# indices to the points used # indices to the points used
idx = r_[0:Nsig] + max(0, min(i_unknown[0] - overlap, idx = r_[0:num_sig] + max(0, min(i_unknown[0] - overlap,
num_x - Nsig)) num_x - num_sig))
mask_unknown = zeros(num_x, dtype=bool) mask_unknown = zeros(num_x, dtype=bool)
# temporary storage of indices to missing points # temporary storage of indices to missing points
mask_unknown[i_unknown] = True mask_unknown[i_unknown] = True
@ -637,29 +639,29 @@ class CovData1D(PlotData):
x2 = x.copy() x2 = x.copy()
while ns > 0: while ns > 0:
Soo, So1, S11 = self._split_cov(Sigma, t_known, t_unknown) soo, so1, s11 = self._split_cov(sigma, t_known, t_unknown)
if issparse(Soo): if issparse(soo):
So1 = So1.todense() so1 = so1.todense()
S11 = S11.todense() s11 = s11.todense()
S1o_Sooinv = spsolve(Soo + Soo.T, 2 * So1).T s1o_sooinv = spsolve(soo + soo.T, 2 * so1).T
else: else:
Sooinv_So1, _res, _rank, _s = lstsq(Soo + Soo.T, 2 * So1, sooinv_so1, _res, _rank, _s = lstsq(soo + soo.T, 2 * so1,
cond=1e-4) cond=1e-4)
S1o_Sooinv = Sooinv_So1.T s1o_sooinv = sooinv_so1.T
Sigma1o = S11 - S1o_Sooinv.dot(So1) sigma1o = s11 - s1o_sooinv.dot(so1)
if (diag(Sigma1o) < 0).any(): if (diag(sigma1o) < 0).any():
raise ValueError('Failed to converge to a solution') raise ValueError('Failed to converge to a solution')
ix = slice((num_restored), (num_restored + ns)) ix = slice((num_restored), (num_restored + ns))
# standard deviation of the expected surface # standard deviation of the expected surface
mu1o_std[ix] = np.maximum(mu1o_std[ix], sqrt(diag(Sigma1o))) mu1o_std[ix] = np.maximum(mu1o_std[ix], sqrt(diag(sigma1o)))
# expected surface conditioned on the closest known # expected surface conditioned on the closest known
# observations from x # observations from x
mu1o[ix] = S1o_Sooinv.dot(x2[idx[t_known]]) mu1o[ix] = s1o_sooinv.dot(x2[idx[t_known]])
# sample conditioned on the known observations from x # sample conditioned on the known observations from x
mu1os = S1o_Sooinv.dot(x[idx[t_known]]) mu1os = s1o_sooinv.dot(x[idx[t_known]])
sample[ix] = rndnormnd(mu1os, Sigma1o, cases=1) sample[ix] = rndnormnd(mu1os, sigma1o, cases=1)
if idx[-1] == num_x - 1: if idx[-1] == num_x - 1:
ns = 0 # no more points to simulate ns = 0 # no more points to simulate
else: else:
@ -732,7 +734,7 @@ def main():
if __name__ == '__main__': if __name__ == '__main__':
if False: # True: # if False: # True: #
import doctest from wafo.testing import test_docstrings
doctest.testmod() test_docstrings(__file__)
else: else:
main() main()

Write Preview
Loading…
Cancel
Save