jakevdp / nfft

Lightweight non-uniform Fast Fourier Transform in Python
MIT License
201 stars 29 forks source link

memory issue when using multiple NFFT's #10

Open stormyweathers opened 5 years ago

stormyweathers commented 5 years ago

In my application, I don't know what the optimal number of samples to compute is. So I compute the adjoint and then forward NFFT and check my reconstruction error. I repeatedly call this checker for various sampling factors to find the optimal value. My laptop has enough memory to work up to factor=400.

But when I try to loop through factors=np.arange(200,400,5) it segfaults after 205 due to lack of memory.

I've monitored the memory usage with free -m during execution, and the memory is not freed. I've tried using gc.collect() after each reconstruction check, but that didn't change anything.

I included my sample data that produces this error.

[2019-02-18 17-01_composite.txt](https://github.com/jakevdp/nfft/files/2878118/2019-02-18.17-01_composite.txt)

 import numpy as np
 import scipy.special as sp
 import nfft
 import matplotlib.pyplot as plt
 import gc

 def checkReconstruction(freq,mag,phase,factor):
     #factor=350 seems to work best

     bandWidth=freq[-1]*2
     minFreq=freq[0]
     ## Transform data to be in freq range: (-1/2,1/2) using response(-w) = response(w)^* and rescaling
     freq=np.concatenate((np.flip(-freq),freq))/bandWidth
     mag=np.concatenate((np.flip(mag),mag))
     phase=np.concatenate((np.flip(-phase),phase))
     numSamples=factor*int(np.ceil(bandWidth/minFreq))
     time=(np.arange(numSamples)-numSamples//2)/bandWidth

     invResponse=mag**-1 *np.e**(-1j * phase)

     FIR=nfft.nfft_adjoint(freq,invResponse,numSamples,sigma=100,tol=1E-25)
     invResponse_reconstructed = nfft.nfft(freq,FIR)/numSamples

     error=sum( np.abs(invResponse**-1)- np.abs(invResponse_reconstructed**-1)  )
     print("for oversampling factor {0}  accumulated error={1}".format(factor,error))
     return error

 def sweepSamplingFactor(freq,mag,phase,factors):
     errors=np.array([])
     for factor in factors:
         gc.collect()
         errors=np.append(errors,checkReconstruction(freq,mag, phase,factor) )
         print("error is {0} for an oversampling factor {1}".format(errors[-1],factor))
     combined= np.vstack((factors,errors))
     np.savetxt('sweptSamplingFactor.csv',combined.transpose(),delimiter=",")
     minIdx=np.argmin(errors)
     print("The minimal error is {0} for an oversampling factor {1}".format(errors[minIdx],factors[minIdx]))
     return factors[minIdx]

 if __name__ == '__main__':
   filename='2019-02-18 17-01_composite.txt'
   [freq,mag,phase]= np.loadtxt(filename,delimiter=',',dtype=float)
   optFactor=sweepSamplingFactor(freq,mag,phase,np.arange(200,400,5))
stormyweathers commented 5 years ago

This isn't an issue with this library but with how scipy caches FFT data. I had to include a cache destruction function and call it on each loop iteration. Since the size of the FFTs was different in each loop iteration, scipy creates a new FFT plan, but doesn't automatically destroy the old one

 import scipy.fftpack._fftpack as sff 

 def destroy_fft_cache():
     sff.destroy_zfft_cache()
     sff.destroy_drfft_cache()
     sff.destroy_zfftnd_cache()