bugGNU Octave - Bugs: bug #61199, Allow 'char' input sets to...

 
 

bug #61199: Allow 'char' input sets to nchoosek for Matlab compatibility

Submitter:  None
Submitted:  Tue 21 Sep 2021 08:51:51 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Matlab Compatibility
Status:  Fixed Assigned to:  None
Originator Name:  Originator Email:  -email is unavailable-
Open/Closed:  * Closed Release:  * dev
Operating System:  * Any Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Wed 29 Sep 2021 06:13:37 PM UTC, comment #41: 

Marking bug as fixed and closing report.

Rik <rik5>
Group administrator
Wed 29 Sep 2021 06:12:31 PM UTC, comment #40: 

I think that was just an oversight of mine.  I was using the test case "nchoosek (63,19)" and hence always going through the branch of code where K was odd.  I made a small change here (http://hg.savannah.gnu.org/hgweb/octave/rev/750de16fd35c) so that the branch where K is even uses the same strategy as when K is odd.  Specifically, I reverse the order of the second range.  As an example, take K = 4.  Previously the denominator (1:K) was split into pairs of 1:2 and 3:4 and the pairwise multiplication was


[1 2] .* [3 4] = [3 8]


where the first entry is odd in the result is odd.  Now I compute


[1 2] .* [4 3] = [4 6]


and the result only contains even numbers which can safely be divided by 2.

Rik <rik5>
Group administrator
Wed 29 Sep 2021 05:56:34 PM UTC, comment #39: 

Thanks, Rik. Your latest patch 750de16fd35c seems to fix it.

Anonymous
Wed 29 Sep 2021 01:44:10 PM UTC, comment #38: 

Rik, sorry, I just now did an older test. The divide-by-two technique fails for nchoosek (56,28). I think it's better to leave numer and denom as they are without dividing by two.

Anonymous
Wed 29 Sep 2021 01:08:33 AM UTC, comment #37: 

Thank you, Rik. Patch looks good to me.

Anonymous
Wed 29 Sep 2021 12:18:52 AM UTC, comment #36: 

I did some more benchmarking and checked in the final changeset here: http://hg.savannah.gnu.org/hgweb/octave/rev/f5aba4d7e651.

Benchmark times (in microseconds)

Original code: 113
New code w/gcd loop: 484
New code w/gcd loop & pairwise multiplication: 371
New code w/gcd loop, pairwise multiplication, factor 2 removal: 339

I was able to get another 8.5% improvement by realizing that pairwise multiplication of an odd and an even number will always produce an even number.  Hence, one can take out a factor of 2 immediately from the numerator and denominator.

Overall, I'm happy with the result.  It is more accurate and no human is going to notice the difference between 0.1 millisecond and 0.3 milliseconds given that reaction times are on the order of 750 milliseconds.

I'm ready to close this bug if everyone else is.

Rik <rik5>
Group administrator
Sat 25 Sep 2021 12:49:26 PM UTC, comment #35: 

Rik, thank you for the profile information. I experimented with primes(k) and factor() and different loop structures, which all turned out to be slower than calling gcd() with the existing loops. This is likely because gcd() is compiled while the others are interpreted. Replacing max() with find (..., 1) was also not an improvement.

But it turns out that making num and denom have fewer elements by pairwise multiplication does make it faster by some 15% to 20%. Even though the individual elements are larger and each call to gcd() takes a little longer, the number of elements is now about half what it was, and the number of invocations of gcd() is reduced correspondingly. On my machine it falls from 280-285 microseconds to 230 microseconds. Still slower than the original code, so it's accuracy at the expense of speed.

Patch is now as follows. Please test it and profile it.

  if (n == 1 && isnumeric (v))
    k = min (k, v-k);

    # make a list of integers for numerator and denominator,
    # then filter out common factors, and multiply what remains.
    # for a 15% to 20% performance boost, we multiply values pairwise so there are fewer elements.
    if (rem(k,2)) # k is odd
      num = [(v-k+1:v-(k+1)/2) .* (v-(k-1)/2:v-1), v];
      denom = [(1:(k-1)/2) .* ((k+1)/2:k-1), k];
    else # k is even
      num = (v-k+1:v-k/2) .* (v-k/2+1:v);
      denom = (1:k/2) .* (k/2+1:k);
    end

    # remove common factors from num and denom
    while (~isempty(denom))
      for i = length(denom):-1:1
        tmp = gcd (denom(i), num);
        [m,f] = max(tmp);
        denom(i) /= m;
        num(f) /= m;
      endfor
      denom = denom(denom > 1);
      num = num(num > 1);
    endwhile

    C = prod (num);
    if (C > flintmax)
      warning ("nchoosek: possible loss of precision");
    end


Anonymous
Sat 25 Sep 2021 05:02:51 AM UTC, comment #34: 

The individual calls to nchoosek don't take very long so I put them in a loop to benchmark it.


N = 1e4;

for i = 1:N
  y = nchoosek (63,19);
endfor


Results are an average of 380 microseconds for the patched code and 85 microseconds for the original code.  Even 380 microseconds is pretty acceptable, but if you see room for improvement that would be a bonus.  I ran the profiler on the patched code and I got this


  bm_nchoosek: 1 calls, 5.043 total, 0.166 self
    nchoosek: 10000 calls, 4.877 total, 3.862 self
      1) gcd: 200000 calls, 0.437 total, 0.437 self
      2) max: 200000 calls, 0.387 total, 0.387 self
      3) binary >: 50000 calls, 0.036 total, 0.036 self
      4) prefix !: 50000 calls, 0.026 total, 0.026 self
      5) min: 10000 calls, 0.019 total, 0.019 self
      6) isempty: 30000 calls, 0.014 total, 0.014 self
      7) prod: 10000 calls, 0.013 total, 0.013 self
      8) isscalar: 20000 calls, 0.009 total, 0.009 self
      9) length: 20000 calls, 0.009 total, 0.009 self
      10) fix: 20000 calls, 0.007 total, 0.007 self
      11) isvector: 10000 calls, 0.007 total, 0.007 self
      12) isnumeric: 20000 calls, 0.005 total, 0.005 self
      13) binary !=: 20000 calls, 0.005 total, 0.005 self
      14) prefix -: 20000 calls, 0.004 total, 0.004 self
      15) flintmax: 10000 calls, 0.004 total, 0.004 self
      16) binary ==: 20000 calls, 0.004 total, 0.004 self
      17) nargin: 10000 calls, 0.004 total, 0.004 self
      18) binary <: 20000 calls, 0.004 total, 0.004 self
      19) binary -: 20000 calls, 0.004 total, 0.004 self
      20) isreal: 10000 calls, 0.003 total, 0.003 self
      21) iscomplex: 10000 calls, 0.003 total, 0.003 self
      22) numel: 10000 calls, 0.003 total, 0.003 self
      23) binary >=: 10000 calls, 0.003 total, 0.003 self
      24) binary +: 10000 calls, 0.002 total, 0.002 self


