python - Pythonでの関数の適応プロット

Question

簡単な質問：[0,1]のある時点で鋭いピークを持つと思われる関数f（t）があります。自然なアイデアは、この関数の適応サンプリングを使用して、優れた「適応」プロットを取得することです。Python + matplotlib + numpy +でそれを高速に行うにはどうすればよいですか？[0,1]の任意のtについてf（t）を計算できます。

Mathematicaにはこのオプションがあるようですが、Pythonにはオプションがありますか？

score 7 · Accepted Answer

私が見つけたものを見てください：1D関数の適応サンプリング、scipy-central.orgからのリンク。

コードは次のとおりです。

# License: Creative Commons Zero (almost public domain) http://scpyce.org/cc0

import numpy as np

def sample_function(func, points, tol=0.05, min_points=16, max_level=16,
                    sample_transform=None):
    """
    Sample a 1D function to given tolerance by adaptive subdivision.

    The result of sampling is a set of points that, if plotted,
    produces a smooth curve with also sharp features of the function
    resolved.

    Parameters
    ----------
    func : callable
        Function func(x) of a single argument. It is assumed to be vectorized.
    points : array-like, 1D
        Initial points to sample, sorted in ascending order.
        These will determine also the bounds of sampling.
    tol : float, optional
        Tolerance to sample to. The condition is roughly that the total
        length of the curve on the (x, y) plane is computed up to this
        tolerance.
    min_point : int, optional
        Minimum number of points to sample.
    max_level : int, optional
        Maximum subdivision depth.
    sample_transform : callable, optional
        Function w = g(x, y). The x-samples are generated so that w
        is sampled.

    Returns
    -------
    x : ndarray
        X-coordinates
    y : ndarray
        Corresponding values of func(x)

    Notes
    -----
    This routine is useful in computing functions that are expensive
    to compute, and have sharp features --- it makes more sense to
    adaptively dedicate more sampling points for the sharp features
    than the smooth parts.

    Examples
    --------
    >>> def func(x):
    ...     '''Function with a sharp peak on a smooth background'''
    ...     a = 0.001
    ...     return x + a**2/(a**2 + x**2)
    ...
    >>> x, y = sample_function(func, [-1, 1], tol=1e-3)

    >>> import matplotlib.pyplot as plt
    >>> xx = np.linspace(-1, 1, 12000)
    >>> plt.plot(xx, func(xx), '-', x, y[0], '.')
    >>> plt.show()

    """
    return _sample_function(func, points, values=None, mask=None, depth=0,
                            tol=tol, min_points=min_points, max_level=max_level,
                            sample_transform=sample_transform)

def _sample_function(func, points, values=None, mask=None, tol=0.05,
                     depth=0, min_points=16, max_level=16,
                     sample_transform=None):
    points = np.asarray(points)

    if values is None:
        values = np.atleast_2d(func(points))

    if mask is None:
        mask = Ellipsis

    if depth > max_level:
        # recursion limit
        return points, values

    x_a = points[...,:-1][...,mask]
    x_b = points[...,1:][...,mask]

    x_c = .5*(x_a + x_b)
    y_c = np.atleast_2d(func(x_c))

    x_2 = np.r_[points, x_c]
    y_2 = np.r_['-1', values, y_c]
    j = np.argsort(x_2)

    x_2 = x_2[...,j]
    y_2 = y_2[...,j]

    # -- Determine the intervals at which refinement is necessary

    if len(x_2) < min_points:
        mask = np.ones([len(x_2)-1], dtype=bool)
    else:
        # represent the data as a path in N dimensions (scaled to unit box)
        if sample_transform is not None:
            y_2_val = sample_transform(x_2, y_2)
        else:
            y_2_val = y_2

        p = np.r_['0',
                  x_2[None,:],
                  y_2_val.real.reshape(-1, y_2_val.shape[-1]),
                  y_2_val.imag.reshape(-1, y_2_val.shape[-1])
                  ]

        sz = (p.shape[0]-1)//2

        xscale = x_2.ptp(axis=-1)
        yscale = abs(y_2_val.ptp(axis=-1)).ravel()

        p[0] /= xscale
        p[1:sz+1] /= yscale[:,None]
        p[sz+1:]  /= yscale[:,None]

        # compute the length of each line segment in the path
        dp = np.diff(p, axis=-1)
        s = np.sqrt((dp**2).sum(axis=0))
        s_tot = s.sum()

        # compute the angle between consecutive line segments
        dp /= s
        dcos = np.arccos(np.clip((dp[:,1:] * dp[:,:-1]).sum(axis=0), -1, 1))

        # determine where to subdivide: the condition is roughly that
        # the total length of the path (in the scaled data) is computed
        # to accuracy `tol`
        dp_piece = dcos * .5*(s[1:] + s[:-1])
        mask = (dp_piece > tol * s_tot)

        mask = np.r_[mask, False]
        mask[1:] |= mask[:-1].copy()


    # -- Refine, if necessary

    if mask.any():
        return _sample_function(func, x_2, y_2, mask, tol=tol, depth=depth+1,
                                min_points=min_points, max_level=max_level,
                                sample_transform=sample_transform)
    else:
        return x_2, y_2

score 2 · Accepted Answer

https://github.com/python-adaptive/adaptiveは、これとはるかに多くのことを行う試みのようです。

アダプティブ

数学関数の適応並列サンプリングのためのツール。

adaptiveは、適応型並列関数の評価を簡単にするために設計されたオープンソースのPythonライブラリです。adaptive関数にその境界を指定するだけで、パラメーター空間の「最良の」ポイントで評価されます。ほんの数行のコードで、コンピューティングクラスター上の関数を評価し、データが返されるときにデータをライブプロットし、適応サンプリングアルゴリズムを微調整できます。

このプロジェクトは、この質問に対する別の回答（または少なくとも関連プロジェクト）のコードにも触発されました。

クレジット

..。

〜adaptive.Learner2Dのインスピレーションとして役立ったAdaptiveTriSamplingスクリプト（SciPy Centralがダウンしたため、オンラインでは利用できなくなりました）のPauliVirtanen。

score -3 · Accepted Answer

プロットの目的では、適応サンプリングはあまり必要ありません。画面解像度以上でサンプリングしてみませんか？

POINTS=1920

from pylab import *
x = arange(0,1,1.0/POINTS)
y = sin(3.14*x)
plot(x,y)
axes().set_aspect('equal') ## optional aspect-ratio control
show()

任意のサンプリング密度が必要な場合、または関数がベクトル化されたアプローチと互換性がない場合は、x、y配列をポイントごとに作成できます。中間値は、関数によって線形補間されplot()ます。

POINTS=1980
from pylab import *
ax,ay = [],[]
for x in linspace(0.0,POINTS,POINTS):
   if randint(0,50)==0 or x==0 or abs(x-POINTS)<1e-9:
      y = math.sin(4*2*math.pi*x/POINTS)
      ax.append(x), ay.append(y)
# print(ax,ay)
plot(ax,ay)

python - Pythonでの関数の適応プロット

3 に答える 3

Related

Reference