python - resample matlab

enteros espaciados logarítmicamente

Question

Digamos que tengo un vector de 10,000 pts y quiero tomar una porción de solo 100 puntos espaciados logarítmicamente. Quiero que una función me dé valores enteros para los índices. Aquí hay una solución simple que simplemente usa alrededor de + logspace, y luego deshacerse de los duplicados.

def genLogSpace( array_size, num ):
    lspace = around(logspace(0,log10(array_size),num)).astype(uint64)
    return array(sorted(set(lspace.tolist())))-1

ls=genLogspace(1e4,100)

print ls.size
>>84
print ls
array([   0,    1,    2,    3,    4,    5,    6,    7,    8,    9,   10,
         11,   13,   14,   15,   17,   19,   21,   23,   25,   27,   30,
         33,   37,   40,   44,   49,   54,   59,   65,   71,   78,   86,
         94,  104,  114,  125,  137,  151,  166,  182,  200,  220,  241,
        265,  291,  319,  350,  384,  422,  463,  508,  558,  613,  672,
        738,  810,  889,  976, 1071, 1176, 1291, 1416, 1555, 1706, 1873,
       2056, 2256, 2476, 2718, 2983, 3274, 3593, 3943, 4328, 4750, 5213,
       5721, 6279, 6892, 7564, 8301, 9111, 9999], dtype=uint64)

Note que había 16 duplicados, así que ahora solo tengo 84 puntos.

¿Alguien tiene una solución que asegure eficientemente que el número de muestras de salida sea num? Para este ejemplo específico, los valores de entrada para los números 121 y 122 dan 100 puntos de salida.

Answer 1

El enfoque en la respuesta de Avaris de generar sus puntos espaciados por el registro directamente, es definitivamente el camino a seguir. Pero pensé que sería interesante ver cómo elegir el valor apropiado para pasar al logspace de logspace para obtener lo que desea.

Los valores en la matriz generados por el logspace(0, k, n) de logspace(0, k, n) son los números 10 ^{ik / ( n −1)} para 0 ≤ i < n :

>>> numpy.logspace(0, 2, 10)
array([   1.        ,    1.66810054,    2.7825594 ,    4.64158883,
          7.74263683,   12.91549665,   21.5443469 ,   35.93813664,
         59.94842503,  100.        ])
>>> [10 ** (i * 2 / 9.0) for i in xrange(10)]
[1.0, 1.6681005372000588, 2.7825594022071245, 4.641588833612778,
 7.742636826811269, 12.91549665014884, 21.544346900318832,
 35.938136638046274, 59.94842503189409, 100.0]

Esta secuencia consiste en un segmento inicial en el que los valores están más cerca que la unidad espaciada (por lo que puede haber duplicados cuando se redondean al número entero más cercano), seguida de un segmento donde los valores son más amplios que la unidad espaciada y no hay duplicados

>>> '' ''.join(''{:.2f}''.format(10 ** (i * 2 / 19.0)) for i in xrange(20))
''1.00 1.27 1.62 2.07 2.64 3.36 4.28 5.46 6.95 8.86 11.29 14.38 18.33 23.36
 29.76 37.93 48.33 61.58 78.48 100.00''
>>> [int(0.5 + 10 ** (i * 2 / 19.0)) for i in xrange(20)]
[1, 1, 2, 2, 3, 3, 4, 5, 7, 9, 11, 14, 18, 23, 30, 38, 48, 62, 78, 100]

El espaciado entre los valores es s ( i ) = 10 ^iK - 10 ^{( i −1) K} , donde K = k / ( n - 1). Sea m el valor más pequeño tal que s ( m ) ≥ 1. ( m = 7 en el ejemplo anterior). Luego, cuando se eliminan los duplicados, quedan exactamente ½½ +10 ^{( m −1) K} ⌋ + n - m restante números.

Un poco de álgebra encuentra:

m = ⌈ - log (1 - 10 ^{- K} ) / K log 10

Vamos a comprobar eso.

from math import ceil, floor, log

def logspace_size(k, n):
    """
    Return the number of distinct integers we''ll get if we round
    `numpy.logspace(0, k, n)` to the nearest integers and remove
    duplicates.

    >>> logspace_size(4, 100)
    84
    >>> logspace_size(4, 121)
    100
    >>> from numpy import around, logspace
    >>> all(logspace_size(k, n) == len(set(around(logspace(0, k, n))))
    ...     for k in xrange(1,10) for n in xrange(2,100))
    True
    """
    K = float(k) / (n - 1)
    m = int(ceil(- log(1 - 10 ** -K) / (K * log(10))))
    if m < n:
        return int(0.5 + 10 ** ((m - 1) * K)) + n - m
    else:
        return int(0.5 + 10 ** ((n - 1) * K))

