diff --git a/wafo/SpecData1D.mm b/wafo/SpecData1D.mm
deleted file mode 100644
index f70b0c8..0000000
--- a/wafo/SpecData1D.mm
+++ /dev/null
@@ -1,69 +0,0 @@
-
diff --git a/wafo/autumn.gif b/wafo/autumn.gif
deleted file mode 100644
index fa351dd..0000000
Binary files a/wafo/autumn.gif and /dev/null differ
diff --git a/wafo/bitwise.py b/wafo/bitwise.py
deleted file mode 100644
index beb72fd..0000000
--- a/wafo/bitwise.py
+++ /dev/null
@@ -1,97 +0,0 @@
-"""
-Module extending the bitoperator capabilites of numpy
-"""
-
-from numpy import (bitwise_and, bitwise_or,
- bitwise_not, binary_repr, # @UnusedImport
- bitwise_xor, where, arange) # @UnusedImport
-__all__ = ['bitwise_and', 'bitwise_or', 'bitwise_not', 'binary_repr',
- 'bitwise_xor', 'getbit', 'setbit', 'getbits', 'setbits']
-
-
-def getbit(i, bit):
- """
- Get bit at specified position
-
- Parameters
- ----------
- i : array-like of uints, longs
- value to
- bit : array-like of ints or longs
- bit position between 0 and the number of bits in the uint class.
-
- Examples
- --------
- >>> import numpy as np
- >>> binary_repr(13)
- '1101'
- >>> getbit(13,np.arange(3,-1,-1))
- array([1, 1, 0, 1])
- >>> getbit(5, np.r_[0:4])
- array([1, 0, 1, 0])
- """
- return bitwise_and(i, 1 << bit) >> bit
-
-
-def getbits(i, numbits=8):
- """
- Returns bits of i in a list
- """
- return getbit(i, arange(0, numbits))
-
-
-def setbit(i, bit, value=1):
- """
- Set bit at specified position
-
- Parameters
- ----------
- i : array-like of uints, longs
- value to
- bit : array-like of ints or longs
- bit position between 0 and the number of bits in the uint class.
- value : array-like of 0 or 1
- value to set the bit to.
-
- Examples
- --------
- Set bit fifth bit in the five bit binary binary representation of 9 (01001)
- yields 25 (11001)
- >>> setbit(9,4)
- array(25)
- """
- val1 = 1 << bit
- val0 = bitwise_not(val1)
- return where((value == 0) & (i == i) & (bit == bit), bitwise_and(i, val0),
- bitwise_or(i, val1))
-
-
-def setbits(bitlist):
- """
- Set bits of val to values in bitlist
-
- Example
- -------
- >>> setbits([1,1])
- 3
- >>> setbits([1,0])
- 1
- """
- val = 0
- for i, j in enumerate(bitlist):
- val |= j << i
- return val
-
-if __name__ == '__main__':
- from wafo.testing import test_docstrings
- test_docstrings(__file__)
-
-# t = set(np.arange(8),1,1)
-# t=get(0x84,np.arange(0,8))
-# t=getbyte(0x84)
-# t=get(0x84,[0, 1, 2, 3, 4, 5, 6, 7])
-# t=get(0x20, 6)
-# bit = [0 for i in range(8)]
-# bit[7]=1
-# t = setbits(bit)
-# print(hex(t))
diff --git a/wafo/conftest.py b/wafo/conftest.py
deleted file mode 100644
index 6bef122..0000000
--- a/wafo/conftest.py
+++ /dev/null
@@ -1,3 +0,0 @@
-# TODO Fix doctests in fig.py
-collect_ignore = ["fig.py", "MSO.py", "MSPPT.py", "powerpoint.py",
- "win32_utils.py"]
diff --git a/wafo/definitions.py b/wafo/definitions.py
deleted file mode 100644
index d75e648..0000000
--- a/wafo/definitions.py
+++ /dev/null
@@ -1,310 +0,0 @@
-"""
-WAFO defintions and numenclature
-
- crossings :
- cycle_pairs :
- turning_points :
- wave_amplitudes :
- wave_periods :
- waves :
-
-Examples
---------
-In order to view the documentation do the following in an ipython window:
-
-import wafo.definitions as wd
-wd.crossings()
-
-or
-
-wd.crossings?
-
-
-"""
-
-
-def wave_amplitudes():
- r"""
- Wave amplitudes and heights definitions and nomenclature
-
- Definition of wave amplitudes and wave heights
- ---------------------------------------------
-
- <----- Direction of wave propagation
-
-
- |..............c_..........|
- | /| \ |
- Hd | _/ | \ | Hu
- M | / | \ |
- / \ | M / Ac | \_ | c_
- F \ | / \m/ | \ | / \
- ------d----|---u------------------d---|---u----d------ level v
- \ | /| \ | / \L
- \_ | / | At \_|_/
- \|/..| t
- t
-
- Parameters
- ----------
- Ac : crest amplitude
- At : trough amplitude
- Hd : wave height as defined for down crossing waves
- Hu : wave height as defined for up crossing waves
-
- See also
- --------
- waves, crossings, turning_points
- """
- print(wave_amplitudes.__doc__)
-
-
-def crossings():
- r"""
- Level v crossing definitions and nomenclature
-
- Definition of level v crossings
- -------------------------------
- M
- . . M M
- . . . . . .
- F d . . L
- -----------------------u-------d-------o----------------- level v
- . . . . u
- . m
- m
-
- Let the letters 'm', 'M', 'F', 'L','d' and 'u' in the
- figure above denote local minimum, maximum, first value, last
- value, down- and up-crossing, respectively. The remaining
- sampled values are indicated with a '.'. Values that are identical
- with v, but do not cross the level is indicated with the letter 'o'.
- We have a level up-crossing at index, k, if
-
- x(k) < v and v < x(k+1)
- or if
- x(k) == v and v < x(k+1) and x(r) < v for some di < r <= k-1
-
- where di is the index to the previous downcrossing.
- Similarly there is a level down-crossing at index, k, if
-
- x(k) > v and v > x(k+1)
- or if
- x(k) == v and v > x(k+1) and x(r) > v for some ui < r <= k-1
-
- where ui is the index to the previous upcrossing.
-
- The first (F) value is a up crossing if x(1) = v and x(2) > v.
- Similarly, it is a down crossing if x(1) = v and x(2) < v.
-
- See also
- --------
- wave_periods, waves, turning_points, findcross, findtp
- """
- print(crossings.__doc__)
-
-
-def cycle_pairs():
- r"""
- Cycle pairs definitions and numenclature
-
- Definition of Max2min and min2Max cycle pair
- --------------------------------------------
- A min2Max cycle pair (mM) is defined as the pair of a minimum
- and the following Maximum. Similarly a Max2min cycle pair (Mm)
- is defined as the pair of a Maximum and the following minimum.
- (all turning points possibly rainflowfiltered before pairing into cycles.)
-
- See also
- --------
- turning_points
- """
- print(cycle_pairs.__doc__)
-
-
-def wave_periods():
- r"""
- Wave periods (lengths) definitions and nomenclature
-
- Definition of wave periods (lengths)
- ------------------------------------
-
-
- <----- Direction of wave propagation
-
- <-------Tu--------->
- : :
- <---Tc-----> :
- : : : <------Tcc---->
- M : c : : : :
- / \ : M / \_ : : c_ c
- F \ :/ \m/ \: :/ \ / \
- ------d--------u----------d-------u----d--------u---d-------- level v
- \ / \ / :\_ _/: :\_ L
- \_ / \_t_/ : \t_/ : : \m/
- \t/ : : : :
- : : <---Tt---> :
- <--------Ttt-------> : :
- <-----Td----->
- Tu = Up crossing period
- Td = Down crossing period
- Tc = Crest period, i.e., period between up crossing and
- the next down crossing
- Tt = Trough period, i.e., period between down crossing and
- the next up crossing
- Ttt = Trough2trough period
- Tcc = Crest2crest period
-
-
- <----- Direction of wave propagation
-
- <--Tcf-> Tuc
- : : <-Tcb-> <->
- M : c : : : :
- / \ : M / \_ c_ : : c
- F \ :/ \m/ \ / \___: :/ \
- ------d---------u----------d---------u-------d--------u---d------ level v
- :\_ / \ __/: \_ _/ \_ L
- : \_ / \_t_/ : \t_/ \m/
- : \t/ : :
- : : : :
- <-Ttf-> <-Ttb->
-
-
- Tcf = Crest front period, i.e., period between up crossing and crest
- Tcb = Crest back period, i.e., period between crest and down crossing
- Ttf = Trough front period, i.e., period between down crossing and trough
- Ttb = Trough back period, i.e., period between trough and up crossing
- Also note that Tcf and Ttf can also be abbreviated by their crossing
- marker, e.g. Tuc (u2c) and Tdt (d2t), respectively. Similar applies
- to all the other wave periods and wave lengths.
-
- (The nomenclature for wave length is similar, just substitute T and
- period with L and length, respectively)
-
- <----- Direction of wave propagation
-
- <--TMm-->
- <-TmM-> : :
- M : : M :
- / \ : M /:\_ : M_ M
- F \ : / \m/ : \ : /: \ / \
- \ : / : \ : / : \ / \
- \ : / : \ : / : \_ _/ \_ L
- \_ : / : \_m_/ : \m_/ \m/
- \m/ : : : :
- <-----TMM-----> <----Tmm----->
-
-
- TmM = Period between minimum and the following Maximum
- TMm = Period between Maximum and the following minimum
- TMM = Period between Maximum and the following Maximum
- Tmm = Period between minimum and the following minimum
-
- See also
- --------
- waves,
- wave_amplitudes,
- crossings,
- turning_points
- """
- print(wave_periods.__doc__)
-
-
-def turning_points():
- r"""
- Turning points definitions and numenclature
-
- Definition of turningpoints
- ---------------------------
- <----- Direction of wave propagation
-
- M M
- / \ .... M /:\_ M_ M
- F \ | / \m/ : \ /: \ / \
- \ h | / : \ / : \ / \
- \ | / : \ / : \_ _/ \_ L
- \_ | / : \_m_/ : \m_/ \m/
- \m/ : : : :
- <------Mw-----> <-----mw----->
-
- Local minimum or maximum are indicated with the
- letters 'm' or 'M'. Turning points in this connection are all
- local max (M) and min (m) and the last (L) value and the
- first (F) value if the first local extremum is a max.
-
- (This choice is made in order to get the exact up-crossing intensity
- from rfc by mm2lc(tp2mm(rfc)) )
-
-
- See also
- --------
- waves,
- crossings,
- cycle_pairs
- findtp
-
- """
- print(turning_points.__doc__)
-
-
-def waves():
- r"""
- Wave definitions and nomenclature
-
- Definition of trough and crest
- ------------------------------
- A trough (t) is defined as the global minimum between a
- level v down-crossing (d) and the next up-crossing (u)
- and a crest (c) is defined as the global maximum between a
- level v up-crossing and the following down-crossing.
-
- Definition of down- and up -crossing waves
- ------------------------------------------
- A level v-down-crossing wave (dw) is a wave from a
- down-crossing to the following down-crossing.
- Similarly, a level v-up-crossing wave (uw) is a wave from an up-crossing
- to the next up-crossing.
-
- Definition of trough and crest waves
- ------------------------------------
- A trough-to-trough wave (tw) is a wave from a trough (t) to the
- following trough. The crest-to-crest wave (cw) is defined similarly.
-
-
- Definition of min2min and Max2Max wave
- --------------------------------------
- A min2min wave (mw) is defined starting from a minimum (m) and
- ending in the following minimum.
- Similarly a Max2Max wave (Mw) is thus a wave from a maximum (M)
- to the next maximum (all waves optionally rainflow filtered).
-
- <----- Direction of wave propagation
-
-
- <------Mw-----> <----mw---->
- M : : c :
- / \ M : / \_ : c_ c
- F \ / \m/ \ : /: \ /:\
- ------d--------u----------d-------u----d--------u---d------ level v
- \ /: \ : /: : :\_ _/ : :\_ L
- \_ / : \_t_/ : : : \t_/ : : \m/
- \t/ <-------uw---------> : <-----dw----->
- : : : :
- <--------tw--------> <------cw----->
-
- (F=first value and L=last value).
-
- See also
- --------
- turning_points,
- crossings,
- wave_periods
- findtc,
- findcross
- """
- print(waves.__doc__)
-
-if __name__ == '__main__':
- import doctest
- doctest.testmod(optionflags=doctest.NORMALIZE_WHITESPACE)
diff --git a/wafo/demos.py b/wafo/demos.py
deleted file mode 100644
index f1dfbcc..0000000
--- a/wafo/demos.py
+++ /dev/null
@@ -1,146 +0,0 @@
-"""
-Created on 20. jan. 2011
-
-@author: pab
-"""
-import numpy as np
-from numpy import exp, meshgrid
-__all__ = ['peaks', 'humps', 'magic']
-
-
-def _magic_odd_order(n):
- ix = np.arange(n) + 1
- J, I = np.meshgrid(ix, ix)
- A = np.mod(I + J - (n + 3) / 2, n)
- B = np.mod(I + 2 * J - 2, n)
- M = n * A + B + 1
- return M
-
-
-def _magic_doubly_even_order(n):
- M = np.arange(1, n * n + 1).reshape(n, n)
- ix = np.mod(np.arange(n) + 1, 4) // 2
- J, I = np.meshgrid(ix, ix)
- iz = np.flatnonzero(I == J)
- M.put(iz, n * n + 1 - M.flat[iz])
- return M
-
-
-def _magic_even_order(n):
- p = n // 2
- M0 = magic(p)
- M = np.hstack((np.vstack((M0, M0 + 3 * p * p)),
- np.vstack((M0 + 2 * p * p, M0 + p * p))))
- if n > 2:
- k = (n - 2) // 4
- jvec = np.hstack((np.arange(k), np.arange(n - k + 1, n)))
- for i in range(p):
- for j in jvec:
- temp = M[i][j]
- M[i][j] = M[i + p][j]
- M[i + p][j] = temp
-
- i = k
- j = 0
- temp = M[i][j]
- M[i][j] = M[i + p][j]
- M[i + p][j] = temp
- j = i
- temp = M[i + p][j]
- M[i + p][j] = M[i][j]
- M[i][j] = temp
- return M
-
-
-def magic(n):
- """
- Return magic square for n of any orders > 2.
-
- A magic square has the property that the sum of every row and column,
- as well as both diagonals, is the same number.
-
- Examples
- --------
- >>> np.allclose(magic(3),
- ... [[8, 1, 6],
- ... [3, 5, 7],
- ... [4, 9, 2]])
- True
-
- >>> np.allclose(magic(4),
- ... [[16, 2, 3, 13],
- ... [ 5, 11, 10, 8],
- ... [ 9, 7, 6, 12],
- ... [ 4, 14, 15, 1]])
- True
-
- >>> np.allclose(magic(6),
- ... [[35, 1, 6, 26, 19, 24],
- ... [ 3, 32, 7, 21, 23, 25],
- ... [31, 9, 2, 22, 27, 20],
- ... [ 8, 28, 33, 17, 10, 15],
- ... [30, 5, 34, 12, 14, 16],
- ... [ 4, 36, 29, 13, 18, 11]])
- True
- """
- if (n < 3):
- raise ValueError('n must be greater than 2.')
-
- if np.mod(n, 2) == 1:
- return _magic_odd_order(n)
- elif np.mod(n, 4) == 0:
- return _magic_doubly_even_order(n)
- return _magic_even_order(n)
-
-
-def peaks(x=None, y=None, n=51):
- """
- Return the "well" known MatLab (R) peaks function
- evaluated in the [-3,3] x,y range
-
- Example
- -------
- >>> import matplotlib.pyplot as plt
- >>> x,y,z = peaks()
-
- h = plt.contourf(x,y,z)
-
- """
- if x is None:
- x = np.linspace(-3, 3, n)
- if y is None:
- y = np.linspace(-3, 3, n)
-
- [x1, y1] = meshgrid(x, y)
-
- z = (3 * (1 - x1) ** 2 * exp(-(x1 ** 2) - (y1 + 1) ** 2) -
- 10 * (x1 / 5 - x1 ** 3 - y1 ** 5) * exp(-x1 ** 2 - y1 ** 2) -
- 1. / 3 * exp(-(x1 + 1) ** 2 - y1 ** 2))
-
- return x1, y1, z
-
-
-def humps(x=None):
- """
- Computes a function that has three roots, and some humps.
-
- Example
- -------
- >>> import matplotlib.pyplot as plt
- >>> x = np.linspace(0,1)
- >>> y = humps(x)
-
- h = plt.plot(x,y)
- """
- if x is None:
- y = np.linspace(0, 1)
- else:
- y = np.asarray(x)
-
- return 1.0 / ((y - 0.3) ** 2 + 0.01) + 1.0 / ((y - 0.9) ** 2 + 0.04) + \
- 2 * y - 5.2
-
-
-if __name__ == '__main__':
- from wafo.testing import test_docstrings
- test_docstrings(__file__)
diff --git a/wafo/f2py_tools.py b/wafo/f2py_tools.py
deleted file mode 100644
index 74b0e97..0000000
--- a/wafo/f2py_tools.py
+++ /dev/null
@@ -1,65 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-f2py c_library.pyf c_functions.c -c
-
-See also http://www.scipy.org/Cookbook/CompilingExtensionsOnWindowsWithMinGW
-"""
-import os
-import sys
-
-
-def which(program):
- """
- Return filepath to program if it exists
-
- In order to test if a certain executable exists, it will search for the
- program name in the environment variables.
- If program is a full path to an executable, it will check it exists
-
- Copied from:
- http://stackoverflow.com/questions/377017/test-if-executable-exists-in-python/
- It is supposed to mimic the UNIX command "which"
- """
-
- def is_exe(fpath):
- return os.path.exists(fpath) and os.access(fpath, os.X_OK)
-
- fpath, unused_fname = os.path.split(program)
- if fpath:
- if is_exe(program):
- return program
- else:
- for path in os.environ["PATH"].split(os.pathsep):
- exe_file = os.path.join(path, program)
- if is_exe(exe_file):
- return exe_file
-
- return None
-
-
-def f2py_call_str():
- '''Return which f2py callable is in the path regardless of platform'''
-
- # define possible options:
- # on Arch Linux, python and f2py are the calls corresponding to python 3
- # and python2/f2py2 for python 2
- # other Linux versions might still use python/f2py for python 2
-
- if os.path.basename(sys.executable).endswith('2'):
- options = ('f2py2', 'f2py2.6', 'f2py2.7',)
- else: # on Windows and other Linux using python/f2py
- options = ('f2py.exe', 'f2py.bat', 'f2py', 'f2py2.6', 'f2py2.7',
- 'f2py.py',)
- for k in options:
- if which(k):
- # Found the f2py path, no need to look further
- f2py_call = k
- f2py_path = which(k)
- break
-
- try:
- print('found f2py in:', f2py_path)
- return f2py_call
- except NameError:
- raise UserWarning('Couldn\'t locate f2py. '
- 'Should be part of NumPy installation.')
diff --git a/wafo/fig.py b/wafo/fig.py
deleted file mode 100644
index f700103..0000000
--- a/wafo/fig.py
+++ /dev/null
@@ -1,908 +0,0 @@
-# /usr/bin/env python
-'''
-Module FIG
-------------
-Module for manipulating windows/figures created using
-pylab or enthought.mayavi.mlab on the windows platform.
-
-Figure manipulation involves
-maximization, minimization, hiding, closing, stacking or tiling.
-
-It is assumed that the figures are uniquely numbered in the following way:
-Figure 1
-Figure 2
-....
-or
-TVTK scene 1
-TVTK scene 2
-TVTK scene 3
-...
-
-Example
--------
->>> import pylab as p
->>> import wafo.fig as fig
->>> for ix in range(6):
-... f = p.figure(ix)
->>> fig.stack('all')
->>> fig.stack(1,2)
->>> fig.hide(1)
->>> fig.restore(1)
->>> fig.tile()
->>> fig.pile()
->>> fig.maximize(4)
->>> fig.close('all')
-'''
-
-from __future__ import absolute_import, division, print_function
-# import win32api
-import win32gui
-import win32con
-import wx
-import numpy
-
-from win32gui import (EnumWindows, MoveWindow, GetWindowRect, FindWindow,
- ShowWindow, BringWindowToTop)
-
-__all__ = ['close', 'cycle', 'hide', 'keep', 'maximize', 'minimize', 'pile',
- 'restore', 'stack', 'tile', 'find_all_figure_numbers', 'set_size']
-
-# Figure format strings to recognize in window title
-FIGURE_TITLE_FORMATS = ('Figure', 'TVTK Scene', 'Chaco Plot Window: Figure')
-_SCREENSIZE = None
-
-
-class CycleDialog(wx.Dialog):
-
- def _get_buttons(self):
- hbox = wx.BoxSizer(wx.HORIZONTAL)
- buttons = ['Forward', 'Back', 'Cancel']
- callbacks = [self.on_forward, self.on_backward, self.on_cancel]
- for button, callback in zip(buttons, callbacks):
- button = wx.Button(self, -1, button, size=(70, 30))
- self.Bind(wx.EVT_BUTTON, callback, button)
- hbox.Add(button, 1, wx.ALIGN_CENTER)
- return hbox
-
- def _get_message(self):
- label = ('Press back or forward to display previous or next figure(s),'
- ' respectively. Press cancel to quit.')
- message = wx.StaticText(self, label=label, size=(240, 25))
- return message
-
- def __init__(self, parent, interval=None, title='Cycle dialog'):
- super(CycleDialog, self).__init__(parent, title=title, size=(260, 130))
- if isinstance(interval, (float, int)):
- self.interval_milli_sec = interval * 1000
- else:
- self.interval_milli_sec = 30
-
- self.timer = wx.Timer(self)
- self.Bind(wx.EVT_TIMER, self.on_forward, self.timer)
-
- vbox = wx.BoxSizer(wx.VERTICAL)
- vbox.Add(self._get_message(), 0, wx.ALIGN_CENTER | wx.TOP, 20)
- vbox.Add(self._get_buttons(), 1, wx.ALIGN_CENTER | wx.TOP | wx.BOTTOM, 10)
- self.SetSizer(vbox)
-
- def ShowModal(self, *args, **kwargs):
- self.timer.Start(self.interval_milli_sec, oneShot=True)
- return super(CycleDialog, self).ShowModal(*args, **kwargs)
-
- def on_forward(self, evt):
- self.EndModal(wx.ID_FORWARD)
-
- def on_backward(self, evt):
- self.EndModal(wx.ID_BACKWARD)
-
- def on_cancel(self, evt):
- self.EndModal(wx.ID_CANCEL)
-
-
-def _get_cycle_dialog(parent=None, interval=None):
- app = wx.GetApp()
- if not app:
- app = wx.App(redirect=False)
- frame = wx.Frame(None)
- app.SetTopWindow(frame)
- dlg = CycleDialog(parent, interval)
- return dlg
-
-
-def get_window_position_and_size(window_handle):
- pos = GetWindowRect(window_handle)
- return pos[0], pos[1], pos[2] - pos[0], pos[3] - pos[1]
-
-
-def get_screen_position_and_size(window_handles):
- """Return screen position; X, Y and size; width, height.
-
- Parameters
- ----------
- window_handles: list of handles to open window figures
- (Note: only needed the first time)
-
- Returns
- --------
- X : coordinate of the left side of the screen.
- Y : coordinate of the top of the screen.
- width : screen horizontal size
- height : screen vertical size
-
- """
- # pylint: disable=global-statement
- global _SCREENSIZE
- if _SCREENSIZE is None:
- window_handle = window_handles[0]
- pos = get_window_position_and_size(window_handle)
- _show_windows((window_handle,), win32con.SW_SHOWMAXIMIZED)
- _SCREENSIZE = get_window_position_and_size(window_handle)
- MoveWindow(window_handle, pos[0], pos[1], pos[2], pos[3], 1)
- return _SCREENSIZE
-
-
-def _get_screen_size(wnds):
- screen_width, screen_height = get_screen_position_and_size(wnds)[2:4]
- return screen_width, screen_height
-
-
-def _windowEnumerationHandler(handle, result_list):
- """Pass to win32gui.EnumWindows() to generate list of window handle, window
- text tuples."""
- # pylint: disable=no-member
- if win32gui.IsWindowVisible(handle):
- result_list.append((handle, win32gui.GetWindowText(handle)))
-
-
-def _find_window_handles_and_titles(wantedTitle=None):
- """Return list of window handle and window title tuples.
-
- Parameter
- ---------
- wantedTitle:
-
- """
- handles_n_titles = []
- EnumWindows(_windowEnumerationHandler, handles_n_titles)
- if wantedTitle is None:
- return handles_n_titles
- else:
- return [(handle, title)
- for handle, title in handles_n_titles
- if title.startswith(wantedTitle)]
-
-
-def find_figure_handles(*figure_numbers):
- """Find figure handles from figure numbers."""
- wnd_handles = []
- for figure_number in _parse_figure_numbers(*figure_numbers):
- for format_ in FIGURE_TITLE_FORMATS:
- winTitle = format_ + ' %d' % figure_number
- handle = FindWindow(None, winTitle)
- if not handle == 0:
- wnd_handles.append(handle)
- return wnd_handles
-
-
-def find_all_figure_numbers():
- """Return list of all figure numbers.
-
- Example
- -------
- >>> import fig
- >>> import pylab as p
- >>> for ix in range(5):
- ... f = p.figure(ix)
- ... p.draw()
-
- fig.find_all_figure_numbers()
- [0, 1, 2, 3, 4]
-
- >>> fig.close()
-
- """
- figure_numbers = []
- for wantedTitle in FIGURE_TITLE_FORMATS:
- handles_n_titles = _find_window_handles_and_titles(wantedTitle)
- for _handle, title in handles_n_titles:
- try:
- number = int(title.split()[-1])
- figure_numbers.append(number)
- except (TypeError, ValueError):
- pass
- # pylint: disable=no-member
- return numpy.unique(figure_numbers).tolist()
-
-
-def _parse_figure_numbers(*args):
- figure_numbers = []
- for arg in args:
- if isinstance(arg, (list, tuple, set)):
- for val in arg:
- figure_numbers.append(int(val))
- elif isinstance(arg, int):
- figure_numbers.append(arg)
- elif arg == 'all':
- figure_numbers = find_all_figure_numbers()
- break
- else:
- raise TypeError('Only integers arguments accepted!')
-
- if len(figure_numbers) == 0:
- figure_numbers = find_all_figure_numbers()
- return figure_numbers
-
-
-def _show_figure(figure_numbers, command):
- """Sets the specified figure's show state.
-
- Parameters
- ----------
- figure_numbers: list of figure numbers
- command: one of following commands:
- SW_FORCEMINIMIZE:
- Minimizes a window, even if the thread that owns the window is not
- responding. This flag should only be used when minimizing windows
- from a different thread.
- SW_HIDE:
- Hides the window and activates another window.
- SW_MAXIMIZE:
- Maximizes the specified window.
- SW_MINIMIZE:
- Minimizes the specified window and activates the next top-level window
- in the Z order.
- SW_RESTORE:
- Activates and displays the window. If the window is minimized or
- maximized, the system restores it to its original size and position.
- An application should specify this flag when restoring a minimized
- window.
- SW_SHOW:
- Activates the window and displays it in its current size and position.
- SW_SHOWDEFAULT:
- Sets the show state based on the SW_ value specified in the STARTUPINFO
- structure passed to the CreateProcess function by the program that
- started the application.
- SW_SHOWMAXIMIZED:
- Activates the window and displays it as a maximized window.
- SW_SHOWMINIMIZED:
- Activates the window and displays it as a minimized window.
- SW_SHOWMINNOACTIVE:
- Displays the window as a minimized window. This value is similar to
- SW_SHOWMINIMIZED, except the window is not activated.
- SW_SHOWNA:
- Displays the window in its current size and position. This value is
- similar to SW_SHOW, except the window is not activated.
- SW_SHOWNOACTIVATE:
- Displays a window in its most recent size and position. This value is
- similar to SW_SHOWNORMAL, except the window is not actived.
- SW_SHOWNORMAL:
- Activates and displays a window. If the window is minimized or
- maximized, the system restores it to its original size and position.
- An application should specify this flag when displaying the window for
- the first time.
-
- """
- for number in _parse_figure_numbers(*figure_numbers):
- for format_ in FIGURE_TITLE_FORMATS:
- title = format_ + ' %d' % number
- handle = FindWindow(None, title)
- if not handle == 0:
- BringWindowToTop(handle)
- ShowWindow(handle, command)
-
-
-def _show_windows(handles, command, redraw_now=False):
- """Sets the specified window's show state.
-
- Parameters
- ----------
- handles: list of window handles
- command: one of following commands:
- SW_FORCEMINIMIZE:
- Minimizes a window, even if the thread that owns the window is not
- responding. This flag should only be used when minimizing windows
- from a different thread.
- SW_HIDE:
- Hides the window and activates another window.
- SW_MAXIMIZE:
- Maximizes the specified window.
- SW_MINIMIZE:
- Minimizes the specified window and activates the next top-level window
- in the Z order.
- SW_RESTORE:
- Activates and displays the window. If the window is minimized or
- maximized, the system restores it to its original size and position.
- An application should specify this flag when restoring a minimized
- window.
- SW_SHOW:
- Activates the window and displays it in its current size and position.
- SW_SHOWDEFAULT:
- Sets the show state based on the SW_ value specified in the STARTUPINFO
- structure passed to the CreateProcess function by the program that
- started the application.
- SW_SHOWMAXIMIZED:
- Activates the window and displays it as a maximized window.
- SW_SHOWMINIMIZED:
- Activates the window and displays it as a minimized window.
- SW_SHOWMINNOACTIVE:
- Displays the window as a minimized window. This value is similar to
- SW_SHOWMINIMIZED, except the window is not activated.
- SW_SHOWNA:
- Displays the window in its current size and position. This value is
- similar to SW_SHOW, except the window is not activated.
- SW_SHOWNOACTIVATE:
- Displays a window in its most recent size and position. This value is
- similar to SW_SHOWNORMAL, except the window is not actived.
- SW_SHOWNORMAL:
- Activates and displays a window. If the window is minimized or
- maximized, the system restores it to its original size and position.
- An application should specify this flag when displaying the window for
- the first time.
-
- redraw_now :
-
- """
- # pylint: disable=no-member
- for handle in handles:
- if not handle == 0:
- BringWindowToTop(handle)
- ShowWindow(handle, command)
- if redraw_now:
- rect = GetWindowRect(handle)
- win32gui.RedrawWindow(handle, rect, None, win32con.RDW_UPDATENOW)
-
-
-def keep(*figure_numbers):
- """Keeps figure windows of your choice and closes the rest.
-
- Parameters
- ----------
- figure_numbers : list of integers specifying which figures to keep.
-
- Example:
- --------
- # keep only figures 1,2,3,5 and 7
- >>> import pylab as p
- >>> import wafo.fig as fig
- >>> for ix in range(10):
- ... f = p.figure(ix)
- >>> fig.keep( range(1,4), 5, 7)
-
- or
- fig.keep([range(1,4), 5, 7])
- >>> fig.close()
-
- See also
- --------
- fig.close
-
- """
- figs2keep = []
- for fig in figure_numbers:
- if isinstance(fig, (list, tuple, set)):
- for val in fig:
- figs2keep.append(int(val))
- elif isinstance(fig, int):
- figs2keep.append(fig)
- else:
- raise TypeError('Only integers arguments accepted!')
-
- if len(figs2keep) > 0:
- allfigs = set(find_all_figure_numbers())
- figs2delete = allfigs.difference(figs2keep)
- close(figs2delete)
-
-
-def close(*figure_numbers):
- """ Close figure window(s)
-
- Parameters
- ----------
- figure_numbers : list of integers or string
- specifying which figures to close (default 'all').
-
- Examples
- --------
- >>> import pylab as p
- >>> import wafo.fig as fig
- >>> for ix in range(5):
- ... f = p.figure(ix)
- >>> fig.close(3,4) # close figure 3 and 4
- >>> fig.close('all') # close all remaining figures
-
- or even simpler
- fig.close() # close all remaining figures
-
- See also
- --------
- fig.keep
-
- """
- # pylint: disable=no-member
- for handle in find_figure_handles(*figure_numbers):
- if win32gui.SendMessage(handle, win32con.WM_CLOSE, 0, 0):
- win32gui.SendMessage(handle, win32con.WM_DESTROY, 0, 0)
-
-
-def restore(*figure_numbers):
- """Restore figures window size and position to its default value.
-
- Parameters
- ----------
- figure_numbers : list of integers or string
- specifying which figures to restor (default 'all').
-
- Description
- -----------
- RESTORE Activates and displays the window. If the window is minimized
- or maximized, the system restores it to its original size and position.
-
- Examples
- ---------
- >>> import pylab as p
- >>> import wafo.fig as fig
- >>> for ix in range(5):
- ... f = p.figure(ix)
- >>> fig.restore('all') #Restores all figures
- >>> fig.restore() #same as restore('all')
- >>> fig.restore(p.gcf().number) #Restores the current figure
- >>> fig.restore(3) #Restores figure 3
- >>> fig.restore([2, 4]) #Restores figures 2 and 4
-
- or alternatively
- fig.restore(2, 4)
- >>> fig.close()
-
- See also
- --------
- fig.close,
- fig.keep
-
- """
- SW_RESTORE = win32con.SW_RESTORE
- # SW_RESTORE = win32con.SW_SHOWDEFAULT
- # SW_RESTORE = win32con.SW_SHOWNORMAL
- _show_figure(figure_numbers, SW_RESTORE)
-
-
-def hide(*figure_numbers):
- """hide figure(s) window.
-
- Parameters
- ----------
- figure_numbers : list of integers or string
- specifying which figures to hide (default 'all').
-
- Examples:
- --------
- >>> import wafo.fig as fig
- >>> import pylab as p
- >>> for ix in range(5):
- ... f = p.figure(ix)
- >>> fig.hide('all') #hides all unhidden figures
- >>> fig.hide() #same as hide('all')
- >>> fig.hide(p.gcf().number) #hides the current figure
- >>> fig.hide(3) #hides figure 3
- >>> fig.hide([2, 4]) #hides figures 2 and 4
-
- or alternatively
- fig.hide(2, 4)
- >>> fig.restore(list(range(5)))
- >>> fig.close()
-
- See also
- --------
- fig.cycle,
- fig.keep,
- fig.restore
-
- """
- _show_figure(figure_numbers, win32con.SW_HIDE)
-
-
-def minimize(*figure_numbers):
- """Minimize figure(s) window size.
-
- Parameters
- ----------
- figure_numbers : list of integers or string
- specifying which figures to minimize (default 'all').
-
- Examples:
- ---------
- >>> import wafo.fig as fig
- >>> import pylab as p
- >>> for ix in range(5):
- ... f = p.figure(ix)
- >>> fig.minimize('all') #Minimizes all unhidden figures
- >>> fig.minimize() #same as minimize('all')
- >>> fig.minimize(p.gcf().number) #Minimizes the current figure
- >>> fig.minimize(3) #Minimizes figure 3
- >>> fig.minimize([2, 4]) #Minimizes figures 2 and 4
-
- or alternatively
- fig.minimize(2, 4)
- >>> fig.close()
-
- See also
- --------
- fig.cycle,
- fig.keep,
- fig.restore
-
- """
- _show_figure(figure_numbers, win32con.SW_SHOWMINIMIZED)
-
-
-def maximize(*figure_numbers):
- """Maximize figure(s) window size.
-
- Parameters
- ----------
- figure_numbers : list of integers or string
- specifying which figures to maximize (default 'all').
-
- Examples:
- ---------
- >>> import pylab as p
- >>> import wafo.fig as fig
- >>> for ix in range(5):
- ... f = p.figure(ix)
- >>> fig.maximize('all') #Maximizes all unhidden figures
- >>> fig.maximize() #same as maximize('all')
- >>> fig.maximize(p.gcf().number) #Maximizes the current figure
- >>> fig.maximize(3) #Maximizes figure 3
- >>> fig.maximize([2, 4]) #Maximizes figures 2 and 4
-
- or alternatively
- fig.maximize(2, 4)
- >>> fig.close()
-
- See also
- --------
- fig.cycle,
- fig.keep,
- fig.restore
-
- """
- _show_figure(figure_numbers, win32con.SW_SHOWMAXIMIZED)
-
-
-def pile(*figure_numbers, **kwds):
- """Pile figure windows.
-
- Parameters
- ----------
- figure_numbers : list of integers or string
- specifying which figures to pile (default 'all').
- kwds : dict with the following keys
- position :
- width :
- height :
-
- Description
- -------------
- PILE piles all open figure windows on top of eachother
- with complete overlap. PILE(FIGS) can be used to specify which
- figures that should be piled. Figures are not sorted when specified.
-
- Example:
- --------
- >>> import pylab as p
- >>> import wafo.fig as fig
- >>> for ix in range(7):
- ... f = p.figure(ix)
- >>> fig.pile() # pile all open figures
- >>> fig.pile(range(1,4), 5, 7) # pile figure 1,2,3,5 and 7
- >>> fig.close()
-
- See also
- --------
- fig.cycle, fig.keep, fig.maximize, fig.restore,
- fig.stack, fig.tile
-
- """
- wnds = find_figure_handles(*figure_numbers)
- numfigs = len(wnds)
- if numfigs > 0:
- screen_width, screen_height = _get_screen_size(wnds)
- pos = kwds.get(
- 'position', (int(screen_width / 5), int(screen_height / 4)))
- width = kwds.get('width', int(screen_width / 2.5))
- height = kwds.get('height', int(screen_height / 2))
-
- for wnd in wnds:
- MoveWindow(wnd, pos[0], pos[1], width, height, 1)
- BringWindowToTop(wnd)
-
-
-def set_size(*figure_numbers, **kwds):
- """Set size for figure windows.
-
- Parameters
- ----------
- figure_numbers : list of integers or string
- specifying which figures to pile (default 'all').
- kwds : dict with the following keys
- width :
- height :
-
- Description
- -------------
- Set size sets the size of all open figure windows. SET_SIZE(FIGS)
- can be used to specify which figures that should be resized.
- Figures are not sorted when specified.
-
- Example:
- --------
- >>> import pylab as p
- >>> import fig
- >>> for ix in range(7):
- ... f = p.figure(ix)
- >>> fig.set_size(7, width=150, height=100)
- >>> fig.set_size(range(1,4), 5,width=250, height=170)
- >>> fig.close()
-
- See also
- --------
- fig.cycle, fig.keep, fig.maximize, fig.restore,
- fig.stack, fig.tile
-
- """
- handles = find_figure_handles(*figure_numbers)
- numfigs = len(handles)
- if numfigs > 0:
- screen_width, screen_height = _get_screen_size(handles)
- width = kwds.get('width', int(screen_width / 2.5))
- height = kwds.get('height', int(screen_height / 2))
- new_pos = kwds.get('position', None)
- pos = new_pos
- for handle in handles:
- if not new_pos:
- pos = get_window_position_and_size(handle)
- MoveWindow(handle, pos[0], pos[1], width, height, 1)
- BringWindowToTop(handle)
-
-
-def stack(*figure_numbers, **kwds):
- """Stack figure windows.
-
- Parameters
- ----------
- figure_numbers : list of integers or string
- specifying which figures to stack (default 'all').
- kwds : dict with the following keys
- figs_per_stack :
- number of figures per stack (default depends on screenheight)
-
- Description
- -----------
- STACK stacks all open figure windows on top of eachother
- with maximum overlap. STACK(FIGS) can be used to specify which
- figures that should be stacked. Figures are not sorted when specified.
-
- Example:
- --------
- >>> import pylab as p
- >>> import wafo.fig as fig
- >>> for ix in range(7):
- ... f = p.figure(ix)
- >>> fig.stack() # stack all open figures
- >>> fig.stack(range(1,4), 5, 7) # stack figure 1,2,3,5 and 7
- >>> fig.close()
-
- See also
- --------
- fig.cycle, fig.keep, fig.maximize, fig.restore,
- fig.pile, fig.tile
-
- """
- wnds = find_figure_handles(*figure_numbers)
- numfigs = len(wnds)
- if numfigs > 0:
- screenpos = get_screen_position_and_size(wnds)
- y_step = 25
- x_step = border = 5
-
- figs_per_stack = kwds.get(
- 'figs_per_stack',
- int(numpy.fix(0.7 * (screenpos[3] - border) / y_step)))
-
- for iy in range(numfigs):
- pos = get_window_position_and_size(wnds[iy])
- # print('[x, y, w, h] = ', pos)
- ix = iy % figs_per_stack
- ypos = int(screenpos[1] + ix * y_step + border)
- xpos = int(screenpos[0] + ix * x_step + border)
- MoveWindow(wnds[iy], xpos, ypos, pos[2], pos[3], 1)
- BringWindowToTop(wnds[iy])
-
-
-def tile(*figure_numbers, **kwds):
- """Tile figure windows.
-
- Parameters
- ----------
- figure_numbers : list of integers or string
- specifying which figures to tile (default 'all').
- kwds : dict with key pairs
- specifying how many pairs of figures that are tiled at a time
-
- Description
- -----------
- TILE places all open figure windows around on the screen with no
- overlap. TILE(FIGS) can be used to specify which figures that
- should be tiled. Figures are not sorted when specified.
-
- Example:
- --------
- >>> import pylab as p
- >>> import wafo.fig as fig
- >>> for ix in range(7):
- ... f = p.figure(ix)
- >>> fig.tile() # tile all open figures
- >>> fig.tile( range(1,4), 5, 7) # tile figure 1,2,3,5 and 7
- >>> fig.tile(range(1,11), pairs=2) # tile figure 1 to 10 two at a time
- >>> fig.tile(range(1,11), pairs=3) # tile figure 1 to 10 three at a time
- >>> fig.close()
-
- See also
- --------
- fig.cycle, fig.keep, fig.maximize, fig.minimize
- fig.restore, fig.pile, fig.stack
-
- """
- wnds = find_figure_handles(*figure_numbers)
-
- nfigs = len(wnds)
- # Number of windows.
-
- if nfigs > 0:
- nfigspertile = kwds.get('pairs', nfigs)
-
- ceil = numpy.ceil
- sqrt = numpy.sqrt
- maximum = numpy.maximum
-
- nlayers = int(ceil(nfigs / nfigspertile))
-
- # Number of figures horisontally.
- nh = maximum(int(ceil(sqrt(nfigspertile))), 2)
- # Number of figures vertically.
- nv = maximum(int(ceil(nfigspertile / nh)), 2)
-
- screenpos = get_screen_position_and_size(wnds)
- screen_width, screen_heigth = screenpos[2:4]
-
- hspc = 10 # Horisontal space.
- topspc = 20 # Space above top figure.
- medspc = 10 # Space between figures.
- botspc = 20 # Space below bottom figure.
-
- figwid = (screen_width - (nh + 1) * hspc) / nh
- fighgt = (screen_heigth - (topspc + botspc) - (nv - 1) * medspc) / nv
-
- figwid = int(numpy.round(figwid))
- fighgt = int(numpy.round(fighgt))
-
- idx = 0
- for unused_ix in range(nlayers):
- for row in range(nv):
- figtop = int(screenpos[1] + topspc + row * (fighgt + medspc))
- for col in range(nh):
- if (row) * nh + col < nfigspertile:
- if idx < nfigs:
- figlft = int(
- screenpos[0] + (col + 1) * hspc + col * figwid)
- fighnd = wnds[idx]
- MoveWindow(fighnd, figlft, figtop, figwid, fighgt,
- 1)
- # Set position.
- BringWindowToTop(fighnd)
- idx += 1
-
-
-class _CycleGenerator(object):
-
- """Cycle through figure windows.
-
- Parameters
- ----------
- figure_numbers : list of integers or string
- specifying which figures to cycle through (default 'all').
- kwds : dict with the following keys
- pairs : number of figures to cycle in pairs (default 1)
- maximize: If True maximize figure when viewing (default False)
- interval : pause interval in seconds
-
- Description
- -----------
- CYCLE brings up all open figure in ascending order and pauses after
- each figure. Press escape to quit cycling, backspace to display previous
- figure(s) and press any other key to display next figure(s)
- When done, the figures are sorted in ascending order.
-
- CYCLE(maximize=True) does the same thing, except figures are maximized.
- CYCLE(pairs=2) cycle through all figures in pairs of 2.
-
- Examples:
- >>> import pylab as p
- >>> import wafo.fig as fig
- >>> for ix in range(4):
- ... f = p.figure(ix)
-
- fig.cycle(range(3), interval=1) # Cycle trough figure 0 to 2
-
- # Cycle trough figure 0 to 2 with figures maximized
- fig.cycle(range(3), maximize=True, interval=1)
- fig.cycle(interval=1) # Cycle through all figures one at a time
- fig.tile(pairs=2, interval=1)
- fig.cycle(pairs=2, interval=2) # Cycle through all figures two at a time
-
- fig.cycle(pairs=2) # Manually cycle through all figures two at a time
- >>> fig.close()
-
- See also
- --------
- fig.keep, fig.maximize, fig.restore, fig.pile,
- fig.stack, fig.tile
-
- """
- escape_key = chr(27)
- backspace_key = chr(8)
-
- def __init__(self, **kwds):
- self.dialog = None
- maximize = kwds.get('maximize', False)
- pairs = kwds.get('pairs', 1)
- self.interval = kwds.get('interval', 'user_defined')
- self.nfigspercycle = 1
- if maximize:
- self.command = win32con.SW_SHOWMAXIMIZED
- else:
- self.command = win32con.SW_SHOWNORMAL
- if pairs is not None:
- self.nfigspercycle = pairs
-
- def _set_options(self, kwds):
- self.__init__(**kwds)
-
- def _iterate(self, handles):
- i = 0
- numfigs = len(handles)
- self.dialog = _get_cycle_dialog(parent=None, interval=self.interval)
- while 0 <= i and i < numfigs:
- iu = min(i + self.nfigspercycle, numfigs)
- yield handles[i:iu]
- i = self.next_index(i)
- self.dialog.Destroy()
- raise StopIteration
-
- def next_index(self, i):
- result = self.dialog.ShowModal()
- if result == wx.ID_FORWARD:
- i += self.nfigspercycle
- elif result == wx.ID_BACKWARD:
- i -= self.nfigspercycle
- else:
- i = -1
- return i
-
- def __call__(self, *figure_numbers, **kwds):
- handles = find_figure_handles(*figure_numbers)
- numfigs = len(handles)
- if numfigs > 0:
- self._set_options(kwds)
- for handle in self._iterate(handles):
- _show_windows(handle, self.command, redraw_now=True)
-
- _show_windows(handles, win32con.SW_SHOWNORMAL)
-
-cycle = _CycleGenerator()
-
-
-if __name__ == '__main__':
- from wafo.testing import test_docstrings
- import matplotlib
- matplotlib.interactive(True)
- test_docstrings(__file__)
diff --git a/wafo/gaussian.py b/wafo/gaussian.py
deleted file mode 100644
index d2cc947..0000000
--- a/wafo/gaussian.py
+++ /dev/null
@@ -1,1040 +0,0 @@
-from __future__ import absolute_import
-from numpy import (r_, minimum, maximum, atleast_1d, atleast_2d, mod, ones,
- floor, random, eye, nonzero, where, repeat, sqrt, exp, inf,
- diag, zeros, sin, arcsin, nan)
-from numpy import triu
-from scipy.special import ndtr as cdfnorm, ndtri as invnorm
-from scipy.special import erfc
-import warnings
-import numpy as np
-
-try:
- from wafo import mvn # @UnresolvedImport
-except ImportError:
- warnings.warn('mvn not found. Check its compilation.')
- mvn = None
-try:
- from wafo import mvnprdmod # @UnresolvedImport
-except ImportError:
- warnings.warn('mvnprdmod not found. Check its compilation.')
- mvnprdmod = None
-try:
- from wafo import rindmod # @UnresolvedImport
-except ImportError:
- warnings.warn('rindmod not found. Check its compilation.')
- rindmod = None
-
-
-__all__ = ['Rind', 'rindmod', 'mvnprdmod', 'mvn', 'cdflomax', 'prbnormtndpc',
- 'prbnormndpc', 'prbnormnd', 'cdfnorm2d', 'prbnorm2d', 'cdfnorm',
- 'invnorm']
-
-
-class Rind(object):
-
- """
- RIND Computes multivariate normal expectations
-
- Parameters
- ----------
- S : array-like, shape Ntdc x Ntdc
- Covariance matrix of X=[Xt,Xd,Xc] (Ntdc = Nt+Nd+Nc)
- m : array-like, size Ntdc
- expectation of X=[Xt,Xd,Xc]
- Blo, Bup : array-like, shape Mb x Nb
- Lower and upper barriers used to compute the integration limits,
- Hlo and Hup, respectively.
- indI : array-like, length Ni
- vector of indices to the different barriers in the indicator function.
- (NB! restriction indI(1)=-1, indI(NI)=Nt+Nd, Ni = Nb+1)
- (default indI = 0:Nt+Nd)
- xc : values to condition on (default xc = zeros(0,1)), size Nc x Nx
- Nt : size of Xt (default Nt = Ntdc - Nc)
-
- Returns
- -------
- val: ndarray, size Nx
- expectation/density as explained below
- err, terr : ndarray, size Nx
- estimated sampling error and estimated truncation error, respectively.
- (err is with 99 confidence level)
-
- Notes
- -----
- RIND computes multivariate normal expectations, i.e.,
- E[Jacobian*Indicator|Condition ]*f_{Xc}(xc(:,ix))
- where
- "Indicator" = I{ Hlo(i) < X(i) < Hup(i), i = 1:N_t+N_d }
- "Jacobian" = J(X(Nt+1),...,X(Nt+Nd+Nc)), special case is
- "Jacobian" = |X(Nt+1)*...*X(Nt+Nd)|=|Xd(1)*Xd(2)..Xd(Nd)|
- "condition" = Xc=xc(:,ix), ix=1,...,Nx.
- X = [Xt, Xd, Xc], a stochastic vector of Multivariate Gaussian
- variables where Xt,Xd and Xc have the length Nt,Nd and Nc, respectively
- (Recommended limitations Nx,Nt<=100, Nd<=6 and Nc<=10)
-
- Multivariate probability is computed if Nd = 0.
-
- If Mb>> import wafo.gaussian as wg
- >>> import numpy as np
- >>> n = 5
- >>> Blo =-np.inf; Bup=-1.2; indI=np.array([-1, n-1], dtype=int) # Barriers
- >>> m = np.zeros(n); rho = 0.3;
- >>> Sc =(np.ones((n,n))-np.eye(n))*rho+np.eye(n)
- >>> rind = wg.Rind()
- >>> E0, err0, terr0 = rind(Sc,m,Blo,Bup,indI) # exact prob. 0.001946
-
- >>> A = np.repeat(Blo,n); B = np.repeat(Bup,n) # Integration limits
- >>> E1 = rind(np.triu(Sc),m,A,B) #same as E0
-
- Compute expectation E( abs(X1*X2*...*X5) )
- >>> xc = np.zeros((0,1))
- >>> infinity = 37
- >>> dev = np.sqrt(np.diag(Sc)) # std
- >>> ind = np.nonzero(indI[1:])[0]
- >>> Bup, Blo = np.atleast_2d(Bup,Blo)
- >>> Bup[0,ind] = np.minimum(Bup[0,ind] , infinity*dev[indI[ind+1]])
- >>> Blo[0,ind] = np.maximum(Blo[0,ind] ,-infinity*dev[indI[ind+1]])
- >>> val, err, terr = rind(Sc,m,Blo,Bup,indI, xc, nt=0)
- >>> np.allclose(val, 0.05494076, rtol=4e-2)
- True
- >>> err[0] < 5e-3, terr[0] < 1e-7
- (True, True)
-
- Compute expectation E( X1^{+}*X2^{+} ) with random
- correlation coefficient,Cov(X1,X2) = rho2.
- >>> m2 = [0, 0]; rho2 = np.random.rand(1)
- >>> Sc2 = [[1, rho2], [rho2 ,1]]
- >>> Blo2 = 0; Bup2 = np.inf; indI2 = [-1, 1]
- >>> rind2 = wg.Rind(method=1)
- >>> def g2(x):
- ... return (x*(np.pi/2+np.arcsin(x))+np.sqrt(1-x**2))/(2*np.pi)
- >>> E2 = g2(rho2) # exact value
- >>> E3 = rind(Sc2,m2,Blo2,Bup2,indI2,nt=0)
- >>> E4 = rind2(Sc2,m2,Blo2,Bup2,indI2,nt=0)
- >>> E5 = rind2(Sc2,m2,Blo2,Bup2,indI2,nt=0,abseps=1e-4)
-
- See also
- --------
- prbnormnd, prbnormndpc
-
- References
- ----------
- Podgorski et al. (2000)
- "Exact distributions for apparent waves in irregular seas"
- Ocean Engineering, Vol 27, no 1, pp979-1016.
-
- P. A. Brodtkorb (2004),
- Numerical evaluation of multinormal expectations
- In Lund university report series
- and in the Dr.Ing thesis:
- The probability of Occurrence of dangerous Wave Situations at Sea.
- Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU,
- Trondheim, Norway.
-
- Per A. Brodtkorb (2006)
- "Evaluating Nearly Singular Multinormal Expectations with Application to
- Wave Distributions",
- Methodology And Computing In Applied Probability, Volume 8, Number 1,
- pp. 65-91(27)
- """
-
- def __init__(self, **kwds):
- """
- Parameters
- ----------
- method : integer, optional
- defining the integration method
- 0 Integrate by Gauss-Legendre quadrature (Podgorski et al. 1999)
- 1 Integrate by SADAPT for Ndim<9 and by KRBVRC otherwise
- 2 Integrate by SADAPT for Ndim<20 and by KRBVRC otherwise
- 3 Integrate by KRBVRC by Genz (1993) (Fast Ndim<101) (default)
- 4 Integrate by KROBOV by Genz (1992) (Fast Ndim<101)
- 5 Integrate by RCRUDE by Genz (1992) (Slow Ndim<1001)
- 6 Integrate by SOBNIED (Fast Ndim<1041)
- 7 Integrate by DKBVRC by Genz (2003) (Fast Ndim<1001)
-
- xcscale : real scalar, optional
- scales the conditinal probability density, i.e.,
- f_{Xc} = exp(-0.5*Xc*inv(Sxc)*Xc + XcScale) (default XcScale=0)
- abseps, releps : real scalars, optional
- absolute and relative error tolerance.
- (default abseps=0, releps=1e-3)
- coveps : real scalar, optional
- error tolerance in Cholesky factorization (default 1e-13)
- maxpts, minpts : scalar integers, optional
- maximum and minimum number of function values allowed. The
- parameter, maxpts, can be used to limit the time. A sensible
- strategy is to start with MAXPTS = 1000*N, and then increase MAXPTS
- if ERROR is too large.
- (Only for METHOD~=0) (default maxpts=40000, minpts=0)
- seed : scalar integer, optional
- seed to the random generator used in the integrations
- (Only for METHOD~=0)(default floor(rand*1e9))
- nit : scalar integer, optional
- maximum number of Xt variables to integrate. This parameter can be
- used to limit the time. If NIT is less than the rank of the
- covariance matrix, the returned result is a upper bound for the
- true value of the integral. (default 1000)
- xcutoff : real scalar, optional
- cut off value where the marginal normal distribution is truncated.
- (Depends on requested accuracy. A value between 4 and 5 is
- reasonable.)
- xsplit : real scalar
- parameter controlling performance of quadrature integration:
- if Hup>=xCutOff AND Hlo<-XSPLIT OR
- Hup>=XSPLIT AND Hlo<=-xCutOff then
- do a different integration to increase speed
- in rind2 and rindnit. This give slightly different results
- if XSPILT>=xCutOff => do the same integration always
- (Only for METHOD==0)(default XSPLIT = 1.5)
- quadno : scalar integer
- Quadrature formulae number used in integration of Xd variables.
- This number implicitly determines number of nodes
- used. (Only for METHOD==0)
- speed : scalar integer
- defines accuracy of calculations by choosing different parameters,
- possible values: 1,2...,9 (9 fastest, default []).
- If not speed is None the parameters, ABSEPS, RELEPS, COVEPS,
- XCUTOFF, MAXPTS and QUADNO will be set according to
- INITOPTIONS.
- nc1c2 : scalar integer, optional
- number of times to use the regression equation to restrict
- integration area. Nc1c2 = 1,2 is recommended. (default 2)
- (note: works only for method >0)
- """
- self.method = 3
- self.xcscale = 0
- self.abseps = 0
- self.releps = 1e-3,
- self.coveps = 1e-10
- self.maxpts = 40000
- self.minpts = 0
- self.seed = None
- self.nit = 1000,
- self.xcutoff = None
- self.xsplit = 1.5
- self.quadno = 2
- self.speed = None
- self.nc1c2 = 2
-
- self.__dict__.update(**kwds)
- self.initialize(self.speed)
- self.set_constants()
-
- def initialize(self, speed=None):
- """
- Initializes member variables according to speed.
-
- Parameter
- ---------
- speed : scalar integer
- defining accuracy of calculations.
- Valid numbers: 1,2,...,10
- (1=slowest and most accurate,10=fastest, but less accuracy)
-
-
- Member variables initialized according to speed:
- -----------------------------------------------
- speed : Integer defining accuracy of calculations.
- abseps : Absolute error tolerance.
- releps : Relative error tolerance.
- covep : Error tolerance in Cholesky factorization.
- xcutoff : Truncation limit of the normal CDF
- maxpts : Maximum number of function values allowed.
- quadno : Quadrature formulae used in integration of Xd(i)
- implicitly determining # nodes
- """
- if speed is None:
- return
- self.speed = min(max(speed, 1), 13)
-
- self.maxpts = 10000
- self.quadno = r_[1:4] + (10 - min(speed, 9)) + (speed == 1)
- if speed in (11, 12, 13):
- self.abseps = 1e-1
- elif speed == 10:
- self.abseps = 1e-2
- elif speed in (7, 8, 9):
- self.abseps = 1e-2
- elif speed in (4, 5, 6):
- self.maxpts = 20000
- self.abseps = 1e-3
- elif speed in (1, 2, 3):
- self.maxpts = 30000
- self.abseps = 1e-4
-
- if speed < 12:
- tmp = max(abs(11 - abs(speed)), 1)
- expon = mod(tmp + 1, 3) + 1
- self.coveps = self.abseps * ((1.0e-1) ** expon)
- elif speed < 13:
- self.coveps = 0.1
- else:
- self.coveps = 0.5
-
- self.releps = min(self.abseps, 1.0e-2)
-
- if self.method == 0:
- # This gives approximately the same accuracy as when using
- # RINDDND and RINDNIT
- # xCutOff= MIN(MAX(xCutOff+0.5d0,4.d0),5.d0)
- self.abseps = self.abseps * 1.0e-1
- trunc_error = 0.05 * max(0, self.abseps)
- self.xcutoff = max(min(abs(invnorm(trunc_error)), 7), 1.2)
- self.abseps = max(self.abseps - trunc_error, 0)
-
- def set_constants(self):
- if self.xcutoff is None:
- trunc_error = 0.1 * self.abseps
- self.nc1c2 = max(1, self.nc1c2)
- xcut = abs(invnorm(trunc_error / (self.nc1c2 * 2)))
- self.xcutoff = max(min(xcut, 8.5), 1.2)
- # self.abseps = max(self.abseps- truncError,0);
- # self.releps = max(self.releps- truncError,0);
-
- if self.method > 0:
- names = ['method', 'xcscale', 'abseps', 'releps', 'coveps',
- 'maxpts', 'minpts', 'nit', 'xcutoff', 'nc1c2', 'quadno',
- 'xsplit']
-
- constants = [getattr(self, name) for name in names]
- constants[0] = mod(constants[0], 10)
- rindmod.set_constants(*constants) # @UndefinedVariable
-
- def __call__(self, cov, m, ab, bb, indI=None, xc=None, nt=None, **kwds):
- if any(kwds):
- self.__dict__.update(**kwds)
- self.set_constants()
- if xc is None:
- xc = zeros((0, 1))
-
- big, b_lo, b_up, xc = atleast_2d(cov, ab, bb, xc)
- b_lo = b_lo.copy()
- b_up = b_up.copy()
-
- Ntdc = big.shape[0]
- Nc = xc.shape[0]
- if nt is None:
- nt = Ntdc - Nc
-
- unused_Mb, Nb = b_lo.shape
- Nd = Ntdc - nt - Nc
- Ntd = nt + Nd
-
- if indI is None:
- if Nb != Ntd:
- raise ValueError('Inconsistent size of b_lo and b_up')
- indI = r_[-1:Ntd]
-
- Ex, indI = atleast_1d(m, indI)
- if self.seed is None:
- seed = int(floor(random.rand(1) * 1e10)) # @UndefinedVariable
- else:
- seed = int(self.seed)
-
- # INFIN = INTEGER, array of integration limits flags: size 1 x Nb
- # if INFIN(I) < 0, Ith limits are (-infinity, infinity);
- # if INFIN(I) = 0, Ith limits are (-infinity, Hup(I)];
- # if INFIN(I) = 1, Ith limits are [Hlo(I), infinity);
- # if INFIN(I) = 2, Ith limits are [Hlo(I), Hup(I)].
- infinity = 37
- dev = sqrt(diag(big)) # std
- ind = nonzero(indI[1:] > -1)[0]
- infin = repeat(2, len(indI) - 1)
- infin[ind] = (2 - (b_up[0, ind] > infinity * dev[indI[ind + 1]]) -
- 2 * (b_lo[0, ind] < -infinity * dev[indI[ind + 1]]))
-
- b_up[0, ind] = minimum(b_up[0, ind], infinity * dev[indI[ind + 1]])
- b_lo[0, ind] = maximum(b_lo[0, ind], -infinity * dev[indI[ind + 1]])
- ind2 = indI + 1
-
- return rindmod.rind(big, Ex, xc, nt, ind2, b_lo, b_up, infin, seed)
-
-
-def test_rind():
- """ Small test function
- """
- n = 5
- Blo = -inf
- Bup = -1.2
- indI = [-1, n - 1] # Barriers
- # A = np.repeat(Blo, n)
- # B = np.repeat(Bup, n) # Integration limits
- m = zeros(n)
- rho = 0.3
- Sc = (ones((n, n)) - eye(n)) * rho + eye(n)
- rind = Rind()
- E0 = rind(Sc, m, Blo, Bup, indI) # exact prob. 0.001946 A)
- print(E0)
-
- A = repeat(Blo, n)
- B = repeat(Bup, n) # Integration limits
- _E1 = rind(triu(Sc), m, A, B) # same as E0
-
- xc = zeros((0, 1))
- infinity = 37
- dev = sqrt(diag(Sc)) # std
- ind = nonzero(indI[1:])[0]
- Bup, Blo = atleast_2d(Bup, Blo)
- Bup[0, ind] = minimum(Bup[0, ind], infinity * dev[indI[ind + 1]])
- Blo[0, ind] = maximum(Blo[0, ind], -infinity * dev[indI[ind + 1]])
- _E3 = rind(Sc, m, Blo, Bup, indI, xc, nt=1)
-
-
-def cdflomax(x, alpha, m0):
- """
- Return CDF for local maxima for a zero-mean Gaussian process
-
- Parameters
- ----------
- x : array-like
- evaluation points
- alpha, m0 : real scalars
- irregularity factor and zero-order spectral moment (variance of the
- process), respectively.
-
- Returns
- -------
- prb : ndarray
- distribution function evaluated at x
-
- Notes
- -----
- The cdf is calculated from an explicit expression involving the
- standard-normal cdf. This relation is sometimes written as a convolution
-
- M = sqrt(m0)*( sqrt(1-a^2)*Z + a*R )
-
- where M denotes local maximum, Z is a standard normal r.v.,
- R is a standard Rayleigh r.v., and "=" means equality in distribution.
-
- Note that all local maxima of the process are considered, not
- only crests of waves.
-
- Example
- -------
- >>> import pylab
- >>> import wafo.gaussian as wg
- >>> import wafo.spectrum.models as wsm
- >>> import wafo.objects as wo
- >>> import wafo.stats as ws
- >>> S = wsm.Jonswap(Hm0=10).tospecdata();
- >>> xs = S.sim(10000)
- >>> ts = wo.mat2timeseries(xs)
- >>> tp = ts.turning_points()
- >>> mM = tp.cycle_pairs()
- >>> m0 = S.moment(1)[0]
- >>> alpha = S.characteristic('alpha')[0]
- >>> x = np.linspace(-10,10,200);
- >>> mcdf = ws.edf(mM.data)
-
- h = mcdf.plot(), pylab.plot(x,wg.cdflomax(x,alpha,m0))
-
- See also
- --------
- spec2mom, spec2bw
- """
- c1 = 1.0 / (sqrt(1 - alpha ** 2)) * x / sqrt(m0)
- c2 = alpha * c1
- return cdfnorm(c1) - alpha * exp(-x ** 2 / 2 / m0) * cdfnorm(c2)
-
-
-def prbnormtndpc(rho, a, b, d=None, df=0, abseps=1e-4, ierc=0, hnc=0.24):
- """
- Return Multivariate normal or T probability with product correlation.
-
- Parameters
- ----------
- rho : array-like
- vector of coefficients defining the correlation coefficient by:
- correlation(I,J) = rho[i]*rho[j]) for J!=I
- where -1 < rho[i] < 1
- a,b : array-like
- vector of lower and upper integration limits, respectively.
- Note: any values greater the 37 in magnitude, are considered as
- infinite values.
- d : array-like
- vector of means (default zeros(size(rho)))
- df = Degrees of freedom, NDF<=0 gives normal probabilities (default)
- abseps = absolute error tolerance. (default 1e-4)
- ierc = 1 if strict error control based on fourth derivative
- 0 if error control based on halving the intervals (default)
- hnc = start interval width of simpson rule (default 0.24)
-
- Returns
- -------
- value = estimated value for the integral
- bound = bound on the error of the approximation
- inform = INTEGER, termination status parameter:
- 0, if normal completion with ERROR < EPS;
- 1, if N > 1000 or N < 1.
- 2, IF any abs(rho)>=1
- 4, if ANY(b(I)<=A(i))
- 5, if number of terms exceeds maximum number of evaluation points
- 6, if fault accurs in normal subroutines
- 7, if subintervals are too narrow or too many
- 8, if bounds exceeds abseps
-
- PRBNORMTNDPC calculates multivariate normal or student T probability
- with product correlation structure for rectangular regions.
- The accuracy is as best around single precision, i.e., about 1e-7.
-
- Example:
- --------
- >>> import wafo.gaussian as wg
- >>> rho2 = np.random.rand(2)
- >>> a2 = np.zeros(2)
- >>> b2 = np.repeat(np.inf,2)
- >>> [val2,err2, ift2] = wg.prbnormtndpc(rho2,a2,b2)
- >>> def g2(x):
- ... return 0.25+np.arcsin(x[0]*x[1])/(2*np.pi)
- >>> E2 = g2(rho2) # exact value
- >>> np.abs(E2-val2)>> rho3 = np.random.rand(3)
- >>> a3 = np.zeros(3)
- >>> b3 = np.repeat(inf,3)
- >>> [val3, err3, ift3] = wg.prbnormtndpc(rho3,a3,b3)
- >>> def g3(x):
- ... return 0.5-sum(np.sort(np.arccos([x[0]*x[1],
- ... x[0]*x[2],x[1]*x[2]])))/(4*np.pi)
- >>> E3 = g3(rho3) # Exact value
- >>> np.abs(E3-val3) < 5 * err2
- True
-
-
- See also
- --------
- prbnormndpc, prbnormnd, Rind
-
- Reference
- ---------
- Charles Dunnett (1989)
- "Multivariate normal probability integrals with product correlation
- structure", Applied statistics, Vol 38,No3, (Algorithm AS 251)
- """
-
- if d is None:
- d = zeros(len(rho))
- # Make sure integration limits are finite
- aa = np.clip(a - d, -100, 100)
- bb = np.clip(b - d, -100, 100)
-
- return mvnprdmod.prbnormtndpc(rho, aa, bb, df, abseps, ierc, hnc)
-
-
-def prbnormndpc(rho, a, b, abserr=1e-4, relerr=1e-4, usesimpson=True,
- usebreakpoints=False):
- """
- Return Multivariate Normal probabilities with product correlation
-
- Parameters
- ----------
- rho = vector defining the correlation structure, i.e.,
- corr(Xi,Xj) = rho(i)*rho(j) for i~=j
- = 1 for i==j
- -1 <= rho <= 1
- a,b = lower and upper integration limits respectively.
- tol = requested absolute tolerance
-
- Returns
- -------
- value = value of integral
- error = estimated absolute error
-
- PRBNORMNDPC calculates multivariate normal probability
- with product correlation structure for rectangular regions.
- The accuracy is up to almost double precision, i.e., about 1e-14.
-
- Example:
- -------
- >>> import wafo.gaussian as wg
- >>> rho2 = np.random.rand(2)
- >>> a2 = np.zeros(2)
- >>> b2 = np.repeat(np.inf,2)
- >>> [val2,err2, ift2] = wg.prbnormndpc(rho2,a2,b2)
- >>> g2 = lambda x : 0.25+np.arcsin(x[0]*x[1])/(2*np.pi)
- >>> E2 = g2(rho2) #% exact value
- >>> np.abs(E2-val2)>> rho3 = np.random.rand(3)
- >>> a3 = np.zeros(3)
- >>> b3 = np.repeat(inf,3)
- >>> [val3,err3, ift3] = wg.prbnormndpc(rho3,a3,b3)
- >>> def g3(x):
- ... return 0.5-sum(np.sort(np.arccos([x[0]*x[1],
- ... x[0]*x[2],x[1]*x[2]])))/(4*np.pi)
- >>> E3 = g3(rho3) # Exact value
- >>> np.abs(E3-val3) 0:
- warnings.warn('Abnormal termination ier = %d\n\n%s' %
- (ier, _ERRORMESSAGE[ier]))
- return val, err, ier
-
-_ERRORMESSAGE = {}
-_ERRORMESSAGE[0] = ''
-_ERRORMESSAGE[1] = """
- Maximum number of subdivisions allowed has been achieved. one can allow
- more subdivisions by increasing the value of limit (and taking the
- according dimension adjustments into account). however, if this yields
- no improvement it is advised to analyze the integrand in order to
- determine the integration difficulties. if the position of a local
- difficulty can be determined (i.e. singularity discontinuity within
- the interval), it should be supplied to the routine as an element of
- the vector points. If necessary an appropriate special-purpose
- integrator must be used, which is designed for handling the type of
- difficulty involved.
- """
-_ERRORMESSAGE[2] = """
- the occurrence of roundoff error is detected, which prevents the requested
- tolerance from being achieved. The error may be under-estimated."""
-
-_ERRORMESSAGE[3] = """
- Extremely bad integrand behaviour occurs at some points of the integration
- interval."""
-_ERRORMESSAGE[4] = """
- The algorithm does not converge. Roundoff error is detected in the
- extrapolation table. It is presumed that the requested tolerance cannot be
- achieved, and that the returned result is the best which can be obtained.
- """
-_ERRORMESSAGE[5] = """
- The integral is probably divergent, or slowly convergent.
- It must be noted that divergence can occur with any other value of ier>0.
- """
-_ERRORMESSAGE[6] = """the input is invalid because:
- 1) npts2 < 2
- 2) break points are specified outside the integration range
- 3) (epsabs<=0 and epsrel EPS and MAXPTS
- function vaules used; increase MAXPTS to
- decrease ERROR;
- if INFORM = 2, N > NMAX or N < 1. where NMAX depends on the
- integration method
- Example
- -------
- Compute the probability that X1<0,X2<0,X3<0,X4<0,X5<0,
- Xi are zero-mean Gaussian variables with variances one
- and correlations Cov(X(i),X(j))=0.3:
- indI=[0 5], and barriers B_lo=[-inf 0], B_lo=[0 inf]
- gives H_lo = [-inf -inf -inf -inf -inf] H_lo = [0 0 0 0 0]
-
- >>> Et = 0.001946 # # exact prob.
- >>> n = 5; nt = n
- >>> Blo =-np.inf; Bup=0; indI=[-1, n-1] # Barriers
- >>> m = 1.2*np.ones(n); rho = 0.3;
- >>> Sc =(np.ones((n,n))-np.eye(n))*rho+np.eye(n)
- >>> rind = Rind()
- >>> E0, err0, terr0 = rind(Sc,m,Blo,Bup,indI, nt=nt)
-
- >>> A = np.repeat(Blo,n)
- >>> B = np.repeat(Bup,n)-m
- >>> val, err, inform = prbnormnd(Sc,A,B)
- >>> np.allclose([val, err, inform],
- ... [0.0019456719705212067, 1.0059406844578488e-05, 0])
- True
-
- >>> np.allclose(np.abs(val-Et) < err0+terr0, True)
- True
- >>> 'val = %2.5f' % val
- 'val = 0.00195'
-
- See also
- --------
- prbnormndpc, Rind
- """
-
- m, n = correl.shape
- Na = len(a)
- Nb = len(b)
- if (m != n or m != Na or m != Nb):
- raise ValueError('Size of input is inconsistent!')
-
- if maxpts is None:
- maxpts = 1000 * n
-
- maxpts = max(round(maxpts), 10 * n)
-
-# % array of correlation coefficients; the correlation
-# % coefficient in row I column J of the correlation matrix
-# % should be stored in CORREL( J + ((I-2)*(I-1))/2 ), for J < I.
-# % The correlation matrix must be positive semidefinite.
-
- D = np.diag(correl)
- if (any(D != 1)):
- raise ValueError('This is not a correlation matrix')
-
- # Make sure integration limits are finite
- A = np.clip(a, -100, 100)
- B = np.clip(b, -100, 100)
- ix = np.where(np.triu(np.ones((m, m)), 1) != 0)
- L = correl[ix].ravel() # % return only off diagonal elements
-
- infinity = 37
- infin = np.repeat(2, n) - (B > infinity) - 2 * (A < -infinity)
-
- err, val, inform = mvn.mvndst(A, B, infin, L, maxpts, abseps, releps)
-
- return val, err, inform
-
- # CALL the mexroutine
-# t0 = clock;
-# if ((method==0) && (n<=100)),
-# %NMAX = 100
-# [value, err,inform] = mexmvnprb(L,A,B,abseps,releps,maxpts);
-# elseif ( (method<0) || ((method<=0) && (n>100)) ),
-# % NMAX = 500
-# [value, err,inform] = mexmvnprb2(L,A,B,abseps,releps,maxpts);
-# else
-# [value, err,inform] = mexGenzMvnPrb(L,A,B,abseps,releps,maxpts,method);
-# end
-# exTime = etime(clock,t0);
-# '
-
-# gauss legendre points and weights, n = 6
-_W6 = [0.1713244923791705e+00, 0.3607615730481384e+00, 0.4679139345726904e+00]
-_X6 = [-0.9324695142031522e+00, -
- 0.6612093864662647e+00, -0.2386191860831970e+00]
-# gauss legendre points and weights, n = 12
-_W12 = [0.4717533638651177e-01, 0.1069393259953183e+00, 0.1600783285433464e+00,
- 0.2031674267230659e+00, 0.2334925365383547e+00, 0.2491470458134029e+00]
-_X12 = [-0.9815606342467191e+00, -0.9041172563704750e+00,
- -0.7699026741943050e+00,
- - 0.5873179542866171e+00, -0.3678314989981802e+00,
- -0.1252334085114692e+00]
-# gauss legendre points and weights, n = 20
-_W20 = [0.1761400713915212e-01, 0.4060142980038694e-01,
- 0.6267204833410906e-01, 0.8327674157670475e-01,
- 0.1019301198172404e+00, 0.1181945319615184e+00,
- 0.1316886384491766e+00, 0.1420961093183821e+00,
- 0.1491729864726037e+00, 0.1527533871307259e+00]
-_X20 = [-0.9931285991850949e+00, -0.9639719272779138e+00,
- - 0.9122344282513259e+00, -0.8391169718222188e+00,
- - 0.7463319064601508e+00, -0.6360536807265150e+00,
- - 0.5108670019508271e+00, -0.3737060887154196e+00,
- - 0.2277858511416451e+00, -0.7652652113349733e-01]
-
-
-def cdfnorm2d(b1, b2, r):
- """
- Returnc Bivariate Normal cumulative distribution function
-
- Parameters
- ----------
-
- b1, b2 : array-like
- upper integration limits
- r : real scalar
- correlation coefficient (-1 <= r <= 1).
-
- Returns
- -------
- bvn : ndarray
- distribution function evaluated at b1, b2.
-
- Notes
- -----
- CDFNORM2D computes bivariate normal probabilities, i.e., the probability
- Prob(X1 <= B1 and X2 <= B2) with an absolute error less than 1e-15.
-
- This function is based on the method described by Drezner, z and
- G.O. Wesolowsky, (1989), with major modifications for double precision,
- and for |r| close to 1.
-
- Example
- -------
- >>> import wafo.gaussian as wg
- >>> x = np.linspace(-5,5,20)
- >>> [B1,B2] = np.meshgrid(x, x)
- >>> r = 0.3;
- >>> F = wg.cdfnorm2d(B1,B2,r)
-
- surf(x,x,F)
-
- See also
- --------
- cdfnorm
-
- Reference
- ---------
- Drezner, z and g.o. Wesolowsky, (1989),
- "On the computation of the bivariate normal integral",
- Journal of statist. comput. simul. 35, pp. 101-107,
- """
- # Translated into Python
- # Per A. Brodtkorb
- #
- # Original code
- # by alan genz
- # department of mathematics
- # washington state university
- # pullman, wa 99164-3113
- # email : alangenz@wsu.edu
-
- b1, b2, r = np.broadcast_arrays(b1, b2, r)
- cshape = b1.shape
- h, k, r = -b1.ravel(), -b2.ravel(), r.ravel()
-
- bvn = where(abs(r) > 1, nan, 0.0)
-
- two = 2.e0
- twopi = 6.283185307179586e0
-
- hk = h * k
-
- k0, = nonzero(abs(r) < 0.925e0)
- if len(k0) > 0:
- hs = (h[k0] ** 2 + k[k0] ** 2) / two
- asr = arcsin(r[k0])
- k1, = nonzero(r[k0] >= 0.75)
- if len(k1) > 0:
- k01 = k0[k1]
- for i in range(10):
- for sign in - 1, 1:
- sn = sin(asr[k1] * (sign * _X20[i] + 1) / 2)
- bvn[k01] = bvn[k01] + _W20[i] * \
- exp((sn * hk[k01] - hs[k1]) / (1 - sn * sn))
-
- k1, = nonzero((0.3 <= r[k0]) & (r[k0] < 0.75))
- if len(k1) > 0:
- k01 = k0[k1]
- for i in range(6):
- for sign in - 1, 1:
- sn = sin(asr[k1] * (sign * _X12[i] + 1) / 2)
- bvn[k01] = bvn[k01] + _W12[i] * \
- exp((sn * hk[k01] - hs[k1]) / (1 - sn * sn))
-
- k1, = nonzero(r[k0] < 0.3)
- if len(k1) > 0:
- k01 = k0[k1]
- for i in range(3):
- for sign in - 1, 1:
- sn = sin(asr[k1] * (sign * _X6[i] + 1) / 2)
- bvn[k01] = bvn[k01] + _W6[i] * \
- exp((sn * hk[k01] - hs[k1]) / (1 - sn * sn))
-
- bvn[k0] *= asr / (two * twopi)
- bvn[k0] += fi(-h[k0]) * fi(-k[k0])
-
- k1, = nonzero((0.925 <= abs(r)) & (abs(r) <= 1))
- if len(k1) > 0:
- k2, = nonzero(r[k1] < 0)
- if len(k2) > 0:
- k12 = k1[k2]
- k[k12] = -k[k12]
- hk[k12] = -hk[k12]
-
- k3, = nonzero(abs(r[k1]) < 1)
- if len(k3) > 0:
- k13 = k1[k3]
- a2 = (1 - r[k13]) * (1 + r[k13])
- a = sqrt(a2)
- b = abs(h[k13] - k[k13])
- bs = b * b
- c = (4.e0 - hk[k13]) / 8.e0
- d = (12.e0 - hk[k13]) / 16.e0
- asr = -(bs / a2 + hk[k13]) / 2.e0
- k4, = nonzero(asr > -100.e0)
- if len(k4) > 0:
- bvn[k13[k4]] = (a[k4] * exp(asr[k4]) *
- (1 - c[k4] * (bs[k4] - a2[k4]) *
- (1 - d[k4] * bs[k4] / 5) / 3 +
- c[k4] * d[k4] * a2[k4] ** 2 / 5))
-
- k5, = nonzero(hk[k13] < 100.e0)
- if len(k5) > 0:
- # b = sqrt(bs);
- k135 = k13[k5]
- bvn[k135] = bvn[k135] - exp(-hk[k135] / 2) * sqrt(twopi) * \
- fi(-b[k5] / a[k5]) * b[k5] * \
- (1 - c[k5] * bs[k5] * (1 - d[k5] * bs[k5] / 5) / 3)
-
- a /= two
- for i in range(10):
- for sign in - 1, 1:
- xs = (a * (sign * _X20[i] + 1)) ** 2
- rs = sqrt(1 - xs)
- asr = -(bs / xs + hk[k13]) / 2
- k6, = nonzero(asr > -100.e0)
- if len(k6) > 0:
- k136 = k13[k6]
- bvn[k136] += (a[k6] * _W20[i] * exp(asr[k6]) *
- (exp(-hk[k136] * (1 - rs[k6]) /
- (2 * (1 + rs[k6]))) / rs[k6] -
- (1 + c[k6] * xs[k6] *
- (1 + d[k6] * xs[k6]))))
-
- bvn[k3] = -bvn[k3] / twopi
-
- k7, = nonzero(r[k1] > 0)
- if len(k7):
- k17 = k1[k7]
- bvn[k17] += fi(-np.maximum(h[k17], k[k17]))
-
- k8, = nonzero(r[k1] < 0)
- if len(k8) > 0:
- k18 = k1[k8]
- bvn[k18] = -bvn[k18] + np.maximum(0, fi(-h[k18]) - fi(-k[k18]))
-
- bvn.shape = cshape
- return bvn
-
-
-def fi(x):
- return 0.5 * (erfc((-x) / sqrt(2)))
-
-
-def prbnorm2d(a, b, r):
- """
- Returns Bivariate Normal probability
-
- Parameters
- ---------
- a, b : array-like, size==2
- vector with lower and upper integration limits, respectively.
- r : real scalar
- correlation coefficient
-
- Returns
- -------
- prb : real scalar
- computed probability Prob(A[0] <= X1 <= B[0] and A[1] <= X2 <= B[1])
- with an absolute error less than 1e-15.
-
- Example
- -------
- >>> import wafo.gaussian as wg
- >>> a = [-1, -2]
- >>> b = [1, 1]
- >>> r = 0.3
- >>> np.allclose(wg.prbnorm2d(a,b,r), 0.56659121350428077)
- True
-
- See also
- --------
- cdfnorm2d,
- cdfnorm,
- prbnormndpc
- """
- infinity = 37
- lower = np.asarray(a)
- upper = np.asarray(b)
- if np.all((lower <= -infinity) & (infinity <= upper)):
- return 1.0
- if (lower >= upper).any():
- return 0.0
- correl = r
- infin = np.repeat(2, 2) - (upper > infinity) - 2 * (lower < -infinity)
-
- if np.all(infin == 2):
- return (bvd(lower[0], lower[1], correl) -
- bvd(upper[0], lower[1], correl) -
- bvd(lower[0], upper[1], correl) +
- bvd(upper[0], upper[1], correl))
- elif (infin[0] == 2 and infin[1] == 1):
- return (bvd(lower[0], lower[1], correl) -
- bvd(upper[0], lower[1], correl))
- elif (infin[0] == 1 and infin[1] == 2):
- return (bvd(lower[0], lower[1], correl) -
- bvd(lower[0], upper[1], correl))
- elif (infin[0] == 2 and infin[1] == 0):
- return (bvd(-upper[0], -upper[1], correl) -
- bvd(-lower[0], -upper[1], correl))
- elif (infin[0] == 0 and infin[1] == 2):
- return (bvd(-upper[0], -upper[1], correl) -
- bvd(-upper[0], -lower[1], correl))
- elif (infin[0] == 1 and infin[1] == 0):
- return bvd(lower[0], -upper[1], -correl)
- elif (infin[0] == 0 and infin[1] == 1):
- return bvd(-upper[0], lower[1], -correl)
- elif (infin[0] == 1 and infin[1] == 1):
- return bvd(lower[0], lower[1], correl)
- elif (infin[0] == 0 and infin[1] == 0):
- return bvd(-upper[0], -upper[1], correl)
- return 1
-
-
-def bvd(lo, up, r):
- return cdfnorm2d(-lo, -up, r)
-
-
-if __name__ == '__main__':
- from wafo.testing import test_docstrings
- test_docstrings(__file__)
-# if __name__ == '__main__':
-# if False: #True: #
-# test_rind()
-# else:
-# import doctest
-# doctest.testmod()
diff --git a/wafo/integrate.py b/wafo/integrate.py
deleted file mode 100644
index 779fbcf..0000000
--- a/wafo/integrate.py
+++ /dev/null
@@ -1,1476 +0,0 @@
-from __future__ import absolute_import, division, print_function
-import warnings
-import numpy as np
-from numpy import pi, sqrt, ones, zeros
-from scipy import integrate as intg
-import scipy.special.orthogonal as ort
-from scipy import special as sp
-
-from scipy.integrate import simps, trapz
-from wafo.plotbackend import plotbackend as plt
-from wafo.demos import humps
-from numdifftools.extrapolation import dea3
-# from wafo.dctpack import dct
-from collections import defaultdict
-# from pychebfun import Chebfun
-
-_EPS = np.finfo(float).eps
-_NODES_AND_WEIGHTS = defaultdict(list)
-
-__all__ = ['dea3', 'clencurt', 'romberg',
- 'h_roots', 'j_roots', 'la_roots', 'p_roots', 'qrule',
- 'gaussq', 'richardson', 'quadgr', 'qdemo']
-
-
-def _assert(cond, msg):
- if not cond:
- raise ValueError(msg)
-
-
-def _assert_warn(cond, msg):
- if not cond:
- warnings.warn(msg)
-
-
-def clencurt(fun, a, b, n=5, trace=False):
- """
- Numerical evaluation of an integral, Clenshaw-Curtis method.
-
- Parameters
- ----------
- fun : callable
- a, b : array-like
- Lower and upper integration limit, respectively.
- n : integer
- defines number of evaluation points (default 5)
-
- Returns
- -------
- q_val = evaluated integral
- tol = Estimate of the approximation error
-
- Notes
- -----
- CLENCURT approximates the integral of f(x) from a to b
- using an 2*n+1 points Clenshaw-Curtis formula.
- The error estimate is usually a conservative estimate of the
- approximation error.
- The integral is exact for polynomials of degree 2*n or less.
-
- Example
- -------
- >>> import numpy as np
- >>> val, err = clencurt(np.exp, 0, 2)
- >>> np.allclose(val, np.expm1(2)), err[0] < 1e-10
- (True, True)
-
-
- See also
- --------
- simpson,
- gaussq
-
- References
- ----------
- [1] Goodwin, E.T. (1961),
- "Modern Computing Methods",
- 2nd edition, New yourk: Philosophical Library, pp. 78--79
-
- [2] Clenshaw, C.W. and Curtis, A.R. (1960),
- Numerische Matematik, Vol. 2, pp. 197--205
- """
- # make sure n_2 is even
- n_2 = 2 * int(n)
- a, b = np.atleast_1d(a, b)
- a_shape = a.shape
- a = a.ravel()
- b = b.ravel()
-
- a_size = np.prod(a_shape)
-
- s = np.c_[0:n_2 + 1:1]
- s_2 = np.c_[0:n_2 + 1:2]
- x = np.cos(np.pi * s / n_2) * (b - a) / 2. + (b + a) / 2
-
- if hasattr(fun, '__call__'):
- f = fun(x)
- else:
- x_0 = np.flipud(fun[:, 0])
- n_2 = len(x_0) - 1
- _assert(abs(x - x_0) <= 1e-8,
- 'Input vector x must equal cos(pi*s/n_2)*(b-a)/2+(b+a)/2')
-
- f = np.flipud(fun[:, 1::])
-
- if trace:
- plt.plot(x, f, '+')
-
- # using a Gauss-Lobatto variant, i.e., first and last
- # term f(a) and f(b) is multiplied with 0.5
- f[0, :] = f[0, :] / 2
- f[n_2, :] = f[n_2, :] / 2
-
- # x = cos(pi*0:n_2/n_2)
- # f = f(x)
- #
- # N+1
- # c(k) = (2/N) sum f''(n)*cos(pi*(2*k-2)*(n-1)/N), 1 <= k <= N/2+1.
- # n=1
- n = n_2 // 2
- fft = np.fft.fft
- tmp = np.real(fft(f[:n_2, :], axis=0))
- c = 2 / n_2 * (tmp[0:n + 1, :] + np.cos(np.pi * s_2) * f[n_2, :])
- c[0, :] = c[0, :] / 2
- c[n, :] = c[n, :] / 2
-
- c = c[0:n + 1, :] / ((s_2 - 1) * (s_2 + 1))
- q_val = (a - b) * np.sum(c, axis=0)
-
- abserr = (b - a) * np.abs(c[n, :])
-
- if a_size > 1:
- abserr = np.reshape(abserr, a_shape)
- q_val = np.reshape(q_val, a_shape)
- return q_val, abserr
-
-
-def romberg(fun, a, b, releps=1e-3, abseps=1e-3):
- """
- Numerical integration with the Romberg method
-
- Parameters
- ----------
- fun : callable
- function to integrate
- a, b : real scalars
- lower and upper integration limits, respectively.
- releps, abseps : scalar, optional
- requested relative and absolute error, respectively.
-
- Returns
- -------
- Q : scalar
- value of integral
- abserr : scalar
- estimated absolute error of integral
-
- ROMBERG approximates the integral of F(X) from A to B
- using Romberg's method of integration. The function F
- must return a vector of output values if a vector of input values is given.
-
-
- Example
- -------
- >>> import numpy as np
- >>> [q,err] = romberg(np.sqrt,0,10,0,1e-4)
- >>> np.allclose(q, 21.08185107)
- True
- >>> err[0] < 1e-4
- True
- """
- h = b - a
- h_min = 1.0e-9
- # Max size of extrapolation table
- table_limit = max(min(np.round(np.log2(h / h_min)), 30), 3)
-
- rom = zeros((2, table_limit))
-
- rom[0, 0] = h * (fun(a) + fun(b)) / 2
- ipower = 1
- f_p = ones(table_limit) * 4
-
- # q_val1 = 0
- q_val2 = 0.
- q_val4 = rom[0, 0]
- abserr = q_val4
- # epstab = zeros(1,decdigs+7)
- # newflg = 1
- # [res,abserr,epstab,newflg] = dea(newflg,q_val4,abserr,epstab)
- two = 1
- one = 0
- converged = False
- for i in range(1, table_limit):
- h *= 0.5
- u_n5 = np.sum(fun(a + np.arange(1, 2 * ipower, 2) * h)) * h
-
- # trapezoidal approximations
- # T2n = 0.5 * (Tn + Un) = 0.5*Tn + u_n5
- rom[two, 0] = 0.5 * rom[one, 0] + u_n5
-
- f_p[i] = 4 * f_p[i - 1]
- # Richardson extrapolation
- for k in range(i):
- rom[two, k + 1] = (rom[two, k] +
- (rom[two, k] - rom[one, k]) / (f_p[k] - 1))
-
- q_val1 = q_val2
- q_val2 = q_val4
- q_val4 = rom[two, i]
-
- if 2 <= i:
- res, abserr = dea3(q_val1, q_val2, q_val4)
- # q_val4 = res
- converged = abserr <= max(abseps, releps * abs(res))
- if converged:
- break
-
- # rom(1,1:i) = rom(2,1:i)
- two = one
- one = (one + 1) % 2
- ipower *= 2
- _assert(converged, "Integral did not converge to the required accuracy!")
- return res, abserr
-
-
-def _h_roots_newton(n, releps=3e-14, max_iter=10):
- # pim4=0.7511255444649425
- pim4 = np.pi ** (-1. / 4)
-
- # The roots are symmetric about the origin, so we have to
- # find only half of them.
- m = int(np.fix((n + 1) / 2))
-
- # Initial approximations to the roots go into z.
- anu = 2.0 * n + 1
- rhs = np.arange(3, 4 * m, 4) * np.pi / anu
- theta = _get_theta(rhs)
- z = sqrt(anu) * np.cos(theta)
-
- p = zeros((3, len(z)))
- k_0 = 0
- k_p1 = 1
-
- for _i in range(max_iter):
- # Newtons method carried out simultaneously on the roots.
- p[k_0, :] = 0
- p[k_p1, :] = pim4
-
- for j in range(1, n + 1):
- # Loop up the recurrence relation to get the Hermite
- # polynomials evaluated at z.
- k_m1 = k_0
- k_0 = k_p1
- k_p1 = np.mod(k_p1 + 1, 3)
-
- p[k_p1, :] = (z * sqrt(2 / j) * p[k_0, :] -
- sqrt((j - 1) / j) * p[k_m1, :])
-
- # p now contains the desired Hermite polynomials.
- # We next compute p_deriv, the derivatives,
- # by the relation (4.5.21) using p2, the polynomials
- # of one lower order.
-
- p_deriv = sqrt(2 * n) * p[k_0, :]
- d_z = p[k_p1, :] / p_deriv
-
- z = z - d_z # Newtons formula.
-
- converged = not np.any(abs(d_z) > releps)
- if converged:
- break
-
- _assert_warn(converged, 'Newton iteration did not converge!')
- weights = 2. / p_deriv ** 2
- return _expand_roots(z, weights, n, m)
-
-
-def h_roots(n, method='newton'):
- """
- Returns the roots (x) of the nth order Hermite polynomial,
- H_n(x), and weights (w) to use in Gaussian Quadrature over
- [-inf,inf] with weighting function exp(-x**2).
-
- Parameters
- ----------
- n : integer
- number of roots
- method : 'newton' or 'eigenvalue'
- uses Newton Raphson to find zeros of the Hermite polynomial (Fast)
- or eigenvalue of the jacobi matrix (Slow) to obtain the nodes and
- weights, respectively.
-
- Returns
- -------
- x : ndarray
- roots
- w : ndarray
- weights
-
- Example
- -------
- >>> import numpy as np
- >>> x, w = h_roots(10)
- >>> np.allclose(np.sum(x*w), -5.2516042729766621e-19)
- True
-
- See also
- --------
- qrule, gaussq
-
- References
- ----------
- [1] Golub, G. H. and Welsch, J. H. (1969)
- 'Calculation of Gaussian Quadrature Rules'
- Mathematics of Computation, vol 23,page 221-230,
-
- [2]. Stroud and Secrest (1966), 'gaussian quadrature formulas',
- prentice-hall, Englewood cliffs, n.j.
- """
-
- if not method.startswith('n'):
- return ort.h_roots(n)
- return _h_roots_newton(n)
-
-
-def _j_roots_newton(n, alpha, beta, releps=3e-14, max_iter=10):
- # Initial approximations to the roots go into z.
- alfbet = alpha + beta
-
- z = np.cos(np.pi * (np.arange(1, n + 1) - 0.25 + 0.5 * alpha) /
- (n + 0.5 * (alfbet + 1)))
-
- p = zeros((3, len(z)))
- k_0 = 0
- k_p1 = 1
- for _i in range(max_iter):
- # Newton's method carried out simultaneously on the roots.
- tmp = 2 + alfbet
- p[k_0, :] = 1
- p[k_p1, :] = (alpha - beta + tmp * z) / 2
-
- for j in range(2, n + 1):
- # Loop up the recurrence relation to get the Jacobi
- # polynomials evaluated at z.
- k_m1 = k_0
- k_0 = k_p1
- k_p1 = np.mod(k_p1 + 1, 3)
-
- a = 2. * j * (j + alfbet) * tmp
- tmp = tmp + 2
- c = 2 * (j - 1 + alpha) * (j - 1 + beta) * tmp
- b = (tmp - 1) * (alpha ** 2 - beta ** 2 + tmp * (tmp - 2) * z)
-
- p[k_p1, :] = (b * p[k_0, :] - c * p[k_m1, :]) / a
-
- # p now contains the desired Jacobi polynomials.
- # We next compute p_deriv, the derivatives with a standard
- # relation involving the polynomials of one lower order.
-
- p_deriv = ((n * (alpha - beta - tmp * z) * p[k_p1, :] +
- 2 * (n + alpha) * (n + beta) * p[k_0, :]) /
- (tmp * (1 - z ** 2)))
- d_z = p[k_p1, :] / p_deriv
- z = z - d_z # Newton's formula.
-
- converged = not any(abs(d_z) > releps * abs(z))
- if converged:
- break
-
- _assert_warn(converged, 'too many iterations in jrule')
-
- x = z # Store the root and the weight.
- f = (sp.gammaln(alpha + n) + sp.gammaln(beta + n) -
- sp.gammaln(n + 1) - sp.gammaln(alpha + beta + n + 1))
- weights = (np.exp(f) * tmp * 2 ** alfbet / (p_deriv * p[k_0, :]))
- return x, weights
-
-
-def j_roots(n, alpha, beta, method='newton'):
- """
- Returns the roots of the nth order Jacobi polynomial, P^(alpha,beta)_n(x)
- and weights (w) to use in Gaussian Quadrature over [-1,1] with weighting
- function (1-x)**alpha (1+x)**beta with alpha,beta > -1.
-
- Parameters
- ----------
- n : integer
- number of roots
- alpha,beta : scalars
- defining shape of Jacobi polynomial
- method : 'newton' or 'eigenvalue'
- uses Newton Raphson to find zeros of the Hermite polynomial (Fast)
- or eigenvalue of the jacobi matrix (Slow) to obtain the nodes and
- weights, respectively.
-
- Returns
- -------
- x : ndarray
- roots
- w : ndarray
- weights
-
-
- Example
- --------
- >>> [x,w]= j_roots(10,0,0)
- >>> sum(x*w)
- 2.7755575615628914e-16
-
- See also
- --------
- qrule, gaussq
-
-
- Reference
- ---------
- [1] Golub, G. H. and Welsch, J. H. (1969)
- 'Calculation of Gaussian Quadrature Rules'
- Mathematics of Computation, vol 23,page 221-230,
-
- [2]. Stroud and Secrest (1966), 'gaussian quadrature formulas',
- prentice-hall, Englewood cliffs, n.j.
- """
- _assert((-1 < alpha) & (-1 < beta),
- 'alpha and beta must be greater than -1')
- if not method.startswith('n'):
- return ort.j_roots(n, alpha, beta)
- return _j_roots_newton(n, alpha, beta)
-
-
-def la_roots(n, alpha=0, method='newton'):
- """
- Returns the roots (x) of the nth order generalized (associated) Laguerre
- polynomial, L^(alpha)_n(x), and weights (w) to use in Gaussian quadrature
- over [0,inf] with weighting function exp(-x) x**alpha with alpha > -1.
-
- Parameters
- ----------
- n : integer
- number of roots
- method : 'newton' or 'eigenvalue'
- uses Newton Raphson to find zeros of the Laguerre polynomial (Fast)
- or eigenvalue of the jacobi matrix (Slow) to obtain the nodes and
- weights, respectively.
-
- Returns
- -------
- x : ndarray
- roots
- w : ndarray
- weights
-
- Example
- -------
- >>> import numpy as np
- >>> [x,w] = h_roots(10)
- >>> np.allclose(np.sum(x*w) < 1e-16, True)
- True
-
- See also
- --------
- qrule, gaussq
-
- References
- ----------
- [1] Golub, G. H. and Welsch, J. H. (1969)
- 'Calculation of Gaussian Quadrature Rules'
- Mathematics of Computation, vol 23,page 221-230,
-
- [2]. Stroud and Secrest (1966), 'gaussian quadrature formulas',
- prentice-hall, Englewood cliffs, n.j.
- """
- _assert(-1 < alpha, 'alpha must be greater than -1')
-
- if not method.startswith('n'):
- return ort.la_roots(n, alpha)
- return _la_roots_newton(n, alpha)
-
-
-def _get_theta(rhs):
- r_3 = rhs ** (1. / 3)
- r_2 = r_3 ** 2
- c = [9.084064e-01, 5.214976e-02, 2.579930e-03, 3.986126e-03]
- theta = r_3 * (c[0] + r_2 * (c[1] + r_2 * (c[2] + r_2 * c[3])))
- return theta
-
-
-def _la_roots_newton(n, alpha, releps=3e-14, max_iter=10):
-
- # Initial approximations to the roots go into z.
- anu = 4.0 * n + 2.0 * alpha + 2.0
- rhs = np.arange(4 * n - 1, 2, -4) * np.pi / anu
- theta = _get_theta(rhs)
- z = anu * np.cos(theta) ** 2
-
- d_z = zeros(len(z))
- p = zeros((3, len(z)))
- p_previous = zeros((1, len(z)))
- p_deriv = zeros((1, len(z)))
- k_0 = 0
- k_p1 = 1
- k = slice(len(z))
- for _i in range(max_iter):
- # Newton's method carried out simultaneously on the roots.
- p[k_0, k] = 0.
- p[k_p1, k] = 1.
-
- for j in range(1, n + 1):
- # Loop up the recurrence relation to get the Laguerre
- # polynomials evaluated at z.
- km1 = k_0
- k_0 = k_p1
- k_p1 = np.mod(k_p1 + 1, 3)
-
- p[k_p1, k] = ((2 * j - 1 + alpha - z[k]) * p[k_0, k] -
- (j - 1 + alpha) * p[km1, k]) / j
- # end
- # p now contains the desired Laguerre polynomials.
- # We next compute p_deriv, the derivatives with a standard
- # relation involving the polynomials of one lower order.
-
- p_previous[k] = p[k_0, k]
- p_deriv[k] = (n * p[k_p1, k] - (n + alpha) * p_previous[k]) / z[k]
-
- d_z[k] = p[k_p1, k] / p_deriv[k]
- z[k] = z[k] - d_z[k] # Newton?s formula.
- # k = find((abs(d_z) > releps.*z))
-
- converged = not np.any(abs(d_z) > releps)
- if converged:
- break
-
- _assert_warn(converged, 'too many iterations!')
-
- nodes = z
- weights = -np.exp(sp.gammaln(alpha + n) -
- sp.gammaln(n)) / (p_deriv * n * p_previous)
- return nodes, weights
-
-
-def _p_roots_newton_start(n):
- m = int(np.fix((n + 1) / 2))
- t = (np.pi / (4 * n + 2)) * np.arange(3, 4 * m, 4)
- a = 1 - (1 - 1 / n) / (8 * n * n)
- x = a * np.cos(t)
- return m, x
-
-
-def _p_roots_newton(n):
- """
- Algorithm given by Davis and Rabinowitz in 'Methods
- of Numerical Integration', page 365, Academic Press, 1975.
- """
- m, x = _p_roots_newton_start(n)
-
- e_1 = n * (n + 1)
- for _j in range(2):
- p_km1 = 1
- p_k = x
- for k in range(2, n + 1):
- t_1 = x * p_k
- p_kp1 = t_1 - p_km1 - (t_1 - p_km1) / k + t_1
- p_km1 = p_k
- p_k = p_kp1
-
- den = 1. - x * x
- d_1 = n * (p_km1 - x * p_k)
- d_pn = d_1 / den
- d_2pn = (2. * x * d_pn - e_1 * p_k) / den
- d_3pn = (4. * x * d_2pn + (2 - e_1) * d_pn) / den
- d_4pn = (6. * x * d_3pn + (6 - e_1) * d_2pn) / den
- u = p_k / d_pn
- v = d_2pn / d_pn
- h = -u * (1 + (.5 * u) * (v + u * (v * v - u * d_3pn / (3 * d_pn))))
- p = p_k + h * (d_pn + (.5 * h) * (d_2pn + (h / 3) *
- (d_3pn + .25 * h * d_4pn)))
- d_p = d_pn + h * (d_2pn + (.5 * h) * (d_3pn + h * d_4pn / 3))
- h = h - p / d_p
- x = x + h
-
- nodes = -x - h
- f_x = d_1 - h * e_1 * (p_k + (h / 2) * (d_pn + (h / 3) *
- (d_2pn + (h / 4) *
- (d_3pn + (.2 * h) * d_4pn))))
- weights = 2 * (1 - nodes ** 2) / (f_x ** 2)
- return _expand_roots(nodes, weights, n, m)
-
-
-def _p_roots_newton1(n, releps=1e-15, max_iter=100):
- m, x = _p_roots_newton_start(n)
- # Compute the zeros of the N+1 Legendre Polynomial
- # using the recursion relation and the Newton-Raphson method
-
- # Legendre-Gauss Polynomials
- p = zeros((3, m))
-
- # Derivative of LGP
- p_deriv = zeros((m,))
- d_x = zeros((m,))
-
- # Compute the zeros of the N+1 Legendre Polynomial
- # using the recursion relation and the Newton-Raphson method
-
- # Iterate until new points are uniformly within epsilon of old
- # points
- k = slice(m)
- k_0 = 0
- k_p1 = 1
- for _ix in range(max_iter):
- p[k_0, k] = 1
- p[k_p1, k] = x[k]
-
- for j in range(2, n + 1):
- k_m1 = k_0
- k_0 = k_p1
- k_p1 = np.mod(k_0 + 1, 3)
- p[k_p1, k] = ((2 * j - 1) * x[k] * p[k_0, k] -
- (j - 1) * p[k_m1, k]) / j
-
- p_deriv[k] = n * (p[k_0, k] - x[k] * p[k_p1, k]) / (1 - x[k] ** 2)
-
- d_x[k] = p[k_p1, k] / p_deriv[k]
- x[k] = x[k] - d_x[k]
- k, = np.nonzero((abs(d_x) > releps * np.abs(x)))
- converged = len(k) == 0
- if converged:
- break
-
- _assert(converged, 'Too many iterations!')
-
- nodes = -x
- weights = 2. / ((1 - nodes ** 2) * (p_deriv ** 2))
- return _expand_roots(nodes, weights, n, m)
-
-
-def _expand_roots(x, w, n, m):
- if (m + m) > n:
- x[m - 1] = 0.0
- if not (m + m) == n:
- m = m - 1
- x = np.hstack((x, -x[m - 1::-1]))
- w = np.hstack((w, w[m - 1::-1]))
- return x, w
-
-
-def p_roots(n, method='newton', a=-1, b=1):
- """
- Returns the roots (x) of the nth order Legendre polynomial, P_n(x),
- and weights to use in Gaussian Quadrature over [-1,1] with weighting
- function 1.
-
- Parameters
- ----------
- n : integer
- number of roots
- method : 'newton' or 'eigenvalue'
- uses Newton Raphson to find zeros of the Hermite polynomial (Fast)
- or eigenvalue of the jacobi matrix (Slow) to obtain the nodes and
- weights, respectively.
-
- Returns
- -------
- nodes : ndarray
- roots
- weights : ndarray
- weights
-
-
- Example
- -------
- Integral of exp(x) from a = 0 to b = 3 is: exp(3)-exp(0)=
- >>> import numpy as np
- >>> nodes, weights = p_roots(11, a=0, b=3)
- >>> np.allclose(np.sum(np.exp(nodes) * weights), 19.085536923187668)
- True
- >>> nodes, weights = p_roots(11, method='newton1', a=0, b=3)
- >>> np.allclose(np.sum(np.exp(nodes) * weights), 19.085536923187668)
- True
- >>> nodes, weights = p_roots(11, method='eigenvalue', a=0, b=3)
- >>> np.allclose(np.sum(np.exp(nodes) * weights), 19.085536923187668)
- True
-
- See also
- --------
- quadg.
-
-
- References
- ----------
- [1] Davis and Rabinowitz (1975) 'Methods of Numerical Integration',
- page 365, Academic Press.
-
- [2] Golub, G. H. and Welsch, J. H. (1969)
- 'Calculation of Gaussian Quadrature Rules'
- Mathematics of Computation, vol 23,page 221-230,
-
- [3] Stroud and Secrest (1966), 'gaussian quadrature formulas',
- prentice-hall, Englewood cliffs, n.j.
- """
-
- if not method.startswith('n'):
- nodes, weights = ort.p_roots(n)
- else:
- if method.endswith('1'):
- nodes, weights = _p_roots_newton1(n)
- else:
- nodes, weights = _p_roots_newton(n)
-
- if (a != -1) | (b != 1):
- # Linear map from[-1,1] to [a,b]
- d_h = (b - a) / 2
- nodes = d_h * (nodes + 1) + a
- weights = weights * d_h
-
- return nodes, weights
-
-
-def q5_roots(n):
- """
- 5 : p(x) = 1/sqrt((x-a)*(b-x)), a =-1, b = 1 Chebyshev 1'st kind
- """
- j = np.arange(1, n + 1)
- weights = ones(n) * np.pi / n
- nodes = np.cos((2 * j - 1) * np.pi / (2 * n))
- return nodes, weights
-
-
-def q6_roots(n):
- """
- 6 : p(x) = sqrt((x-a)*(b-x)), a =-1, b = 1 Chebyshev 2'nd kind
- """
- j = np.arange(1, n + 1)
- x_j = j * np.pi / (n + 1)
- weights = np.pi / (n + 1) * np.sin(x_j) ** 2
- nodes = np.cos(x_j)
- return nodes, weights
-
-
-def q7_roots(n):
- """
- 7 : p(x) = sqrt((x-a)/(b-x)), a = 0, b = 1
- """
- j = np.arange(1, n + 1)
- x_j = (j - 0.5) * pi / (2 * n + 1)
- nodes = np.cos(x_j) ** 2
- weights = 2 * np.pi * nodes / (2 * n + 1)
- return nodes, weights
-
-
-def q8_roots(n):
- """
- 8 : p(x) = 1/sqrt(b-x), a = 0, b = 1
- """
- nodes_1, weights_1 = p_roots(2 * n)
- k, = np.where(0 <= nodes_1)
- weights = 2 * weights_1[k]
- nodes = 1 - nodes_1[k] ** 2
- return nodes, weights
-
-
-def q9_roots(n):
- """
- 9 : p(x) = sqrt(b-x), a = 0, b = 1
- """
- nodes_1, weights_1 = p_roots(2 * n + 1)
- k, = np.where(0 < nodes_1)
- weights = 2 * nodes_1[k] ** 2 * weights_1[k]
- nodes = 1 - nodes_1[k] ** 2
- return nodes, weights
-
-
-def qrule(n, wfun=1, alpha=0, beta=0):
- """
- Return nodes and weights for Gaussian quadratures.
-
- Parameters
- ----------
- n : integer
- number of base points
- wfun : integer
- defining the weight function, p(x). (default wfun = 1)
- 1 : p(x) = 1 a =-1, b = 1 Gauss-Legendre
- 2 : p(x) = exp(-x^2) a =-inf, b = inf Hermite
- 3 : p(x) = x^alpha*exp(-x) a = 0, b = inf Laguerre
- 4 : p(x) = (x-a)^alpha*(b-x)^beta a =-1, b = 1 Jacobi
- 5 : p(x) = 1/sqrt((x-a)*(b-x)), a =-1, b = 1 Chebyshev 1'st kind
- 6 : p(x) = sqrt((x-a)*(b-x)), a =-1, b = 1 Chebyshev 2'nd kind
- 7 : p(x) = sqrt((x-a)/(b-x)), a = 0, b = 1
- 8 : p(x) = 1/sqrt(b-x), a = 0, b = 1
- 9 : p(x) = sqrt(b-x), a = 0, b = 1
-
- Returns
- -------
- bp = base points (abscissas)
- wf = weight factors
-
- The Gaussian Quadrature integrates a (2n-1)th order
- polynomial exactly and the integral is of the form
- b n
- Int ( p(x)* F(x) ) dx = Sum ( wf_j* F( bp_j ) )
- a j=1
- where p(x) is the weight function.
- For Jacobi and Laguerre: alpha, beta >-1 (default alpha=beta=0)
-
- Examples:
- ---------
- >>> import numpy as np
-
- # integral of x^2 from a = -1 to b = 1
- >>> bp, wf = qrule(10)
- >>> np.allclose(sum(bp**2*wf), 0.66666666666666641)
- True
-
- # integral of exp(-x**2)*x**2 from a = -inf to b = inf
- >>> bp, wf = qrule(10,2)
- >>> np.allclose(sum(bp ** 2 * wf), 0.88622692545275772)
- True
-
- # integral of (x+1)*(1-x)**2 from a = -1 to b = 1
- >>> bp, wf = qrule(10,4,1,2)
- >>> np.allclose((bp * wf).sum(), 0.26666666666666755)
- True
-
- See also
- --------
- gaussq
-
- Reference
- ---------
- Abromowitz and Stegun (1954)
- (for method 5 to 9)
- """
-
- if wfun == 3: # Generalized Laguerre
- return la_roots(n, alpha)
- if wfun == 4: # Gauss-Jacobi
- return j_roots(n, alpha, beta)
-
- _assert(0 < wfun < 10, 'unknown weight function')
-
- root_fun = [None, p_roots, h_roots, la_roots, j_roots, q5_roots, q6_roots,
- q7_roots, q8_roots, q9_roots][wfun]
- return root_fun(n)
-
-
-class _Gaussq(object):
- """
- Numerically evaluate integral, Gauss quadrature.
-
- Parameters
- ----------
- fun : callable
- a,b : array-like
- lower and upper integration limits, respectively.
- releps, abseps : real scalars, optional
- relative and absolute tolerance, respectively.
- (default releps=abseps=1e-3).
- wfun : scalar integer, optional
- defining the weight function, p(x). (default wfun = 1)
- 1 : p(x) = 1 a =-1, b = 1 Gauss-Legendre
- 2 : p(x) = exp(-x^2) a =-inf, b = inf Hermite
- 3 : p(x) = x^alpha*exp(-x) a = 0, b = inf Laguerre
- 4 : p(x) = (x-a)^alpha*(b-x)^beta a =-1, b = 1 Jacobi
- 5 : p(x) = 1/sqrt((x-a)*(b-x)), a =-1, b = 1 Chebyshev 1'st kind
- 6 : p(x) = sqrt((x-a)*(b-x)), a =-1, b = 1 Chebyshev 2'nd kind
- 7 : p(x) = sqrt((x-a)/(b-x)), a = 0, b = 1
- 8 : p(x) = 1/sqrt(b-x), a = 0, b = 1
- 9 : p(x) = sqrt(b-x), a = 0, b = 1
- trace : bool, optional
- If non-zero a point plot of the integrand (default False).
- gn : scalar integer
- number of base points to start the integration with (default 2).
- alpha, beta : real scalars, optional
- Shape parameters of Laguerre or Jacobi weight function
- (alpha,beta>-1) (default alpha=beta=0)
-
- Returns
- -------
- val : ndarray
- evaluated integral
- err : ndarray
- error estimate, absolute tolerance abs(int-intold)
-
- Notes
- -----
- GAUSSQ numerically evaluate integral using a Gauss quadrature.
- The Quadrature integrates a (2m-1)th order polynomial exactly and the
- integral is of the form
- b
- Int (p(x)* Fun(x)) dx
- a
- GAUSSQ is vectorized to accept integration limits A, B and
- coefficients P1,P2,...Pn, as matrices or scalars and the
- result is the common size of A, B and P1,P2,...,Pn.
-
- Examples
- ---------
- integration of x**2 from 0 to 2 and from 1 to 4
-
- >>> import numpy as np
- >>> A = [0, 1]
- >>> B = [2, 4]
- >>> fun = lambda x: x**2
- >>> val1, err1 = gaussq(fun,A,B)
- >>> np.allclose(val1, [ 2.6666667, 21. ])
- True
- >>> np.allclose(err1, [ 1.7763568e-15, 1.0658141e-14])
- True
-
- Integration of x^2*exp(-x) from zero to infinity:
- >>> fun2 = lambda x : np.ones(np.shape(x))
- >>> val2, err2 = gaussq(fun2, 0, np.inf, wfun=3, alpha=2)
- >>> val3, err3 = gaussq(lambda x: x**2,0, np.inf, wfun=3, alpha=0)
- >>> np.allclose(val2, 2), err2[0] < 1e-14
- (True, True)
- >>> np.allclose(val3, 2), err3[0] < 1e-14
- (True, True)
-
- Integrate humps from 0 to 2 and from 1 to 4
- >>> val4, err4 = gaussq(humps, A, B, trace=True)
-
- See also
- --------
- qrule
- gaussq2d
- """
-
- @staticmethod
- def _get_dx(wfun, jacob, alpha, beta):
- def fun1(x):
- return x
- if wfun == 4:
- d_x = jacob ** (alpha + beta + 1)
- else:
- d_x = [None, fun1, fun1, fun1, None, lambda x: ones(np.shape(x)),
- lambda x: x ** 2, fun1, sqrt,
- lambda x: sqrt(x) ** 3][wfun](jacob)
- return d_x.ravel()
-
- @staticmethod
- def _nodes_and_weights(num_nodes, wfun, alpha, beta):
- global _NODES_AND_WEIGHTS
- name = 'wfun{:d}_{:d}_{:g}_{:g}'.format(wfun, num_nodes, alpha, beta)
- nodes_and_weights = _NODES_AND_WEIGHTS[name]
- if len(nodes_and_weights) == 0:
- nodes_and_weights.extend(qrule(num_nodes, wfun, alpha, beta))
- nodes, weights = nodes_and_weights
- return nodes, weights
-
- def _initialize_trace(self, max_iter):
- if self.trace:
- self.x_trace = [0] * max_iter
- self.y_trace = [0] * max_iter
-
- def _plot_trace(self, x, y):
- if self.trace:
- self.x_trace.append(x.ravel())
- self.y_trace.append(y.ravel())
- hfig = plt.plot(x, y, 'r.')
- plt.setp(hfig, 'color', 'b')
-
- def _plot_final_trace(self):
- if self.trace > 0:
- plt.clf()
- plt.plot(np.hstack(self.x_trace), np.hstack(self.y_trace), '+')
-
- @staticmethod
- def _get_jacob(wfun, a, b):
- if wfun in [2, 3]:
- jacob = ones((np.size(a), 1))
- else:
- jacob = (b - a) * 0.5
- if wfun in [7, 8, 9]:
- jacob *= 2
- return jacob
-
- @staticmethod
- def _warn_msg(k, a_shape):
- n = len(k)
- if n > 1:
- if n == np.prod(a_shape):
- msg = 'All integrals did not converge'
- else:
- msg = '%d integrals did not converge' % (n, )
- return msg + '--singularities likely!'
- return 'Integral did not converge--singularity likely!'
-
- @staticmethod
- def _initialize(wfun, a, b, args):
- args = np.broadcast_arrays(*np.atleast_1d(a, b, *args))
- a_shape = args[0].shape
- args = [np.reshape(x, (-1, 1)) for x in args]
- a_out, b_out = args[:2]
- args = args[2:]
- if wfun in [2, 3]:
- a_out = zeros((a_out.size, 1))
- return a_out, b_out, args, a_shape
-
- @staticmethod
- def _revert_nans_with_old(val, val_old):
- if any(np.isnan(val)):
- val[np.isnan(val)] = val_old[np.isnan(val)]
-
- @staticmethod
- def _update_error(i, abserr, val, val_old, k):
- if i > 1:
- abserr[k] = abs(val_old[k] - val[k]) # absolute tolerance
-
- def __call__(self, fun, a, b, releps=1e-3, abseps=1e-3, alpha=0, beta=0,
- wfun=1, trace=False, args=(), max_iter=11):
- self.trace = trace
- num_nodes = 2
-
- a_0, b_0, args, a_shape = self._initialize(wfun, a, b, args)
-
- jacob = self._get_jacob(wfun, a_0, b_0)
- shift = int(wfun in [1, 4, 5, 6])
- d_x = self._get_dx(wfun, jacob, alpha, beta)
-
- self._initialize_trace(max_iter)
-
- # Break out of the iteration loop for three reasons:
- # 1) the last update is very small (compared to int and to releps)
- # 2) There are more than 11 iterations. This should NEVER happen.
- dtype = np.result_type(fun((a_0 + b_0) * 0.5, *args))
- n_k = np.prod(a_shape) # # of integrals we have to compute
- k = np.arange(n_k)
- opt = (n_k, dtype)
- val, val_old, abserr = zeros(*opt), np.nan * ones(*opt), 1e100 * ones(*opt)
- nodes_and_weights = self._nodes_and_weights
- for i in range(max_iter):
- x_n, weights = nodes_and_weights(num_nodes, wfun, alpha, beta)
- x = (x_n + shift) * jacob[k, :] + a_0[k, :]
-
- params = [xi[k, :] for xi in args]
- y = fun(x, *params)
- self._plot_trace(x, y)
- val[k] = np.sum(weights * y, axis=1) * d_x[k] # do the integration
- self._revert_nans_with_old(val, val_old)
- self._update_error(i, abserr, val, val_old, k)
-
- k, = np.where(abserr > np.maximum(abs(releps * val), abseps))
- converged = len(k) == 0
- if converged:
- break
- val_old[k] = val[k]
- num_nodes *= 2 # double the # of basepoints and weights
-
- _assert_warn(converged, self._warn_msg(k, a_shape))
-
- # make sure int is the same size as the integration limits
- val.shape = a_shape
- abserr.shape = a_shape
-
- self._plot_final_trace()
- return val, abserr
-
-
-gaussq = _Gaussq()
-
-
-def richardson(q_val, k):
- # license BSD
- # Richardson extrapolation with parameter estimation
- c = np.real((q_val[k - 1] - q_val[k - 2]) / (q_val[k] - q_val[k - 1])) - 1.
- # The lower bound 0.07 admits the singularity x.^-0.9
- c = max(c, 0.07)
- return q_val[k] + (q_val[k] - q_val[k - 1]) / c
-
-
-class _Quadgr(object):
- """
- Gauss-Legendre quadrature with Richardson extrapolation.
-
- [q_val,ERR] = QUADGR(FUN,A,B,TOL) approximates the integral of a function
- FUN from A to B with an absolute error tolerance TOL. FUN is a function
- handle and must accept vector arguments. TOL is 1e-6 by default. q_val is
- the integral approximation and ERR is an estimate of the absolute
- error.
-
- QUADGR uses a 12-point Gauss-Legendre quadrature. The error estimate is
- based on successive interval bisection. Richardson extrapolation
- accelerates the convergence for some integrals, especially integrals
- with endpoint singularities.
-
- Examples
- --------
- >>> import numpy as np
- >>> q_val, err = quadgr(np.log,0,1)
- >>> q, err = quadgr(np.exp,0,9999*1j*np.pi)
- >>> np.allclose(q, -2.0000000000122662), err < 1.0e-08
- (True, True)
-
- >>> q, err = quadgr(lambda x: np.sqrt(4-x**2), 0, 2, abseps=1e-12)
- >>> np.allclose(q, 3.1415926535897811), err < 1.0e-12
- (True, True)
-
- >>> q, err = quadgr(lambda x: np.sqrt(4-x**2), 0, 0, abseps=1e-12)
- >>> np.allclose(q, 0), err < 1.0e-12
- (True, True)
-
- >>> q, err = quadgr(lambda x: x**-0.75, 0, 1)
- >>> np.allclose(q, 4), err < 1.e-13
- (True, True)
-
- >>> q, err = quadgr(lambda x: 1./np.sqrt(1-x**2), -1, 1)
- >>> np.allclose(q, 3.141596056985029), err < 1.0e-05
- (True, True)
-
- >>> q, err = quadgr(lambda x: np.exp(-x**2), -np.inf, np.inf, 1e-9)
- >>> np.allclose(q, np.sqrt(np.pi)), err < 1e-9
- (True, True)
-
- >>> q, err = quadgr(lambda x: np.cos(x)*np.exp(-x), 0, np.inf, 1e-9)
- >>> np.allclose(q, 0.5), err < 1e-9
- (True, True)
- >>> q, err = quadgr(lambda x: np.cos(x)*np.exp(-x), np.inf, 0, 1e-9)
- >>> np.allclose(q, -0.5), err < 1e-9
- (True, True)
- >>> q, err = quadgr(lambda x: np.cos(x)*np.exp(x), -np.inf, 0, 1e-9)
- >>> np.allclose(q, 0.5), err < 1e-9
- (True, True)
-
- See also
- --------
- QUAD,
- QUADGK
- """
- # Author: jonas.lundgren@saabgroup.com, 2009. license BSD
- # Order limits (required if infinite limits)
-
- def _change_variable_and_integrate(self, fun, a, b, abseps, max_iter):
- isreal = np.isreal(a) & np.isreal(b) & ~np.isnan(a) & ~np.isnan(b)
- _assert(isreal, 'Infinite intervals must be real.')
- integrate = self._integrate
- # Change of variable
- if np.isfinite(a) & np.isinf(b): # a to inf
- val, err = integrate(lambda t: fun(a + t / (1 - t)) / (1 - t) ** 2,
- 0, 1, abseps, max_iter)
- elif np.isinf(a) & np.isfinite(b): # -inf to b
- val, err = integrate(lambda t: fun(b + t / (1 + t)) / (1 + t) ** 2,
- -1, 0, abseps, max_iter)
- else: # -inf to inf
- val1, err1 = integrate(lambda t: fun(t / (1 - t)) / (1 - t) ** 2,
- 0, 1, abseps / 2, max_iter)
- val2, err2 = integrate(lambda t: fun(t / (1 + t)) / (1 + t) ** 2,
- -1, 0, abseps / 2, max_iter)
- val = val1 + val2
- err = err1 + err2
- return val, err
-
- @staticmethod
- def _nodes_and_weights():
- # Gauss-Legendre quadrature (12-point)
- x = np.asarray(
- [0.12523340851146894, 0.36783149899818018, 0.58731795428661748,
- 0.76990267419430469, 0.9041172563704748, 0.98156063424671924])
- w = np.asarray(
- [0.24914704581340288, 0.23349253653835478, 0.20316742672306584,
- 0.16007832854334636, 0.10693932599531818, 0.047175336386511842])
- nodes = np.hstack((x, -x))
- weights = np.hstack((w, w))
- return nodes, weights
-
- @staticmethod
- def _get_best_estimate(k, vals0, vals1, vals2):
- if k >= 6:
- q_v = np.hstack((vals0[k], vals1[k], vals2[k]))
- q_w = np.hstack((vals0[k - 1], vals1[k - 1], vals2[k - 1]))
- elif k >= 4:
- q_v = np.hstack((vals0[k], vals1[k]))
- q_w = np.hstack((vals0[k - 1], vals1[k - 1]))
- else:
- q_v = np.atleast_1d(vals0[k])
- q_w = vals0[k - 1]
- # Estimate absolute error
- errors = np.atleast_1d(abs(q_v - q_w))
- j = errors.argmin()
- err = errors[j]
- q_val = q_v[j]
-# if k >= 2: # and not iscomplex:
-# _val, err1 = dea3(vals0[k - 2], vals0[k - 1], vals0[k])
- return q_val, err
-
- def _extrapolate(self, k, val0, val1, val2):
- # Richardson extrapolation
- if k >= 5:
- val1[k] = richardson(val0, k)
- val2[k] = richardson(val1, k)
- elif k >= 3:
- val1[k] = richardson(val0, k)
- q_val, err = self._get_best_estimate(k, val0, val1, val2)
- return q_val, err
-
- def _integrate(self, fun, a, b, abseps, max_iter):
- dtype = np.result_type(fun((a + b) / 2), fun((a + b) / 4))
-
- # Initiate vectors
- val0 = zeros(max_iter, dtype=dtype) # Quadrature
- val1 = zeros(max_iter, dtype=dtype) # First Richardson extrapolation
- val2 = zeros(max_iter, dtype=dtype) # Second Richardson extrapolation
-
- x_n, weights = self._nodes_and_weights()
- n = len(x_n)
- # One interval
- d_x = (b - a) / 2 # Half interval length
- x = (a + b) / 2 + d_x * x_n # Nodes
- # Quadrature
- val0[0] = d_x * np.sum(weights * fun(x), axis=0)
-
- # Successive bisection of intervals
- for k in range(1, max_iter):
-
- # Interval bisection
- d_x = d_x / 2
- x = np.hstack([x + a, x + b]) / 2
- # Quadrature
- val0[k] = np.sum(np.sum(np.reshape(fun(x), (-1, n)), axis=0) *
- weights, axis=0) * d_x
-
- q_val, err = self._extrapolate(k, val0, val1, val2)
-
- converged = (err < abseps) | ~np.isfinite(q_val)
- if converged:
- break
-
- _assert_warn(converged, 'Max number of iterations reached without '
- 'convergence.')
- _assert_warn(np.isfinite(q_val),
- 'Integral approximation is Infinite or NaN.')
-
- # The error estimate should not be zero
- err = err + 2 * np.finfo(q_val).eps
- return q_val, err
-
- @staticmethod
- def _order_limits(a, b):
- if np.real(a) > np.real(b):
- return b, a, True
- return a, b, False
-
- def __call__(self, fun, a, b, abseps=1e-5, max_iter=17):
- a = np.asarray(a)
- b = np.asarray(b)
- if a == b:
- q_val = b - a
- err = np.abs(b - a)
- return q_val, err
-
- a, b, reverse = self._order_limits(a, b)
-
- improper_integral = np.isinf(a) | np.isinf(b)
- if improper_integral: # Infinite limits
- q_val, err = self._change_variable_and_integrate(fun, a, b, abseps,
- max_iter)
- else:
- q_val, err = self._integrate(fun, a, b, abseps, max_iter)
-
- # Reverse direction
- if reverse:
- q_val = -q_val
-
- return q_val, err
-
-
-quadgr = _Quadgr()
-
-
-def boole(y, x):
- a, b = x[0], x[-1]
- n = len(x)
- h = (b - a) / (n - 1)
- return (2 * h / 45) * (7 * (y[0] + y[-1]) + 12 * np.sum(y[2:n - 1:4]) +
- 32 * np.sum(y[1:n - 1:2]) +
- 14 * np.sum(y[4:n - 3:4]))
-
-
-def _plot_error(neval, err_dic, plot_error):
- if plot_error:
- plt.figure(0)
- for name in err_dic:
- plt.loglog(neval, err_dic[name], label=name)
-
- plt.xlabel('number of function evaluations')
- plt.ylabel('error')
- plt.legend()
-
-
-def _print_headers(formats_h, headers, names):
- print(''.join(fi % t for (fi, t) in zip(formats_h,
- ['ftn'] + names)))
- print(' '.join(headers))
-
-
-def _stack_values_and_errors(neval, vals_dic, err_dic, names):
- data = [neval]
- for name in names:
- data.append(vals_dic[name])
- data.append(err_dic[name])
-
- data = np.vstack(tuple(data)).T
- return data
-
-
-def _print_data(formats, data):
- for row in data:
- print(''.join(fi % t for (fi, t) in zip(formats, row.tolist())))
-
-
-def _print_values_and_errors(neval, vals_dic, err_dic):
- names = sorted(vals_dic.keys())
- num_cols = 2
- formats = ['%4.0f, '] + ['%10.10f, '] * num_cols * 2
- formats[-1] = formats[-1].split(',')[0]
- formats_h = ['%4s, '] + ['%20s, '] * num_cols
- formats_h[-1] = formats_h[-1].split(',')[0]
- headers = ['evals'] + ['%12s %12s' % ('approx', 'error')] * num_cols
- while len(names) > 0:
- names_c = names[:num_cols]
- _print_headers(formats_h, headers, names_c)
- data = _stack_values_and_errors(neval, vals_dic, err_dic, names_c)
- _print_data(formats, data)
- names = names[num_cols:]
-
-
-def _display(neval, vals_dic, err_dic, plot_error):
- # display results
- _print_values_and_errors(neval, vals_dic, err_dic)
- _plot_error(neval, err_dic, plot_error)
-
-
-def chebychev(y, x, n=None):
- if n is None:
- n = len(y)
- c_k = np.polynomial.chebyshev.chebfit(x, y, deg=min(n - 1, 36))
- c_ki = np.polynomial.chebyshev.chebint(c_k)
- q = np.polynomial.chebyshev.chebval(x[-1], c_ki)
- return q
-
-
-def qdemo(f, a, b, kmax=9, plot_error=False):
- """
- Compares different quadrature rules.
-
- Parameters
- ----------
- f : callable
- function
- a,b : scalars
- lower and upper integration limits
-
- Details
- -------
- qdemo(f,a,b) computes and compares various approximations to
- the integral of f from a to b. Three approximations are used,
- the composite trapezoid, Simpson's, and Boole's rules, all with
- equal length subintervals.
- In a case like qdemo(exp,0,3) one can see the expected
- convergence rates for each of the three methods.
- In a case like qdemo(sqrt,0,3), the convergence rate is limited
- not by the method, but by the singularity of the integrand.
-
- Example
- -------
- >>> import numpy as np
- >>> qdemo(np.exp,0,3, plot_error=True)
- true value = 19.08553692
- ftn, Boole, Chebychev
- evals approx error approx error
- 3, 19.4008539142, 0.3153169910, 19.5061466023, 0.4206096791
- 5, 19.0910191534, 0.0054822302, 19.0910191534, 0.0054822302
- 9, 19.0856414320, 0.0001045088, 19.0855374134, 0.0000004902
- 17, 19.0855386464, 0.0000017232, 19.0855369232, 0.0000000000
- 33, 19.0855369505, 0.0000000273, 19.0855369232, 0.0000000000
- 65, 19.0855369236, 0.0000000004, 19.0855369232, 0.0000000000
- 129, 19.0855369232, 0.0000000000, 19.0855369232, 0.0000000000
- 257, 19.0855369232, 0.0000000000, 19.0855369232, 0.0000000000
- 513, 19.0855369232, 0.0000000000, 19.0855369232, 0.0000000000
- ftn, Clenshaw-Curtis, Gauss-Legendre
- evals approx error approx error
- 3, 19.5061466023, 0.4206096791, 19.0803304585, 0.0052064647
- 5, 19.0834145766, 0.0021223465, 19.0855365951, 0.0000003281
- 9, 19.0855369150, 0.0000000082, 19.0855369232, 0.0000000000
- 17, 19.0855369232, 0.0000000000, 19.0855369232, 0.0000000000
- 33, 19.0855369232, 0.0000000000, 19.0855369232, 0.0000000000
- 65, 19.0855369232, 0.0000000000, 19.0855369232, 0.0000000000
- 129, 19.0855369232, 0.0000000000, 19.0855369232, 0.0000000000
- 257, 19.0855369232, 0.0000000000, 19.0855369232, 0.0000000000
- 513, 19.0855369232, 0.0000000000, 19.0855369232, 0.0000000000
- ftn, Simps, Trapz
- evals approx error approx error
- 3, 19.5061466023, 0.4206096791, 22.5366862979, 3.4511493747
- 5, 19.1169646189, 0.0314276957, 19.9718950387, 0.8863581155
- 9, 19.0875991312, 0.0020622080, 19.3086731081, 0.2231361849
- 17, 19.0856674267, 0.0001305035, 19.1414188470, 0.0558819239
- 33, 19.0855451052, 0.0000081821, 19.0995135407, 0.0139766175
- 65, 19.0855374350, 0.0000005118, 19.0890314614, 0.0034945382
- 129, 19.0855369552, 0.0000000320, 19.0864105817, 0.0008736585
- 257, 19.0855369252, 0.0000000020, 19.0857553393, 0.0002184161
- 513, 19.0855369233, 0.0000000001, 19.0855915273, 0.0000546041
- """
- true_val, _tol = intg.quad(f, a, b)
- print('true value = %12.8f' % (true_val,))
- neval = zeros(kmax, dtype=int)
- vals_dic = {}
- err_dic = {}
-
- # try various approximations
- methods = [trapz, simps, boole, chebychev]
-
- for k in range(kmax):
- n = 2 ** (k + 1) + 1
- neval[k] = n
- x = np.linspace(a, b, n)
- y = f(x)
- for method in methods:
- name = method.__name__.title()
- q = method(y, x)
- vals_dic.setdefault(name, []).append(q)
- err_dic.setdefault(name, []).append(abs(q - true_val))
-
- name = 'Clenshaw-Curtis'
- q = clencurt(f, a, b, (n - 1) // 2)[0]
- vals_dic.setdefault(name, []).append(q[0])
- err_dic.setdefault(name, []).append(abs(q[0] - true_val))
-
- name = 'Gauss-Legendre' # quadrature
- q = intg.fixed_quad(f, a, b, n=n)[0]
- vals_dic.setdefault(name, []).append(q)
- err_dic.setdefault(name, []).append(abs(q - true_val))
-
- _display(neval, vals_dic, err_dic, plot_error)
-
-
-def main():
- # val, err = clencurt(np.exp, 0, 2)
- # valt = np.exp(2) - np.exp(0)
- # [Q, err] = quadgr(lambda x: x ** 2, 1, 4, 1e-9)
- # [Q, err] = quadgr(humps, 1, 4, 1e-9)
- #
- # [x, w] = h_roots(11, 'newton')
- # sum(w)
- # [x2, w2] = la_roots(11, 1, 't')
- #
- # from scitools import numpyutils as npu #@UnresolvedImport
- # fun = npu.wrap2callable('x**2')
- # p0 = fun(0)
- # A = [0, 1, 1]; B = [2, 4, 3]
- # area, err = gaussq(fun, A, B)
- #
- # fun = npu.wrap2callable('x**2')
- # [val1, err1] = gaussq(fun, A, B)
- #
- #
- # Integration of x^2*exp(-x) from zero to infinity:
- # fun2 = npu.wrap2callable('1')
- # [val2, err2] = gaussq(fun2, 0, np.inf, wfun=3, alpha=2)
- # [val2, err2] = gaussq(lambda x: x ** 2, 0, np.inf, wfun=3, alpha=0)
- #
- # Integrate humps from 0 to 2 and from 1 to 4
- # [val3, err3] = gaussq(humps, A, B)
- #
- # [x, w] = p_roots(11, 'newton', 1, 3)
- # y = np.sum(x ** 2 * w)
-
- x = np.linspace(0, np.pi / 2)
- _q0 = np.trapz(humps(x), x)
- [q, err] = romberg(humps, 0, np.pi / 2, 1e-4)
- print(q, err)
-
-
-if __name__ == '__main__':
- from wafo.testing import test_docstrings
- test_docstrings(__file__)
- # qdemo(np.exp, 0, 3, plot_error=True)
- # plt.show('hold')
- # main()
diff --git a/wafo/integrate_oscillating.py b/wafo/integrate_oscillating.py
deleted file mode 100644
index 81217ef..0000000
--- a/wafo/integrate_oscillating.py
+++ /dev/null
@@ -1,538 +0,0 @@
-"""
-Created on 20. aug. 2015
-
-@author: pab
-"""
-from __future__ import division
-from collections import namedtuple
-import warnings
-import numdifftools as nd
-import numdifftools.nd_algopy as nda
-from numdifftools.extrapolation import dea3
-from numdifftools.limits import Limit
-import numpy as np
-from numpy import linalg
-from numpy.polynomial.chebyshev import chebval, Chebyshev
-from numpy.polynomial import polynomial
-from wafo.misc import piecewise, findcross, ecross
-
-_FINFO = np.finfo(float)
-EPS = _FINFO.eps
-_EPS = EPS
-_TINY = _FINFO.tiny
-_HUGE = _FINFO.max
-
-
-def _assert(cond, msg):
- if not cond:
- raise ValueError(msg)
-
-
-def _assert_warn(cond, msg):
- if not cond:
- warnings.warn(msg)
-
-
-class PolyBasis(object):
- @staticmethod
- def _derivative(c, m):
- return polynomial.polyder(c, m)
-
- @staticmethod
- def eval(t, c):
- return polynomial.polyval(t, c)
-
- @staticmethod
- def _coefficients(k):
- c = np.zeros(k + 1)
- c[k] = 1
- return c
-
- def derivative(self, t, k, n=1):
- c = self._coefficients(k)
- d_c = self._derivative(c, m=n)
- return self.eval(t, d_c)
-
- def __call__(self, t, k):
- return t**k
-
-
-poly_basis = PolyBasis()
-
-
-class ChebyshevBasis(PolyBasis):
-
- @staticmethod
- def _derivative(c, m):
- cheb = Chebyshev(c)
- dcheb = cheb.deriv(m=m)
- return dcheb.coef
-
- @staticmethod
- def eval(t, c):
- return chebval(t, c)
-
- def __call__(self, t, k):
- c = self._coefficients(k)
- return self.eval(t, c)
-
-
-chebyshev_basis = ChebyshevBasis()
-
-
-def richardson(q_val, k):
- # license BSD
- # Richardson extrapolation with parameter estimation
- c = np.real((q_val[k - 1] - q_val[k - 2]) / (q_val[k] - q_val[k - 1])) - 1.
- # The lower bound 0.07 admits the singularity x.^-0.9
- c = max(c, 0.07)
- return q_val[k] + (q_val[k] - q_val[k - 1]) / c
-
-
-def evans_webster_weights(omega, g, d_g, x, basis, *args, **kwds):
-
- def psi(t, k):
- return d_g(t, *args, **kwds) * basis(t, k)
-
- j_w = 1j * omega
- n = len(x)
- a_matrix = np.zeros((n, n), dtype=complex)
- rhs = np.zeros((n,), dtype=complex)
-
- dbasis = basis.derivative
- lim_g = Limit(g)
- b_1 = np.exp(j_w * lim_g(1, *args, **kwds))
- if np.isnan(b_1):
- b_1 = 0.0
- a_1 = np.exp(j_w * lim_g(-1, *args, **kwds))
- if np.isnan(a_1):
- a_1 = 0.0
-
- lim_psi = Limit(psi)
- for k in range(n):
- rhs[k] = basis(1, k) * b_1 - basis(-1, k) * a_1
- a_matrix[k] = (dbasis(x, k, n=1) + j_w * lim_psi(x, k))
-
- solution = linalg.lstsq(a_matrix, rhs)
- return solution[0]
-
-
-def osc_weights(omega, g, d_g, x, basis, a_b, *args, **kwds):
- def _g(t):
- return g(scale * t + offset, *args, **kwds)
-
- def _d_g(t):
- return scale * d_g(scale * t + offset, *args, **kwds)
-
- w = []
-
- for a, b in zip(a_b[::2], a_b[1::2]):
- scale = (b - a) / 2
- offset = (a + b) / 2
-
- w.append(evans_webster_weights(omega, _g, _d_g, x, basis))
-
- return np.asarray(w).ravel()
-
-
-class _Integrator(object):
- info = namedtuple('info', ['error_estimate', 'n'])
-
- def __init__(self, f, g, dg=None, a=-1, b=1, basis=chebyshev_basis, s=1,
- precision=10, endpoints=True, full_output=False):
- self.f = f
- self.g = g
- self.dg = nd.Derivative(g) if dg is None else dg
- self.basis = basis
- self.a = a
- self.b = b
- self.s = s
- self.endpoints = endpoints
- self.precision = precision
- self.full_output = full_output
-
-
-class QuadOsc(_Integrator):
- def __init__(self, f, g, dg=None, a=-1, b=1, basis=chebyshev_basis, s=15,
- precision=10, endpoints=False, full_output=False, maxiter=17):
- self.maxiter = maxiter
- super(QuadOsc, self).__init__(f, g, dg=dg, a=a, b=b, basis=basis, s=s,
- precision=precision, endpoints=endpoints,
- full_output=full_output)
-
- @staticmethod
- def _change_interval_to_0_1(f, g, d_g, a, _b):
- def f_01(t, *args, **kwds):
- den = 1 - t
- return f(a + t / den, *args, **kwds) / den ** 2
-
- def g_01(t, *args, **kwds):
- return g(a + t / (1 - t), *args, **kwds)
-
- def d_g_01(t, *args, **kwds):
- den = 1 - t
- return d_g(a + t / den, *args, **kwds) / den ** 2
- return f_01, g_01, d_g_01, 0., 1.
-
- @staticmethod
- def _change_interval_to_m1_0(f, g, d_g, _a, b):
- def f_m10(t, *args, **kwds):
- den = 1 + t
- return f(b + t / den, *args, **kwds) / den ** 2
-
- def g_m10(t, *args, **kwds):
- return g(b + t / (1 + t), *args, **kwds)
-
- def d_g_m10(t, *args, **kwds):
- den = 1 + t
- return d_g(b + t / den, *args, **kwds) / den ** 2
- return f_m10, g_m10, d_g_m10, -1.0, 0.0
-
- @staticmethod
- def _change_interval_to_m1_1(f, g, d_g, _a, _b):
- def f_m11(t, *args, **kwds):
- den = (1 - t**2)
- return f(t / den, *args, **kwds) * (1 + t**2) / den ** 2
-
- def g_m11(t, *args, **kwds):
- den = (1 - t**2)
- return g(t / den, *args, **kwds)
-
- def d_g_m11(t, *args, **kwds):
- den = (1 - t**2)
- return d_g(t / den, *args, **kwds) * (1 + t**2) / den ** 2
- return f_m11, g_m11, d_g_m11, -1., 1.
-
- def _get_functions(self):
- a, b = self.a, self.b
- reverse = np.real(a) > np.real(b)
- if reverse:
- a, b = b, a
- f, g, dg = self.f, self.g, self.dg
-
- if a == b:
- pass
- elif np.isinf(a) | np.isinf(b):
- # Check real limits
- if ~np.isreal(a) | ~np.isreal(b) | np.isnan(a) | np.isnan(b):
- raise ValueError('Infinite intervals must be real.')
- # Change of variable
- if np.isfinite(a) & np.isinf(b):
- f, g, dg, a, b = self._change_interval_to_0_1(f, g, dg, a, b)
- elif np.isinf(a) & np.isfinite(b):
- f, g, dg, a, b = self._change_interval_to_m1_0(f, g, dg, a, b)
- else: # -inf to inf
- f, g, dg, a, b = self._change_interval_to_m1_1(f, g, dg, a, b)
-
- return f, g, dg, a, b, reverse
-
- def __call__(self, omega, *args, **kwds):
- f, g, dg, a, b, reverse = self._get_functions()
-
- val, err = self._quad_osc(f, g, dg, a, b, omega, *args, **kwds)
-
- if reverse:
- val = -val
- if self.full_output:
- return val, err
- return val
-
- @staticmethod
- def _get_best_estimate(k, q_0, q_1, q_2):
- if k >= 5:
- q_v = np.hstack((q_0[k], q_1[k], q_2[k]))
- q_w = np.hstack((q_0[k - 1], q_1[k - 1], q_2[k - 1]))
- elif k >= 3:
- q_v = np.hstack((q_0[k], q_1[k]))
- q_w = np.hstack((q_0[k - 1], q_1[k - 1]))
- else:
- q_v = np.atleast_1d(q_0[k])
- q_w = q_0[k - 1]
- errors = np.atleast_1d(abs(q_v - q_w))
- j = np.nanargmin(errors)
- return q_v[j], errors[j]
-
- def _extrapolate(self, k, q_0, q_1, q_2):
- if k >= 4:
- q_1[k] = dea3(q_0[k - 2], q_0[k - 1], q_0[k])[0]
- q_2[k] = dea3(q_1[k - 2], q_1[k - 1], q_1[k])[0]
- elif k >= 2:
- q_1[k] = dea3(q_0[k - 2], q_0[k - 1], q_0[k])[0]
- # # Richardson extrapolation
- # if k >= 4:
- # q_1[k] = richardson(q_0, k)
- # q_2[k] = richardson(q_1, k)
- # elif k >= 2:
- # q_1[k] = richardson(q_0, k)
- q, err = self._get_best_estimate(k, q_0, q_1, q_2)
- return q, err
-
- def _quad_osc(self, f, g, dg, a, b, omega, *args, **kwds):
- if a == b:
- q_val = b - a
- err = np.abs(b - a)
- return q_val, err
-
- abseps = 10**-self.precision
- max_iter = self.maxiter
- basis = self.basis
- if self.endpoints:
- x_n = chebyshev_extrema(self.s)
- else:
- x_n = chebyshev_roots(self.s)
- # x_n = tanh_sinh_open_nodes(self.s)
-
- # One interval
- hh = (b - a) / 2
- x = (a + b) / 2 + hh * x_n # Nodes
-
- dtype = complex
- val0 = np.zeros((max_iter, 1), dtype=dtype) # Quadrature
- val1 = np.zeros((max_iter, 1), dtype=dtype) # First extrapolation
- val2 = np.zeros((max_iter, 1), dtype=dtype) # Second extrapolation
-
- lim_f = Limit(f)
- a_b = np.hstack([a, b])
- wq = osc_weights(omega, g, dg, x_n, basis, a_b, *args, **kwds)
- val0[0] = hh * np.sum(wq * lim_f(x, *args, **kwds))
-
- # Successive bisection of intervals
- nq = len(x_n)
- n = nq
- for k in range(1, max_iter):
- n += nq * 2**k
-
- hh = hh / 2
- x = np.hstack([x + a, x + b]) / 2
- a_b = np.hstack([a_b + a, a_b + b]) / 2
- wq = osc_weights(omega, g, dg, x_n, basis, a_b, *args, **kwds)
-
- val0[k] = hh * np.sum(wq * lim_f(x, *args, **kwds))
-
- q_val, err = self._extrapolate(k, val0, val1, val2)
-
- converged = (err <= abseps) | ~np.isfinite(q_val)
- if converged:
- break
- _assert_warn(converged, 'Max number of iterations reached '
- 'without convergence.')
- _assert_warn(np.isfinite(q_val),
- 'Integral approximation is Infinite or NaN.')
-
- # The error estimate should not be zero
- err += 2 * np.finfo(q_val).eps
- return q_val, self.info(err, n)
-
-
-def adaptive_levin_points(m, delta):
- m_1 = m - 1
- prm = 0.5
- while prm * m_1 / delta >= 1:
- delta = 2 * delta
- k = np.arange(m)
- x = piecewise([k < prm * m_1, k == np.ceil(prm * m_1)],
- [-1 + k / delta, 0 * k, 1 - (m_1 - k) / delta])
- return x
-
-
-def open_levin_points(m, delta):
- return adaptive_levin_points(m + 2, delta)[1:-1]
-
-
-def chebyshev_extrema(m, delta=None):
- k = np.arange(m)
- x = np.cos(k * np.pi / (m - 1))
- return x
-
-
-def tanh_sinh_nodes(m, delta=None, tol=_EPS):
- tmax = np.arcsinh(np.arctanh(1 - _EPS) * 2 / np.pi)
- # tmax = 3.18
- m_1 = int(np.floor(-np.log2(tmax / max(m - 1, 1)))) - 1
- h = 2.0**-m_1
- t = np.arange((m + 1) // 2 + 1) * h
- x = np.tanh(np.pi / 2 * np.sinh(t))
- k = np.flatnonzero(np.abs(x - 1) <= 10 * tol)
- y = x[:k[0] + 1] if len(k) else x
- return np.hstack((-y[:0:-1], y))
-
-
-def tanh_sinh_open_nodes(m, delta=None, tol=_EPS):
- return tanh_sinh_nodes(m + 1, delta, tol)[1:-1]
-
-
-def chebyshev_roots(m, delta=None):
- k = np.arange(1, 2 * m, 2) * 0.5
- x = np.cos(k * np.pi / m)
- return x
-
-
-class AdaptiveLevin(_Integrator):
- """Return integral for the Levin-type and adaptive Levin-type methods"""
-
- @staticmethod
- def _a_levin(omega, f, g, d_g, x, s, basis, *args, **kwds):
-
- def psi(t, k):
- return d_g(t, *args, **kwds) * basis(t, k)
-
- j_w = 1j * omega
- nu = np.ones((len(x),), dtype=int)
- nu[0] = nu[-1] = s
- S = np.cumsum(np.hstack((nu, 0)))
- S[-1] = 0
- n = int(S[-2])
- a_matrix = np.zeros((n, n), dtype=complex)
- rhs = np.zeros((n,))
- dff = Limit(nda.Derivative(f))
- d_psi = Limit(nda.Derivative(psi))
- dbasis = basis.derivative
- for r, t in enumerate(x):
- for j in range(S[r - 1], S[r]):
- order = ((j - S[r - 1]) % nu[r]) # derivative order
- dff.fun.n = order
- rhs[j] = dff(t, *args, **kwds)
- d_psi.fun.n = order
- for k in range(n):
- a_matrix[j, k] = (dbasis(t, k, n=order + 1) +
- j_w * d_psi(t, k))
- k1 = np.flatnonzero(1 - np.isfinite(rhs))
- if k1.size > 0: # Remove singularities
- warnings.warn('Singularities detected! ')
- a_matrix[k1] = 0
- rhs[k1] = 0
- solution = linalg.lstsq(a_matrix, rhs)
- v = basis.eval([-1, 1], solution[0])
-
- lim_g = Limit(g)
- g_b = np.exp(j_w * lim_g(1, *args, **kwds))
- if np.isnan(g_b):
- g_b = 0
- g_a = np.exp(j_w * lim_g(-1, *args, **kwds))
- if np.isnan(g_a):
- g_a = 0
- return v[1] * g_b - v[0] * g_a
-
- def _get_integration_limits(self, omega, args, kwds):
- a, b = self.a, self.b
- M = 30
- ab = [a]
- scale = (b - a) / 2
- n = 30
- x = np.linspace(a, b, n + 1)
- dg_x = np.asarray([scale * omega * self.dg(xi, *args, **kwds)
- for xi in x])
- i10 = findcross(dg_x, M)
- i1 = findcross(dg_x, 1)
- i0 = findcross(dg_x, 0)
- im1 = findcross(dg_x, -1)
- im10 = findcross(dg_x, -M)
- x10 = ecross(x, dg_x, i10, M) if len(i10) else ()
- x1 = ecross(x, dg_x, i1, 1) if len(i1) else ()
- x0 = ecross(x, dg_x, i0, 0) if len(i0) else ()
- xm1 = ecross(x, dg_x, im1, -1) if len(im1) else ()
- xm10 = ecross(x, dg_x, im10, -M) if len(im10) else ()
-
- for i in np.unique(np.hstack((x10, x1, x0, xm1, xm10))):
- if x[0] < i < x[n]:
- ab.append(i)
- ab.append(b)
- return ab
-
- def __call__(self, omega, *args, **kwds):
- ab = self._get_integration_limits(omega, args, kwds)
- s = self.s
- val = 0
- n = 0
- err = 0
- for ai, bi in zip(ab[:-1], ab[1:]):
- vali, infoi = self._QaL(s, ai, bi, omega, *args, **kwds)
- val += vali
- err += infoi.error_estimate
- n += infoi.n
- if self.full_output:
- info = self.info(err, n)
- return val, info
- return val
-
- @staticmethod
- def _get_num_points(s, prec, betam):
- return 1 if s > 1 else int(prec / max(np.log10(betam + 1), 1) + 1)
-
- def _QaL(self, s, a, b, omega, *args, **kwds):
- """if s>1,the integral is computed by Q_s^L"""
- scale = (b - a) / 2
- offset = (a + b) / 2
- prec = self.precision # desired precision
-
- def ff(t, *args, **kwds):
- return scale * self.f(scale * t + offset, *args, **kwds)
-
- def gg(t, *args, **kwds):
- return self.g(scale * t + offset, *args, **kwds)
-
- def dgg(t, *args, **kwds):
- return scale * self.dg(scale * t + offset, *args, **kwds)
- dg_a = abs(omega * dgg(-1, *args, **kwds))
- dg_b = abs(omega * dgg(1, *args, **kwds))
- g_a = abs(omega * gg(-1, *args, **kwds))
- g_b = abs(omega * gg(1, *args, **kwds))
- delta, alpha = min(dg_a, dg_b), min(g_a, g_b)
-
- betam = delta # * scale
- if self.endpoints:
- if delta < 10 or alpha <= 10 or s > 1:
- points = chebyshev_extrema
- else:
- points = adaptive_levin_points
- elif delta < 10 or alpha <= 10 or s > 1:
- points = chebyshev_roots
- else:
- points = open_levin_points # tanh_sinh_open_nodes
-
- m = self._get_num_points(s, prec, betam)
- abseps = 10 * 10.0**-prec
- num_collocation_point_list = m * 2**np.arange(1, 5) + 1
- basis = self.basis
-
- q_val = 1e+300
- num_function_evaluations = 0
- n = 0
- for num_collocation_points in num_collocation_point_list:
- n_old = n
- q_old = q_val
- x = points(num_collocation_points, betam)
- n = len(x)
- if n > n_old:
- q_val = self._a_levin(omega, ff, gg, dgg, x, s, basis, *args,
- **kwds)
- num_function_evaluations += n
- err = np.abs(q_val - q_old)
- if err <= abseps:
- break
- info = self.info(err, num_function_evaluations)
- return q_val, info
-
-
-class EvansWebster(AdaptiveLevin):
- """Return integral for the Evans Webster method"""
-
- def __init__(self, f, g, dg=None, a=-1, b=1, basis=chebyshev_basis, s=8,
- precision=10, endpoints=False, full_output=False):
- super(EvansWebster,
- self).__init__(f, g, dg=dg, a=a, b=b, basis=basis, s=s,
- precision=precision, endpoints=endpoints,
- full_output=full_output)
-
- def _a_levin(self, omega, ff, gg, dgg, x, s, basis, *args, **kwds):
- w = evans_webster_weights(omega, gg, dgg, x, basis, *args, **kwds)
-
- f = Limit(ff)(x, *args, **kwds)
- return np.sum(f * w)
-
- def _get_num_points(self, s, prec, betam):
- return 8 if s > 1 else int(prec / max(np.log10(betam + 1), 1) + 1)
-
-
-if __name__ == '__main__':
- tanh_sinh_nodes(16)
diff --git a/wafo/padua.py b/wafo/padua.py
deleted file mode 100644
index a372248..0000000
--- a/wafo/padua.py
+++ /dev/null
@@ -1,532 +0,0 @@
-'''
-
-All the software contained in this library is protected by copyright.
-Permission to use, copy, modify, and distribute this software for any
-purpose without fee is hereby granted, provided that this entire notice
-is included in all copies of any software which is or includes a copy
-or modification of this software and in all copies of the supporting
-documentation for such software.
-
-THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
-WARRANTY. IN NO EVENT, NEITHER THE AUTHORS, NOR THE PUBLISHER, NOR ANY
-MEMBER OF THE EDITORIAL BOARD OF THE JOURNAL "NUMERICAL ALGORITHMS",
-NOR ITS EDITOR-IN-CHIEF, BE LIABLE FOR ANY ERROR IN THE SOFTWARE, ANY
-MISUSE OF IT OR ANY DAMAGE ARISING OUT OF ITS USE. THE ENTIRE RISK OF
-USING THE SOFTWARE LIES WITH THE PARTY DOING SO.
-
-ANY USE OF THE SOFTWARE CONSTITUTES ACCEPTANCE OF THE TERMS OF THE
-ABOVE STATEMENT.
-
-
-AUTHORS:
-Per A Brodtkorb
-Python code Based on matlab code written by:
-
-Marco Caliari
-University of Verona, Italy
-E-mail: marco.caliari@univr.it
-
-Stefano de Marchi, Alvise Sommariva, Marco Vianello
-University of Padua, Italy
-E-mail: demarchi@math.unipd.it, alvise@math.unipd.it,
- marcov@math.unipd.it
-
-Reference
----------
-Padua2DM: fast interpolation and cubature at the Padua points in Matlab/Octave
-NUMERICAL ALGORITHMS, 56 (2011), PP. 45-60
-
-
-Padua module
-------------
-In polynomial interpolation of two variables, the Padua points are the first
-known example (and up to now the only one) of a unisolvent point set
-(that is, the interpolating polynomial is unique) with minimal growth of their
-Lebesgue constant, proven to be O(log2 n).
-This module provides all the functions needed to perform interpolation and
-cubature at the Padua points, together with the functions and the demos used
-in the paper.
-
-pdint.m : main function for interpolation at the Padua points
-pdcub.m : main function for cubature at the Padua points
-pdpts.m : function for the computation of the Padua points
-padua_fit.m : function for the computation of the interpolation
- coefficients by FFT (recommended)
-pdcfsMM.m : function for the computation of the interpolation
- coefficients by matrix multiplications
-pdval.m : function for the evaluation of the interpolation
- polynomial
-pdwtsFFT.m : function for the computation of the cubature
- weights by FFT
-pdwtsMM.m : function for the computation of the cubature
- weights by matrix multiplications (recommended)
-funct.m : function containing some test functions
-demo_pdint.m : demo script for pdint
-demo_cputime_pdint.m : demo script for the computation of CPU time for
- interpolation
-demo_errors_pdint.m : demo script for the comparison of interpolation with
- coefficients computed by FFT or by matrix
- multiplications
-demo_pdcub : demo script for pdcub
-demo_cputime_pdcub.m : demo script for the computation of CPU time for
- cubature
-demo_errors_pdcub.m : demo script for the comparison of cubature with
- weights computed by FFT or by matrix multiplications
-demo_errors_pdcub_gl.m : demo script for the comparison of different cubature
- formulas
-cubature_square.m : function for the computation of some cubature
- formulas for the square
-omelyan_solovyan_rule.m : function for the computation of Omelyan-Solovyan
- cubature points and weights
-Contents.m : Contents file for Matlab
-
-
-'''
-from __future__ import absolute_import, division
-import numpy as np
-from numpy.fft import fft
-from wafo.dctpack import dct
-from wafo.polynomial import map_from_interval, map_to_interval
-# from scipy.fftpack.realtransforms import dct
-
-
-class _ExampleFunctions(object):
- '''
- Computes test function in the points (x, y)
-
- Parameters
- ----------
- x,y : array-like
- evaluate the function in the points (x,y)
- i : scalar int (default 0)
- defining which test function to use. Options are
- 0: franke
- 1: half_sphere
- 2: poly_degree20
- 3: exp_fun1
- 4: exp_fun100
- 5: cos30
- 6: constant
- 7: exp_xy
- 8: runge
- 9: abs_cubed
- 10: gauss
- 11: exp_inv
-
- Returns
- -------
- z : array-like
- value of the function in the points (x,y)
- '''
- @staticmethod
- def franke(x, y):
- '''Franke function.
-
- The value of the definite integral on the square [-1,1] x [-1,1],
- obtained using a Padua Points cubature formula of degree 500, is
- 2.1547794245591083e+000 with an estimated absolute error of 8.88e-016.
-
- The value of the definite integral on the square [0,1] x [0,1],
- obtained using a Padua Points cubature formula of degree 500, is
- 4.06969589491556e-01 with an estimated absolute error of 8.88e-016.
-
- Maple: 0.40696958949155611906
- '''
- def _exp(x, y, loc, scale, p2=2):
- return np.exp(- (x - loc[0])**2 / scale[0] - (y - loc[1])**p2 / scale[1])
- # exp = np.exp
- x9, y9 = 9. * x, 9. * y
- return (3. / 4 * _exp(x9, y9, [2, 2], [4, 4]) +
- 3. / 4 * _exp(x9, y9, [-1, -1], [49, 10], p2=1) +
- 1. / 2 * _exp(x9, y9, [7, 3], [4, 4]) -
- 1. / 5 * _exp(x9, y9, [4, 7], [1, 1]))
-
- @staticmethod
- def half_sphere(x, y):
- '''The value of the definite integral on the square [-1,1] x [-1,1],
- obtained using a Padua Points cubature formula of degree 2000, is
- 3.9129044444568244e+000 with an estimated absolute error of 3.22e-010.
- '''
- return ((x - 0.5)**2 + (y - 0.5)**2)**(1. / 2)
-
- @staticmethod
- def poly_degree20(x, y):
- ''''Bivariate polynomial having moderate degree.
- The value of the definite integral on the square [-1,1] x
- [-1,1], up to machine precision, is 18157.16017316017 (see ref. 6).
- The value of the definite integral on the square [-1,1] x [-1,1],
- obtained using a Padua Points cubature formula of degree 500,
- is 1.8157160173160162e+004.
-
- 2D modification of an example by L.N.Trefethen (see ref. 7), f(x)=x^20.
- '''
- return (x + y)**20
-
- @staticmethod
- def exp_fun1(x, y):
- ''' The value of the definite integral on the square [-1,1] x [-1,1],
- obtained using a Padua Points cubature formula of degree 2000,
- is 2.1234596326670683e+001 with an estimated absolute error of
- 7.11e-015.
- '''
- return np.exp((x - 0.5)**2 + (y - 0.5)**2)
-
- @staticmethod
- def exp_fun100(x, y):
- '''The value of the definite integral on the square [-1,1] x [-1,1],
- obtained using a Padua Points cubature formula of degree 2000,
- is 3.1415926535849605e-002 with an estimated absolute error of
- 3.47e-017.
- '''
- return np.exp(-100 * ((x - 0.5)**2 + (y - 0.5)**2))
-
- @staticmethod
- def cos30(x, y):
- ''' The value of the definite integral on the square [-1,1] x [-1,1],
- obtained using a Padua Points cubature formula of degree 500,
- is 4.3386955120336568e-003 with an estimated absolute error of
- 2.95e-017.
- '''
- return np.cos(30 * (x + y))
-
- @staticmethod
- def constant(x, y):
- '''Constant.
- To test interpolation and cubature at degree 0.
- The value of the definite integral on the square [-1,1] x [-1,1]
- is 4.
- '''
- return np.ones(np.shape(x + y))
-
- @staticmethod
- def exp_xy(x, y):
- '''The value of the definite integral on the square [-1,1] x [-1,1]
- is up to machine precision is 5.524391382167263 (see ref. 6).
- 2. The value of the definite integral on the square [-1,1] x [-1,1],
- obtained using a Padua Points cubature formula of degree 500,
- is 5.5243913821672628e+000 with an estimated absolute error of
- 0.00e+000.
- 2D modification of an example by L.N.Trefethen (see ref. 7),
- f(x)=exp(x).
- '''
- return np.exp(x + y)
-
- @staticmethod
- def runge(x, y):
- ''' Bivariate Runge function: as 1D complex function is analytic
- in a neighborhood of [-1; 1] but not throughout the complex plane.
-
- The value of the definite integral on the square [-1,1] x [-1,1],
- up to machine precision, is 0.597388947274307 (see ref. 6).
- The value of the definite integral on the square [-1,1] x [-1,1],
- obtained using a Padua Points cubature formula of degree 500,
- is 5.9738894727430725e-001 with an estimated absolute error of
- 0.00e+000.
-
- 2D modification of an example by L.N.Trefethen (see ref. 7),
- f(x)=1/(1+16*x^2).
- '''
- return 1. / (1 + 16 * (x**2 + y**2))
-
- @staticmethod
- def abs_cubed(x, y):
- '''Low regular function.
- The value of the definite integral on the square [-1,1] x [-1,1],
- up to machine precision, is 2.508723139534059 (see ref. 6).
- The value of the definite integral on the square [-1,1] x [-1,1],
- obtained using a Padua Points cubature formula of degree 500,
- is 2.5087231395340579e+000 with an estimated absolute error of
- 0.00e+000.
-
- 2D modification of an example by L.N.Trefethen (see ref. 7),
- f(x)=abs(x)^3.
- '''
- return (x**2 + y**2)**(3 / 2)
-
- @staticmethod
- def gauss(x, y):
- '''Bivariate gaussian: smooth function.
- The value of the definite integral on the square [-1,1] x [-1,1],
- up to machine precision, is 2.230985141404135 (see ref. 6).
- The value of the definite integral on the square [-1,1] x [-1,1],
- obtained using a Padua Points cubature formula of degree 500,
- is 2.2309851414041333e+000 with an estimated absolute error of
- 2.66e-015.
-
- 2D modification of an example by L.N.Trefethen (see ref. 7),
- f(x)=exp(-x^2).
- '''
- return np.exp(-x**2 - y**2)
-
- @staticmethod
- def exp_inv(x, y):
- '''Bivariate example stemming from a 1D C-infinity function.
- The value of the definite integral on the square [-1,1] x [-1,1],
- up to machine precision, is 0.853358758654305 (see ref. 6).
- The value of the definite integral on the square [-1,1] x [-1,1],
- obtained using a Padua Points cubature formula of degree 2000,
- is 8.5335875865430544e-001 with an estimated absolute error of
- 3.11e-015.
-
- 2D modification of an example by L.N.Trefethen (see ref. 7),
- f(x)=exp(-1/x^2).
- '''
- arg_z = (x**2 + y**2)
- # Avoid cases in which "arg_z=0", setting only in those instances
- # "arg_z=eps".
- arg_z = arg_z + (1 - np.abs(np.sign(arg_z))) * 1.e-100
- arg_z = 1. / arg_z
- return np.exp(-arg_z)
-
- def __call__(self, x, y, i=0):
- s = self
- test_function = [s.franke, s.half_sphere, s.poly_degree20, s.exp_fun1,
- s.exp_fun100, s.cos30, s.constant, s.exp_xy, s.runge,
- s.abs_cubed, s.gauss, s.exp_inv]
- return test_function[i](x, y)
-
-
-example_functions = _ExampleFunctions()
-
-
-def _find_m(n):
- ix = np.r_[1:(n + 1) * (n + 2):2]
- if np.mod(n, 2) == 0:
- n2 = n // 2
- offset = np.array([[0, 1] * n2 + [0, ]] * (n2 + 1))
- ix = ix - offset.ravel(order='F')
- return ix
-
-
-def padua_points(n, domain=(-1, 1, -1, 1)):
- ''' Return Padua points
-
- Parameters
- ----------
- n : scalar integer
- interpolation degree
- domain : vector [a,b,c,d]
- defining the rectangle [a,b] x [c,d]. (default domain = (-1,1,-1,1))
-
- Returns
- -------
- pad : array of shape (2 x (n+1)*(n+2)/2) such that
- (pad[0,:], pad[1,: ]) defines the Padua points in the domain
- rectangle [a,b] x [c,d].
- or
- X1,Y1,X2,Y2 : arrays
- Two subgrids X1,Y1 and X2,Y2 defining the Padua points
- -------------------------------------------------------------------------------
- '''
- a, b, c, d = domain
- t0 = [np.pi] if n == 0 else np.linspace(0, np.pi, n + 1)
- t1 = np.linspace(0, np.pi, n + 2)
- zn = map_to_interval(np.cos(t0), a, b)
- zn1 = map_to_interval(np.cos(t1), c, d)
-
- Pad1, Pad2 = np.meshgrid(zn, zn1)
- ix = _find_m(n)
- return np.vstack((Pad1.ravel(order='F')[ix],
- Pad2.ravel(order='F')[ix]))
-
-
-def error_estimate(C0f):
- ''' Return interpolation error estimate from Padua coefficients
- '''
- n = C0f.shape[1]
- C0f2 = np.fliplr(C0f)
- errest = sum(np.abs(np.diag(C0f2)))
- if (n >= 1):
- errest = errest + sum(np.abs(np.diag(C0f2, -1)))
- if (n >= 2):
- errest = errest + sum(np.abs(np.diag(C0f2, -2)))
- return 2 * errest
-
-
-def padua_fit(Pad, fun, *args):
- '''
- Computes the Chebyshevs coefficients
-
- so that f(x, y) can be approximated by:
-
- f(x, y) = sum cjk*Tjk(x, y)
-
- Parameters
- ----------
- Pad : array-like
- Padua points, as computed with padua_points function.
- fun : function to be interpolated in the form
- fun(x, y, *args), where *args are optional arguments for fun.
-
- Returns
- -------
- coefficents: coefficient matrix
- abs_err : interpolation error estimate
-
- Example
- ------
- >>> import numpy as np
- >>> import wafo.padua as wp
- >>> domain = [0, 1, 0, 1]
- >>> a, b, c, d = domain
- >>> points = wp.padua_points(21, domain)
- >>> C0f, abs_error = wp.padua_fit(points, wp.example_functions.franke)
- >>> x1 = np.linspace(a, b, 100)
- >>> x2 = np.linspace(c, d, 101)
- >>> val = wp.padua_val(x1, x2, C0f, domain, use_meshgrid=True)
- >>> X1, X2 = np.meshgrid(x1, x2)
- >>> true_val = wp.example_functions.franke(X1, X2)
-
- >>> np.allclose(val, true_val, atol=10*abs_error)
- True
- >>> np.allclose(np.abs(val-true_val).max(), 0.0073174614275738296)
- True
- >>> np.allclose(abs_error, 0.0022486904061664046)
- True
-
- import matplotlib.pyplot as plt
- plt.contour(x1, x2, val)
-
- '''
-
- N = np.shape(Pad)[1]
- # recover the degree n from N = (n+1)(n+2)/2
- n = int(round(-3 + np.sqrt(1 + 8 * N)) / 2)
- C0f = fun(Pad[0], Pad[1], *args)
- if (n > 0):
- ix = _find_m(n)
- GfT = np.zeros((n + 2) * (n + 1))
- GfT[ix] = C0f * 2 / (n * (n + 1))
- GfT = GfT.reshape(n + 1, n + 2)
- GfT = GfT.T
- GfT[0] = GfT[0] / 2
- GfT[n + 1] = GfT[n + 1] / 2
- GfT[:, 0] = GfT[:, 0] / 2
- GfT[:, n] = GfT[:, n] / 2
- Gf = GfT.T
- # compute the interpolation coefficient matrix C0f by FFT
- Gfhat = np.real(fft(Gf, 2 * n, axis=0))
- Gfhathat = np.real(fft(Gfhat[:n + 1, :], 2 * (n + 1), axis=1))
- C0f = 2 * Gfhathat[:, 0:n + 1]
- C0f[0] = C0f[0, :] / np.sqrt(2)
- C0f[:, 0] = C0f[:, 0] / np.sqrt(2)
- C0f = np.fliplr(np.triu(np.fliplr(C0f)))
- C0f[n] = C0f[n] / 2
-
- return C0f, error_estimate(C0f)
-
-
-def paduavals2coefs(f):
- m = len(f)
- n = int(round(-1.5 + np.sqrt(.25 + 2 * m)))
- x = padua_points(n)
- idx = _find_m(n)
- w = 0 * x[0] + 1. / (n * (n + 1))
- idx1 = np.all(np.abs(x) == 1, axis=0)
- w[idx1] = .5 * w[idx1]
- idx2 = np.all(np.abs(x) != 1, axis=0)
- w[idx2] = 2 * w[idx2]
-
- G = np.zeros(idx.max() + 1)
- G[idx] = 4 * w * f
-
- use_dct = 100 < n
- if use_dct:
- C = np.rot90(dct(dct(G.T).T)) # , axis=1)
- else:
- t1 = np.r_[0:n + 1].reshape(-1, 1)
- Tn1 = np.cos(t1 * t1.T * np.pi / n)
- t2 = np.r_[0:n + 2].reshape(-1, 1)
- Tn2 = np.cos(t2 * t2.T * np.pi / (n + 1))
- C = np.dot(Tn2, np.dot(G, Tn1))
-
- C[0] = .5 * C[0]
- C[:, 1] = .5 * C[:, 1]
- C[0, -1] = .5 * C[0, -1]
- del C[-1]
-
- # Take upper-left triangular part:
- return np.fliplr(np.triu(np.fliplr(C)))
- # C = triu(C(:,end:-1:1));
- # C = C(:,end:-1:1);
-
-
-# TODO: padua_fit2 does not work correctly yet.
-def padua_fit2(Pad, fun, *args):
- # N = np.shape(Pad)[1]
- # recover the degree n from N = (n+1)(n+2)/2
- # _n = int(round(-3 + np.sqrt(1 + 8 * N)) / 2)
- C0f = fun(Pad[0], Pad[1], *args)
- return paduavals2coefs(C0f)
-
-
-def _compute_moments(n):
- k = np.r_[0:n:2]
- mom = 2 * np.sqrt(2) / (1 - k ** 2)
- mom[0] = 2
- return mom
-
-
-def padua_cubature(coefficients, domain=(-1, 1, -1, 1)):
- '''
- Compute the integral through the coefficient matrix.
- '''
- n = coefficients.shape[1]
- mom = _compute_moments(n)
- M1, M2 = np.meshgrid(mom, mom)
- M = M1 * M2
- C0fM = coefficients[0:n:2, 0:n:2] * M
- a, b, c, d = domain
- integral = (b - a) * (d - c) * C0fM.sum() / 4
- return integral
-
-
-def padua_val(X, Y, coefficients, domain=(-1, 1, -1, 1), use_meshgrid=False):
- '''
- Evaluate polynomial in padua form at X, Y.
-
- Evaluate the interpolation polynomial defined through its coefficient
- matrix coefficients at the target points X(:,1),X(:,2) or at the
- meshgrid(X1,X2)
-
- Parameters
- ----------
- X, Y: array-like
- evaluation points.
- coefficients : array-like
- coefficient matrix
- domain : a vector [a,b,c,d]
- defining the rectangle [a,b] x [c,d]
- use_meshgrid: bool
- If True interpolate at the points meshgrid(X, Y)
-
- Returns
- -------
- fxy : array-like
- evaluation of the interpolation polynomial at the target points
- '''
- def transform(tn, x, a, b):
- xn = map_from_interval(x, a, b).clip(min=-1, max=1).reshape(1, -1)
- tx = np.cos(tn * np.arccos(xn)) * np.sqrt(2)
- tx[0] = 1
- return tx
-
- X, Y = np.atleast_1d(X, Y)
- original_shape = X.shape
- a, b, c, d = domain
- n = np.shape(coefficients)[1]
-
- tn = np.r_[0:n][:, None]
- tx1 = transform(tn, X.ravel(), a, b)
- tx2 = transform(tn, Y.ravel(), c, d)
-
- if use_meshgrid: # eval on meshgrid points
- return np.dot(tx1.T, np.dot(coefficients, tx2)).T
- # scattered points
- val = np.sum(np.dot(tx1.T, coefficients) * tx2.T, axis=1)
- return val.reshape(original_shape)
-
-
-if __name__ == '__main__':
- from wafo.testing import test_docstrings
- test_docstrings(__file__)
diff --git a/wafo/setup.py b/wafo/setup.py
deleted file mode 100644
index a535100..0000000
--- a/wafo/setup.py
+++ /dev/null
@@ -1,18 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-Created on Sun Oct 25 14:55:34 2015
-
-@author: dave
-"""
-
-
-def configuration(parent_package='', top_path=None):
- from numpy.distutils.misc_util import Configuration
- config = Configuration('wafo', parent_package, top_path)
- config.add_subpackage('source')
- config.make_config_py()
- return config
-
-if __name__ == "__main__":
- from numpy.distutils.core import setup
- setup(**configuration(top_path='').todict())
diff --git a/wafo/wavemodels.py b/wafo/wavemodels.py
deleted file mode 100644
index 36a77c1..0000000
--- a/wafo/wavemodels.py
+++ /dev/null
@@ -1,268 +0,0 @@
-'''
-Created on 13. mar. 2018
-
-@author: pab
-'''
-import numpy as np
-from numpy import pi, sqrt
-import wafo.transform.estimation as te
-import wafo.transform as wt
-from wafo.containers import PlotData
-from wafo.kdetools.kernels import qlevels
-from wafo.misc import tranproc
-import warnings
-
-
-
-def _set_default_t_h_g(t, h, g, m0, m2):
- if g is None:
- y = np.linspace(-5, 5)
- x = sqrt(m0) * y + 0
- g = wt.TrData(y, x)
- if t is None:
- tt1 = 2 * pi * sqrt(m0 / m2)
- t = np.linspace(0, 1.7 * tt1, 51)
- if h is None:
- px = g.gauss2dat([0, 4.])
- px = abs(px[1] - px[0])
- h = np.linspace(0, 1.3 * px, 41)
- return h, t, g
-
-def lh83pdf(t=None, h=None, mom=None, g=None):
- """
- LH83PDF Longuet-Higgins (1983) approximation of the density (Tc,Ac)
- in a stationary Gaussian transform process X(t) where
- Y(t) = g(X(t)) (Y zero-mean Gaussian, X non-Gaussian).
-
- CALL: f = lh83pdf(t,h,[m0,m1,m2],g);
-
- f = density of wave characteristics of half-wavelength
- in a stationary Gaussian transformed process X(t),
- where Y(t) = g(X(t)) (Y zero-mean Gaussian)
- t,h = vectors of periods and amplitudes, respectively.
- default depending on the spectral moments
- m0,m1,m2 = the 0'th,1'st and 2'nd moment of the spectral density
- with angular frequency.
- g = space transformation, Y(t)=g(X(t)), default: g is identity
- transformation, i.e. X(t) = Y(t) is Gaussian,
- The transformation, g, can be estimated using lc2tr
- or dat2tr or given apriori by ochi.
-
- Example
- -------
- >>> import wafo.spectrum.models as sm
- >>> Sj = sm.Jonswap()
- >>> w = np.linspace(0,4,256)
- >>> S = Sj.tospecdata(w) #Make spectrum object from numerical values
- >>> S = sm.SpecData1D(Sj(w),w) # Alternatively do it manually
- >>> mom, mom_txt = S.moment(nr=2, even=False)
- >>> f = lh83pdf(mom=mom)
- >>> f.plot()
-
- See also
- --------
- cav76pdf, lc2tr, dat2tr
-
- References
- ----------
- Longuet-Higgins, M.S. (1983)
- "On the joint distribution wave periods and amplitudes in a
- random wave field", Proc. R. Soc. A389, pp 24--258
-
- Longuet-Higgins, M.S. (1975)
- "On the joint distribution wave periods and amplitudes of sea waves",
- J. geophys. Res. 80, pp 2688--2694
- """
-
- # tested on: matlab 5.3
- # History:
- # Revised pab 01.04.2001
- # - Added example
- # - Better automatic scaling for h,t
- # revised by IR 18.06.2000, fixing transformation and transposing t and h to fit simpson req.
- # revised by pab 28.09.1999
- # made more efficient calculation of f
- # by Igor Rychlik
-
- m0, m1, m2 = mom
- h, t, g = _set_default_t_h_g(t, h, g, m0, m2)
-
- L0 = m0
- L1 = m1 / (2 * pi)
- L2 = m2 / (2 * pi)**2
- eps2 = sqrt((L2 * L0) / (L1**2) - 1)
-
- if np.any(~np.isreal(eps2)):
- raise ValueError('input moments are not correct')
-
- const = 4 / sqrt(pi) / eps2 / (1 + 1 / sqrt(1 + eps2**2))
-
- a = len(h)
- b = len(t)
- der = np.ones((a, 1))
-
- h_lh = g.dat2gauss(h.ravel(), der.ravel())
-
- der = abs(h_lh[1]) # abs(h_lh[:, 1])
- h_lh = h_lh[0]
-
- # Normalization + transformation of t and h ???????
- # Without any transformation
-
- t_lh = t / (L0 / L1)
- #h_lh = h_lh/sqrt(2*L0)
- h_lh = h_lh / sqrt(2)
- t_lh = 2 * t_lh
-
- # Computation of the distribution
- T, H = np.meshgrid(t_lh[1:b], h_lh)
- f_th = np.zeros((a, b))
- tmp = const * der[:, None] * (H / T)**2 * np.exp(-H**2. *
- (1 + ((1 - 1. / T) / eps2)**2)) / ((L0 / L1) * sqrt(2) / 2)
- f_th[:, 1:b] = tmp
-
- f = PlotData(f_th, (t, h),
- xlab='Tc', ylab='Ac',
- title='Joint density of (Tc,Ac) - Longuet-Higgins (1983)',
- plot_kwds=dict(plotflag=1))
-
- return _add_contour_levels(f)
-
-
-def cav76pdf(t=None, h=None, mom=None, g=None):
- """
- CAV76PDF Cavanie et al. (1976) approximation of the density (Tc,Ac)
- in a stationary Gaussian transform process X(t) where
- Y(t) = g(X(t)) (Y zero-mean Gaussian, X non-Gaussian).
-
- CALL: f = cav76pdf(t,h,[m0,m2,m4],g);
-
- f = density of wave characteristics of half-wavelength
- in a stationary Gaussian transformed process X(t),
- where Y(t) = g(X(t)) (Y zero-mean Gaussian)
- t,h = vectors of periods and amplitudes, respectively.
- default depending on the spectral moments
- m0,m2,m4 = the 0'th, 2'nd and 4'th moment of the spectral density
- with angular frequency.
- g = space transformation, Y(t)=g(X(t)), default: g is identity
- transformation, i.e. X(t) = Y(t) is Gaussian,
- The transformation, g, can be estimated using lc2tr
- or dat2tr or given a priori by ochi.
- [] = default values are used.
-
- Example
- -------
- >>> import wafo.spectrum.models as sm
- >>> Sj = sm.Jonswap()
- >>> w = np.linspace(0,4,256)
- >>> S = Sj.tospecdata(w) #Make spectrum object from numerical values
- >>> S = sm.SpecData1D(Sj(w),w) # Alternatively do it manually
- >>> mom, mom_txt = S.moment(nr=4, even=True)
- >>> f = cav76pdf(mom=mom)
- >>> f.plot()
-
- See also
- --------
- lh83pdf, lc2tr, dat2tr
-
- References
- ----------
- Cavanie, A., Arhan, M. and Ezraty, R. (1976)
- "A statistical relationship between individual heights and periods of
- storm waves".
- In Proceedings Conference on Behaviour of Offshore Structures,
- Trondheim, pp. 354--360
- Norwegian Institute of Technology, Trondheim, Norway
-
- Lindgren, G. and Rychlik, I. (1982)
- Wave Characteristics Distributions for Gaussian Waves --
- Wave-lenght, Amplitude and Steepness, Ocean Engng vol 9, pp. 411-432.
- """
- # tested on: matlab 5.3 NB! note
- # History:
- # revised pab 04.11.2000
- # - fixed xlabels i.e. f.labx={'Tc','Ac'}
- # revised by IR 4 X 2000. fixed transform and normalisation
- # using Lindgren & Rychlik (1982) paper.
- # At the end of the function there is a text with derivation of the density.
- #
- # revised by jr 21.02.2000
- # - Introduced cell array for f.x for use with pdfplot
- # by pab 28.09.1999
-
- m0, m2, m4 = mom
- h, t, g = _set_default_t_h_g(t, h, g, m0, m2)
-
- eps4 = 1.0 - m2**2 / (m0 * m4)
- alfa = m2 / sqrt(m0 * m4)
- if np.any(~np.isreal(eps4)):
- raise ValueError('input moments are not correct')
-
- a = len(h)
- b = len(t)
- der = np.ones((a, 1))
-
- h_lh = g.dat2gauss(h.ravel(), der.ravel())
- der = abs(h_lh[1])
- h_lh = h_lh[0]
-
- # Normalization + transformation of t and h
-
- pos = 2 / (1 + alfa) # inverse of a fraction of positive maxima
- cons = 2 * pi**4 * pos / sqrt(2 * pi) / m4 / sqrt((1 - alfa**2))
- # Tm=2*pi*sqrt(m0/m2)/alpha; #mean period between positive maxima
-
- t_lh = t
- h_lh = sqrt(m0) * h_lh
-
- # Computation of the distribution
- T, H = np.meshgrid(t_lh[1:b], h_lh)
- f_th = np.zeros((a, b))
-
- f_th[:, 1:b] = cons * der[:, None] * (H**2 / (T**5)) * np.exp(-0.5 * (
- H / T**2)**2. * ((T**2 - pi**2 * m2 / m4)**2 / (m0 * (1 - alfa**2)) + pi**4 / m4))
- f = PlotData(f_th, (t, h),
- xlab='Tc', ylab='Ac',
- title='Joint density of (Tc,Ac) - Cavanie et al. (1976)',
- plot_kwds=dict(plotflag=1))
- return _add_contour_levels(f)
-
-def _add_contour_levels(f):
- p_levels = np.r_[10:90:20, 95, 99, 99.9]
- try:
- c_levels = qlevels(f.data, p=p_levels)
- f.clevels = c_levels
- f.plevels = p_levels
- except ValueError as e:
- msg = "Could not calculate contour levels!. ({})".format(str(e))
- warnings.warn(msg)
- return f
-
- # Let U,Z be the height and second derivative (curvature) at a local maximum in a Gaussian proces
- # with spectral moments m0,m2,m4. The conditional density ($U>0$) has the following form
- #$$
- # f(z,u)=c \frac{1}{\sqrt{2\pi}}\frac{1}{\sqrt{m0(1-\alpha^2)}}\exp(-0.5\left(\frac{u-z(m2/m4)}
- # {\sqrt{m0(1-\alpha^2)}}\right)^2)\frac{|z|}{m4}\exp(-0.5z^2/m4), \quad z<0,
- #$$
- # where $c=2/(1+\alpha)$, $\alpha=m2/\sqrt{m0\cdot m4}$.
- #
- # The cavanie approximation is based on the model $X(t)=U \cos(\pi t/T)$, consequently
- # we have $U=H$ and by twice differentiation $Z=-U(\pi^2/T)^2\cos(0)$. The variable change has Jacobian
- # $2\pi^2 H/T^3$ giving the final formula for the density of $T,H$
- #$$
- # f(t,h)=c \frac{2\pi^4}{\sqrt{2\pi}}\frac{1}{m4\sqrt{m0(1-\alpha^2)}}\frac{h^2}{t^5}
- # \exp(-0.5\frac{h^2}{t^4}\left(\left(\frac{t^2-\pi^2(m2/m4)}
- # {\sqrt{m0(1-\alpha^2)}}\right)^2+\frac{\pi^4}{m4}\right)).
- #$$
- #
- #
-
-
-def test_docstrings():
- import doctest
- print('Testing docstrings in %s' % __file__)
- doctest.testmod(optionflags=doctest.NORMALIZE_WHITESPACE)
-
-
-if __name__ == '__main__':
- test_docstrings()
diff --git a/wafo/win32_utils.py b/wafo/win32_utils.py
deleted file mode 100644
index ace64e4..0000000
--- a/wafo/win32_utils.py
+++ /dev/null
@@ -1,258 +0,0 @@
-from __future__ import division
-from numpy import round
-from threading import Thread
-from time import sleep
-from win32gui import (InitCommonControls, CallWindowProc, CreateWindowEx,
- CreateWindow, SetWindowLong, SendMessage, ShowWindow,
- PumpWaitingMessages, PostQuitMessage, DestroyWindow,
- MessageBox, EnumWindows, GetClassName)
-from win32api import GetModuleHandle, GetSystemMetrics # @UnresolvedImport
-from win32api import SetWindowLong as api_SetWindowLong # @UnresolvedImport
-from commctrl import (TOOLTIPS_CLASS, TTM_GETDELAYTIME, TTM_SETDELAYTIME,
- TTDT_INITIAL, TTDT_AUTOPOP)
-import win32con
-
-WM_USER = win32con.WM_USER
-PBM_SETRANGE = (WM_USER + 1)
-PBM_SETPOS = (WM_USER + 2)
-PBM_DELTAPOS = (WM_USER + 3)
-PBM_SETSTEP = (WM_USER + 4)
-PBM_STEPIT = (WM_USER + 5)
-PBM_SETRANGE32 = (WM_USER + 6)
-PBM_GETRANGE = (WM_USER + 7)
-PBM_GETPOS = (WM_USER + 8)
-PBM_SETBARCOLOR = (WM_USER + 9)
-PBM_SETMARQUEE = (WM_USER + 10)
-PBM_GETSTEP = (WM_USER + 13)
-PBM_GETBKCOLOR = (WM_USER + 14)
-PBM_GETBARCOLOR = (WM_USER + 15)
-PBM_SETSTATE = (WM_USER + 16)
-PBM_GETSTATE = (WM_USER + 17)
-PBS_SMOOTH = 0x01
-PBS_VERTICAL = 0x04
-PBS_MARQUEE = 0x08
-PBS_SMOOTHREVERSE = 0x10
-PBST_NORMAL = 1
-PBST_ERROR = 2
-PBST_PAUSED = 3
-WC_DIALOG = 32770
-WM_SETTEXT = win32con.WM_SETTEXT
-
-
-def MAKELPARAM(a, b):
- return (a & 0xffff) | ((b & 0xffff) << 16)
-
-
-def _get_tooltip_handles(hwnd, resultList):
- '''
- Adds a window handle to resultList if its class-name is 'tooltips_class32',
- i.e. the window is a tooltip.
- '''
- if GetClassName(hwnd) == TOOLTIPS_CLASS:
- resultList.append(hwnd)
-
-
-def set_max_pop_delay_on_tooltip(tooltip):
- '''
- Sets maximum auto-pop delay (delay before hiding) on an instance of
- wx.ToolTip.
- NOTE: The tooltip's SetDelay method is used just to identify the correct
- tooltip.
- '''
- test_value = 12345
- # Set initial delay just to identify tooltip.
- tooltip.SetDelay(test_value)
- handles = []
- EnumWindows(_get_tooltip_handles, handles)
- for hwnd in handles:
- if SendMessage(hwnd, TTM_GETDELAYTIME, TTDT_INITIAL) == test_value:
- SendMessage(hwnd, TTM_SETDELAYTIME, TTDT_AUTOPOP, 32767)
- tooltip.SetDelay(500) # Restore default value
-
-
-class Waitbar(Thread):
-
- def __init__(self, title='Waitbar', can_abort=True, max_val=100):
- Thread.__init__(self) # Initialize thread
- self.title = title
- self.can_abort = can_abort
- self.max_val = max_val
- InitCommonControls()
- self.hinst = GetModuleHandle(None)
- self.started = False
- self.position = 0
- self.do_update = False
- self.start() # Run the thread
- while not self.started:
- sleep(0.1) # Wait until the dialog is ready
-
- def DlgProc(self, hwnd, uMsg, wParam, lParam):
- if uMsg == win32con.WM_DESTROY:
- api_SetWindowLong(self.dialog,
- win32con.GWL_WNDPROC,
- self.oldWndProc)
- if uMsg == win32con.WM_CLOSE:
- self.started = False
- if uMsg == win32con.WM_COMMAND and self.can_abort:
- self.started = False
- return CallWindowProc(self.oldWndProc, hwnd, uMsg, wParam, lParam)
-
- def BuildWindow(self):
- width = 400
- height = 100
- self.dialog = CreateWindowEx(
- win32con.WS_EX_TOPMOST,
- WC_DIALOG,
- self.title + ' (0%)',
- win32con.WS_VISIBLE | win32con.WS_OVERLAPPEDWINDOW,
- int(round(
- GetSystemMetrics(win32con.SM_CXSCREEN) * .5 - width * .5)),
- int(round(
- GetSystemMetrics(win32con.SM_CYSCREEN) * .5 - height * .5)),
- width,
- height,
- 0,
- 0,
- self.hinst,
- None)
- self.progbar = CreateWindow(
- # win32con.WS_EX_DLGMODALFRAME,
- 'msctls_progress32',
- '',
- win32con.WS_VISIBLE | win32con.WS_CHILD,
- 10,
- 10,
- width - 30,
- 20,
- self.dialog,
- 0,
- 0,
- None)
- if self.can_abort:
- self.button = CreateWindow(
- # win32con.WS_EX_DLGMODALFRAME,
- 'BUTTON',
- 'Cancel',
- win32con.WS_VISIBLE | win32con.WS_CHILD | win32con.BS_PUSHBUTTON, # @IgnorePep8
- int(width / 2.75),
- 40,
- 100,
- 20,
- self.dialog,
- 0,
- 0,
- None)
- self.oldWndProc = SetWindowLong(
- self.dialog,
- win32con.GWL_WNDPROC,
- self.DlgProc)
- SendMessage(self.progbar, PBM_SETRANGE, 0, MAKELPARAM(0, self.max_val))
-# win32gui.SendMessage(self.progbar, PBM_SETSTEP, 0, 10)
-# win32gui.SendMessage(self.progbar, PBM_SETMARQUEE, 0, 0)
- ShowWindow(self.progbar, win32con.SW_SHOW)
-
- def run(self):
- self.BuildWindow()
- self.started = True
- while self.started:
- PumpWaitingMessages()
- if self.do_update:
- SendMessage(self.progbar, PBM_SETPOS,
- int(self.position % self.max_val), 0)
- percentage = int(round(100.0 * self.position / self.max_val))
- SendMessage(self.dialog, WM_SETTEXT, 0,
- self.title + ' (%d%%)' % percentage)
- # SendMessage(self.progbar, PBM_STEPIT, 0, 0)
- self.do_update = False
- sleep(0.1)
- PostQuitMessage(0)
- DestroyWindow(self.dialog)
-
- def update(self, pos):
- if self.started:
- if not self.do_update:
- self.position = pos
- self.do_update = True
- return True
- return False
-
- def close(self):
- self.started = False
-
-
-# class Waitbar2(Dialog, Thread):
-# def __init__(self, title='Waitbar'):
-# template = [[title, (0, 0, 215, 36),
-# (win32con.DS_MODALFRAME | win32con.WS_POPUP |
-# win32con.WS_VISIBLE | win32con.WS_CAPTION |
-# win32con.WS_SYSMENU | win32con.DS_SETFONT |
-# win32con.WS_GROUP | win32con.WS_EX_TOPMOST),
-# | win32con.DS_SYSMODAL),
-# None, (8, "MS Sans Serif")], ]
-# Dialog.__init__(self, id=template)
-# Thread.__init__(self) # Initialize thread
-# self.started = False
-# self.start() # Run the thread
-# while not self.started:
-# sleep(0.1) # Wait until the dialog is ready
-#
-# def OnInitDialog(self):
-# rc = Dialog.OnInitDialog(self)
-# self.pbar = CreateProgressCtrl()
-# self.pbar.CreateWindow (win32con.WS_CHILD | win32con.WS_VISIBLE,
-# (10, 10, 310, 24), self, 1001)
-# self.started = True
-# return rc
-#
-# def run(self):
-# self.DoModal()
-#
-# def update(self, pos):
-# self.pbar.SetPos(int(pos))
-#
-# def close(self):
-# self.OnCancel()
-
-
-class WarnDlg(Thread):
-
- def __init__(self, message='', title='Warning!'):
- Thread.__init__(self) # Initialize thread
- self.title = title
- self.message = message
- self.start() # Run the thread
-
- def run(self):
- # MessageBox(self.message, self.title, win32con.MB_ICONWARNING)
- MessageBox(0, self.message, self.title,
- win32con.MB_ICONWARNING | win32con.MB_SYSTEMMODAL)
-
-
-class ErrorDlg(Thread):
-
- def __init__(self, message='', title='Error!', blocking=False):
- Thread.__init__(self) # Initialize thread
- self.title = title
- self.message = message
- if blocking:
- self.run() # Run without threading
- else:
- self.start() # Run in thread
-
- def run(self):
- # MessageBox(self.message, self.title, win32con.MB_ICONERROR)
- MessageBox(0, self.message, self.title,
- win32con.MB_ICONERROR | win32con.MB_SYSTEMMODAL)
-
-
-if __name__ == '__main__':
- WarnDlg('This is an example of a warning', 'Warning!')
- ErrorDlg('This is an example of an error message')
- wb = Waitbar('Waitbar example')
-# wb2 = Waitbar2('Waitbar example')
- for i in range(20):
- print(wb.update(i * 5))
-# wb2.update(i)
- sleep(0.1)
- wb.close()
-# wb2.close()