7
$\begingroup$

I am looking for an implementation of MATLAB's numgrid function, in particular the "B" mode.

For example, I want to get the matrix corresponding to the following matrix in MATLAB:

numgrid('B', 80)
$\endgroup$
7
  • 6
    $\begingroup$ 1. More information is needed. The linked page doesn't involve enough details about numgrid (and it's a bit surprising to me there doesn't seem to be a stand-alone document page for numgrid! ), according to this page, 'B' stands for exterior of a "Butterfly", then how is the butterfly defined? 2. Somewhat related: mathematica.stackexchange.com/a/43039/1871 3. Are you aware that it's not good idea to use For and AppendTo in Mathematica?: mathematica.stackexchange.com/q/134609/1871 $\endgroup$
    – xzczd
    Jul 8, 2022 at 6:18
  • 4
    $\begingroup$ The laziness is astounding. You include a couple of links I'm supposed to follow to figure out what this function does (I did not follow the links) when you should have described plainly what it does. Then you include a screenshot from some textbook that is - I guess - related to your problem, with a couple of pieces highlighted that in no way describe your problem. $\endgroup$
    – Jason B.
    Jul 8, 2022 at 13:06
  • $\begingroup$ @JasonB. sorry if this wasted your time. I sometimes start questions as placeholders for problems whose solutions would benefit the wider community, and come back later to fill in the results. The community can be quite good at mind reading :) $\endgroup$ Jul 8, 2022 at 17:41
  • $\begingroup$ Well, I don't think it's good idea. Surely the community is kind. Finally there come 3 answers for an (currently) unclear question! But kindness is non-renewable. $\endgroup$
    – xzczd
    Jul 8, 2022 at 18:07
  • 4
    $\begingroup$ It would be great if someone rewrote the question now that it has answers, future readers would surely appreciate clarity. $\endgroup$
    – Jason B.
    Jul 9, 2022 at 0:02

4 Answers 4

21
$\begingroup$

