Skip to main content

memory - Plotting large datasets


I have a lot of data (e.g., 64x8192), to plot and it consumes too much memory (more than 2 GB) causing Mathematica to return "no memory available".


For example,


lst = Table[Sin[x^2 + y^2], {x, 2^7}, {y, 2^13}];
ListDensityPlot[lst]

Is there any way to cope with it?



Answer



Scroll to the end for the latest edit.


If you really want all the points in a large data set, I'd use a more basic display method:



Let's say the matrix of data is m, then do this (assuming the entries m[[i, j]] are real numbers):


Image[Rescale[m]]

This will make one pixel for each data point without trying to do any interpolation as is the default in density plots.


Edit


As suggested in the comment, here is a comparison of the timings:


(* A test matrix with random entries between 0 and 1: *)
m = RandomReal[1, {300, 300}];
Timing@Image[m]


fast


Timing@ListDensityPlot[m, MaxPlotPoints -> 1000]

slow


So the ListDensityPlot is several orders of magnitude slower than my Image suggestion. This conclusion seems to be almost independent of the choice for MaxPlotPoints.


With the Image approach, it is no problem to create extremely large images (say $5000\times 5000$). I would not advise using ListDensityPlot for anything like this.


If you want a colored output, just use Colorize as follows:


m = RandomReal[1, {300, 300}];

Timing[plot = Colorize[Image[m], ColorFunction -> "LakeColors"]]


colorized


Edit 2


The final question would be, how to re-create the frame and ticks that you get from ListDensityPlot and its relatives. Let's assume that the range of x and y values corresponding to the matrix dimensions of m are $x \in [0, 2]$ and $y\in[0, 3]$, respectively. Then you get the standard display for the plot we just created as follows:


xmin = 0;
ymin = 0;
xmax = 2;
ymax = 3;
Graphics[Inset[
  Show[plot, AspectRatio -> Full], 

  {xmin, ymin}, 
  {0, 0}, 
  {xmax - xmin, ymax - ymin}], 
 PlotRange -> {{xmin, xmax}, {ymin, ymax}}, 
 Frame -> True]

framed plot


The trick here is to align the bottom left corner {0,0} of the Image to the lower coordinate limits of the PlotRange I chose for the Graphics, and then (using the fourth argument of Inset) to stretch the image to fit exactly into that PlotRange. The content is only stretchable this way if you show it with the option AspectRatio -> Full.


Edit 3


For completeness, here is the example matrix given in the question, constructed in a slightly modified way to make its entries lie between 0 and 1. Here, Rescale turns out to be very slow when applied to an array of exact numbers, so I speed things up by adding N[] during the construction of the matrix. Thanks to MichaelE2 for pointing this out in the comment. The matrix construction uses Outer and automatically threads Sin over lists in its argument.



In this construction, the x coordinate is represented by Range[2^7] and the y coordinate by Range[2^13] - if you switch these then the plot will have the axes switched. The Reverse below is needed to have the vertical axis pointing up:


s = Sin[Outer[Plus, N[Reverse[Range[2^13]]^2], N[Range[2^7]^2]] ];
plot = Colorize[Image[s], ColorFunction -> "LakeColors"];

xmin = 1;
ymin = 1;
xmax = 2^7;
ymax = 2^13;
Graphics[Inset[
  Show[plot, AspectRatio -> Full], {xmin, ymin}, {0, 0}, {xmax - xmin, 

   ymax - ymin}], PlotRange -> {{xmin, xmax}, {ymin, ymax}}, 
 Frame -> True, AspectRatio -> 1]

sin 2d


The point of this example is to show how much more detail we get with this approach, compared to the slower ListDensityPlot.


Alternative: ArrayPlot


I mentioned using Rescale earlier because the arguments to Image have to be in the unit interval to be mapped onto the grayscale range from white to black. This grayscale can later be colorized too. It's important to keep this fact in mind: matrix entries that fall outside the interval $[0, 1]$ will not lie in the range where the grayscale varies and will therefore appear at the same extreme grayscale (or color afte using Colorize).


Once we add the rescaling and the frame ticks, the Image approach de-facto converges to the built-in function ArrayPlot. This is another easy-to-use way of getting a fast plot of a large matrix. You would use it like this: ArrayPlot[m, ColorFunction -> "LakeColors", DataReversed -> True]. I just tried this with the last example matrix, and it works nicely.


