@@ -41,10 +41,10 @@ require Exporter;
4141our @ISA =qw( Exporter)
4242our @EXPORT_OK =qw( vectors2Dlist tile RtoG GtoR
4343 HProd MHProd EProd VSProd SProd
44- corner_rotate mvN top_slice linearCombineIt lentzCF any_complex tensor
44+ corner_rotate mvN top_slice linearCombine lentzCF any_complex tensor
4545 make_haydock make_greenp
4646 cartesian_product dummyN triangle_coords incarnate_as
47- wave_operator apply_longitudinal_projection make_dyads
47+ apply_longitudinal_projection make_dyads
4848 cgtsv lu_decomp lu_solve convert_units upper_sqrt
4949)
5050use PDL::LiteF;
@@ -58,12 +58,16 @@ require List::Util;
5858
5959
6060sub top_slice :lvalue {
61+ # slices an ndarray for a given last dimension
6162 my ($pdl , $index ) = @_ ;
6263 my $slice_arg = join ' ,' map ' :' $pdl -> ndims-1)), $index ;
6364 $pdl -> slice($slice_arg );
6465}
6566
6667sub dummyN {
68+ # Inserts a block of $how_many dummy dimensions of size $dim_size
69+ # (default 1) starting at dimension $where (default 0) into an
70+ # ndarray $pdl.
6771 my ($pdl , $how_many , $where , $dim_size ) = @_ ;
6872 return $pdl if $how_many <= 0;
6973 my $ndims =$pdl -> ndims;
@@ -76,35 +80,21 @@ sub dummyN {
7680 $pdl -> slice($slice_arg );
7781}
7882
79- sub linearCombineIt { # complex linear combination of states
80- my ($coefficients , $states , $thread_dims )=@_ ;
81- $thread_dims //= 0;
83+ sub linearCombine { # linear combination of states
84+ my ($coefficients , $states )=@_ ;
8285 $coefficients =dummyN($coefficients , $states -> ndims-$coefficients -> ndims);
83- ($coefficients *$states )-> mv(-$thread_dims - 1,0)-> sumover;
86+ ($coefficients *$states )-> mv(-1,0)-> sumover;
8487}
8588
86- sub any_complex {
89+ sub any_complex { # test if an ndarray is any kind of complex
8790 grep ref $_ && ($_ -> isnull || !$_ -> type-> real), @_ ;
8891}
8992
90- sub wave_operator {
91- my ($green , $nd ) = @_ ;
92- lu_solve([lu_decomp($green )], r2C(PDL::MatrixOps::identity($nd )));
93- }
94-
95- sub cartesian_product {
96- my ($s1 , $s2 ) = @_ ;
97- my $ndims_target = List::Util::max(2, map $_ -> ndims, $s1 , $s2 );
98- $_ = dummyN($_ , $ndims_target -$_ -> ndims) for grep $_ -> ndims < $ndims_target , $s1 , $s2 ;
99- my ($nd1 , $nd2 ) = map $_ -> ndims, $s1 , $s2 ;
100- my @dims = $s1 -> dims;
101- $dims [-2] += $s2 -> dim(-2); # X1+X2
102- $dims [-1] *= $s2 -> dim(-1); # m*n
103- my $res = zeroes(@dims );
104- my ($res_mv , $s1_mv , $s2_mv ) = map mvN($_ , 0, $nd1 -3, -1), $res , $s1 , $s2 ; # now work with first 2 dims
105- $res -> slice(' 0:' $s1_mv -> dim(-2)-1)) .= $s1_mv -> dummy(2, $s2_mv -> dim(-1))-> clump(1,2);
106- $res -> slice($s1_mv -> dim(-2).' :-1' $s2_mv -> dummy(1, $s1_mv -> dim(-1))-> clump(1,2);
107- $res ;
93+ sub cartesian_product { # all pairs of elements of $s1 and $s2
94+ my ($s1 , $s2 ) = @_ ;
95+ my ($a1 , $a2 )=map {$_ -> ndims==1?$_ -> dummy(0):$_ } ($s1 , $s2 );
96+ my $r =append($a1 -> dummy(2), $a2 -> dummy(1));
97+ $r -> reshape($r -> dim(0), $r -> dim(1)*$r -> dim(2));
10898}
10999
110100sub triangle_coords {
@@ -507,7 +497,7 @@ ndarray.
507497Adds C<$how_many > (no-op if <= 0) dummy dimensions of size C<$dim_size >
508498(default 1) in the C<$which_dim > (default 0) position.
509499
510- =item * $r =linearCombineIt ($c , $it , $thread_dims )
500+ =item * $r =linearCombine ($c , $it , $thread_dims )
511501
512502Complex linear combination of states. $c is a 'complex'
513503ndarray and $it is an ndarray of states from a L<Photonic::Roles::Haydock>.
@@ -609,21 +599,15 @@ True if any of the args are a complex PDL.
609599
610600=item * cartesian_product
611601
612- Given two ndarrays a(Z,x1,m), b(Z,x2,n), return c(Z,x1+x2,m*n), with
613- each row from C<b> appended to all rows from C<a>. C<Z> can be empty
614- but must be compatible; the shorter one will be "dummied up" from the
615- zero end as necessary.
602+ Given two ndarrays a(x1,m), b(x2,n), return c(x1+x2,m*n), with
603+ each row from C<b> appended to all rows from C<a>. If one of them is
604+ 1D, a dummy first index is added.
616605
617606=item * tensor
618607
619608Given a complex PDL, an LU decomposition array-ref as returned by
620609L</lu_decomp>, and the size of the tensor, returns the tensor.
621610
622- =item * wave_operator
623-
624- Given a Green tensor and number of dimension in the geometry, returns
625- a wave operator.
626-
627611=item * make_haydock
628612
629613=item * make_greenp
0 commit comments