Items 1 (gcd) and 2 (max) are the things to concentrate on.

Rik <rik5>
Group administrator
Sat 25 Sep 2021 02:36:25 AM UTC, comment #33: 

This might not be a performance issue in a practical sense but I observe that the gcd code is basically O(n^3) because of the while loop, for loop and max inside it. Please do keep this report open for a few days. I'll try to speed up the gcd code so that it's a bit faster than what it is now. Right now the unpatched code takes 60 microseconds while the gcd takes 280 microseconds. This probably doesn't affect anything of importance but I might as well make it faster while the changes are happening rather than after everything is settled.

Anonymous
Sat 25 Sep 2021 02:07:33 AM UTC, comment #32: 

Matlab gives same result as unpatched Octave:


nchoosek (63,19)
ans = 6131164307078482


So, it looks like a good idea to use the gcd() code and we will be better than Matlab.

Rik <rik5>
Group administrator
Sat 25 Sep 2021 12:37:40 AM UTC, comment #31: 

@jwe: Until I read your comment, I had not even realized we even had a preview button! I had looked for it years ago next to the submit button and not finding it there it never occurred to me that it's above the comment box. Speaking of Rik's comment that an old function like nchoosek is being heavily improved suddenly, I find myself suddenly enlightened about the preview button after years of using Octave.

Anonymous
Sat 25 Sep 2021 12:29:47 AM UTC, comment #30: 

The comments are getting messed up because of mismatched verbatim tags.  In one of the comments for this report, the end verbatim tag had + and - instead of - and - and in the other case, the second - was missing.  The comment processor used by the savannah trackers is unfortunately very unforgiving.  The preview button can help but it's also easy to forget to preview before posting.

