How to combine coordinates and pixel values

[deeptho]

I have been experimenting with your code, with nice results.

However, one thing I am really struggling with is how to implement
operations that need both the coordinates of pixels and the pixel values
themselves. Although my use case needs something more complex, here
is a simple (artificial) program in pseuode code

double input[100]={1,3,100,50, ...} ; //an array containing input data, e.g. grey values
double output[100]; //an array for storing output

double fn(int x, double v)
{
return x/(100.)+v;
}

for(int x =0; x<100; ++x)
  output[i] = fn(x, input[x]);

How is one supposed to handle this in zimt?
E.g., transform does not work in this case as it operates either on coordinates or on values,
but ot both.

Maybe a solution could be to store a reference to the input array in the
unary functor, use the coordinates version of transform and then index "input"
in the functor? This seems quite involved and perhaps not efficient. Also my real code
is m ore complex and uses a 2D linspace, but the solution would not allow to use it (as the functor
needs to receive the coordinates to index the stored references array.

What is added to v in fn is of course the result of a linspace. I could store the result
in a temporary array, combine it with input using a join, but that seems inefficient as well.

What seems to be missing is some functionality to allow, e.g., linspace or get_crd to be passed in combination with array views as the input to transform or process. In essence all of these produce "data at some coordinates", one by generating the data and the other by loading it from disk

Maybe I am missing something?

[kfjahnke]

No, you're right. You need a different strategy.
zimt's functional composition facilities are not made to handle binary functors - I've kept it deliberately simple. You can think of zimt's functor chaining mechanism as a sort of lambda calculus - it's always a single input and a single output, à la void ( const in_t & , out_t & ). To implement n-ary functors, you have to use a different strategy. You set up a wrapper functor. This outer functor can then call one or more inner functors, combine their results and write them to the output. The outer functor is coded as to take a single argument. To use your example, you would code the outer functor as taking int and passing double, and inside the outer functor you'd use the ints to gather doubles from 'input'.
Let's look at your concrete example. it uses 'x' in two different ways - as a factor (denominator in x/100) and as an index (input[x]). Common to both uses is x itself, so this is what you'd feed as input to the functor. Inside the functor you build the result and pass it out. zimt has an overload of zimt::transform which takes no input array - instead the 'act' functor receives canononical indices. This overload exists because such indices are such a common starting point.

Let's code some. I assume your artificial example should read

for(int x =0; x<100; ++x)
  output[x] = fn(x, input[x]);

(it says output[i] in your post)

Here's how you'd get the functionality with zimt:

// compile with e.g. g++ binary_f.cc
// note that zimt headers must be reachable

#include <zimt/zimt.h>

// next we define the fn functor, taking int input and producing
// double. note that we don't use a linspace input: we need the
// ints coming in as indices to gather from 'v', so if we had linspace
// input, we'd have to 'calculate backwards' to regain them.

struct fn_t
: public zimt::unary_functor < int , double , 4 >
{
  const double * p_input ; // pointer to your 'input' array, ex c'tor
  
  fn_t ( const double * _p_input )
  : p_input ( _p_input )
  { }

  // zimt functors need an 'eval' function. this is also most conveniently
  // coded as a template. incoming, we'll have canonical indices, outgoing
  // the result of 'fn'

  template < typename in_t , typename out_t >
  void eval ( const in_t & in , out_t & out )
  {
    std::cout << "eval - in:       " << in << std::endl ;

    // let's have a variable of the output's type:

    out_t v ;

    // we gather from p_input, using 'in' as gather indices

    v.gather ( p_input , in ) ;

    std::cout << "eval - gathered: " << v << std::endl ;

    // and we form fn(x,v)

    out = in / 100.0 + v ;

    std::cout << "eval - out:      " << out
              << std::endl << std::endl ;
  }
} ;

// finally, with two zimt arrays, we can do this:

int main ( int argc , char * argv[] )
{
  zimt::array_t < 1 , double > input ( 10 ) ;
  zimt::array_t < 1 , double > output ( 10 ) ;

  // let's initialize the input array to something recognizable

  for ( std::size_t i = 0 ; i < 10 ; i++ )
  {
    input[i] = i *.0001 ;
  }

  // now the 'zimt magic'. first we set up the fn functor with a pointer
  // to the 'input' data

  fn_t fn ( input.data() ) ;

  // now we call the zimt::transform overload feeding canonical indices
  // to fn, storing the result in the zimt::array 'output'

  zimt::transform ( fn , output ) ;

  // bingo:

  for ( std::size_t i = 0 ; i < 10 ; i++ )
  {
    std::cout << "output[" << i << "] = " << output[i] << std::endl ;
  }
}

and on my system, I get this output:

$ ./a.out
eval - in:       { 0, 1, 2, 3 }
eval - gathered: { 0, 0.0001, 0.0002, 0.0003 }
eval - out:      { 0, 0.0101, 0.0202, 0.0303 }

eval - in:       { 4, 5, 6, 7 }
eval - gathered: { 0.0004, 0.0005, 0.0006, 0.0007 }
eval - out:      { 0.0404, 0.0505, 0.0606, 0.0707 }

eval - in:       { 8, 9, 6, 7 }
eval - gathered: { 0.0008, 0.0009, 0.0006, 0.0007 }
eval - out:      { 0.0808, 0.0909, 0.0606, 0.0707 }

output[0] = 0
output[1] = 0.0101
output[2] = 0.0202
output[3] = 0.0303
output[4] = 0.0404
output[5] = 0.0505
output[6] = 0.0606
output[7] = 0.0707
output[8] = 0.0808
output[9] = 0.0909

The 'unary functors only' paradigm does at times require some rearrangements in the operation, and figuring out 'x' as the 'common ancestor' of both operands leads us to an implementation as the one above - a zimt::linspace object as input to zimt::process doesn't work here, because we would have to regain the canonical indices needed for gathering from 'input'. Of course there are other ways of coding the task at hand, but this strategy is typical.

You note that using an intermediate array is inefficient, and indeed functional programming can avoid intermediate mass storage most of the time, resulting in code which 'happens' only in registers, apart from necessary read operations from input data and storage to output data. This is maximally efficient.

I hope this helps you solve your problem - don't hesitate to ask if something's not clear.

Out of curiosity: what's your intended use case?

[deeptho]

Thanks for the quick answer. I find it useful as I have learned something.
The code also compiles and runs on my computer.

In the mean time, I tried a different approach: I made a version "process2" of "process" that accepts two getters.
This also required some minor changes vs_adapter, where I added an "eval2" function, which then calls
a "unary functor" which of course is now really a binary functor.

I then call this with linspace as one getter and a loader as the second getter. This way, at least there is no need
for linspace to save to memory first.

The use case is related to time of flight cameras:

1. I start from 9 images, which are threated as three groups of three.
2. For each pixel in each group I compute a complex number, which is stored as a (real, imag) pair as complex numbers
   are not supported by zimt (but that is a small detail). This output can be seen as 6 images
   (3 for the real parts and and 3 for the imaginary parts). All processing is "unary functor"/
3. Each of those 6 images needs to be filtered spatially (using convolve), which means that the 6 input images and
   the 6 output images are stored in memory. There is some efficiency here: except to run the convolve, I
   have no need for keeping the input and output images; the output images will be immediately processed and then discarded. I understand it may be quite difficult to eliminate this inefficiency, so I am prepared to accept it.
4. The 6 output images are now treated as 3 "complex-valued" images, which need to be phase corrected based with the correction depending on the pixel location. This is where I got stuck and opened the ticket and ended up with converting
   the unary functor into a binary one.
5. Now the 6 results are treated as 3 complex numbers. They are combined into a single complex numer
   I have already implemented steps 1-4 with an output to memory. I can also see how to run step 5 by reading from that memory in which step 4 saves.

However, I would also like to eliminate this needless "storing to memory" between 4 and 5, which should be
much easier to solve than the inefficiency in step 3, as there is only "pixel per pixel" processing involved in steps 4 and 5.
Again: I am not sure how to implement this cleanly.... Suggestions are welcome.

I have already implemented all of this without zimt, directly in C++. Using the standard clang++ optimisations,
and carefully written code and multiprocessing, but no explicit vectorised code. For example, one optimisation is
to store only a few lines of the output image in step 3, which saves a bit of memory and perhaps of perhaps not even some computation time. The main optimisations are of course related to paying attention to memory accesses.

The detailed documentation in your code proved a great help, but the learning curve still was/is quite steep. Although I
program in C++ a lot I also find c++ code horribly complex in general. I think your type-erasure strategy at least helps to easy some of that pain. There were still some horrible compiler messages from time to time.

I already see quite high performance without zimt (200 frames per second on standard intel pc and 25 fps on some arm cpu; each frame is 1000x1000 pixels) , but I am wondering how much more zimt can squeeze out, compared to the compiler itself. I am also wondering how easy it would be to vectorize some more complicated (but still regular) code
from a programming point of view. E.g., stereo processing algorithms, which have data dependent branches.

PS. My modifications are currently just a hack: undocumented, not cleanly named etc. I can provide them if you would like to inspect them, but they are not yet of good quality. The system of provind more than one getter could also be easily generalized to providing an array of getters.

[kfjahnke]

complex numbers are not supported by zimt

They are not supported explicitly, but many things you do with complex numbers can be done in just the same way using e.g. a zimt::xel_t < float , 2 >, which can be used with 'normal' arithmetic like any other xel_t - the channels are linearly independent. Of course then you need to code complex-specific operations 'on foot'.

you might consider bunching processing of several images into one pipeline - there is e.g. the join_t to pull data from a set of N arrays of fundamentals into a pipeline processing xel_t of N elements. And when it comes to processing several images at the same time, the code I gave above might serve as an ispiration: start out feeding coordinates, then gather from several images. so instead of just one pointer in fn_t, you pass several, gather to vectors v1, v2, v3 and then you can combine them any way you want.

convolution does of course throw a spanner in the works of a single pixel pipeline - the specialized convolution code in zimt for separable convolutions processes entire arrays only. but it's good code, processing the arrays with multithreaded SIMD code. where convolutions are involved, it's often possible to move the convolution to earler or later in the process. If that's an option, you might convolve the input images in-place and then start the pipeline with the convolution results.

the learning curve still was/is quite steep

Sorry for that. I also took a good while to get into that kind of programming. I started out using 'generic programming' techniques using U. Koethe's excellent vigra library, which isn't evolving much anymore these days. vigra's TinyVector was my template for zimt::xel_t, and vigra's arrays (which in turn are similar to NumPy ndarrays) were the inspiration for zimt::array_t. When I wrote my vspline library, I used vigra for *xel and array data. But the vigra design is quite old and predates a lot of modern C++, and I tried to come up with a leaner and more modern design, omitting some of vigra's complexity along the way. Since vigra is now no longer evolving and my attempts to interest vigra for SIMD processing failed, I wrote zimt. There is even a compatibility layer in zimt to integrate with vigra code. The SIMD processing and efficient multithreading also started out in vspline, and convolution, IIR filtering and b-spline processing in zimt are actually vspline code ported to zimt. If you want to see my code 'in action', you might be interested in two of my projects: the lux image and panorama viewer which still uses vigra/vspline and envutil. There you can see a lot of code working pixel pipelines and handling imagery which might inspire you.

I think your type-erasure strategy at least helps to easy some of that pain. There were still some horrible compiler messages from time to time.

Yes, generic code produces lots of errors. mostly all it takes is understanding and fixing the first one, and then the whole rat's tail disappears. but it's definitely a pain to debug. 'grokking' is really helpful. I also set up code to grok get_t and put_t in zimt, not only the 'act' functor.

wondering how much more zimt can squeeze out

This also depends on the backend. For lux and envutil, I now favour google's highway, and I found that using xel_t with twice as many lanes as the hardware registers contain often makes the code even a tad faster. I think my 'bushfire' strategy of multithreading has next to no overhead and with SIMD and good optimization I'm quite confident zimt code will perform well. The biggest gains are to be had by avoiding intermediate arrays and setting up a single pipeline from beginning to end. In lux and envutil I have complex pipelines: I start out with image coordinates, convert to 3D polar (rays), rotate some with quaternions, project to a plane again, interpolate and postprocess to output data format. All of that is in a single pipeline, so I can run at 60 fps fullHD if the interpolators aren't too expensive and the geometric transforms don't need too many transcendental functions.

My modifications are currently just a hack

No problem. That's just how I work: hack my path with the machete, later come back with the bulldozer ;)
Feel free to point me to some of your work, and if you feel you have code which fits well into zimt, I'll take a look!

[deeptho]

I will upload my code on github soon.

complex numbers are not supported by zimt

They are not supported explicitly, but many things you do with complex numbers can be done in just the same way using e.g. a zimt::xel_t < float , 2 >, which can be used with 'normal' arithmetic like any other xel_t - the channels are linearly independent. Of course then you need to code complex-specific operations 'on foot'.

This is what I am doing now and it works well.

you might consider bunching processing of several images into one pipeline - there is e.g. the join_t to pull data from a set of N arrays of fundamentals into a pipeline processing xel_t of N elements. And when it comes to processing several images at the same time, the code I gave above might serve as an ispiration: start out feeding coordinates, then gather from several images. so instead of just one pointer in fn_t, you pass several, gather to vectors v1, v2, v3 and then you can combine them any way you want.

This is what I had in mind. The details do require modifying some of your code rather than using it "as is",
but it remains reasonably simple.

convolution does of course throw a spanner in the works of a single pixel pipeline - the specialized convolution code in zimt for separable convolutions processes entire arrays only. but it's good code, processing the arrays with multithreaded SIMD code. where convolutions are involved, it's often possible to move the convolution to earler or later in the process. If that's an option, you might convolve the input images in-place and then start the pipeline with the convolution results.

I think your approach is indeed efficient. My old approach relies on some redundant computations instead.

the learning curve still was/is quite steep

Sorry for that. I also took a good while to get into that kind of programming. I started out using 'generic programming'

Well, I think Bjarne Stroestrup is the one who should be a it sorry (despite his acievements): c++ has grown too
complex, with too many rules, too many exceptions an still many simple things that are missing. Debugability is an
issue as well. I do not have a better alternative for it, so I stick with it

techniques using U. Koethe's excellent vigra library, which isn't evolving much anymore these days. vigra's TinyVector was my template for zimt::xel_t, and vigra's arrays (which in turn are similar to NumPy ndarrays) were the inspiration for zimt::array_t. When I wrote my vspline library, I used vigra for *xel and array data. But the vigra design is quite old and predates a lot of modern C++, and I tried to come up with a leaner and more modern design, omitting some of vigra's complexity along the way. Since vigra is now no longer evolving and my attempts to interest vigra for SIMD processing failed, I wrote zimt. There is even a compatibility layer in zimt to integrate with vigra code. The SIMD processing and efficient multithreading also started out in vspline, and convolution, IIR filtering and b-spline processing in zimt are actually vspline code ported to zimt. If you want to see my code 'in action', you might be interested in two of my projects: the lux image and panorama viewer which still uses vigra/vspline and envutil. There you can see a lot of code working pixel pipelines and handling imagery which might inspire you.

I have read about your approach. I think it is a very sound one. I have opted to only use the hwy variant for now.

I think your type-erasure strategy at least helps to easy some of that pain. There were still some horrible compiler messages from time to time.

Yes, generic code produces lots of errors. mostly all it takes is understanding and fixing the first one, and then the whole rat's tail disappears. but it's definitely a pain to debug. 'grokking' is really helpful. I also set up code to grok get_t and put_t in zimt, not only the 'act' functor.

wondering how much more zimt can squeeze out

This also depends on the backend. For lux and envutil, I now favour google's highway, and I found that using xel_t with twice as many lanes as the hardware registers contain often makes the code even a tad faster. I think my 'bushfire' strategy of multithreading has next to no overhead and with SIMD and good optimization I'm quite confident zimt code will perform well. The biggest gains are to be had by avoiding intermediate arrays and setting up a single pipeline from beginning to end. In lux and envutil I have complex pipelines: I start out with image coordinates, convert to 3D polar (rays), rotate some with quaternions, project to a plane again, interpolate and postprocess to output data format. All of that is in a single pipeline, so I can run at 60 fps fullHD if the interpolators aren't too expensive and the geometric transforms don't need too many transcendental functions.

Sound advice indeed.

My modifications are currently just a hack

No problem. That's just how I work: hack my path with the machete, later come back with the bulldozer ;) Feel free to point me to some of your work, and if you feel you have code which fits well into zimt, I'll take a look!

I wlll upload them soon, so that you can have a look.

[kfjahnke]

I dug into my code base and noticed that I have code which may fit your need in my envutil project, or at least serve as an inspiration. In envutil, my code is often ahead of what's in zimt - it's where I try out new stuff, or implement things which are needed specifically for envutil Once the code matures, I update the zimt library, but envutil's zimt code has two new getter classes: fusion_t and zip_t. I plain forgot I had them, I should have told you right away! Please have a look at get.h in the zimt folder in the envutil repo. I use the code in envutil, but it may not be fully 'library grade'.

[deeptho]

I will have a look at envutil

In the mean time I have working code, although not well tested yet but at least
it seems to produce correct results for the most common cases.

There are a few places in my code with conditional computations. So it took a while to figure out how to do that well (with "masks" and with"allof" to detect the fast path). Other things that took some time was simple things such as rounding,
modulo operations.

My code has two versions:
A: one which outputs some extra data needlessly (for debugging)
B: one which does not and therefore is by nature much faster

Unfortunately version B still crashes with zimt (work in progress, probably something simple).

Version A works but still saves data to memory in between steps 4 and 5. Also, due to limitations in
zimt (only one data type possible in the putter for process) the zimt version saves some of that data
as double rather than int So zimt is a bit penalized.

So in summary: version A is not the most representative, but as it is working, I can run some timings:

These are the results on my modern desktop:
phase_diff_from_images: fps=147.6 Call took 6773us (32 threads 1000 runs)
ZIMT: phase_diff_from_images: fps=186.0 Call took 5376us (32 threads 1000 runs)

(note that the reported number of threads for zimt is incorrect. I think you use double that by default,
so it is probably 64)

So the speedup is not yet very large (yet): only 25%

On my laptop, which is less powerful, the difference between pure c++ code and zimst is negligible. I get around 38
fps in each case (c++ versus zimt). Of course this is before any meta parameter optimisation.

Still on my todo list
1 fix version B and time that. This version will have far fewer memory transfers and there might be larger
speedups with zimt as te effect of optimizing with SIMD may be larger
2. remove the needless memory re-reading between steps 4 and 5. This will make the zimt version faster
3. experiment with some parameters (lane size and such)

Any tips for the parameter optimisation?

Learnings for zimt: it may be useful to add some more variety in both put and get and the functors:
-coordinates, and multiple arrays with multiple data types as inputs
-multiple arrays with multiple data types as outputs.

[kfjahnke]

I will have a look at envutil

The link is specifically to the get.h header in envutil's zimt folder, which has more get_t classes than the library version. The last two classes in that header look like they might be helpful for your problem.

So the speedup is not yet very large (yet): only 25%

At least you do get some speedup. Whether an explicit SIMD implementation actually runs faster than what the compiler outputs by way of optimization can't be foretold. The gain tends to increase if computational load is higher, but having to work with masking to emulate conditionals in conventional code is better avoided. It looks as if you're using double precision. Are you sure you need it? 32 bit floats are often perfectly good enough, and the speedup is large. And if you mix types of different sizes this is often bad for performance - so if you have integers and floating point, best use int32 and float32 together. long ints don't usually verctorize well.

Learnings for zimt: it may be useful to add some more variety in both put and get and the functors:
-coordinates, and multiple arrays with multiple data types as inputs
-multiple arrays with multiple data types as outputs.

One way of achieving that in a 'generic' way is to work with tuples as arguments, but the resulting code might become quite involved. when it comes to input, I think you pass in coordinates as ints, because it feels natural. But I guess you have to convert them to floating point somewhere down the line to do calculations with them. You might be better off feeding float coodinates and float values throughout and just use a single larger xel_t in the pipeline which contains both - e.g. v[0] and v[1] for x, y and v[2] and v[3] for re and im.

zimt gives a default value for the lane count, which I set because it seemed to work well for me, but you can easily experiment with different values. If, for example, your code produces strong register pressure you might be better off with the lane count as large as one register - e.g. 8 for float/int on AVX2. Any small multiple of the SIMD register's size could be good.

You seem to have many cores. zimt grabs many threads then, and this may turn out too much of a good thing. Did you notice the optional 'bill' argument to zimt's processing functions? You can set bill.njobs to different values than the default to experiment with fewer or more threads. I often even use njobs=1 for debugging, because with N threads running concurrently it's sometimes hard to figure out what wen wrong. Your large core count seems to indicate you have a system with P and E cores. The E cores may not contribute too much to computation. memory traffic may also become an issue with too many threads.

[deeptho]

The link is specifically to the get.h header in envutil's zimt folder, which has more get_t classes than the library version. The last two classes in that header look like they might be helpful for your problem.

Thanks!

At least you do get some speedup. Whether an explicit SIMD implementation actually runs faster than what the compiler outputs by way of optimization can't be foretold. The gain tends to increase if computational load is higher, but having to work with masking to emulate conditionals in conventional code is better avoided. It looks as if you're using double precision. Are you sure you need it? 32 bit floats are often perfectly good enough, and the speedup is large. And if you mix types of different sizes this is often bad for performance - so if you have integers and floating point, best use int32 and float32 together. long ints don't usually verctorize well.

The conditional code is only run a very small fraction of the time (less than 1%) and it has to be there.
About double versus single precision: I already tried that with my own code and there was not much difference.
I will check this of course.

One way of achieving that in a 'generic' way is to work with tuples as arguments, but the resulting code might become quite involved. when it comes to input, I think you pass in coordinates as ints, because it feels natural. But I guess you have to convert them to floating point somewhere down the line to do calculations with them. You might be better off feeding float coodinates and float values throughout and just use a single larger xel_t in the pipeline which contains both - e.g. v[0] and v[1] for x, y and v[2] and v[3] for re and im.

That is why I currently pass in linspace. In order to integrate even more (to avoid some intermediate storage),
I need to pass in three different linspaces (or one coordinate set, but then I need more multiplications)

zimt gives a default value for the lane count, which I set because it seemed to work well for me, but you can easily experiment with different values. If, for example, your code produces strong register pressure you might be better off with the lane count as large as one register - e.g. 8 for float/int on AVX2. Any small multiple of the SIMD register's size could be good.

I tried that already. There are observable differences in speed, but not major ones.

You seem to have many cores. zimt grabs many threads then, and this may turn out too much of a good thing. Did you notice the optional 'bill' argument to zimt's processing functions? You can set bill.njobs to different values than the default to experiment with fewer or more threads. I often even use njobs=1 for debugging, because with N threads running concurrently it's sometimes hard to figure out what wen wrong. Your large core count seems to indicate you have a system with P and E cores. The E cores may not contribute too much to computation. memory traffic may also become an issue with too many threads.

I also reduced the number of jobs to 32. It had a marginal effect.
With minimal_output=True I skip writing 5 large extra arrays, but the zimt version still needs
one large array extra.

In the mean time, I made the code that produces a "minimal output" work. Disappointingly, now zimt is slower than
my own code. I guess this all comes down to still using one extra intermediary storage step.

Here are some typical results. With pre_filter_uniform=-1 I skip the 3x3 filter
with pre_filter_uniform=3 I don't

------------------------
minimal_output=True pre_filter_uniform=-1
ZIMT: phase_diff_from_images: fps=228.6 Call took 4374us (32 threads 1000 runs)
phase_diff_from_images: fps=318.3 Call took 3142us (32 threads 1000 runs)

-------------------------
minimal_output=False pre_filter_uniform=-1
ZIMT: phase_diff_from_images: fps=186.6 Call took 5360us (32 threads 1000 runs)
phase_diff_from_images: fps=147.6 Call took 6777us (32 threads 1000 runs)
-------------------------
minimal_output=True pre_filter_uniform=3
ZIMT: phase_diff_from_images: fps=227.1 Call took 4403us (32 threads 1000 runs)
phase_diff_from_images: fps=205.1 Call took 4875us (32 threads 1000 runs)
-------------------------
minimal_output=False pre_filter_uniform=3
ZIMT: phase_diff_from_images: fps=186.6 Call took 5359us (32 threads 1000 runs)
phase_diff_from_images: fps=122.7 Call took 8152us (32 threads 1000 runs)

So currently my own code (c++ only) often is on par (tiny changes in code sometimes produce
suprising speedups or slowdowns) with zimt, or even faster
when filtering is not needed. With filtering on zimt is the winner.

There are reasons for this

• ineffficiencies by trying to use current zimt interfaces : I think I need a "process" with multiple
  getters of multiple types (e.g. 3 linspaces and one join_t all stored in one container "getters")
• my own still not perfect use of zimt. E.g. vstorer may be more efficient for the intermediary storage
• my own spatial filter code does not use any memory but instead it performs redundant
  computations (and quite a few of them).
  Basically it is a sliding window filter that computes in each window the average of f(data(i,j)).
  f(.) does require several additions/multiplications, but I do not cache the results, leading to 9x redundant computations
  when the filter is on. Similar filters often occur in image processing.
  The zmt version first stores the f() results in a large image (6 million elements)
  and then runs your separable filter (a modified version as all weights are equal)
  I need to figure out some hybrid that combines the best of both versions.

[deeptho]

Regarding envutil:

I think the "synopsis" getter is indeed something that could be of use to me.

[kfjahnke]

Some of your speed test results are puzzling. Did you notice that the result of your zimt-based tests are the same for using and not using the sliding average filter? That seems odd - it should make some difference, and for the tests of your own code without zimt, the differences are significant.
Just in case you don't - make sure you measure the relevant part of program runs only - so, grab a high-precision time point just before and after the workhorse code happens, then form the difference. And I assume you average over a few tens of runs, rather than just using 'time' on the CL. Do try runs with less threads than you have cores as well, down to single-threading the code.
Also, would you share your precise compiler invocation - especially for the zimt builds - with all flags you pass? I can have a look and see if something's amiss or can be improved.

[deeptho]

Some of your speed test results are puzzling. Did you notice that the result of your zimt-based tests are the same for using and not using the sliding average filter? That seems odd - it should make some difference, and for the tests of your own code without zimt, the differences are significant.

I explained why: my own code performs 9x redundant input computations when applying the filter.
It is probably not the filter but those needless computations that show up in the results. The filter itself
has little effect in zimt, probably because it is is negligible compared to memory access, which happens
with and without the filter.

Just in case you don't - make sure you measure the relevant part of program runs only - so, grab a high-precision time point just before and after the workhorse code happens, then form the difference. And I assume you average over a few tens of runs, rather than just using 'time' on the CL. Do try runs with less

Of course: the times is started before 1000 runs and then stopped.

threads than you have cores as well, down to single-threading the code. Also, would you share your precise compiler invocation - especially for the zimt builds - with all flags you pass? I can have a look and see if something's amiss or can be improved.

Sure, that can be done, but it may take a while as I know have to focus on some other things first.

The parameters are basically copied from "examples.sh", except for the c++ version (to be
compatible with the rest of the code). It is possible that CMake has added a few

/usr/lib64/ccache/clang++ -Dzimt_EXPORTS -I/xxx/zimt/include --std=c++23 -Wall -Wextra -fno-omit-frame-pointer -fvisibility=default -O3 -fPIC -Wno-unused-parameter -Wno-deprecated-copy -fPIC -g -O3 -DUSE_HWY --std=c++23 -msse2 -mssse3 -msse4.1 -msse4.2 -mpclmul -maes -mavx -mavx2 -mbmi -mbmi2 -mfma -mf16c -MD -MT zimt/CMakeFiles/zimt.dir/phasors_zimt.cc.o -MF CMakeFiles/zimt.dir/phasors_zimt.cc.o.d -o CMakeFiles/zimt.dir/phasors_zimt.cc.o -c /xxx/zimt/phasors_zimt.cc

[deeptho]

A small update, with more or less the same code: I have now tested the single precision version of the code,
included for my own code. In my own code (not zimt) I discovered a few places where computations were still
done in double precision, even when they were supposed to run in single precision. After correcting this, there
is now a clear difference between single/double precision.

I also made a 2nd version of my code (pure c++) to avoid some
needless recomputations of the input to the sliding window filter, at the expense of storing them temporarily
in memory.

So the comparison is

1. zimt, with a few modifications like a variant of 'process'
2. my existing c++ code . This just recomputes all inputs needed for the sliding filter instead of storing the inputs somewhere. This overhead is acceptable as the most useful filter size is only 3x3
3. new c++ code, where the filter implementation does not exploit separability of the kernel (so still some needless computations), but inputs are computed only once and cached in memory

---------------
minimal_output=True pre_filter_uniform=3 with_multipath_correction=False
ZIMT:  fps=413.1 Call took 2421us (32 threads 1000 runs)
old c++: fps=234.2 Call took 4269us (32 threads 1000 runs)
new c++: fps=271.4 Call took 3684us (32 threads 1000 runs)
-------------------------
minimal_output=True pre_filter_uniform=5 with_multipath_correction=False
ZIMT: fps=410.3 Call took 2437us (32 threads 1000 runs)
old c++: fps=116.7 Call took 8569us (32 threads 1000 runs)
new c++: fps=233.5 Call took 4282us (32 threads 1000 runs)

Although 3x3 is the needed filter size, I also included 5x5 to show how complexity changes.
ZIMT achieves over 400 fps and (as expected) its complexity does not depend very strongly on the filter size.
On the other hand my own code benefits from the caching (new versus old version) and already achieves 66% of zimt speed, without explicit use of SIMD.

There are still many inefficiencies in my zimt implementation, but I have not had time to experiment further
with that.

For instance the current processing results in:

1. pack multi-image data into simd vectors
2. perform some "per pixel location" processing;
3. unpack the simd vectors and store them as multiple images
4. perform the convolution (simplified to a sliding window sum), reading from images and writing to images
5. pack multi-image data into simd vectors
6. perform some more "per pixel location" processing
7. unpacks simd vectors and store as images.

Would it be possible to eliminate steps 3 and 5? These steps are used to glue process/convolve/process

Alternatively, I am thinking of just directly making one big "process" like function and then runs code
similar to my own c++ program, in which "convolve" is not a separate step, but rather a "dumb" combination
of operations on neighbouring pixel values. That requires some deviations of zimt's functor design, which then no longer acts on single pixels. Or I could just not use functors and directly use some of the lower level zimt concepts.

This would be mostly a learning exercise. 400 fps is already nice.

[deeptho]

Here are some timings for several choices of num_lanes and for the number of threads:

• 32 lanes is the optimal. believe avx512 can store 16 floats in one register, but there are two processing heads So 32 lanes would keep these busy?
• performance starts to saturate at 16 threads. I think this is because the number of image lines processed by one thread then drops below 64, which may not be a coincidence. In the float part of my code this corresponds to 128 "lines of floats". Not sure if that is related and what the explanation might be. Also not all cores on this computer are equal
• there are some weird results: 2 threads is sometimes slower(!) than 1 thread, e.g. for 32 lanes
  and for 16 lanes. This is hard to explain.

The bottom graph is a zoom of the top one.

[kfjahnke]

Would it be possible to eliminate steps 3 and 5? These steps are used to glue process/convolve/process

I think it's expensive to split the pipeline into several images. If you do have to store to run the convolution, why not store to an array with larger pixels - e.g. an array of zimt_xel_t < float , 3 > instead of three arrays of float. The convolution code is dimension- and channel-count-agnostic, you can convolve the whole thing in one go and reload into a three- or six-channel pixel pipeline.

"convolve" is not a separate step, but rather a "dumb" combination of operations on neighbouring pixel values

I suggested you might start out with discrete coordinates instead of using specific getters. Transforming such discrete coordinates into the equivalent of what a linspace provides is simple - just one multiplication and one addition, which otherwise happen in the linspace getter anyway. If you have the discrete coordinates, you can run the convolution 'explicitly' as you go along, if you make sure that your algorithm stays within the bounds of the array. to guarantee that, you can either hold your data in slightly oversized arrays and extend the data into the 'frame' outside (a single-pixel frame is enough for 3X3 convolution) or you have to look at the coordinates and do special-casing for marginal coordinates. In any case you need some index manipulations and arithmetic, and some of this is bound to be redundant re. the optimized convolution code, but saving the store-and-reload may well be more expensive than that, you'll have to give it a shot. Without a 'safety margin', the index mathematics are quite simple masking operations, for the 3X3 averaging filter, with the central ccordinate given as as zimt::xel-t < int , 2 > 'c' and the array of px_t 'a' being sized wXh:

void eval ( const I & c , O & px )
{
auto cc = c ;
px_t sum = 0.0f ;

// first row
cc [ 0 ] = c [ 0 ] - 1 ;
cc [ 1 ] = c [ 1 ] - 1 ;
// guard against y < 0 or x < 0
cc [ 1 ] ( cc [ 1 ] < 0 ) = 0 ;
cc [ 0 ] ( cc [ 0 ] < 0 ) = 0 ;
sum += a [ cc ] ;
cc [ 0 ] = c [ 0 ] ;
sum += a [ cc ] ;
cc [ 0 ] += 1 ;
// guard against x >= w
cc [ 0 ] ( cc [ 0 ] >= w ) = w - 1 ;
sum += a [ cc ] ;

// second row
cc [ 0 ] = c [ 0 ] - 1 ;
cc [ 1 ] = c [ 1 ] ;
// guard against x < 0 ( y is in-range)
cc [ 0 ] ( cc [ 0 ] < 0 ) = 0 ;
sum += a [ cc ] ;
cc [ 0 ] = c [ 0 ] ;
sum += a [ cc ] ;
cc [ 0 ] += 1 ;
// guard against x >= w
cc [ 0 ] ( cc [ 0 ] >= w ) = w - 1 ;
sum += a [ cc ] ;

// third row
cc [ 0 ] = c [ 0 ] - 1 ;
cc [ 1 ] = c [ 1 ] + 1 ;
// guard against y >= h or x < 0
cc [ 1 ] ( cc [ 1 ] >= h ) = h - 1 ;
cc [ 0 ] ( cc [ 0 ] < 0 ) = 0 ;
sum += a [ cc ] ;
cc [ 0 ] = c [ 0 ] ;
sum += a [ cc ] ;
cc [ 0 ] += 1 ;
cc [ 0 ] ( cc [ 0 ] >= w ) = w - 1 ;
sum += a [ cc ] ;

px = sum * 1 / 9.0f ;
}

This is untested, but I think you'll get my drift. You'd have to set up an 'act' type functor taking the view in the c'tor to have a, w and h handy in 'eval'. The code looks like a lot, but since it will run as multithreaded SIMD, every cycle processes ncores X nlanes, and this may well be much faster than store-reload of entire arrays. I'd be curious to hear which variant wins out!

there are some weird results: 2 threads is sometimes slower(!) than 1 thread, e.g. for 32 lanes and for 16 lanes. This is hard to explain.

single-threading eliminates need for setting up threads in the first place and then for inter-thread synchronization. this may well explain your result. My experience is also that twice as many lanes as the platform's register size is optimal - one explanation is that some SIMD ops need set-up, and the second one following up can be processed faster. too many lanes lead to too much register pressure and performance drops again. two times the 'native' lane count is often the 'sweet spot'.

You write about avx512. Looking at your compilation parameters, I fail to see the relevant flags for avx512. maybe you just copied my flags for avx2? using native avx512 code may well boost your code's performance a good deal. You'll need at least -mavx512f, and possibly more, depending on your CPU - sorry I don't have that handy right now. Have to sign off now - other stuff to do - hope this helps you a bit further.

[kfjahnke]

If you build highway from source, it provides 'hwy_list_targets', a small utility which gives you the right flags for your CPU. e.g. on mine I get:

$ ~/src/highway/build/hwy_list_targets
Compiler: Clang 1901
Target attributes: sse2,ssse3,sse4.1,sse4.2,pclmul,aes,avx,avx2,bmi,bmi2,fma,f16c
Config: emu128:0 scalar:0 static:0 all_attain:0 is_test:0
Have: constexpr_lanes:1 runtime_dispatch:1 auxv:1 f16 type:1/ops1 bf16 type:1/ops0
Compiled HWY_TARGETS:   AVX3_SPR AVX3_ZEN4 AVX3_DL AVX3 AVX2 SSE4 SSSE3 SSE2
HWY_ATTAINABLE_TARGETS: AVX3_SPR AVX3_ZEN4 AVX3_DL AVX3 AVX2 SSE4 SSSE3 SSE2 EMU128
HWY_BASELINE_TARGETS:   SSE2 EMU128
HWY_STATIC_TARGET:      SSE2
HWY_BROKEN_TARGETS:     AVX10_2 SVE2_128 SVE_256 SVE2 SVE
HWY_DISABLED_TARGETS:  
CPU supports:           AVX2 SSE4 SSSE3 SSE2 EMU128 SCALAR
CPU: Intel(R) Core(TM) i5-4570 CPU @ 3.20GHz

see the 'target attribute' line, pick those flags.