bugGNU Octave - Bugs: bug #32980, OpenGL plotting backend uses...

 
 

bug #32980: OpenGL plotting backend uses single precision, not double

Submitter:  Michael Godfrey <godfrey>
Submitted:  Sun 03 Apr 2011 01:04:16 AM UTC
   
 
Category:  Plotting with OpenGL Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Confirmed Assigned to:  None
Originator Name:  Michael Godfrey Open/Closed:  * Open
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Fri 13 Oct 2023 01:01:11 PM UTC, comment #32: 

because the choice to limit opengl to single precision is outside the control of the octave developers, and there is not a trivial solution or workaround to expand precision, no there has not been any significant change to this behavior.

The comment #30 workaround idea was promising, but there does not appear to have been any further work in that direction. We are always welcoming of any additional interest, assistance, or code contributions to help improve Octave.

Nicholas Jankowski <nrjank>
Group Member
Fri 13 Oct 2023 09:15:03 AM UTC, comment #31: 

Are there news on this long-standing issue? Moreover it seems there are many bugs related to this one

Fabio <efferre79>
Tue 04 Apr 2023 08:07:38 AM UTC, comment #30: 

Is this the location where the precision is no longer controlled by Octave?

It looks like opengl_renderer in gl-render.cc contains a member of type _graphics_xform which in turn contains can do matrix / affine transformations via the member Matrix m_xform and do (non linear) scaling via the members scaler m_sx, m_sy, m_sz.

This graphics_xform object is used in forward direction to transform the graphic objects before passing to OpenGL (e.g. renderer::draw_line).

It is also (appears to be) used in reverse direction to untransform mouse position in screen coordinates back to data coordinates (particularly important for semilog and loglog plots).

Can this parameter (which is already a present feature and well tested) be used to solve this issue?

IMO, the transform should be handled by an axis as the axis object will need to apply the same transform to all its children if one of the children needs transformation.

Should a plot* function or axis-limit-setting-function detect large values (or large dynamic range or large axis limits) and set this graphic_transform object to non default values to shift and scale the data to within the range of single precision?

Is it worth pursuing along this line?

octavecontrib
Fri 04 Oct 2019 07:22:09 AM UTC, comment #29: 

@Roland: No, there is no patch and the only known workaround is to subtract large offsets from your data when it is option.

Pantxo Diribarne <pantxo>
Group Member
Thu 03 Oct 2019 07:51:38 AM UTC, comment #28: 

Hi there,

I reported yesterday bug #56983, which seems to be a duplicate of this issue.

I've read the comments here, but it's not clear to me if there is a patch that fixes (at least partially) the bug?

Roland Baudin <roland65>
Thu 25 Oct 2018 08:04:09 AM UTC, comment #27: 

A possible Workarount is to subtract the Year 2000 from the date like in the following code:


  %Get file and path
  [file,path]= uigetfile ();
  %Read file
  impDaten=csv2cell(strcat(path,file));

  %Read Zime and Date
  zeit = impDaten(2:end,2);
  Zeit=datenum(cell2mat(zeit),'dd/mm/yyyy HH:MM:SS');
  Zeit=Zeit-730485; %Big dates makes truble. Use the year 00xx.

  %in zeit is the original time / date
  %in Zeit is the substracted time / date.
  %in the plots only display the date in the following format:
  datestr(Zeit(1),'dd.mm.yy') %without time
  datestr(Zeit(1),'dd.mm.yy HH:MM:SS.FFF') %with time


M Kater <melody8801>
Mon 04 Jun 2018 12:39:36 PM UTC, comment #26: 

attached the plot...


Lars Kindermann <larskindermann>
Mon 04 Jun 2018 12:37:00 PM UTC, comment #25: 

I really want to bring this topic up again. At first glance it seems like a rare corner case problem, occurring only when working with pretty ill scaled data.

But have a look at this fairly simple example (plot example is attached):


#the last 24 hours in one minute steps
t = now-1:(1/1440):now;

#some example timeseries
x = -cos(2*pi*t) + rand(size(t))/10;

subplot(1,2,1)
plot(t,x)
datetick x
title(datestr(t(end)))