John W. Eaton <jwe>
Group administrator
Sat 25 Sep 2021 12:25:19 AM UTC, comment #29: 

Case where the result is off by 7 AND incorrectly gives no precision warning at all even though it's off by 7: "nchoosek (63,19)"

Unpatched Octave, incorrect answer with no warning:

>> nchoosek (63,19)
ans = 6131164307078482


Patched Octave:

>> nchoosek (63,19)
ans = 6131164307078475


Please compare Matlab as well. Expecting it will be correct since the result is less than flintmax.

Anonymous
Fri 24 Sep 2021 11:51:35 PM UTC, comment #28: 

Please also try nchoosek (100,13) in both Matlab and unpatched Octave. I do not have Matlab but in unpatched Octave it's off by 3 along with a warning. Patched Octave gives the right answer as shown in Comment #23. Not sure what Matlab gives for that.

Anonymous
Fri 24 Sep 2021 11:43:40 PM UTC, comment #27: 

What in the world is going on with the comments being cut off? Let me try again. This is what I get with unpatched Octave. Observe the warning, which doesn't show up in Matlab.


octave:4> format long g compact
octave:5>
octave:5> nchoosek (56,28)
warning: nchoosek: possible loss of precision
warning: called from
    nchoosek at line 118 column 7

ans = 7648690600760440
octave:6> ver
----------------------------------------------------------------------
GNU Octave Version: 7.0.0 (hg id: 54774a713d7c)


Anonymous
Fri 24 Sep 2021 11:41:10 PM UTC, comment #26: 

@Rik:

First, your email with the cut-off portion is here:

https://lists.gnu.org/archive/html/octave-bug-tracker/2021-09/msg00238.html

Second, this is what I get for unpatched Octave:


octave:4> format long g compact
octave:5>
octave:5> nchoosek (56,28)
warning: nchoosek: possible loss of precision
warning: called from
    nchoosek at line 118 column 7

ans = 7648690600760440
octave:6> ver
----------------------------------------------------------------------
GNU Octave Version: 7.0.0 (hg id: 54774a713d7c)
-verbatim

The result is accurate but I get the warning about loss of precision. I've included version information.

With the patch, there's no warning or discrepancy up to flintmax.

Do you get the same result?

Anonymous
Fri 24 Sep 2021 10:09:47 PM UTC, comment #25: 

Some of my response was cut off.  If someone receives it in e-mail it would be nice to re-post the entire response to this bug report.

In order to stay compatible with Matlab I think we have to return doubles.  However, we could switch to the the algorithm that gives more exact results when the return value is close to flintmax.

Rik <rik5>
Group administrator
Fri 24 Sep 2021 10:05:04 PM UTC, comment #24: 

@Anonymous: Matlab gives no warnings for nchoosek (54,27) and nchoosek (56,28).  The results are the same for Matlab and for unpatched Octave, so it is not clear that the improved code with gcd() is necessary for these cases.

In Matlab, nchoosek (58,29) does produce a warning.  Specifically, "Warning: Result may not be exact.  Coefficient is greater than 9.007199e15 and is only accurate to 15 digits.

The result is the same in either Matlab or Octave


format long
z = nchoosek (58,29)
z = 3.006726649954106e+16
+verbatim-

With the patched code that uses gcd() the answer differs by 2 in the last position

+verbatim+
z = 3.006726649954104e+16


There are some differences.  For Matlab and Octave unpatched code,


nchoosek (56,26)
ans = 6.646448384109071e+15


BUT, for patched code


nchoosek (56,26)
ans = 6646448384109072


There is a difference of 1 in the final position.  This is correct as the answer should be even.  But, is it worth it?  eps (6.646448384109071e+15) is 1 so essentially the current calculation is within the error bar for a numerical approximation.  I don't have much of an opinion here.

Rik <rik5>
Group administrator
Fri 24 Sep 2021 09:42:05 PM UTC, comment #23: 

Could someone please check if Matlab gives integers beyond flintmax, like nchoosek (64,32)? Where does Matlab start issuing warnings about accuracy?

Test for comment #21:


>> format long g

>> nchoosek (100,13)  ## NO PATCH, GETS WARNING, RESULT SLIGHTLY OFF
warning: nchoosek: possible loss of precision
warning: called from
    nchoosek at line 118 column 7

