From 9305c426ffb189a178c946cdf96bfa46ebe0b976 Mon Sep 17 00:00:00 2001
From: Ed J <mohawk2@users.noreply.github.com>
Date: Tue, 7 Jan 2025 04:12:56 +0000
Subject: [PATCH] add Ufunc::partial - #516

---
 Changes          |  2 +-
 lib/PDL/Ufunc.pd | 73 ++++++++++++++++++++++++++++++++++++++++++++++--
 t/ufunc.t        | 12 ++++++++
 3 files changed, 83 insertions(+), 4 deletions(-)

diff --git a/Changes b/Changes
index 36c4ee932..e04eb7d3c 100644
--- a/Changes
+++ b/Changes
@@ -8,7 +8,7 @@
 - now an error to upd_data if datasv set to different size from nbytes
 - add PDL::readonly method (#516) - thanks @wlmb for suggestion
 - rename Ufunc::diffover to Ufunc::numdiff (#516)
-- add Ufunc::diffcentred (#516)
+- add Ufunc::{diffcentred,partial} (#516)
 
 2.098 2025-01-03
 - fix Windows build problems
diff --git a/lib/PDL/Ufunc.pd b/lib/PDL/Ufunc.pd
index 3f69f7b8c..6fd5b2354 100644
--- a/lib/PDL/Ufunc.pd
+++ b/lib/PDL/Ufunc.pd
@@ -9,6 +9,8 @@ pp_addpm({At=>'Top'},<<'EOD');
 use strict;
 use warnings;
 
+=encoding utf8
+
 =head1 NAME
 
 PDL::Ufunc - primitive ufunc operations for pdl
@@ -225,7 +227,7 @@ Unlike L</diff2>, output vector is same length.
 Originally by Maggie J. Xiong.
 
 Compare to L</cumusumover>, which acts as the converse of this.
-See also L</diff2>, L</diffcentred>, L<PDL::Primitive/pchip_chim>.
+See also L</diff2>, L</diffcentred>, L</partial>, L<PDL::Primitive/pchip_chim>.
 EOF
 );
 
@@ -256,7 +258,7 @@ calculating "curl".
 
 By using L<PDL::Slices/xchg> etc. it is possible to use I<any> dimension.
 
-See also L</diff2>, L</numdiff>, L<PDL::Primitive/pchip_chim>.
+See also L</diff2>, L</partial>, L</numdiff>, L<PDL::Primitive/pchip_chim>.
 
 =for usage
 
@@ -269,6 +271,71 @@ A bad value at C<n> means the affected output values at C<n-2>,C<n>
 EOF
 );
 
+pp_add_exported('partial');
+pp_addpm(<<'EOF');
+=head2 partial
+
+=for ref
+
+Take a numerical partial derivative along a given dimension, either
+forward, backward, or centred.
+
+See also L</numdiff>, L</diffcentred>, and L<PDL::Slices/mv>, which
+are currently used to implement this.
+See also L<PDL::Primitive/pchip_chim> and  L</numdiff>.
+
+Can be used to implement divergence and curl calculations (adapted
+from Luis Mochán's work at
+https://sourceforge.net/p/pdl/mailman/message/58843767/):
+
+  use v5.36;
+  use PDL;
+  sub curl ($f) {
+    my ($f0, $f1, $f2) = $f->using(0..2);
+    my $o = {dir=>'c'};
+    pdl(
+      $f2->partial(1,$o) - $f1->partial(2,$o),
+      $f0->partial(2,$o) - $f2->partial(0,$o),
+      $f1->partial(0,$o) - $f0->partial(1,$o),
+    )->mv(-1,0);
+  }
+  sub div ($f) {
+    my ($f0, $f1, $f2) = $f->using(0..2);
+    my $o = {dir=>'c'};
+    $f0->partial(0,$o) + $f1->partial(1,$o) + $f2->partial(2,$o);
+  }
+  sub trim3d ($f) { $f->slice(':,1:-2,1:-2,1:-2') } # adjust if change "dir"
+  my $z=zeroes(5,5,5);
+  my $v=pdl(-$z->yvals, $z->xvals, $z->zvals)->mv(-1,0);
+  say trim3d(curl($v));
+  say div($v);
+
+=for usage
+
+  $pdl->partial(2);           # along dim 2, centred
+  $pdl->partial(2, {d=>'c'}); # along dim 2, centred
+  $pdl->partial(2, {d=>'f'}); # along dim 2, forward
+  $pdl->partial(2, {d=>'b'}); # along dim 2, backward
+
+=cut
+
+my %dirtype2func = (
+  f => \&numdiff,
+  b => sub { $_[0]->slice('-1:0')->numdiff },
+  c => \&diffcentred,
+);
+*partial = \&PDL::partial;
+sub PDL::partial {
+  my ($f, $dim, $opts) = @_;
+  $opts ||= {};
+  my $difftype = $opts->{dir} || $opts->{d} || 'c';
+  my $func = $dirtype2func{$difftype} || barf "partial: unknown 'dir' option '$difftype', only know (@{[sort keys %dirtype2func]})";
+  $f = $f->mv($dim, 0) if $dim;
+  my $ret = $f->$func;
+  $dim ? $ret->mv(0, $dim) : $ret;
+}
+EOF
+
 pp_def('diff2',
   HandleBad => 1,
   Pars => 'a(n); [o]o(nminus1=CALC($SIZE(n) - 1))',
@@ -292,7 +359,7 @@ By using L<PDL::Slices/xchg> etc. it is possible to use I<any> dimension.
 Unlike L</numdiff>, output vector is one shorter.
 Combined with C<slice('-1:0')>, can be used for backward differencing.
 
-See also L</numdiff>, L</diffcentred>, L<PDL::Primitive/pchip_chim>.
+See also L</numdiff>, L</diffcentred>, L</partial>, L<PDL::Primitive/pchip_chim>.
 
 =for usage
 
diff --git a/t/ufunc.t b/t/ufunc.t
index 98b810dd1..bf1b365b1 100644
--- a/t/ufunc.t
+++ b/t/ufunc.t
@@ -239,6 +239,18 @@ subtest 'minimum_n_ind' => sub {
     }
 };
 
+subtest partial => sub {
+  my $a = (sequence(4,4) + 2) ** 2;
+  is_pdl $a->partial(0), pdl('-8 6 8 -6; -16 14 16 -14;
+    -24 22 24 -22; -32 30 32 -30'), "partial(0)";
+  is_pdl $a->partial(0,{d=>'f'}), pdl('4 5 7 9; 36 13 15 17;
+    100 21 23 25; 196 29 31 33'), "partial(0,f)";
+  is_pdl $a->partial(1), pdl('-80 -88 -96 -104; 48 56 64 72;
+    80 88 96 104; -48 -56 -64 -72'), "partial(1)";
+  is_pdl $a->partial(1,{d=>'f'}), pdl('4 9 16 25; 32 40 48 56;
+    64 72 80 88; 96 104 112 120'), "partial(1,f)";
+};
+
 subtest numdiff => sub {
   my $a = sequence(5) + 2;
   is_pdl $a->numdiff, pdl(2, 1, 1, 1, 1), "numdiff";