#get rid of the calendar day
t = t - fix(now);

subplot(1,2,2)
plot(t,x)
datetick x
title(datestr(t(end)))


In fact, serial dates as used by all time related functions in Octave are subject to this problem. Everybody plotting time stamped data with high resolution will be hit. This is daily routine in geophysics, meteorology, astronomy, stock charts, renewable energy, web traffic analysis, medicine, etc. And worse, there is no indicator that the plots are incorrect in fact. 

At the moment there are only two options to get around this, use     gnuplot (but then you are limited to just plain plots, no interactivity) or manually tune your code to get rid of the integer date offset. But then you loose many of the benefits of the datetick() function, like calendar dates on the time axis and any interactive handling or zooming becomes much more complicated.

Having looked at the code, I fear that fixing this fundamentally will be really a lot of work, requiring to offset and scale all plot data before sending to the plot engine and keep track of this for transforming correctly backwards when the plot returns interactive messages, or when adding or modifying data. Probably an additional abstraction layer between the octave plot logic and the graphics api needs to be inserted.

But in the meantime I strongly opt for at least throwing a warning when there is some indication that the scaling problem actually hits, not leaving the user unnoticed with an incorrect plot. This is also not trivial as it will be more a guess than some well defined calculation and should not be triggered too often by insignificant or false positives.

I am willing to spend some time on it, as this is one of the most annoying issues in Octave I am frequently confronted with. Does somebody have a suggestion where in the code such a test should ideally be located? This might be also a good start for inserting the rescaling code.   

Lars Kindermann <larskindermann>
Sun 20 Nov 2016 11:52:56 AM UTC, comment #24: 

As discussed below this is a feature of OpenGL. It is 32bit
and there is no plan to change this even though some graphics
hardware supports 64bit. And, there is no real benefit from
use of 64bit

So, as again discussed below, the change that is required is
for the data to be tested for within 32bit range and rescaled
if not. (As is done in Matlab.)

Michael Godfrey <godfrey>
Group Member
Sat 19 Nov 2016 01:56:31 PM UTC, comment #23: 

This behavior is still present in Octave 4.2.0.

Hartmut <hardy>
Fri 24 Feb 2012 09:20:13 PM UTC, comment #22: 

Yes, you're right!

Now that given a bit of thought, log10(double) can reliably represented as a single.

Ben Abbott <bpabbott>
Group Member
Fri 24 Feb 2012 09:08:41 PM UTC, comment #21: 

Aren't log scales already handled?  For example, the following already works properly for me:


graphics_toolkit fltk
x = y = logspace (-200, 200);
loglog (x, y);


John W. Eaton <jwe>
Group administrator
Fri 24 Feb 2012 08:38:34 PM UTC, comment #20: 

Yes. I think we'll need to avoid m-files to avoid triggering listeners.