ans = 7110542499799197


>> nchoosek (100,13) ## PATCH, NO WARNING, RESULT EXACT
ans = 7110542499799200


Please compare with Matlab for presence of warnings and exactness.

Anonymous
Fri 24 Sep 2021 09:08:16 PM UTC, comment #22: 

@nrjank: For extended documentation on the BIST system see Appendix B Test and Demo Functions in the Octave manual.  It covers the various forms and how to use the conditionals like %!testif, etc.

I find it fascinating that this particular function has been in Octave for years, and is suddenly undergoing a great deal of improvement.

I checked in @nrjank's changeset to allow struct array inputs.  I did change the BIST test because the k==2 code was not getting tested because k==n elseif branch was being hit first.  I simply expanded the struct array to 3 elements in order to get this to work.


diff -r 281ab99062cf scripts/specfun/nchoosek.m
--- a/scripts/specfun/nchoosek.m        Thu Sep 23 21:06:41 2021 -0400
+++ b/scripts/specfun/nchoosek.m        Fri Sep 24 14:00:51 2021 -0700
@@ -228,16 +228,17 @@ endfunction
 %! s.a = [1 2 3];
 %! s.b = [4 5 6];
 %! s(2).a = 1;  # make s a struct array rather than scalar struct
+%! s(3).b = 2;  # make s at least three elements for k == 2 test below
 %! x = nchoosek (s, 0);
 %! assert (size (x), [1, 0]);
 %! assert (isstruct (x));
 %! assert (fieldnames (x), {"a"; "b"});
 %! x = nchoosek (s, 2);
-%! assert (size (x), [1, 2]);
+%! assert (size (x), [3, 2]);
 %! assert (isstruct (x));
 %! assert (fieldnames (x), {"a"; "b"});
-%! x = nchoosek (s, 3);
-%! assert (size (x), [0, 3]);
+%! x = nchoosek (s, 4);
+%! assert (size (x), [0, 4]);
 %! assert (isstruct (x));
 %! assert (fieldnames (x), {"a"; "b"});


Final changeset is here: http://hg.savannah.gnu.org/hgweb/octave/rev/54774a713d7c




Rik <rik5>
Group administrator
Fri 24 Sep 2021 04:17:27 PM UTC, comment #21: 

Change this in the existing code:

  if (n == 1 && isnumeric (v))
    ## Improve precision at next step.
    k = min (k, v-k);
    C = round (prod ((v-k+1:v)./(1:k)));
    if (C*2*k*eps >= 0.5)
      warning ("nchoosek: possible loss of precision");
    endif
  elseif (k == 0)


to this:

  if (n == 1 && isnumeric (v))
    k = min (k, v-k);
    num = (v-k+1):v;
    denom = 2:k; # no need for 1

    # remove common factors from num and denom
    while (~isempty(denom))
      for i = length(denom):-1:1
        tmp = gcd (denom(i), num);
        [m,f] = max(tmp);
        denom(i) /= m;
        num(f) /= m;
      endfor
      denom = denom(denom > 1);
      num = num(num > 1);
    endwhile

    C = prod (num);
    if (C > flintmax)
      warning ("nchoosek: possible loss of precision");
    end
  elseif (k == 0)


in order to achieve this difference:

>> nchoosek (58,29)   ### PATCHED
ans = 30067266499541040

>> nchoosek (58,29)   ### UNPATCHED
warning: nchoosek: possible loss of precision
warning: called from
    nchoosek at line 118 column 7

ans = 3.006726649954106e+16


If it is acceptable to cast into uint64, we can change the end of the patch as follows:

    C = prod (num);
    if (C > flintmax)
      if (C < intmax("uint64")) # can cast up to get precision
        C = prod (uint64(num), "native");
      else # can't get precision from casting up
        warning ("nchoosek: possible loss of precision");
      end
    end


That will go all the way to 2^64-1 but will return a uint64 instead of a double for values bigger than flintmax but less than intmax("uint64").

>> nchoosek (67,33)
ans = 14226520737620288370
>> class(nchoosek (67,33))
ans = uint64


Anonymous
Fri 24 Sep 2021 08:01:40 AM UTC, comment #20: 

Matlab 2021b prerelease (to expire in 6 days):

