in reply to Re^4: first stumbling steps in PDL
in thread first stumbling steps in PDL

Now I have two different solutions to cope with a missing argument to tricpy. Both meet the target, as all these statements are valid and produce the same result: 

tricpy($m, 0, $c);
$c = $m->tricpy(0);
$c = $m->tricpy

[download]
Though I like the first solution more, it raises a documentation issue.

First: 

pp_def(
    'tricpy',
    Pars => 'A(m,n);[o] C(m,n)',
    OtherPars => 'int uplo',
    OtherParsDefaults => {uplo => 0},
    ArgOrder => [qw(A uplo C)],
    Code => '
            ...
    ',
    Doc => <<EOT
=for ref

Copy triangular part to another matrix. If uplo == 0 copy upper triang
+ular part.

=cut

EOT

[download]

Second: 

pp_def(
    'tricpy',
    Pars => 'A(m,n);uplo();[o] C(m,n)',
    PMCode => 'sub PDL::tricpy {
        return PDL::_tricpy_int(@_) if @_ == 3;
        my ($A, $uplo) = @_;
        PDL::_tricpy_int($A, $uplo // 0, my $C = PDL->null);
        $C;
    }',
    Code => '
        ...
    ',
    Doc => <<EOT
=for ref

Copy triangular part to another matrix. If uplo == 0 copy upper triang
+ular part.

=cut

EOT

);

[download]

The generated POD for the first version is misleading, as the signature doesn't conform to the actual argument order: 

  tricpy
      Signature: (A(m,n);[o] C(m,n); int uplo)

    Copy triangular part to another matrix. If uplo == 0 copy upper
    triangular part.

    tricpy does not process bad values. It will set the bad-value flag
+ of
    all output ndarrays if the flag is set for any of the input ndarra
+ys.

[download]

Any suggestions?

Greetings,
-jo

$gryYup$d0ylprbpriprrYpkJl2xyl~rzg??P~5lp2hyl0p$

Replies are listed 'Best First'.
Re^6: first stumbling steps in PDL
by etj (Priest) on Feb 06, 2024 at 17:58 UTC
    I like that you explored both options! You've joined quite a small select band of people who've used PP to make/modify functions in PDL. I do think that's a pity, since with the examples and hopefully the docs these days, more people should give it a go.

    OtherPars

    I thought of the OtherPars option, and was inclining away from it since leaving the flexibility of broadcasting over uplo seemed a bit more powerful. However, the clarity gained from having the calling code say what kind of tricpy is being done, and that not changing per broadcasted item seems more important on balance. To solve the documentation problem, use your example code in a usage section, like in the Ops ones: 
    =for usage

  tricpy($m, 0, $c); # all params
  $c = $m->tricpy(0); # explicitly supply uplo
  $c = $m->tricpy; # uplo defaults to 0

    [download]
    There's another outcome from using OtherPars: the if branches should each contain a broadcastloop for efficiency, since obviously the OtherPars won't change during broadcasting.

    Please do add tests for the new functionality.

    PMCode

    Minor nits for another time (please run with the OtherPars route): the initial shortcut should use  >= 3 in case the user supplied too many arguments and the XS can provide a good error. It looks like you dropped the int from the type, which is probably to be avoided.

    Complex data

    One of the legacy things that, at least for now, PDL::LinearAlgebra has to deal with is the complications of dealing with complex-number data in one of two ways, with the pp_defc function that does some horrible (albeit successful) mangling. I was just wondering whether that would require modifying here. It doesn't, since tricpy is handled in its own way, but until just now (including in the still-latest released version), the "real" version (which gets called for native-complex data) didn't handle native-complex types, and would have been silently coerced to real. I've now fixed that bug.

    The only consequence for this change is simply that in Complex/complex.pd there is a ctricpy which will need the same changes making to it as for the real version.

[reply]
[d/l]
[select]

      Either I don't understand how ctricpy is supposed to work, or it is broken, or both. If a complex matrix is fed to ctricpy, the result, its type and dimensions look strange. OTOH, with a "matrix" representing Re and Im, the result looks good, but it can be achieved using tricpy as well (after reordering dims). 

      #!/usr/bin/perl

use v5.24;
use warnings;
use PDL;
use PDL::LinearAlgebra::Complex;

my $m1 = pdl '[[0, i], [1, 1+i]]';
ctricpy($m1, 0, my $c = null);
say $c->info;
say "complex: $c";

my $m2 = sequence 2, 3, 3;
ctricpy($m2, 0, $c = null);
say $c->info;
say "Re + Im:", $c->reorder(1, 2, 0);
say "tricpy:", $m2->reorder(1, 2, 0)->tricpy(0);

__DATA__
PDL: LDouble D [2,2,1]
complex: 
[
 [
  [0 0]
  [1 1]
 ]
]

PDL: Double D [2,3,3]
Re + Im:
[
 [
  [ 0  2  4]
  [ 0  8 10]
  [ 0  0 16]
 ]
 [
  [ 1  3  5]
  [ 0  9 11]
  [ 0  0 17]
 ]
]

tricpy:
[
 [
  [ 0  2  4]
  [ 0  8 10]
  [ 0  0 16]
 ]
 [
  [ 1  3  5]
  [ 0  9 11]
  [ 0  0 17]
 ]
]

      [download]

      Adding complex as type to the input and output arguments in tricpy's signature makes it complex-aware and it works for real and complex matrices. Thus ctricpy might be superfluous. The result from an accordingly modified tricpy on the complex matrix $m1 is: 

      [
 [  0   i]
 [  0 1+i]
]

      [download]

      Greetings,
      -jo

      $gryYup$d0ylprbpriprrYpkJl2xyl~rzg??P~5lp2hyl0p$
[reply]
[d/l]
[select]
        It's not broken. The "legacy" stuff I referred to in passing above I now need to explain a bit more. PDL in about 2000 gained a PDL::Complex class, which has a badly-hidden 0th dimension length 2, which forms pairs of real, imaginary. There are also a number of methods that provide some of the operations that real-valued ndarrays, that know about having the first dimension as real, imaginary. Various non-core modules like PDL::FFTW3 and PDL::LinearAlgebra partly worked around the many limitations of PDL::Complex.

        In 2021 I ran with Ingo Schmid's pull-request that partly added support for C99 "native complex" types in PDL, and after quite a lot of bug-fixing, plus fixing up those two non-core modules to work properly with both "native complex" and PDL::Complex.

        ctricpy is only for working with PDL::Complex, which you can tell by looking at its "signature" with the leading dimension length 2. After I added the foolishly missed-out type-support for "native complex" in the last few days, tricpy is now suitable for use with native complex data, as you have shown above.

        Why does each operation need to "opt in" to support native-complex data? Because if you default to all of them supporting it, stuff breaks, as shown after 2.026_01 was first released: https://sourceforge.net/p/pdl/mailman/pdl-general/?viewmonth=202102

        The current state of affairs with PDL and complex numbers is you should only use native-complex types for complex-number data. Don't touch PDL::Complex at all. Don't use ctricpy. PDL::LinearAlgebra and PDL::FFTW3 work correctly with them (unless they don't, in which case please report bugs).

[reply]
[d/l]
[select]

In Section Seekers of Perl Wisdom