The source code of some of the MATLAB functions is available to its users. In this example, if you type edit numgrid, its file will be open, and to my surprise, it wasn't very complicated, so I decided to rewrite it with Mathematica's built-in functions (I'm sure there is room for improvement).

Update

Thanks to @chyanog's suggestion to remove Thread in building SparseArray, it performs almost 3 times faster.

ClearAll[NumGridB];

NumGridB[n_Integer] :=
 Block[{x = ConstantArray[Subdivide[-1., 1, n - 1], n], y, t, r, temp,
    Indeterminate = 0.},

  y = Transpose[x[[All, -1 ;; 1 ;; -1]]];

  t = Quiet@ArcTan[x, y];
  r = Sqrt[x^2 + y^2];

  temp = SparseArray[
   1 - Unitize[
     UnitStep[
       r - (Sin[2 t] + 0.2*Sin[8 t])] + (1 -
        UnitStep[Abs[x] - 1]) + (1 - UnitStep[Abs[y] - 1]) - 3]];

  SparseArray[
 SortBy[temp["ExplicitPositions"], Last] ->
  Range[temp["ExplicitLength"]], {n, n}, 0]
  ]

Notes:

  • I used the same variable names as MATLAB, except temp, so you can follow its logic.
  • In Block scope, I used Indeterminate=0. with Quiet because in Mathematica ArcTan[0, 0] returns Indeterminate with a Message, but in MATLAB, atan2 simply return 0 - @chyanog suggests using Arg[x+I*y], but since it's a little slower, I decided not to use it.

Example:

Mathematica:

ArrayPlot[NumGridB[7], Frame -> False]

MATLAB:

imagesc(numgrid('B',7))

comparison

Higher number (n=500):

comparison

Speed

n MATLAB (second - tic;toc;) Mathematica (second - RepeatedTiming)
500 0.006 0.18 0.044
5000 0.55 17.9 4.75

Without touching *Compile, for not large n, its timing is reasonable. If you did improve the speed, feel free to post it here.

Update 2

Other answers show clearly compiling will further boost the performance. I decide to use the wolfram-library-link-rs crate to use Rust for building a LibraryLink library.

You can get the code and compiled library from this GitHub page.

At first glance, it wasn't much faster than others (clearly slower than @chyanog solution!), but if you open Task Manager while running MATLAB numgrid, you'll notice is not a single-threaded command, So this became an opportunity to leverage Rust's power in multithreading to set a new record.

After downloading the .dll file from the Release Section of the GitHub page, execute the following command and you're good to go:

(* Single-threaded *)
NumGridBCompiled =
 LibraryFunctionLoad[
  "C:\\numgrid.dll", "numgrid_b", {Integer},
  LibraryDataType["NumericArray", "UnsignedInteger32", 2]]


(* Multi-threaded *)
NumGridBParallelCompiled =
 LibraryFunctionLoad[
  "C:\\numgrid.dll",
  "numgrid_b_parallel", {Integer},
  LibraryDataType["NumericArray", "UnsignedInteger32", 2]]

Benchmark

Speed

benchmark comparison

  • For "Ben Izd", AbsoluteTiming was used.

Table (all the timings are in seconds):

n MATLAB NumGridBParallelCompiled NumGridBCompiled chyanog xzczd Ben Izd
1,000 0.023 0.0039 0.033 0.016 0.061 0.16
5,000 0.61 0.25 0.84 0.63 1.49 4.72
10,000 2.4 1.47 3.39 2.06 5.98 46.1

Memory

Although I'm not fully aware of Mathematica accuracy on LibraryLink calls, here is the plot:

memory consumption

Table (all the sizes are in bytes):

n MATLAB NumGridBParallelCompiled NumGridBCompiled chyanog xzczd Ben Izd
1,000 7,864,000 4,000,256 4,000,256 24,020,064 16,000,432 95,388,352
5,000 195,764,000 100,000,256 100,000,256 600,085,664 400,000,432 2,390,195,048
10,000 782,712,000 400,000,256 400,000,256 2,400,163,872 1,600,000,432 9,563,573,936

Notes:

  • MATLAB timings were based on multiple tic;toc; and I try a slightly different size and the timings were around the same number.
  • MATLAB memory usage was based on its profiler with profile('-memory','on'), and its Peak Memory (values were in Kb, so it was converted accordingly)
  • MATLAB version is 9.12.0 (2022b update 2), Mathematica 13.1.0, and system is AMD Ryzen 3950X
$\endgroup$
6
  • 4
    $\begingroup$ Thread is unnecessary here and inefficient, you can also use Arg[x+I*y] instead of Arctan, it will not have a warnings. $\endgroup$
    – chyanog
    Jul 8, 2022 at 8:38
  • 1
    $\begingroup$ @chyanog I really appreciate your comment (~3x times faster). Thank you. $\endgroup$
    – Ben Izd
    Jul 8, 2022 at 9:05
  • 2
    $\begingroup$ As to the MATLAB timing, have you taken the influence of memory pre-allocation into account? Are the timings obtained on first run? $\endgroup$
    – xzczd
    Jul 8, 2022 at 9:50
  • 1
    $\begingroup$ Amazing excavation, thanks! $\endgroup$ Jul 8, 2022 at 13:47
  • 1
    $\begingroup$ @xzczd I hope I answered your questions in update 2. $\endgroup$
    – Ben Izd
    Jul 9, 2022 at 6:09
12
$\begingroup$

After searching for a while I found the source code of numgrid here. So, here's my trial with FunctionCompile:

bfnumgrid =
   FunctionCompile@
    Function[Typed[n, "Real64"], 
     With[{atan2 = {y, x} |-> If[x == 0 && y == 0, 0., ArcTan[x, y]], 
           sqrt = Sqrt, sin = Sin}, 
       Module[{t, r, i = 1},
         Table[t = atan2[y, x];
               r = sqrt@(x^2 + y^2); 
               If[r >= (sin@(2*t) + .2*sin@(8*t)) && -1 < x < 1 && -1 < y < 1, 
                  i++, 0], 
               {x, -1, 1, 2/(n - 1)}, 
               {y, 1, -1, -2/(n - 1)}]\[Transpose]]]]; // AbsoluteTiming
(* {9.24076, Null} *)

bfshape[5000]; // AbsoluteTiming
(* {1.57289, Null} *)

Be careful with the subtle differences between MATLAB and Mathematica e.g. atan2 v.s. ArcTan, row-major v.s. column-major, etc. I don't choose the Arg trick shown by chyanog because it turns out to be a little slower compared with the If[…] work-around in my current implementation.

For comparison, Ben lzd's solution takes about 6.55 s on my laptop for n == 5000.

Visualization:

ArrayPlot[bfnumgrid[500], ColorFunction -> "AvocadoColors"]

enter image description here

$\endgroup$
12
$\begingroup$
Clear[cf, numGridB];
cf = Compile[{{n, _Integer}, x, {range, _Real, 1}},
   Table[With[{r = Sqrt[x^2 + y^2], t = Arg[x + I y]},
     Boole[r >= Sin[2 t] + 0.2*Sin[8 t] && Abs@x < 1 && Abs@y < 1]], {y, range}
    ], RuntimeAttributes -> {Listable}, RuntimeOptions -> "Speed", 
       CompilationTarget -> "C"
   ];

numGridB[n_] :=
  With[{range = Range[-1., 1, 2/(n - 1)]},
    Transpose@ArrayReshape[# Accumulate@#, {n, n}] &@Flatten@cf[n, -range, range]
   ];

numGridB[10] // MatrixForm
numGridB[5000]; // RepeatedTiming

$Version

13.1.0 for Microsoft Windows (64-bit) (June 16, 2022)

For GCC (mingw-w64) compiler, there is no significant difference between -O2, -O3, -Ofast flag

{0.840892, Null}

For VC compiler, /O2 flag

{0.596376, Null}

For Clang compiler, same flag as GCC

{0.561526, Null}

Compared with MATLAB
MATLAB screenshot

Another implementation:

Clear[numGridB2];
numGridB2 = Compile[{{n, _Integer}},
   Module[{x, y, t, r, bin, tol = 2.*^-16},
    y = ConstantArray[Range[-1., 1, 2/(n - 1)], n];
    x = -Transpose[y];
    t = ArcTan[x + tol, y];
    r = Sqrt[x^2 + y^2];
    bin = Flatten@BitAnd[UnitStep[r - (Sin[2 t] + 0.2*Sin[8 t])], 
       UnitStep[1 - tol - Abs@x], UnitStep[1 - tol - Abs@y]];
    Transpose@Partition[bin Accumulate@bin, n]
    ]];

numGridB2[10] // Grid
numGridB2[5000]; // RepeatedTiming

{0.980508, Null}

$\endgroup$
2
  • $\begingroup$ So this is faster than matlab in MSVC case? Why gcc is much slower? $\endgroup$ Jul 9, 2022 at 7:35
  • $\begingroup$ What compiler flags were used? Was it compiled with different optimisation flags? $\endgroup$ Jul 9, 2022 at 20:14
7
$\begingroup$

This is a little slower than Ben's version, but I thought I'd share this anyway since the "ExplicitLength" and "ExplicitPositions" properties are not known to earlier versions, and I wanted to show off the use of CoordinateBoundsArray[]:

numGridB[n_Integer?Positive] := Module[{pos, t, temp, x, y},
   {x, y} = Transpose[CoordinateBoundsArray[{{-1., 1.}, {-1., 1.}},
                                            {Into[n - 1], Into[n - 1]}],
                      {3, 2, 1}];
   y = Reverse[y]; t = ArcTan[x + $MachineEpsilon, y];
   temp = SparseArray[1 - Unitize[UnitStep[Sqrt[x^2 + y^2] -
                                           (Sin[2 t] + 0.2*Sin[8 t])] -
                      UnitStep[Abs[x] - 1] - UnitStep[Abs[y] - 1] - 1]];
   pos = temp["NonzeroPositions"];
   SparseArray[SortBy[pos, Last] -> Range[Length[pos]], {n, n}, 0]]

For visualization purposes, let me include this as well:

(* https://tpfto.wordpress.com/2019/04/17/faking-a-birds-colors-and-miscellanea/ *)
birdie[x_?NumericQ] /; 0 <= x <= 1 := 
Blend[{RGBColor[0.242, 0.15, 0.66], RGBColor[0.27, 0.2135, 0.83], 
       RGBColor[0.281, 0.298, 0.939], RGBColor[0.271, 0.385, 0.99], 
       RGBColor[0.202, 0.479, 0.99], RGBColor[0.172, 0.564, 0.94], 
       RGBColor[0.13, 0.639, 0.892], RGBColor[0.048, 0.7025, 0.827], 
       RGBColor[0.07, 0.746, 0.726], RGBColor[0.203, 0.78, 0.61], 
       RGBColor[0.36, 0.8, 0.463], RGBColor[0.58, 0.79, 0.29], 
       RGBColor[0.786, 0.757, 0.16], RGBColor[0.946, 0.7285, 0.217], 
       RGBColor[0.995, 0.79, 0.2], RGBColor[0.961, 0.891, 0.153], 
       RGBColor[0.977, 0.984, 0.08]}, x]

and then, we have

ArrayPlot[numGridB[7], ColorFunction -> birdie, Frame -> False]

numgrid visualization

ArrayPlot[numGridB[500], ColorFunction -> birdie, Frame -> False]

numgrid visualization, higher resolution


Edit 7/10/2022

The stuff above was written a little before I went to bed. Now with the benefit of getting sleep and seeing Ben's comment and chyanog's new answer, I came up with the following:

numGridB[n_Integer?Positive] := Module[{pos, sgn2, temp, zz},
   zz = {1, I} . Transpose[CoordinateBoundsArray[{{-1., 1.}, {1., -1.}},
                                                 {Into[n - 1], Into[n - 1]}], {3, 2, 1}];
   sgn2 = Sign[zz]^2;
   temp = SparseArray[BitOr[1 - UnitStep[Abs[zz] - Im[sgn2 (1 + 0.2 sgn2^3)]],
                            UnitStep[Abs[Re[zz]] - 1], UnitStep[Abs[Im[zz]] - 1]],
                      Automatic, 1];
   pos = temp["NonzeroPositions"];
   SparseArray[SortBy[pos, Last] -> Range[Length[pos]], {n, n}, 0]]

How this version works is an exercise I will leave for the interested.

$\endgroup$
3
  • 1
    $\begingroup$ No matter how hard we push, there is always room for improvement. I really appreciate your answer. $\endgroup$
    – Ben Izd
    Jul 10, 2022 at 7:55
  • 2
    $\begingroup$ (one small point, wouldn't {1., -1.} be better than Reverse[y]?) $\endgroup$
    – Ben Izd
    Jul 10, 2022 at 7:55
  • 1
    $\begingroup$ Thanks @Ben, that certainly streamlined the new version I came up with a little after waking up. ;) $\endgroup$ Jul 10, 2022 at 13:01

Not the answer you're looking for? Browse other questions tagged or ask your own question.