>> nchoosek (56,28)
ans =
   7.6487e+15
>> nchoosek (54,27)
ans =
   1.9469e+15
>>


No warnings whatsoever.

Philip Nienhuis <philipnienhuis>
Group Member
Fri 24 Sep 2021 05:10:04 AM UTC, comment #19: 

Hello,

I made some changes to the binomil coefficient calculation part of nchoosek so that it returns exact integers all the way to flintmax and optionally up to 2^64‐1. Could someone verify if nchoosek (54,27) and nchoosek (56,28) give warnings of loss of precision in Matlab please? Is there is interest in extending the range of exactness for Octave or is it a Matlab compatibility issue?

Anonymous
Fri 24 Sep 2021 02:54:59 AM UTC, comment #18: 

Yeah, we should probably document the preferred style for writing tests.

%!xtest came before the idea of adding bug report numbers to the tests.  Now that tests can be tagged with bug numbers, I think %!xtest is mostly obsolete.

You can use


%!test <NNNNNN>
%!  additional test lines
%!  are usually indented

%!assert <NNNNNN> (condition)
%!assert <NNNNNN> (test_expr, expected_result)


Note that using "%!  assert (...)" in a multi-line test block is just a call to the assert function and is not the same as "%!assert", so you can't use a bug number there.

Prefixing a bug number with "*" (i.e., "<*123456>" indicates that the bug status is "fixed", so any failure will be reported as a regression, not a known bug.  The "update-bug-status" target in the Makefile will search for tests with numbers, get the status from the bug tracker, and then edit the sources to add the "*" if needed.

John W. Eaton <jwe>
Group administrator
Fri 24 Sep 2021 01:09:26 AM UTC, comment #17: 

ok, now I'm looking at repelem and there are cell and struct tests in there.

nchoosek was written before we had repelem.  it calls repelems, which can't handle general arrays.

changing calls to repelems to repelem with appropriate calling syntax passes all tests. 

Added a struct k=2 test that now passes.

patch attached.



(file #51957)

Nicholas Jankowski <nrjank>
Group Member
Fri 24 Sep 2021 12:08:09 AM UTC, comment #16: 

No, all good. since it's in the input definition it is technically documented matlab and a missing feature.  I'll create a bug report and put it on my todo list to take a look.

regarding adding bug number to the tests, can you point me to a good example of that syntax (where does the # appear for a %!test block or a one line %!assert or %!warning?) and should we maybe add it to https://wiki.octave.org/Tests ?

Nicholas Jankowski <nrjank>
Group Member
Thu 23 Sep 2021 09:43:33 PM UTC, comment #15: 

I didn't mean to imply that there is anything wrong with repelem.  Octave is compatible with the documented behavior of repelem and nchoosek so it is our choice whether we want to go any further.

I would also expect that if repelem expanded capability to include any arbitrary array then nchoosek would automatically gain the ability to operate on these additional types.

Rik <rik5>
Group administrator
Thu 23 Sep 2021 08:27:03 PM UTC, comment #14: 

You can just add a bug number to a test or assert and it is essentially the same as xtest.

With the bug number, it will show up as an XFAIL in the test report.  Then if the bug is fixed, it will be PASS.  Then, after running "make update-bug-status" to add "*" to the reports that are marked as "fixed" in the bug tracker, any subsequent failure will be shown as a REGRESSION.

And yes, since most of what nchoosek does is indexing, I would expect that if repelem is fixed to handle arrays of any kind of object, then nchoosek would also work for those objects.

John W. Eaton <jwe>
Group administrator
Thu 23 Sep 2021 07:49:27 PM UTC, comment #13: 

as a coauthor of repelem... i don't think I even knew struct arrays were a thing back then.

looking at repelem input description:

Input array, specified as a matrix or multidimensional array.

Data Types: single | double | int8 | int16 | int32 | int64 | uint8 | uint16 | uint32 | uint64 | logical | char | struct | table | cell | datetime | duration
Complex Number Support: Yes


i guess that could maybe be a separate incompatibility report.

looking at the same input description for nchoosek:

set of all choices, specified as a vector.

Example: [1 2 3 4 5]

Example: [1+1i 2+1i 3+1i 4+1i]

Example: int16([1 2 3 4 5])

Example: [true false true false]

Example: ['abcd']

Data Types: single | double | int8 | int16 | int32 | int64 | uint8 | uint16 | uint32 | uint64 | logical | char
Complex Number Support: Yes


no mention of either cells or structs (or strings either), so this seems to fall under chasing completely undocumented matlab?  That said, the output rules do seem fairly consistent and 'input class agnostic'


assuming repelem was expanded for struct arrays, should that more or less cover the cases? maybe put those in as  xtests pending that change?

Nicholas Jankowski <nrjank>
Group Member
Thu 23 Sep 2021 07:23:58 PM UTC, comment #12: 

I checked in a changeset (http://hg.savannah.gnu.org/hgweb/octave/rev/23a907b2dbd5) before some of the additional tests became available.

I checked in an additional changeset that adds BISTs for exceptional cases of scalar structs and struct arrays.  See http://hg.savannah.gnu.org/hgweb/octave/rev/aaee7b170cb1.

Octave's version only works for struct inputs for k = 0, k = 1, k = n, k > n.  Otherwise the algorithm fails because repelem doesn't seem to work for struct.


Rik <rik5>
Group administrator
Thu 23 Sep 2021 06:38:19 PM UTC, comment #11: 

here are a handful of matlab-compatibility BISTs if we're trying to hit most things discussed here (chars, cells, structs), with output based on Matlab 2021a (ignoring char and string differentiation. it sees "abc" as a single string element, but 'abc' as a 3-element char array.)


##char inputs
%!assert (nchoosek ("A", 0), char (ones (1,0)))
%!assert (class (nchoosek ("A", 0)), "char")
%!assert (nchoosek ("A", 2), char (ones (0,2)))
%!assert (class (nchoosek ("A", 2)), "char")
%!assert (nchoosek ("A", 1), "A")
%!assert (class (nchoosek ("A", 1)), "char")
%!assert (nchoosek (["AB"], 1), ["A"; "B"])
%!assert (nchoosek (["A"; "B"], 1), ["A"; "B"])
%!assert (nchoosek (["AB"], 2), "AB")
%!assert (nchoosek (["A"; "B"], 2), "AB")
%!assert (sort (nchoosek (["A"; "B"; "C"], 2)),["AB"; "AC"; "BC"])

##cell inputs (undocumented)
%!assert (nchoosek ({1,2,3}, 1), {[1]; [2]; [3]})
%!assert (nchoosek ({1,2,3}, 3), {[1], [2], [3]})
%!assert (nchoosek ({1,2,3}, 2), {[1], [2]; [1], [3]; [2], [3]})
%!assert (nchoosek ({[1 4], 2, 3}, 1), {[1 4]; [2]; [3]})
%!assert (nchoosek ({[1 4], 2, "Ab"}, 1), {[1 4]; [2]; "Ab"})
%!assert (nchoosek ({{1}, {2}, {3}}, 1), {{1}; {2}; {3}})
%!assert (nchoosek ({{1}, {2}, {3}}, 3), {{1}, {2}, {3}})

%!test ##struct/ struct array inputs (undocumented)
%! A.a = [1 2 3];
%! A.b = [4 5 6];
%! assert (nchoosek (A, 1), A)
%! assert (nchoosek (A, 2), struct("a",[],"b",[]))
## not quite right, should be 0x2 struct array with
## empty fields a and b (not [] or "", just 'empty')
## but can't recreate that in Octave (or Matlab with struct.empty)
%! A(2).a = [7 8 9];
%! A(2).b = [10 11 12];
%! assert (size (A), [1 2]);
%! B = nchoosek (A, 1);
%! assert (size (B), [2 1]);
%! assert ({A(1),A(2)}, {B(1),B(2)});
%! C = nchoosek (A, 2);
%! assert (size (B), [1 2]);
%! assert ({A(1),A(2)}, {C(1),C(2)});
%! A(3).a = [13 14 15];
%! A(3).b = [16 17 18];
%! D = nchoosek (A, 1);
%! assert (size (D), [3 1]);
%! assert ({A(1),A(2),A(3)}, {D(1),D(2),D(3)});
%! E = nchoosek (A, 3);
%! assert (size (E), [1 3]);
%! assert ({A(1),A(2),A(3)}, {E(1),E(2),E(3)});
%! F = nchoosek (A, 2);
%! assert (size (F), [3 2]);
%! assert ({E(1),E(2),E(3),E(4),E(5),E(6)}, {A(1),A(1),A(2),A(2),A(3),A(3)});

%!error (nchoosek (["Ab"; "Cd"], 1))  ##only scalars or vectors


most of these don't run through the current nchoosek, so couldn't test all those for syntax, etc., but verified 'intent' over in matlab.

based on the empty outputs, I'm guessing matlab preallocates an array with size (m, k), m is calculated from n and k, and if zero it just results in size(0,k) when k>n. Then k=0 gives size(m,0), and m is either calculated or set to 1 in that case.

I've had a heck of a time trying to get empty outputs matlab compatible in a couple of cases. Seems it usually shouldn't matter, but its always possible someone will expect it to be for later empty math consistency.
https://octave.org/doc/v6.3.0/Empty-Matrices.html#Empty-Matrices

Nicholas Jankowski <nrjank>
Group Member
Thu 23 Sep 2021 04:19:09 PM UTC, comment #10: 

I'm attaching my draft changes.

I know the hash is somewhat inconvenient, but it's unique.  The numeric id for the changeset could change even in the savannah repo if that repo had to be reconstructed from another copy at some point, or if we were forced to strip a change from the repo.  Also, even if we don't have to rewrite the savannah repo, users have to know that the number corresponds to that repo only, not their local copy.

(file #51954)

John W. Eaton <jwe>
Group administrator
Thu 23 Sep 2021 03:54:25 PM UTC, comment #9: 

Doesn't matter to me who makes the changes.  I have everything except the new BIST tests so I might as well do it.

This could count as "undocumented Matlab".  I had simply tried some sample inputs and found that Matlab accepted cell arrays and struct arrays.

I was using the changeset number from Octave's Mercurial repo at http://hg.savannah.gnu.org/hgweb/octave, rather than my local one.  Isn't that number a unique identifier?  The hash key is longer and less fun to work with.

Rik <rik5>
Group administrator
Thu 23 Sep 2021 03:42:43 PM UTC, comment #8: 

This page

https://www.mathworks.com/help/matlab/ref/nchoosek.html

documents nchoosek as accepting numeric, logical, or char for the 'v' argument.  Are there other overloaded versions of nchoosek for other types, or is it just that nchoosek works for other types even though the documentation doesn't mention them because, as you say, it is essentially just an indexing operation?

I didn't think to check the cases like


nchoosek ("A", 1)


because I didn't see anything other than the use of "zeros" in the code that I expected would convert away from the original type.

The if/elseif thing was also because I couldn't think of a way to get the proper result in any other way.

I would say that generating the exact output dimensions that Matlab does is important to avoid breaking code that was written for Matlab and expects that behavior.

I can fix these issues unless you are already working on it.

BTW, the numeric changeset number is local and can vary between repos.

John W. Eaton <jwe>
Group administrator
Thu 23 Sep 2021 03:08:53 PM UTC, comment #7: 

I think a different coding strategy is necessary.  Changeset 30204 allows 'char' inputs, but what about cell arrays? or struct arrays?  Matlab accepts both because when input V is a vector this is about manipulating indices, rather than the underlying data object.

This input validation seems to work for the wider class of inputs.


  if (nargin != 2
      || ! isvector (v)
      || ! (isreal (k) && isscalar (k)))


Preserving the class of the input in a null output is simpler with indexing than if/elseif/else tree.  Particularly if the input is expanded to include cell arrays.

Instead of


  elseif (k == 0)
    if (is_sq_string (v))
      C = resize ('', 1, 0);
    elseif (is_dq_string (v))
      C = resize ("", 1, 0);
    else
      C = zeros (1, 0, class (v));
    endif


one can write


  elseif (k == 0)
    C = v([]);


One question is whether we want to support bug-for-bug compatibility with Matlab.  Matlab returns empty objects of size 1x0, but there's no real reason to do that.  They are empty so returning a 0x0 object seems better to me.  However, if Octave wants to preserve that then indexing can be done with an appropriately sized zero vector


  elseif (k == 0)
    C = v(zeros (1,0));


Similarly, this code cane be simplified with null indexing.


  elseif (k > n)
    if (is_sq_string (v))
      C = resize ('', 0, k);
    elseif (is_dq_string (v))
      C = resize ("", 0, k);
    else
      C = zeros (0, k, class (v));
    endif


Replacement:


  elseif (k > n)
    C = v([]);


For the corner case of comment #5, I think we need to check that 'v' is numeric before assuming that the binomial coefficient needs to be calculated.  This should do it.


  if (n == 1 && isnumeric (v))


I would add a BIST test for


nchoosek ({1,2,3}, 2)


to verify it works for cell arrays.

Rik <rik5>
Group administrator
Thu 23 Sep 2021 02:51:43 PM UTC, comment #6: 

I think changing "if (n == 1)" to "if (n == 1 && isnumeric (v))" fixes the problem from Comment #5.

Anonymous
Thu 23 Sep 2021 02:18:46 PM UTC, comment #5: 

The results look correct, except for this edge case where n is a single character, where it converts n to the ASCII value and returns the numeric value nCk as a double:


>> nchoosek ("A", 0)
ans = 1
>> nchoosek ("A", 1)
ans = 65
>> class (nchoosek ("A", 0))
ans = double
>> class (nchoosek ("A", 1))
ans = double


Anonymous
Thu 23 Sep 2021 04:01:20 AM UTC, comment #4: 

I pushed the following changeset on default.  Marking as ready for test.

http://hg.savannah.gnu.org/hgweb/octave/rev/1cd077e9f127

John W. Eaton <jwe>
Group administrator
Wed 22 Sep 2021 09:41:35 PM UTC, comment #3: 

Changed the Summary and Item Group to be about Matlab compatibility.

The changes will need to be a little deeper.  For example, when k = 0 the following code is executed


  elseif (k == 0)
    C = zeros (1,0);


This will return a numeric, rather than char, result for


nchoosek ('A':'C', 0)



Rik <rik5>
Group administrator
Tue 21 Sep 2021 09:43:02 PM UTC, comment #2: 

Sorry! Wrong thread! Was meant for Bug #61201. Please delete Comment #1.

Anonymous
Tue 21 Sep 2021 09:41:38 PM UTC, comment #1: 

Simpler example:

debug_on_error(1)
debug_on_warning(1)

NOVER = NLIM = 10000;
nsucc = -666;  nchg=0;
j = 0;
while ((j < 100) & nsucc);
  j++;

  nsucc = 0;
  k = 0;
  while ((k < NOVER) & (nsucc < NLIM));
    k += 1;
    1;
    2;
    3;
    4;
    5;
    6;
    7;
    8;
    9;
    10;
  endwhile
endwhile

[j k nsucc nchg]


It appears the debug on warning always stops and points at the last line inside the while loop, hence "10;" in the above case.

Anonymous
Tue 21 Sep 2021 08:51:51 PM UTC, original submission:  

The nchoosek function currently does not accept character vectors (strings) though that is perfectly valid mathematically. To achieve that requires only one line of input validation to be changed:

Change this:

  if (nargin != 2
      || ! (isreal (k) && isscalar (k))
      || ! (isnumeric (v) && isvector (v)))
    print_usage ();
  endif


to this, removing "isnumeric(v)":

  if (nargin != 2
      || ! (isreal (k) && isscalar (k))
      || ! isvector (v))
    print_usage ();
  endif


Then this can be achieved:

octave:7> nchoosek ('A':'E', 3)
ans =
ABC
ABD
ABE
ACD
ACE
ADE
BCD
BCE
BDE
CDE


Anonymous

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #51957:  bug61199-nchoosek-repelemfix.patch added by nrjank (2KiB - application/octet-stream)
file #51954:  nchoosek-diffs.txt added by jwe (3KiB - 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 nrjank (Posted a comment)
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  •  

    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 9 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2021-09-29 rik5 StatusReady For Test Fixed
        Open/ClosedOpen Closed
    2021-09-24 nrjank Attached File- Added bug61199-nchoosek-repelemfix.patch, #51957
    2021-09-23 jwe Attached File- Added nchoosek-diffs.txt, #51954
    2021-09-23 jwe StatusConfirmed Ready For Test
    2021-09-22 rik5 Item GroupOther Matlab Compatibility
        StatusNone Confirmed
        Release6.3.0 dev
        SummaryInput validation for nchoosek Allow 'char' input sets to nchoosek for Matlab compatibility

    Back to the top

    Powered by Savane 3.13-4448.
    Corresponding source code