1 #ifndef SIGMATRANSFORM_H
2 #define SIGMATRANSFORM_H
20 namespace SigmaTransform {
23 using cmpx = std::complex<double>;
24 using cxVec = std::vector<cmpx>;
25 template<
size_t N>
using diffFunc = std::function<point<N>(
const point<N>&)>;
26 template<
size_t N>
using winFunc = std::function<cmpx(const point<N>&)>;
27 template<
size_t N>
using actFunc = std::function<point<N>(
const point<N>&,
const point<N>&)>;
28 template<
size_t N>
using mskFunc = std::function<cmpx(point<N>
const&,point<N>
const&)>;
49 const std::vector<
point<N>> &steps=std::vector<
point<N>>(0), actFunc<N> action=minus<N> ,
int const& numThreads = 4 )
63 const std::vector<
point<N>> &steps=std::vector<
point<N>>(0), actFunc<N> action=minus<N> ,
int const& numThreads = 4 )
67 if( !fftw_init_threads() )
68 std::cerr <<
"thread error\n";
127 {
m_numThreads = numThreads ; fftw_plan_with_nthreads(numThreads);
return *
this; }
138 if( steps.empty() ) {
139 throw std::runtime_error(
"steps must not be empty.");
142 if( steps.size() == 1 ) {
144 throw std::runtime_error(
"Sampling rate and diffeomorphism must be set for automatic step-array determination");
147 std::array<std::vector<double>,N> stps;
148 for(
int k = 0 ; k < N ; ++k ) {
156 for(
int k = 0 ; k < N ; ++k ) {
157 maxi[k] = (step[k]>maxi[k])?step[k]:maxi[k];
158 mini[k] = (step[k]<mini[k])?step[k]:mini[k];
163 for(
int k = 0 ; k < N ; ++k ) {
164 stps[k] =
linspace( mini[k] , maxi[k] , steps[0][k] );
214 throw std::runtime_error(
"steps not set");
217 throw std::runtime_error(
"Fs not set");
220 throw std::runtime_error(
"size not set");
278 if( mask.size() !=
m_coeff.size() ) {
279 throw std::runtime_error(
"Size of mask does not match size of coefficients.");
282 std::unique_ptr<std::thread[]> _threads(
new std::thread[
m_numThreads] );
284 int stepsPerThread = ceil( (
double) (
m_steps.size()) / m_numThreads );
285 for(
int k = 0, stepsLeft =
m_steps.size() ; k <
m_numThreads ; ++k, stepsLeft -= stepsPerThread ) {
286 _threads[k] = std::move( std::thread( [
this,&stepsPerThread,&mask,stepsLeft,k]() {
288 int offset = k*stepsPerThread*
m_size.prod(),
289 numval = ((stepsLeft>stepsPerThread)?stepsPerThread:stepsLeft)*
m_size.prod();
291 for(
int i = 0 ; i < numval ; ++i ) {
292 m_coeff[offset+i] *= mask[offset+i];
314 std::unique_ptr<std::thread[]> _threads(
new std::thread[
m_numThreads] );
316 int stepsPerThread = ceil( (
double) (
m_steps.size()) / m_numThreads );
317 for(
int k = 0, stepsLeft =
m_steps.size() ; k <
m_numThreads ; ++k, stepsLeft -= stepsPerThread ) {
318 _threads[k] = std::move( std::thread( [
this,&maskFunc,&spatialDom,&stepsPerThread,stepsLeft,k]() {
320 int stepsOffset = stepsPerThread*k,
321 numval = (stepsLeft>stepsPerThread)?stepsPerThread:stepsLeft;
324 for(
int i=0;i<numval;++i) {
325 for(
auto const& x : spatialDom ) {
326 *coeff++ *= maskFunc( x ,
m_steps[stepsOffset+i] );
352 std::unique_ptr<std::thread[]> _threads(
new std::thread[
m_numThreads] );
354 int stepsPerThread = ceil( (
double) (
m_steps.size()) / m_numThreads );
355 for(
int k = 0, stepsLeft =
m_steps.size() ; k <
m_numThreads ; ++k, stepsLeft -= stepsPerThread ) {
356 _threads[k] = std::move( std::thread( [
this,&stepsPerThread,k,stepsLeft]() {
360 for(
int i = 0 ; i < ((stepsLeft>stepsPerThread)?stepsPerThread:stepsLeft) ; ++i ) {
381 m_threads.insert( std::make_pair<std::string,std::thread>(
"Transform" , std::move( std::thread( [
this,sig,onFinish]() {
384 std::unique_lock<std::mutex> lk(
m_mtx );
401 m_threads.insert( std::make_pair<std::string,std::thread>(
"Inverse" , std::move( std::thread( [
this,onFinish]() {
404 std::unique_lock<std::mutex> lk(
m_mtx );
423 m_threads.insert( std::make_pair<std::string,std::thread>(
"Multiplier" , std::move( std::thread( [
this,sig,mask,onFinish]() {
426 std::unique_lock<std::mutex> lk(
m_mtx );
445 m_threads.insert( std::make_pair<std::string,std::thread>(
"Multiplier" , std::move( std::thread( [
this,sig,maskFunc,onFinish]() {
448 std::unique_lock<std::mutex> lk(
m_mtx );
464 std::vector<std::string> toErase;
466 if( it->second.joinable() ) {
469 toErase.push_back( it->first );
473 for(
auto const& name : toErase ) {
475 if( name.compare( it->first) == 0 ) {
492 cxVec
fft( cxVec
const& in ,
int const& howmany = 1 ) {
494 cxVec out( in.size() );
496 fftN( reinterpret_cast<fftw_complex*> (out.data()) ,
497 reinterpret_cast<fftw_complex*>(const_cast<cmpx*> (in.data())) ,
498 m_size , howmany , FFTW_FORWARD );
500 return std::move( out );
512 fftN( reinterpret_cast<fftw_complex*> (inout.data()) ,
513 reinterpret_cast<fftw_complex*>(const_cast<cmpx*> (inout.data())) ,
514 m_size , howmany , FFTW_FORWARD );
524 cxVec
ifft( cxVec
const& in ,
int const& howmany = 1 ) {
526 cxVec out( in.size() );
528 fftN( reinterpret_cast<fftw_complex*> (out.data()) ,
529 reinterpret_cast<fftw_complex*>(const_cast<cmpx*> (in.data())) ,
530 m_size , howmany , FFTW_BACKWARD );
532 return std::move( out );
544 fftN( reinterpret_cast<fftw_complex*> (inout.data()) ,
545 reinterpret_cast<fftw_complex*>(const_cast<cmpx*> (inout.data())) ,
546 m_size , howmany , FFTW_BACKWARD );
557 std::array<std::vector<double>,N> doms;
558 auto itFs =
m_fs.begin(),itSz =
m_size.begin();
559 for(
auto& d : doms) {
572 std::array<std::vector<double>,N> doms;
573 auto itFs =
m_fs.begin(), itSz =
m_size.begin();
574 for(
auto& d : doms) {
575 d =
linspace( 0 , (*itSz-1) / *itFs , *itSz );
590 for(
int k = 0 ; k < N ; ++k ) {
591 if( step[k]>maxi[k] ) {
595 if( step[k]<mini[k] ) {
604 auto width = (maxi-mini) / num_steps *
m_winWidth;
617 double sigsize =
m_size.prod();
621 cxVec Fsig =
fft( in );
625 std::unique_ptr<std::thread[]> _threads(
new std::thread[
m_numThreads] );
627 int stepsPerThread = ceil( (
double) (
m_steps.size()) / m_numThreads );
628 for(
int k = 0, stepsLeft =
m_steps.size() ; k <
m_numThreads ; ++k, stepsLeft -= stepsPerThread ) {
629 _threads[k] = std::move( std::thread( [
this,&Fsig,&stepsPerThread,stepsLeft,k]() {
630 cxVec::iterator coeff =
m_coeff.begin() + (int)( k*stepsPerThread*
m_size.prod() );
631 for(
int i = 0 ; i < ((stepsLeft>stepsPerThread)?stepsPerThread:stepsLeft) ; ++i ) {
632 for(
auto const& val : Fsig ) {
633 *coeff = conj(*coeff++) * val / ((double)Fsig.size());
660 auto coeff = temp.begin(),
662 for(
int k=0;k<
m_steps.size();++k) {
664 val += (*coeff++) * (*window++);
685 void fftN( fftw_complex *out, fftw_complex *in,
const point<N> &size,
const int &howmany = 1,
const int& DIR = FFTW_FORWARD ) {
687 for(
int k = 0 ; k < N ; ++k ) sz[k] = (
int) size[k];
689 fftw_plan p = fftw_plan_many_dft( N , sz , howmany , in , NULL , 1 , (
int) size.prod() ,
690 out , NULL , 1 , (int) size.prod() ,
691 DIR , FFTW_ESTIMATE );
693 fftw_destroy_plan( p );
722 #endif //SIGMATRANSFORM_H