3

一般化された状態空間のルートを見つける必要があります。つまり、私は次元の離散グリッドを持っていますが、それがgrid=AxBx(...)xXいくつの次元を持っているかは事前にはわかりません（解決策はすべてに適用できるはずですgrid.size）。

二分法を使用して内部f(z) = 0のすべての状態の roots() を見つけたいです。含むと言って、私は知っています。それから私はする必要がありますzgridremainderf(z)f'(z) < 0

  • > 0のz 場合は増加remainder
  • < 0のz場合は減少remainder

historyブログ、形状の行列にはグリッド内のすべての点(grid.shape, T)の以前の値の履歴が含まれており、増加する必要があるとします( > 0 から)。次に、「より大きいもののうち最小のもの」を選択する必要があります。疑似コードでは、次のようになります。zzremainderzAlternativehistory[z, :]z

zAlternative =  hist[z,:][hist[z,:] > z].min()

私は以前にこれを尋ねました。私が与えられた解決策は

b = sort(history[..., :-1], axis=-1)
mask = b > history[..., -1:]
index = argmax(mask, axis=-1)
indices = tuple([arange(j) for j in b.shape[:-1]])
indices = meshgrid(*indices, indexing='ij', sparse=True)
indices.append(index)
indices = tuple(indices)
lowerZ = history[indices]

b = sort(history[..., :-1], axis=-1)
mask = b <= history[..., -1:]
index = argmax(mask, axis=-1)
indices = tuple([arange(j) for j in b.shape[:-1]])
indices = meshgrid(*indices, indexing='ij', sparse=True)
indices.append(index)
indices = tuple(indices)
higherZ = history[indices]

newZ = history[..., -1]
criterion = 0.05
increase = remainder > 0 + criterion
decrease = remainder < 0 - criterion
newZ[increase] = 0.5*(newZ[increase] + higherZ[increase])
newZ[decrease] = 0.5*(newZ[decrease] + lowerZ[decrease])

ただし、このコードは機能しなくなります。それを認めるのは非常に残念ですが、インデックスで起こっている魔法を理解していなかったので、残念ながら助けが必要です.

コードが実際に行うことは、最低と最高をそれぞれ与えることです。つまり、2 つの特定のz値を修正すると、次のようになります。

history[z1] = array([0.3, 0.2, 0.1])
history[z2] = array([0.1, 0.2, 0.3])

higherZ[z1]=0.3とlowerZ[z2] = 0.1、つまり極値を取得します。両方の場合の正しい値は0.2. ここで何がうまくいかないのですか？

必要に応じて、テスト データを生成するために、次のようなものを使用できます。

history = tile(array([0.1, 0.3, 0.2, 0.15, 0.13])[newaxis,newaxis,:], (10, 20, 1))
remainder = -1*ones((10, 20))

2 番目のケースをテストします。

期待される結果

上記の変数を調整してhistory、上向きと下向きの両方のテストケースを提供しました。期待される結果は

lowerZ = 0.1 * ones((10,20))
higherZ = 0.15 * ones((10,20))

zつまり、履歴[z, :]のすべてのポイントについて、次に高い前の値 ( higherZ) と次に小さい前の値 ( lowerZ) です。すべての点zはまったく同じ履歴 ( [0.1, 0.3, 0.2, 0.15, 0.13]) を持っているため、 と の値はすべて同じにlowerZなりhigherZます。もちろん、一般に、それぞれの履歴はz異なるため、2 つのマトリックスにはすべてのグリッド ポイントで異なる値が含まれる可能性があります。

4

2 に答える 2

1

ここに投稿したものを以前の投稿のソリューションと比較したところ、いくつかの違いがあることに気付きました。

より小さいzについて、あなたは言いました

mask = b > history[..., -1:]
index = argmax(mask, axis=-1)

彼らは言った：

mask = b >= a[..., -1:]
index = np.argmax(mask, axis=-1) - 1

より大きなzについて、あなたは言った

mask = b <= history[..., -1:]
index = argmax(mask, axis=-1)

彼らは言った：

mask = b > a[..., -1:]
index = np.argmax(mask, axis=-1)

以前の投稿のソリューションを使用すると、次のようになります。

import numpy as np
history = np.tile(np.array([0.1, 0.3, 0.2, 0.15, 0.13])[np.newaxis,np.newaxis,:], (10, 20, 1))
remainder = -1*np.ones((10, 20))

a = history

# b is a sorted ndarray excluding the most recent observation
# it is sorted along the observation axis
b = np.sort(a[..., :-1], axis=-1)

# mask is a boolean array, comparing the (sorted)
# previous observations to the current observation - [..., -1:]
mask = b > a[..., -1:]

# The next 5 statements build an indexing array.
# True evaluates to one and False evaluates to zero.
# argmax() will return the index of the first True,
# in this case along the last (observations) axis.
# index is an array with the shape of z (2-d for this test data).
# It represents the index of the next greater
# observation for every 'element' of z.
index = np.argmax(mask, axis=-1)

# The next two statements construct arrays of indices
# for every element of z - the first n-1 dimensions of history.
indices = tuple([np.arange(j) for j in b.shape[:-1]])
indices = np.meshgrid(*indices, indexing='ij', sparse=True)

# Adding index to the end of indices (the last dimension of history)
# produces a 'group' of indices that will 'select' a single observation
# for every 'element' of z
indices.append(index)
indices = tuple(indices)
higherZ = b[indices]