One more observation


According to this linked answer, when using ColorFunction, ArrayPlot isn't as fast as the direct creation of a Raster image. The solution proposed there is not Image as I used here, but instead something like this: Graphics[Raster[Rescale[...], ColorFunction -> "LakeColors"]]



To try this with my frame construction, you would have to write it like this (where s is the same as defined in Edit 3):


xmin = 1;
ymin = 1;
xmax = 2^7;
ymax = 2^13;
Graphics[Inset[
  Show[
   Graphics[Raster[s, ColorFunction -> "LakeColors"]], 
   AspectRatio -> Full, 
   PlotRangePadding -> 0], 

  {xmin, ymin}, 
  {0, 0},
  {xmax - xmin, ymax - ymin}
 ], 
 PlotRange -> {{xmin, xmax}, {ymin, ymax}}, 
 Frame -> True, AspectRatio -> 1
]

When I wrap this in Timing or AbsoluteTiming, I indeed get even lower values. However, I don't really believe those values because the actual time it takes for the plot to appear in the notebook is noticeably slower with this last approach. Nevertheless, I think it's worth including this here.


More timing results



To test these display methods with the Mandelbrot set of this linked answer, I did a different timing measurement based on recording the absolute time before and after the image/graphic has been displayed. Here is the code to prepare the pictures:


mandelComp = 
  Compile[{{c, _Complex}}, 
   Module[{num = 1}, 
    FixedPoint[(num++; #^2 + c) &, 0, 99, 
     SameTest -> (Re[#]^2 + Im[#]^2 >= 4 &)]; num], 
   CompilationTarget -> "C", RuntimeAttributes -> {Listable}, 
   Parallelization -> True];

data = ParallelTable[

   a + I b, {a, -.715, -.61, .0001}, {b, -.5, -.4, .0001}];

rc = Rescale[mandelComp[data]];

Now the timing tests (not showing the images that are generated):


a1 = AbsoluteTime[]; Graphics[
 Raster[rc, ColorFunction -> "TemperatureMap"], ImageSize -> 800, 
 PlotRangePadding -> 0]

a2 = AbsoluteTime[]


a2 - a1

(* ==> 6.878487 *)

This is pretty slow. Compare this to the Image approach:


a1 = AbsoluteTime[]; Colorize[Image[rc, ImageSize -> 800], 
 ColorFunction -> "TemperatureMap"]

a2 = AbsoluteTime[]


a2 - a1

(* ==> 1.119959 *)

This is much faster. The ArrayPlot method is twice as slow as the last result, but still faster than the Raster method:


a1 = AbsoluteTime[];
ArrayPlot[rc, ImageSize -> 800, ColorFunction -> "TemperatureMap"]

a2 = AbsoluteTime[]


a2 - a1

(* ==> 2.089126 *)

Therefore, my conclusion is: when in doubt, go with Image.


Edit (19 June 2014)


It seems useful to wrap the conclusion of this answer in a general-purpose function that works similarly to ListDensityPlot but is based on directly creating an Image from the data:


Clear[listPixelPlot] ;
Options[listPixelPlot] = {Frame -> True, FrameLabel -> {}, 

   ColorFunction -> "LakeColors" ,
   AspectRatio -> 1, DataRange -> All , ImageSize -> Automatic, 
   PlotRange -> All, PlotRangeClipping -> True};

listPixelPlot[array_, OptionsPattern[]] := 
 Module [{field, plot, xmin, xmax, ymin, ymax, 
   dataRange = OptionValue[DataRange], 
   plotRange = 
    Replace[OptionValue[PlotRange], 
     Full | All | Automatic -> {Automatic, Automatic}], range},

  field = Reverse@Transpose[array] ;
  {{xmin, xmax}, {ymin, ymax}} = range = If[
     MatrixQ[dataRange] ,
     Reverse[dataRange],
     Map [ {0 , #} &, Dimensions[array] ]
     ] ;
  plot = Colorize[Image[Rescale[N[field] ] ] , 
    ColorFunction -> OptionValue[ColorFunction]] ;
  Graphics[
   Inset[

    Show [plot, AspectRatio -> Full] , {xmin, ymin} , {0, 
     0} , {xmax - xmin, ymax - ymin}],
   PlotRange -> 
    Table[plotRange[[i]] /. Full | All | Automatic -> range[[i]], {i, 
      2}],
   Sequence @@ 
    Map[# ->  OptionValue[#] & , {Frame, FrameLabel, ImageSize, 
      AspectRatio, PlotRangeClipping}]]
  ]


This function takes some of the same options that ListDensityPlot or ArrayPlot accept (see Options[). With the Mandelbrot test data rc from above, you would call this function (and test its timing) as follows:


a1 = AbsoluteTime[];
listPixelPlot[rc, ColorFunction -> "TemperatureMap", ImageSize -> 800]

a2 = AbsoluteTime[];

a2 - a1

(* ==> 1.679611 *)


This is only slightly slower than the fastest result above. The reason is simply that in this function the argument array is always processed by Rescale before converting it to an Image. To make sure Rescale sees a packed array of machine numbers, the core of the function reads Image[Rescale[N[field] ] ]. A lot of the code deals merely with option handling. It can be extended to provide additional options.


Comments

Post a Comment

Popular posts from this blog

plotting - Filling between two spheres in SphericalPlot3D

Manipulate[ SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, Mesh -> None, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], {n, 0, 1}] I cant' seem to be able to make a filling between two spheres. I've already tried the obvious Filling -> {1 -> {2}} but Mathematica doesn't seem to like that option. Is there any easy way around this or ... Answer There is no built-in filling in SphericalPlot3D . One option is to use ParametricPlot3D to draw the surfaces between the two shells: Manipulate[ Show[SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], ParametricPlot3D[{ r {Sin[t] Cos[1.5 Pi], Sin[t] Sin[1.5 Pi], Cos[t]}, r {Sin[t] Cos[0 Pi], Sin[t] Sin[0 Pi], Cos[t]}}, {r, 1, 2 - n}, {t, 0, Pi}, PlotStyle -> Yellow, Mesh -> {2, 15}]], {n, 0, 1}]
Post a Comment
Read more

plotting - Plot 4D data with color as 4th dimension

I have a list of 4D data (x position, y position, amplitude, wavelength). I want to plot x, y, and amplitude on a 3D plot and have the color of the points correspond to the wavelength. I have seen many examples using functions to define color but my wavelength cannot be expressed by an analytic function. Is there a simple way to do this? Answer Here a another possible way to visualize 4D data: data = Flatten[Table[{x, y, x^2 + y^2, Sin[x - y]}, {x, -Pi, Pi,Pi/10}, {y,-Pi,Pi, Pi/10}], 1]; You can use the function Point along with VertexColors . Now the points are places using the first three elements and the color is determined by the fourth. In this case I used Hue, but you can use whatever you prefer. Graphics3D[ Point[data[[All, 1 ;; 3]], VertexColors -> Hue /@ data[[All, 4]]], Axes -> True, BoxRatios -> {1, 1, 1/GoldenRatio}]
Post a Comment
Read more

plotting - Adding a thick curve to a regionplot

Suppose we have the following simple RegionPlot: f[x_] := 1 - x^2 g[x_] := 1 - 0.5 x^2 RegionPlot[{y < f[x], f[x] < y < g[x], y > g[x]}, {x, 0, 2}, {y, 0, 2}] Now I'm trying to change the curve defined by $y=g[x]$ into a thick black curve, while leaving all other boundaries in the plot unchanged. I've tried adding the region $y=g[x]$ and playing with the plotstyle, which didn't work, and I've tried BoundaryStyle, which changed all the boundaries in the plot. Now I'm kinda out of ideas... Any help would be appreciated! Answer With f[x_] := 1 - x^2 g[x_] := 1 - 0.5 x^2 You can use Epilog to add the thick line: RegionPlot[{y < f[x], f[x] < y < g[x], y > g[x]}, {x, 0, 2}, {y, 0, 2}, PlotPoints -> 50, Epilog -> (Plot[g[x], {x, 0, 2}, PlotStyle -> {Black, Thick}][[1]]), PlotStyle -> {Directive[Yellow, Opacity[0.4]], Directive[Pink, Opacity[0.4]],
Post a Comment
Read more