Las pruebas pasan, así que esto me parece bien. Así que todo lo que necesita hacer es encontrar n tal que logspace_size(4, n) == 100 . Puede hacer esto por medio de un corte binario o uno de los métodos scipy.optimize :

>>> f = lambda x, k, n:(logspace_size(k, x) - n)**2
>>> int(round(scipy.optimize.fmin(f, 100, args=(4,100), xtol=0.5, ftol=0.5)[0]))
Optimization terminated successfully.
         Current function value: 0.015625
         Iterations: 8
         Function evaluations: 17
122

Answer 2

Esto es un poco complicado. No siempre se pueden obtener números espaciados logarítmicamente. Como en tu ejemplo, la primera parte es más bien lineal. Si estás de acuerdo con eso, tengo una solución. Pero para la solución, debes entender por qué tienes duplicados.

La escala logarítmica satisface la condición:

s[n+1]/s[n] = constant

Llamemos a esta constante r para la ratio . Para n de estos números entre el rango 1...size , obtendrá:

1, r, r**2, r**3, ..., r**(n-1)=size

Así que esto te da:

r = size ** (1/(n-1))

En su caso, n=100 y size=10000 , r será ~1.0974987654930561 , lo que significa que si comienza con 1 , su próximo número será 1.0974987654930561 que luego se redondeará nuevamente a 1 . Así tus duplicados. Este problema está presente para los pequeños números. Después de un número suficientemente grande, multiplicar con la relación resultará en un entero redondeado diferente.

Teniendo esto en cuenta, lo mejor es agregar números enteros consecutivos hasta cierto punto para que esta multiplicación con la relación ya no sea un problema. Luego puedes continuar con la escala logarítmica. La siguiente función hace eso:

import numpy as np

def gen_log_space(limit, n):
    result = [1]
    if n>1:  # just a check to avoid ZeroDivisionError
        ratio = (float(limit)/result[-1]) ** (1.0/(n-len(result)))
    while len(result)<n:
        next_value = result[-1]*ratio
        if next_value - result[-1] >= 1:
            # safe zone. next_value will be a different integer
            result.append(next_value)
        else:
            # problem! same integer. we need to find next_value by artificially incrementing previous value
            result.append(result[-1]+1)
            # recalculate the ratio so that the remaining values will scale correctly
            ratio = (float(limit)/result[-1]) ** (1.0/(n-len(result)))
    # round, re-adjust to 0 indexing (i.e. minus 1) and return np.uint64 array
    return np.array(list(map(lambda x: round(x)-1, result)), dtype=np.uint64)

Actualización de Python 3: la última línea solía ser return np.array(map(lambda x: round(x)-1, result), dtype=np.uint64) en Python 2

Aquí hay algunos ejemplos usándolo:

In [157]: x = gen_log_space(10000, 100)

In [158]: x.size
Out[158]: 100

In [159]: len(set(x))
Out[159]: 100

In [160]: y = gen_log_space(2000, 50)

In [161]: y.size
Out[161]: 50

In [162]: len(set(y))
Out[162]: 50

In [163]: y
Out[163]:
array([   0,    1,    2,    3,    4,    5,    6,    7,    8,    9,   11,
         13,   14,   17,   19,   22,   25,   29,   33,   38,   43,   49,
         56,   65,   74,   84,   96,  110,  125,  143,  164,  187,  213,
        243,  277,  316,  361,  412,  470,  536,  612,  698,  796,  908,
       1035, 1181, 1347, 1537, 1753, 1999], dtype=uint64)

Y solo para mostrar qué tan logarítmicos son los resultados, aquí hay un gráfico semilog de la salida para x = gen_log_scale(10000, 100) (como puede ver, la parte izquierda no es realmente logarítmica):

Answer 3

Llegué aquí mientras buscaba un método simple para obtener series espaciadas logarítmicamente (con base de 10) en python (omitiendo el uso de numpy). Pero sus soluciones son demasiado complicadas para mis demandas ultra simples.

def logarithmic_decade(numbers_per_decade, offset=10):
    for n in xrange(numbers_per_decade):
        yield offset * 10.0 ** (n / float(numbers_per_decade))

Ya que es generador para obtener la lista tienes que:

numbers = list(logarithmic_decade(5))
print numbers
[10.0, 15.848931924611136, 25.118864315095802, 39.81071705534972, 63.095734448019336]

for p, n in zip(numbers, numbers[1:] + [100]):
    print ''prev = {p:.2f}, next = {n:.2f}, next/prev = {rt:.4f}''.format(p=p, n=n, rt=n / p)

Da siguiente el siguiente resultado:

prev = 10.00, next = 15.85, next/prev = 1.5849
prev = 15.85, next = 25.12, next/prev = 1.5849
prev = 25.12, next = 39.81, next/prev = 1.5849
prev = 39.81, next = 63.10, next/prev = 1.5849
prev = 63.10, next = 100.00, next/prev = 1.5849