bugGNU Octave - Bugs: bug #59850, uniquetol missing, so I...

 
 

bug #59850: uniquetol missing, so I implemented it

Submitted by:  None
Submitted on:  Sun 10 Jan 2021 12:10:34 PM UTC  
 
Category:  Octave Function Severity:  1 - Wish
Priority:  3 - Low Item Group:  Feature Request
Status:  In Progress Assigned to:  None
Originator Name:  Originator Email:  -email is unavailable-
Open/Closed:  Open Release:  dev
Operating System:  Any

Add a New Comment (Rich Markup)
   

You are not logged in

Please log in, so followups can be emailed to you.

 

( Jump to the original submission)

Wed 24 Feb 2021 11:05:56 AM UTC, comment #25: 

My method for creating this patch:

  • Copied the file to a different directory.
  • Updated repository.
  • Striped all local changes.
  • Copied back to directory.
  • Created a new patch.

Hopefully, it works. The "byrows", false should be fast without the x = 1:.045:3 bug.

(file #50913)

Anonymous
Tue 23 Feb 2021 07:38:22 PM UTC, comment #24: 

I added some general tips to the wiki for working with TortoiseHg. [1]
Maybe some of those could be helpful to produce a consolidated patch with the changes that applies on a newer tip.

[1]: https://wiki.octave.org/Mercurial#Tips_for_working_with_TortoiseHg

Markus Mützel <mmuetzel>
Project Member
Wed 10 Feb 2021 09:54:08 PM UTC, comment #23: 

The m file is the latest version. The first patch was the first attempt at a fast version. I modified the code of the file the first patch was created from to make the second patch to fix the bug it introduced. I will need instructions on how to create a patch from the original or current code.

Anonymous
Wed 10 Feb 2021 04:39:26 PM UTC, comment #22: 

This is getting confusing.
Could you please upload a patch with the changes you propose that applies on a current head of the default branch?

Markus Mützel <mmuetzel>
Project Member
Tue 26 Jan 2021 04:23:02 AM UTC, comment #21: 

Yes

Attached the m file.

(file #50787)

Anonymous
Tue 26 Jan 2021 04:16:35 AM UTC, comment #20: 

The latest patch didn't apply cleanly.  Do I need to have 29313.patch applied first and then 29315.patch?

Rik <rik5>
Project Administrator
Tue 26 Jan 2021 04:06:29 AM UTC, comment #19: 

Fixed x = 1:.045:3 issue. Nan and inf becoming more annoying to handle. Loop is only called for examples where the values are creeping past the tolerance so it should still be fast for removing numerical error. Timings for the hundred thousand random points are still around .02 seconds. Patch attached.

(file #50786)

Anonymous
Mon 25 Jan 2021 07:42:01 PM UTC, comment #18: 

It's not quite good enough because it's not Matlab compatible.  For the test vector 'x = 1:.045:3' the first four values from Matlab are

1.0, 1.135, 1.27, 1.405

while your code produces

1.0, 1.135, 1.225, 1.315

The difference between the second and third elements (1.225 - 1.135) is .09 which is below the tolerance of 0.1 so the third value is not really unique.

This does turn out to be a hard problem.  It may be that the only real way to resolve this is to move to a C++ solution where there is no overhead for employing loops.  If we want to stay in m-files there may be an iterative solution where the loop body contains a repeated cumsum over the the differences dx and then a calculation of the idx vector of unique values to keep.  The first time through the loop establishes the initial guess for where the unique values are.  Subsequent times through the loop it zeroes out the the differences where the unique values occur and then recomputes the cumsum and idx.  It keeps going until there is no change in the idx values. 

Rik <rik5>
Project Administrator
Mon 25 Jan 2021 05:12:05 AM UTC, comment #17: 

Is the following good enough?

tol=.1;
x=(1:.045:3)';
## the difference check
dx=diff(x);
csdx=cumsum(dx);
shouldkeep=[true;diff([0;csdx-mod(csdx,tol)])>eps(max(abs(x)))]
x(shouldkeep)

It's not exactly .1 between each point but it is close. I don't imagine exactly subtracting without a loop. I am interested if this is the same as Matlab. Should I update uniquetol with the above suggestion?

I thought of a way to decrease the number of iterations of the loop for the by rows method but it would require sorting for each unique value so it may not be faster. I could also write it in c++ which would hopefully remove much of the loop overhead.

Anonymous
Mon 25 Jan 2021 03:10:22 AM UTC, comment #16: 

Unfortunately my suggestion for a quick algorithm is imperfect.  Consider this example

x = 1:0.045:3;
y = uniquetol (x, 0.1, 'datascale', 1)

With the quick algorithm this returns just the first value, 1.  But it should return 1, 1.135, 1.27, ...

The problem is in these two lines

dx = diff (xs);
idx = (dx <= tol);

In this case, all of the differences are .045 which is less than the tolerance of .1 so the algorithm assumes none of the values besides the first are unique.  In reality, the differences need to be chained.  The first value is 1, the second is 1.045.  The difference is less than .1 so the algorithm should continue.  The next value is 1.09.  The difference to the original value 1 is .09 which is also below .1 so this value is not unique.  Finally the third value is 1.135 and the difference is .135 which means this last one is a unique value.

I've thought about using cumsum() on the vector dx and checking when it exceeds the threshold tol, but this doesn't quite work because I need to reset the sum of differences to 0 every time a new unique value is found.

Frankly, I'm a bit stumped at the moment.

Rik <rik5>
Project Administrator
Sun 24 Jan 2021 09:30:27 AM UTC, comment #15: 

Implemented a fast version of uniquetol for when byrows is false. My timing for Rik's example is around 0.02 seconds which makes the implementation faster than Matlab (assuming our computers are the same).

I made minor changes elsewhere, ia and ic are now always column vectors regardless of the orientation of the input excluding when a cell array is returned then the cell is a row and the elements within are columns. This makes the output consistent with the current documentation. This implementation also naturally fixed the sorting difference between the previous version and Matlab.

I also changed the unique 'see also' to include uniquetol.

I have attached two patches generated from TortoiseHg.

A note about TortoiseHg in the installation process it replaces files in Microsoft Visual C++ 2015-2019 Redistributable (x64) with older versions of those files which resulted in other software I frequently use to stop working. Uninstalling TortoiseHg does not solve the problem. However, going to add remove programs and clicking remove then repair on Microsoft Visual C++ 2015-2019 Redistributable (x64) fixes the issue and does not appear to break TortoiseHg.

(file #50774, file #50775)

Anonymous
Sat 23 Jan 2021 10:31:33 AM UTC, comment #14: 

The texinfo didn't build for me after the recent change.
Hopefully fixed here:
https://hg.savannah.gnu.org/hgweb/octave/rev/875d799ab0b3

Markus Mützel <mmuetzel>
Project Member
Fri 22 Jan 2021 07:50:18 PM UTC, comment #13: 

I took another stab at making the documentation clearer in this changeset: http://hg.savannah.gnu.org/hgweb/octave/rev/2d26113ddf57.  I also changed the example so that it is not a direct copy of the Matlab one.

Rik <rik5>
Project Administrator
Fri 22 Jan 2021 05:47:18 PM UTC, comment #12: 

There already will be some difference between the regular algorithm and the "ByRows" algorithm.  Also, I would note that we don't have freedom to innovate here.  Matlab has defined how "uniquetol" should behave.  In particular, the criteria deciding for rows is

all(abs(u-v) <= tol*max(abs(A),[],1))

The part after "tol" can be overridden by specifying the "DataScale" option.  I wonder if judicious choices about the tolerance and DataScale options will do what you want.

Representing a point in an N-dimensional vector space as a row vector [x1, x2, ..., xN] then the Matlab criterion above says that a vector is nearly identical if for each separate component (x1 - y1, x2 - y2, ...) the absolute difference is less than some value. 

That seems to work reasonably well.  Take an example in 2-D

x = [0, 1];
y = [0, -1];
norm (x) == norm (y) == 1
So, these two vectors look somewhat the same, BUT
abs (1 - (-1)) == 2 will show that there are unique components and so the vectors themselves are unique.

Rik <rik5>
Project Administrator
Fri 22 Jan 2021 04:34:42 PM UTC, comment #11: 

We could also have different algorithms dependent on "byrows".

Markus Mützel <mmuetzel>
Project Member
Fri 22 Jan 2021 07:20:30 AM UTC, comment #10: 

A comment on the sorting suggestion, it is a very good idea for vector data, when byrows is not true. This is where I started. However, when byrows is true there will be cases where it could miss duplicates.

An example is trying to remove duplicated points in a Clenshaw Curtis multidimensional quadrature sparse grid. The points are positioned such that one and two vector norms along with axis coordinates are the same between points which can be a significant distance apart so when you sort the data not all of the points which are almost on top of each other are grouped together after sorting. The diff function will see a big difference with a point that is not on top of another and then jump back to the group that is on top of other points but have a big difference again leaving a duplicate.

To address this issue I tried sorting several different ways but I was still having difficulties ensuring no duplicates remained. I also did a Google search and tried out every uniquetol implementation I could find and all I could find left duplicates.

So for any algorithm, we implement can I ask that it is able to leave no duplicates in the sparse grid data. I have attached the difficult data generator replace my [w,x]=uniquewx(w,x,1e-10); with the uniquetol implementation or check the sizes of the outputs between the implementations to ensure we are not leaving duplicates. d and r can be around 4 and the rule should be Clenshaw Curtis.

Looks like I'm wrong about rows as an input. I've been misreading a call to unique as uniquetol since I started writing the function.

(file #50765)

Anonymous
Fri 22 Jan 2021 07:10:38 AM UTC, comment #9: 

@Steven: I accidentally pushed the wrong changeset using my name instead of yours. I'm sorry. I'll try to be more concentrated next time.

I agree with the points Rik made. (Thank you for the improvements.) But having the function like it is currently is probably better than not having it at all.

It would still be nice if these points could be improved before a final release of Octave 7 (some time in the future).

If you are still interested in improving this function (and I hope you are), it would be nice if you could provide further changes in the form of patches to the Mercurial repository.
Since the first file you provided had CRLF line endings, I'm assuming you are on Windows. Personally, I like to work with TortoiseHg on that platform because it is very intuitive to use.
See also [1]. But while working on .m files only, you can probably skip the whole part about the VM.

[1]: https://wiki.octave.org/Building_on_Microsoft_Windows

Markus Mützel <mmuetzel>
Project Member
Fri 22 Jan 2021 04:14:51 AM UTC, comment #8: 

I did some further cleanup of the function in this changeset: http://hg.savannah.gnu.org/hgweb/octave/rev/83fe13ca9ce3.

uniquetol now validates the inputs more tightly which is in keeping with Matlab.

I think the algorithm still needs some attention.  I note that

1. Matlab returns zeros(0,1) rather than [] for empty input.

This is a small consideration.  We could implement this or not.

2. Matlab returns sorted results.

This is a more important consideration.

uniquetol ([3 2 1])
ans = [1 2 3]

3. Matlab returns just one Inf value when multiple occur.

uniquetol ([Inf 1 Inf])
ans = [1 Inf]

4. Performance is very slow (for loop).

I benchmarked uniquetol with this code

tic; c = uniquetol (rand (1e3,1e2), 0.1); toc

Matlab: 0.03 seconds
Octave: 14.0 seconds

This isn't complete (needs to handle ia, ic outputs as well as Inf), but it points the way to an efficient algorithm.

xs = sort (A(:));
dx = diff (xs);
idx = (dx <= tol);
xs([false; idx]) = [];

This code runs in .03 seconds for me on the benchmark data above.

5. I didn't find "rows" keyword was valid input to uniquetol.

Rik <rik5>
Project Administrator
Fri 22 Jan 2021 12:32:07 AM UTC, comment #7: 

The "rows" keyword is in the third section of Matlab's documentation. It describes two ways. I would prefer to keep the two ways and it keeps consistency with unique.

The error on line 121 can be removed it became redundant when I called unique for non double or single precision.

Anonymous
Thu 21 Jan 2021 05:59:25 PM UTC, comment #6: 

I received your user details privately.

I didn't see the argument "row" mentioned in Matlab's documentation. I also couldn't get Matlab to accept that argument. Removing it also allowed to simplify (and complete) the input check.

I also wasn't able to trigger the error in line 121. Can it be removed?

Anyway, I pushed the new function here:
https://hg.savannah.gnu.org/hgweb/octave/rev/f3272029d42c

Marking as ready for test.

Markus Mützel <mmuetzel>
Project Member
Mon 18 Jan 2021 06:17:26 PM UTC, comment #5: 

Sorry, I reloaded the page and forgot to attach the file again.

I understand that you'd prefer to not post your email address as plain text in the comments.
Could you download the attached patch, exchange the line starting with "# User" with your name and email address, and re-upload the patch?
I could add your name to the list of contributors in the manual after that. Or do you prefer to not be added to that list?

(file #50746)

Markus Mützel <mmuetzel>
Project Member
Mon 18 Jan 2021 06:13:21 AM UTC, comment #4: 

I fixed the problematic sentence. I also fixed a few typos. I changed the numbers and some of the function calls log10 to log etc to be different from the Matlab documentation but we need to test the same behavior as the function does almost the same thing.

I also changed the function to call unique if the input was not numeric with double or single precision, and issued a warning if it does this.

I did not see another file uploaded by you so I modified my copy.

I can provide contact details but I don't want to give them on a public webpage. How do I go about this?

(file #50736)

Anonymous
Sun 17 Jan 2021 05:07:32 PM UTC, comment #3: 

Thank you very much. That looks very good to me!

All tests pass for me.

I made only some minor changes to the docstring and style (see attachment).

Some remarks:

  • I don't understand the last sentence before the example in the docstring. (I'm no native speaker. So that might be the reason.) Could you please re-phrase that sentence? Maybe split it up in two or more sentences.
  • It looks like some of the BISTs you added are copies from the examples in Matlab's documentation. While this is a good check to see if the function is compatible, I'm not sure if we can take those and distribute them under GPL. Could you please remove those or replace them with other tests?

I didn't check yet if the docstring compiles. But I think that this function is almost ready to be added.

Would you like to appear as the author in the changeset? We usually use a real first (and last) name together with a valid email address as the author of a changeset.
Which name and email address could I use?
Would you like to appear in the list of contributors (that appears at the beginning of the manual)?

If you prefer to remain anonymous, I could also push the change with my name.

Markus Mützel <mmuetzel>
Project Member
Sun 17 Jan 2021 12:51:33 PM UTC, comment #2: 

I have attached the new version, following the code style, with documentation and tests. I copied from unique extensively. There are a few differences between unique and uniquetol. For example, inf does not make sense to combine into a single inf as the difference between two inf values could be anything so I treat it like nan. I do not consider data types other than double and single.

(file #50729)

Anonymous
Thu 14 Jan 2021 12:02:16 PM UTC, comment #1: 

Thank you for your contribution.

I haven't tested yet.
But it would help a lot if you could add some documentation (texinfo). Please, take a look at some other .m file that is part of Octave at how that could look like.
Also please add some built-in self test (see e.g. %!test sections in other files) and maybe also some %!demo-s.
Ideally, please provide a patch that adds that file to the default branch of Octave's Mercurial repository. But if you aren't set up for that, the .m file will also do. (Please don't paste the function as plain text in a comment.)
In that file, please add the GPL license that you can find at the top of other .m files as well. That tells us that you are willing (and able) to distribute the function under that license.
Also, please try and follow Octave's coding style:
https://wiki.octave.org/Octave_style_guide

I sincerely hope this doesn't intimidate in any way. We are always happy to see new contributors. So, please don't hesitate to ask if anything is unclear.

Markus Mützel <mmuetzel>
Project Member
Sun 10 Jan 2021 12:10:34 PM UTC, original submission:  

I get the message uniquetol is missing from Octave and I could not find it in a package so I have implemented it. I do not have Matlab so I can not directly compare it but I have run all of the examples on the Matlab uniquetol help page and I received the same output from my implementation as shown on that page.

The Implementation

function [C,IA,IC]=uniquetol(A,varargin)


maxabsA=max(abs(A(:)));
ByRows=false;
OutputAllIndices=false;
DataScale=maxabsA;
sizeA=size(A);
if nargin()==1||~(isnumeric(varargin{1})&&isscalar(varargin{1}))
  if isa(A,'double')
    tol=1e-12;
  elseif isa(A,'single')
    tol=1e-6;
  endif
else
  tol=varargin{1};
end

for k=1:length(varargin)
  if isnumeric(varargin{k})
    continue
  elseif ischar(varargin{k})
    if strcmpi(varargin{k},'rows')
      ByRows=true;
    elseif strcmpi(varargin{k},'byrows')
      ByRows=logical(varargin{k+1});
    elseif strcmpi(varargin{k},'OutputAllIndices')
      OutputAllIndices=logical(varargin{k+1});
    elseif strcmpi(varargin{k},'DataScale')
      DataScale=varargin{k+1}(:).';
      columnsDataScale=columns(DataScale);
      if ~(columnsDataScale==1||columnsDataScale==sizeA(2))
        error("Incorrect size of data scale.");
      endif
    endif
  endif
end
if ~ByRows
  A=A(:);
end
x=A;
points=rows(x);
d=columns(x);
Iall=zeros(points,1);
I=nan(d,1);
IA={};
J=nan(d,1);
j=1;
ii=0;
tolDataScale=tol*DataScale;
for i=1:points
  if any(Iall==i)
    continue
  else
    eq=all(abs(x-x(i,:))<tolDataScale,2);
    sumeq=sum(eq);
    ia=find(eq);
    if OutputAllIndices
      IA{end+1}=ia;
    endif
    Iall(ii+(1:sumeq))=ia;
    I(j)=ia(1);
    J(eq)=j;
    ii+=sumeq;
    j+=1;
  end
end
C=x(I,:);
if sizeA(1)==1&&sizeA(2)>1
  C=C.';
end
if ~OutputAllIndices
  IA=I(1:j-1);
endif
IC=J;
end

Anonymous

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #50913:  29393.patch added by None (5KiB - application/octet-stream)
file #50787:  uniquetol.m added by None (12KiB - text/plain)
file #50786:  29315.patch added by None (2KiB - application/octet-stream)
file #50774:  29313.patch added by None (4KiB - application/octet-stream)
file #50775:  29314.patch added by None (780B - application/octet-stream)
file #50765:  squaresparse1dtomd.m added by None (3KiB - text/plain)
file #50736:  uniquetol.m added by None (8KiB - text/plain)
file #50729:  uniquetol.m added by None (7KiB - text/plain)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by mmuetzel (Posted a comment)
  • -email is unavailable- added by jwe (Updated the item)
  •  

    Do you think this task is very important?
    If so, you can add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.

     

     

     

    Follow 17 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2021-02-24 None Attached File- => Added 29393.patch, #50913
    2021-01-26 None Attached File- => Added uniquetol.m, #50787
    2021-01-26 None Attached File- => Added 29315.patch, #50786
    2021-01-24 None Attached File- => Added 29313.patch, #50774
        Attached File- => Added 29314.patch, #50775
    2021-01-22 None Attached File- => Added squaresparse1dtomd.m, #50765
    2021-01-22 rik5 StatusReady For Test => In Progress
    2021-01-21 mmuetzel StatusPatch Submitted => Ready For Test
    2021-01-18 mmuetzel Attached File- => Added bug59850_uniquetol.patch, #50746
        StatusIn Progress => Patch Submitted
    2021-01-18 None Attached File- => Added uniquetol.m, #50736
    2021-01-17 mmuetzel StatusNone => In Progress
    2021-01-17 None Attached File- => Added uniquetol.m, #50729
    2021-01-11 jwe Severity3 - Normal => 1 - Wish
        Priority5 - Normal => 3 - Low
        Item GroupMatlab Compatibility => Feature Request
        Release6.1.0 => dev

    Back to the top


    Powered by Savane 3.5