mask = b >= a[..., -1:]
# Since b excludes the current observation, we want the
# index just before the next highest observation for lowerZ,
# hence the minus one.
index = np.argmax(mask, axis=-1) - 1
indices = tuple([np.arange(j) for j in b.shape[:-1]])
indices = np.meshgrid(*indices, indexing='ij', sparse=True)
indices.append(index)
indices = tuple(indices)
lowerZ = b[indices]
assert np.all(lowerZ == .1)
assert np.all(higherZ == .15)

うまくいくようです

于 2014-06-10T23:43:36.243 に答える
0

history 現在の観測がhistory[...,-1:]

history これは、 z の各要素の観測値を反復処理しやすくするために、 のストライドを操作することにより、上位および下位の配列を構築します。これは、Numpy を使用した Efficient Overlapping Windows でnumpy.lib.stride_tricks.as_strided見つかった n-dim generalized 関数を使用して達成されます- 最後にそのソースを含めます

history.shape(10,20,x)に対して 200 回の反復を行う単一の python ループがあります。

import numpy as np

history = np.tile(np.array([0.1, 0.3, 0.2, 0.15, 0.13])[np.newaxis,np.newaxis,:], (10, 20, 1))
remainder = -1*np.ones((10, 20))

z_shape = final_shape = history.shape[:-1]
number_of_observations = history.shape[-1]
number_of_elements_in_z = np.product(z_shape)

# manipulate histories to efficiently iterate over
# the observations of each "element" of z
s = sliding_window(history, (1,1,number_of_observations))
# s.shape will be (number_of_elements_in_z, number_of_observations)

# create arrays of the next lower and next higher observation
lowerZ = np.zeros(number_of_elements_in_z)
higherZ = np.zeros(number_of_elements_in_z)
for ndx, observations in enumerate(s):
    current_observation = observations[-1]
    a = np.sort(observations)
    lowerZ[ndx] = a[a < current_observation][-1]
    higherZ[ndx] = a[a > current_observation][0]

assert np.all(lowerZ == .1)
assert np.all(higherZ == .15)
lowerZ = lowerZ.reshape(z_shape)
higherZ = higherZ.reshape(z_shape)

sliding_windowNumpy を使用した効率的なオーバーラップ ウィンドウから

import numpy as np
from numpy.lib.stride_tricks import as_strided as ast
from itertools import product

def norm_shape(shape):
    '''
    Normalize numpy array shapes so they're always expressed as a tuple, 
    even for one-dimensional shapes.

    Parameters
        shape - an int, or a tuple of ints

    Returns
        a shape tuple

    from http://www.johnvinyard.com/blog/?p=268
    '''
    try:
        i = int(shape)
        return (i,)
    except TypeError:
        # shape was not a number
        pass

    try:
        t = tuple(shape)
        return t
    except TypeError:
        # shape was not iterable
        pass

    raise TypeError('shape must be an int, or a tuple of ints')


def sliding_window(a,ws,ss = None,flatten = True):
    '''
    Return a sliding window over a in any number of dimensions

    Parameters:
        a  - an n-dimensional numpy array
        ws - an int (a is 1D) or tuple (a is 2D or greater) representing the size 
             of each dimension of the window
        ss - an int (a is 1D) or tuple (a is 2D or greater) representing the 
             amount to slide the window in each dimension. If not specified, it
             defaults to ws.
        flatten - if True, all slices are flattened, otherwise, there is an 
                  extra dimension for each dimension of the input.

    Returns
        an array containing each n-dimensional window from a

    from http://www.johnvinyard.com/blog/?p=268
    '''

    if None is ss:
        # ss was not provided. the windows will not overlap in any direction.
        ss = ws
    ws = norm_shape(ws)
    ss = norm_shape(ss)

    # convert ws, ss, and a.shape to numpy arrays so that we can do math in every 
    # dimension at once.
    ws = np.array(ws)
    ss = np.array(ss)
    shape = np.array(a.shape)


    # ensure that ws, ss, and a.shape all have the same number of dimensions
    ls = [len(shape),len(ws),len(ss)]
    if 1 != len(set(ls)):
        error_string = 'a.shape, ws and ss must all have the same length. They were{}'
        raise ValueError(error_string.format(str(ls)))

    # ensure that ws is smaller than a in every dimension
    if np.any(ws > shape):
        error_string = 'ws cannot be larger than a in any dimension. a.shape was {} and ws was {}'
        raise ValueError(error_string.format(str(a.shape),str(ws)))

    # how many slices will there be in each dimension?
    newshape = norm_shape(((shape - ws) // ss) + 1)
    # the shape of the strided array will be the number of slices in each dimension
    # plus the shape of the window (tuple addition)
    newshape += norm_shape(ws)
    # the strides tuple will be the array's strides multiplied by step size, plus
    # the array's strides (tuple addition)
    newstrides = norm_shape(np.array(a.strides) * ss) + a.strides
    strided = ast(a,shape = newshape,strides = newstrides)
    if not flatten:
        return strided

    # Collapse strided so that it has one more dimension than the window.  I.e.,
    # the new array is a flat list of slices.
    meat = len(ws) if ws.shape else 0
    firstdim = (np.product(newshape[:-meat]),) if ws.shape else ()
    dim = firstdim + (newshape[-meat:])
    # remove any dimensions with size 1
    dim = filter(lambda i : i != 1,dim)
    return strided.reshape(dim)
于 2014-06-10T23:16:08.033 に答える