My fist inclination was that we'd place everything in _init_fltk_.cc, but there's probably a nifty OO way of doing this ? (I'm a procedural guy).

We'll also need to consider log scale plots. They'll be a bit more tricky, but the same approach will work.

Ben Abbott <bpabbott>
Group Member
Fri 24 Feb 2012 08:14:48 PM UTC, comment #19: 

Mapping the data seems fine.  To get it to work I had to move the mapping of {x,y,z}tick ahead of the for loop because the call to set (h, 'ydata', ydata) updates the ylimits which are then wrong in the ymap function.

How were you thinking of implementing this?  I just reread the discussion and JWE seemed reluctant to put this in axes.m.  You mentioned putting equivalent code in _init_fltk_.cc?




Rik <rik5>
Group administrator
Fri 24 Feb 2012 05:54:14 PM UTC, comment #18: 

Any problem with mapping the {x,y,z}data ?


xtick = get (gca, "xtick");
xticklabel = get (gca, "xticklabel");
ytick = get (gca, "ytick");
yticklabel = get (gca, "yticklabel");
ztick = get (gca, "ztick");
zticklabel = get (gca, "zticklabel");

xmap = @(x) interp1 (xlim (), [-1 1], x, "linear", "extrap");
ymap = @(y) interp1 (ylim (), [-1 1], y, "linear", "extrap");
zmap = @(z) interp1 (zlim (), [-1 1], z, "linear", "extrap");

limits = [-1, 1, -1, 1, -1, 1];

for h = get (gca, "children")
  xdata = xmap (get (h, "xdata"));
  ydata = ymap (get (h, "ydata"));
  zdata = zmap (get (h, "zdata"));
  set (h, "xdata", xdata, "ydata", ydata, "zdata", zdata);
endfor

xtick = xmap (xtick);
ytick = ymap (ytick);
ztick = zmap (ztick);

set (gca, "xtick", xtick, "xticklabel", xticklabel)
set (gca, "ytick", ytick, "yticklabel", yticklabel)
set (gca, "ztick", ztick, "zticklabel", zticklabel)

axis (limits)


Ben Abbott <bpabbott>
Group Member
Fri 24 Feb 2012 05:15:29 PM UTC, comment #17: 

The approach is right but it is slightly more complicated.  The data needs to be shifted AND normalized.  For example, try this simple plot.


x = 1:10;
y = x * 1e40;
h = plot (x,y);


This produces a blank plot for me because 1e40 is beyond what is representable in a variable of class single.  Offsetting in this case doesn't help because even after offsetting the new y-axis range is from -5e40 to +5e40 and both of those values are still not representable.

I think you can get this to work by both shifting by the mean and normalizing by the mean as well.  I just did some command line work and dividing ydata by 1e40, ytick by 1e40, and ylim by 1e40 caused the correct graph to show up.



Rik <rik5>
Group administrator
Fri 24 Feb 2012 03:45:20 PM UTC, comment #16: 

It should be possible to shift (and normalized?) the xdata, ydata, zdata, cdata in _init_fltk_.cc.

A simple example from the command line that illustrates the problem for xdata is below.


graphics_toolkit fltk
close all
figure (1)
xdata = datenum (2012, 02, 22, 18, 00, 00) : 0.01 : now;
ydata = sin (xdata);
ydata(2:2:end) = NaN;
plot (xdata (~isnan (ydata)), ydata (~isnan (ydata)))
datetick ("x", 13)


By offsetting the xdata by its median the problem by be resolved. After running the commands above, add those below.


xtick = get (gca, "xtick");
xticklabel = get (gca, "xticklabel");
ytick = get (gca, "ytick");
yticklabel = get (gca, "yticklabel");
ztick = get (gca, "ztick");
zticklabel = get (gca, "zticklabel");

xoffset = mean (xlim ());
yoffset = mean (ylim ());
zoffset = mean (zlim ());

limits = axis ();
limits = limits - [xoffset*[1 1], yoffset*[1 1], zoffset*[1 1]](1:numel(limits));

for h = get (gca, "children")
  xdata = get (h, "xdata");
  ydata = get (h, "ydata");
  zdata = get (h, "zdata");
  xdata = xdata - xoffset;
  ydata = ydata - yoffset;
  zdata = zdata - zoffset;
  set (h, "xdata", xdata, "ydata", ydata, "zdata", zdata);
endfor

xtick = xtick - xoffset;
ytick = ytick - yoffset;
ztick = ztick - zoffset;
set (gca, "xtick", xtick, "xticklabel", xticklabel)
set (gca, "ytick", ytick, "yticklabel", yticklabel)
set (gca, "ztick", ztick, "zticklabel", zticklabel)

axis (limits)


Is this a proper approach ?

Ben Abbott <bpabbott>
Group Member
Sat 18 Feb 2012 05:35:11 AM UTC, comment #15: 

I suspect that Michael Goffioul knows how to do this.
It would be a good idea to test all data being processed
for plotting to ensure that the plot comes out right.
For OpenGL this means doing a transform that maps all the
data to within reals. gl-render.cc is a candidate place
to do this.   But it would also be useful to
check for very small data value differences so that the
plot scales and tic values can be made meaningful.
Now, an attempt to plot: plot(1:0.000001:0.000009)
does not work well, due to axis values being truncated.
Note that Matlab does not do this right either.


Michael Godfrey <godfrey>
Group Member
Sat 18 Feb 2012 01:26:20 AM UTC, comment #14: 

ping Any new thoughts/progress here?

Ben Abbott <bpabbott>
Group Member
Mon 25 Apr 2011 03:30:53 PM UTC, comment #13: 

John,

The code is specific in the sense that it tests for
the one current backend that restricts plot data to
single.  I looked in the fltk code, but could not find
a way to implement these tests there.  Also, it is
possible that there will be other backends that only
support single precision.

Currently the code uses realmax("double") for all cases
other than fltk, but this is pretty meaningless.  This
test can never, currently, fail.  I will change the code
to be fltk/OpenGL specific.

I thought about the rescaling and this should be possible,
but a lot more work in code that I am not familiar with.
It would be great if someone would do that.

In the meantime, it is not good that users may experience
mysteriously incorrect plots due to overflow of singles
in fltk/OpenGL.


Michael Godfrey <godfrey>
Group Member
Mon 25 Apr 2011 02:52:26 PM UTC, comment #12: 

The code in axes.m should be generic and not have special cases for specific toolkits.  So it seems to me that it would be better to fix this in the part of Octave that works directly with the fltk/opengl toolkit.  Since we have complete control over the axes labels in that code, it seems that all we would need to do is detect the condition there, scale the values to plot and then place axis labels appropriately.  That should allow us to handle the full range of double precision values even if OpenGL doesn't handle that internally.

John W. Eaton <jwe>
Group administrator
Sun 24 Apr 2011 12:27:31 AM UTC, comment #11: 

Here is another update.  I found that all the tests
can be (probably must be) made in axes.m.  So, this
fix only updates axes.m

It can be implemented in 2 forms, using a variable
do_reset.  If do_reset is non-zero, the code changes
the axis limits be, possibly, "in range."  Here
"in range" means that the plot will have correct
axes and display any data that fits in the axes.

It is likely that the code should be used with do_reset = 0.
Then the axes are not changed, but a warning is issued.
The resulting plot will be garbage.

I plan to do more testing, but this is much better than
the last try.

Attached are 2 files:
diff -c on axes.m and the axes.m file.  For testing
just place the axes.m file in the current directory.

Michael


(file #23293, file #23294)

Michael Godfrey <godfrey>
Group Member
Fri 22 Apr 2011 04:46:26 PM UTC, comment #10: 

Thanks Rik.

Since choice 1 would be at best unreliable I agree that
choice 2 is better, and easier.  It may be appropriate
to issue the warning at a value less than realmax("single")
and change the message to read something like "axis
and data values near realmax("single") or realmin("single")
may produce incorrect plots due to overflow/underflow"

And, your elimination of for loops is appropriate, of
course.  I try to do that too, but this code was a quick
test and copied from similar loops in axis.m...
Rounding to zero may be right, too.

Michael

Michael Godfrey <godfrey>
Group Member
Fri 22 Apr 2011 04:21:32 PM UTC, comment #9: 

My inclination is for option 2.  I subscribe to the philosophy that the user has more information, and is more intelligent, than the program and can make a better choice on how to resolve an error situation.  Thus, I think we can just issue a warning and produce an unreadable mess.  The user may then choose to scale their graph, perhaps choose the gnuplot toolkit, change the plot type to semilogy, etc.  What we have right now isn't workable because one can reach an error state, such as having one of the limits set to Inf, that is impossible to escape.

Also, regarding your patch.  I think it would be acceptable, when the value is below realmin, to round down to 0 rather than round up to realmin.

Finally, even though the for loop in the patch is over only a few items, I hate to see them in Octave.  Here is a sample way to avoid the loops by using indexing that you might consider or modify.


# Assuming the axis limits are in vector axlims = [xlims ylims zlims]
axlims(abs (axlims) < realmin) = 0;
axlims(axlims > realmax) = realmax;
axlims(axlims < -realmax) = -realmax;
bad_lims = axlims(1:2:end) >= axlims(2:2:end);
if (any (bad_lims))
  error (...);
endif


Rik <rik5>
Group administrator
Fri 22 Apr 2011 06:43:05 AM UTC, comment #8: 

Later:

The problem with sombrero was that the test in axes.m
was incorrect.  I now have a patch which works for
sombrero.

However, this has revealed a new problem.  The limit
for "single" in the case of the sombrero plot is
approximately 4e+35.  This is a lot smaller than for
the 2D plots where the limit is about 1.3e38.  Somewhere
in the code values must occur which are above
realmax("single") even though the peak of the sombrero
is only about 4e35.

So, there are at least 2 choices:

1. Find max and min values for the axes which will
   allow a "correct" plot in all known cases and
   reset the axes to these values.  This will produce
   a valid plot of all the data that fix within the
   reset axes.  A max/min value around +-e35 may work.

2. Just issue the warning and let the plot appear as
   an unreadable mess.  If this uses the exact
   realmax/min("single") there will be error cases that
   are not caught.

While waiting for comments, I will clean up the code.

Michael

Michael Godfrey <godfrey>
Group Member
Fri 22 Apr 2011 01:38:55 AM UTC, comment #7: 

Rik,

Thanks for trying the patch.  I am not surprised
that it failed the sombrero.  I had only tried 2D
plots mainly to try to check that the code was in the
right places.  It does check all three dimensions, but
something is wrong for 3D.

If you agree that this code is capable of doing the
right tests, I will improve it (also making the now
duplicated code into a function).

Michael

Michael Godfrey <godfrey>
Group Member
Fri 22 Apr 2011 01:11:26 AM UTC, comment #6: 

Unfortunately the patch does not seem to operate correctly on some plots.

Try the following example:


[x,y,z] = sombrero;
surf (x,y,z);


Without the patch, it produces the sombrero figure.  With the patch it produces a tiny fraction of a surface within the axis limits [0 1 0 1 0 1].  I can get it to display the whole figure with axis("auto"), but it indicates that something isn't right.

Rik <rik5>
Group administrator
Wed 20 Apr 2011 11:36:10 PM UTC, comment #5: 

Here is code that appears to address the scale
problem in fltk due to use of single. 

It was claimed that the code could be placed in
axis.m (around line 280) instead of in axes.m
but after putting the code in axis.m I found
that essentially the same code was needed in axes.m.
For example plot([1 10e40]) does not call axis.m
at all.  So, for now, there is duplicated code.
This code should be made into a called function.

I am posting this code for review.  There may well be
a better way to do this.  The code issues a warning and
rescales the plot.  This will prevent various forms
of inappropriate behavior.

After any comments I will try the cleanup.

Michael


(file #23281)

Michael Godfrey <godfrey>
Group Member
Tue 05 Apr 2011 06:30:15 PM UTC, comment #4: 

It appears that a place where the axes can be
checked for validity is after line 47 in axes.m.

Two questions for anyone (John?) who knows this
code:
1. Is this the right place?

2. Is it OK to put code which is conditional on
   the current backend at this point?

I notice that currently the command:
axis([0 Inf 0 200]) is processed without
complaint, but it obviously does not work and
trying to recover by axis([0 200 0 200]) which
is what the axes were before does not bring back
the plot correctly.

So, it seems a good idea to provide an axes validity
check in any case.  Currently the checks for gnuplot
(essentially valid range -max(double) to max(double))
and fltk/OpenGL are different, but both are needed.

Michael Godfrey <godfrey>
Group Member
Tue 05 Apr 2011 03:24:33 AM UTC, comment #3: 

I added a Caution note to the documentation about only single precision being available with the FLTK toolkit.  (http://hg.savannah.gnu.org/hgweb/octave/rev/562e40e3d5db)

Rik <rik5>
Group administrator
Tue 05 Apr 2011 02:44:51 AM UTC, comment #2: 

I checked that Matlab does not have this
problem.  But, I also noticed that Matlab
makes use of opengl.

So, there may be a way to get opengl to
work correctly.

Ideas?

Michael Godfrey <godfrey>
Group Member
Sun 03 Apr 2011 05:39:03 AM UTC, comment #1: 

As Shai pointed out, the problem seems to be
in OpenGL which accepts double, but converts
to single.  This is not really made clear in
the web pages I have read, but it is mentioned that
they the work in the precision available in the
GPU hardware which is typically 32bit.  Newer
hardware provides 64bit, so there is some
indication that a future OpenGL will provide
64bit.  In the meantime, this report could be
dealt with by providing a warning that the plot
cannot be rendered correctly due to data out of
range of representable values. 

The current behavior is not good since it could
lead to plots which are seriously, but not always
obviously, in error.

Michael

Michael Godfrey <godfrey>
Group Member
Sun 03 Apr 2011 01:04:16 AM UTC, original submission:  

A simple example is:

x = [1e38:1e38:1e40];
plot(x)


This problem does not exist for gnuplot, so the error
must be in fltk-dependent code.

Michael Godfrey <godfrey>
Group Member

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #44290:  32BitDateTimePlot.png added by larskindermann (19KiB - image/png - Example of plotting a timeseries on a serial date axis)
file #23293:  axes.m added by godfrey (4KiB - application/vnd.wolfram.mathematica.package)
file #23294:  axes_2.diff added by godfrey (3KiB - application/octet-stream)
file #23281:  fltk_axes.diff added by godfrey (4KiB - application/octet-stream)

 

Depends on the following items: None found

Items that depend on this one

Digest:
   bug dependencies.

 

Carbon-Copy List
  • -email is unavailable- added by nrjank (Posted a comment)
  • -email is unavailable- added by efferre79 (Posted a comment)
  • -email is unavailable- added by octavecontrib (Posted a comment)
  • -email is unavailable- added by pantxo (Posted a comment)
  • -email is unavailable- added by roland65 (Posted a comment)
  • -email is unavailable- added by melody8801 (Posted a comment)
  • -email is unavailable- added by larskindermann (Posted a comment)
  • -email is unavailable- added by hardy (Posted a comment)
  • -email is unavailable- added by jordigh (Updated the item)
  • -email is unavailable- added by bpabbott (Posted a comment)
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by godfrey (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 25 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2023-02-20 pantxo Dependencies- bugs #63820 is dependent
    2022-10-18 mmuetzel Dependencies- bugs #63234 is dependent
    2021-12-10 mmuetzel Dependencies- bugs #61658 is dependent
    2021-08-29 mmuetzel Dependencies- bugs #61087 is dependent
    2021-04-06 mmuetzel Dependencies- bugs #60350 is dependent
    2021-03-18 mmuetzel Dependencies- bugs #60252 is dependent
    2021-03-16 pantxo Dependencies- bugs #60243 is dependent
    2020-12-15 mmuetzel Dependencies- bugs #59690 is dependent
    2020-12-10 mmuetzel Dependencies- bugs #59650 is dependent
    2020-05-06 mmuetzel Dependencies- bugs #58305 is dependent
    2019-10-04 pantxo StatusNone Confirmed
        Operating SystemGNU/Linux Any
    2019-10-03 pantxo Dependencies- bugs #56983 is dependent
    2018-06-04 larskindermann Attached File- Added 32BitDateTimePlot.png, #44290
    2018-03-23 mtmiller Dependencies- bugs #53436 is dependent
    2016-04-26 mtmiller Dependencies- bugs #47788 is dependent
    2015-11-25 mtmiller Dependencies- bugs #46529 is dependent
    2015-08-31 mtmiller Dependencies- bugs #45851 is dependent
    2015-07-16 pantxo Dependencies- bugs #45541 is dependent
    2015-06-07 mtmiller Dependencies- bugs #45271 is dependent
    2015-04-30 mtmiller Dependencies- bugs #44963 is dependent
    2014-12-26 rik5 Dependencies- bugs #33748 is dependent
    2013-03-14 jordigh Summaryfor backend FLTK, real is used for axes. double should be used OpenGL plotting backend uses single precision, not double
    2013-03-14 jordigh Dependencies- bugs #38515 is dependent
    2013-01-10 jordigh Dependencies- bugs #38054 is dependent

    Back to the top

    Powered by Savane 3.13-cf05.
    Corresponding source code