Anonymous avatar Anonymous committed 1ce1073

Moving to Trunk

Comments (0)

Files changed (8)

+This is the MIT/X11 license:
+
+--------------
+Copyright 2002 Shlomi Fish.
+
+Permission is hereby granted, free of charge, to any person obtaining
+a copy of this software and associated documentation files (the
+"Software"), to deal in the Software without restriction, including
+without limitation the rights to use, copy, modify, merge, publish,
+distribute, sublicense, and/or sell copies of the Software, and to
+permit persons to whom the Software is furnished to do so, subject to
+the following conditions:
+
+The above copyright notice and this permission notice shall be
+included in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
+ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+OTHER DEALINGS IN THE SOFTWARE.
+
+COPYING
+MANIFEST
+Makefile.PL
+linear_solve.pd
+mytest.pl
+t/lin_solve.t
+
+# Makefile.PL for PDL::Ops module.
+
+# Use this as a template for the Makefile.PL for
+# any external PDL module.
+
+use PDL::Core::Dev;
+use ExtUtils::MakeMaker;
+
+PDL::Core::Dev->import();
+$package = ["linear_solve.pd",Solve,PDL::Lib::Linear::Solve];
+%hash = pdlpp_stdargs($package);
+#$hash{OBJECT} .= ' additional_Ccode$(OBJ_EXT) ';
+#$hash{clean}->{FILES} .= ' todelete_Ccode$(OBJ_EXT) ';
+$hash{'VERSION_FROM'} = 'linear_solve.pd';
+WriteMakefile(%hash);
+
+sub MY::postamble { pdlpp_postamble($package); }
+
+package MatInv;
+
+use strict;
+
+use PDL;
+
+require Exporter;
+
+use vars qw(@EXPORT @ISA);
+
+@EXPORT=qw(linear_solver);
+
+@ISA=qw(Exporter);
+
+my $eps = 0.00001;
+
+sub linear_solver
+{
+    my $coeffs = shift;
+    my $values = shift;
+
+    my @dims = dims($coeffs);
+    if ((scalar(@dims) != 2) || ($dims[0] != $dims[1]))
+    {
+        die "First argument is not a matrix";
+    }
+    
+    my @v_dims = dims($values);
+    if (@v_dims == 1)
+    {
+        $values = transpose($values);
+        @v_dims = dims($values);
+    }
+    if ($v_dims[1] != $dims[1])
+    {
+        die "Second argument does not match in length";
+    }
+
+    $coeffs = $coeffs->copy();
+
+    my $mr   = $dims[0];
+    my $mc   = $mr;
+    my $f;
+    my $try;
+
+    ROW: 
+    for(my $i = 0; $i < $mr; $i++) {
+        $try=$i;
+        # make diagonal element nonzero if possible
+        while(abs($coeffs->at($i,$try)) < $eps)
+        {
+            $try++;
+            if ($try >= $mr)
+            {
+                last ROW;
+            }
+        }
+        my ($temp, $mat);
+
+        foreach $mat ($coeffs, $values)
+        {
+            $temp = $mat->slice(":,$i")->copy();
+            $mat->slice(":,$i") .= $mat->slice(":,$try");
+            $mat->slice(":,$try") .= $temp;
+        }
+
+        $f = (1 / ($coeffs->at($i,$i)));
+        foreach $mat ($coeffs, $values)
+        {
+            $mat->slice(":,$i") *= $f;
+        }
+        
+        # subtract multiple of designated row from other rows
+
+        my @matrixes = 
+            (
+                [$coeffs, $coeffs->slice(":,$i")], 
+                [$values, $values->slice(":,$i")]
+            );
+            
+        for(my $j = 0; $j < $mr; $j++) {
+            next if $i == $j;
+            $f = $coeffs->at($i,$j);
+            foreach $mat (@matrixes)
+            {
+                $mat->[0]->slice(":,$j") -= $mat->[1] * $f;
+            }
+        }
+    }
+
+    return ($coeffs, $values);
+}
+
+1;
+
+$VERSION = 0.01;
+
+pp_addpm({ At => 'Top'}, << 'EOD');
+
+=head1 NAME
+
+PDL::Lib::Linear::Solve -- Linear Equations Set Solution or Matrix Inversion
+
+=head1 SYNOPSIS
+
+    use PDL;
+    use PDL::Lib::Linear::Solve;
+
+    my $coeffs = pdl([[1,1],[0,1]]);
+    my $values = pdl([[5],[2]]);
+    my ($new_coeffs, $new_values) = linear_solve($coeffs, $values);
+
+=head1 DESCRIPTION
+
+This package provides routines for solving a set of linear equations and
+inverting a matrix.
+
+=cut
+EOD
+
+pp_addpm({At => 'Bot'}, << 'EOD');
+
+=head1 AUTHOR
+
+Copyright (C) 2003 Shlomi Fish. There is no warranty. You are allowed
+to redistribute this software / documentation under certain
+conditions. For details, see the file COPYING in the PDL
+distribution. If this file is separated from the PDL distribution,
+the copyright notice should be included in the file.
+
+=cut
+
+EOD
+
+pp_def(
+    'linear_solve',
+    'Pars' => 'in_coeffs(incols,inrows); in_values(valcols,inrows); [o]coeffs(incols,inrows) ; [o]values(valcols,inrows);',
+    Code => << 'EOF' ,
+int mr, mc, i,j,col,try, valcols_size;
+$GENERIC() temp;
+$GENERIC() f;
+mr = $SIZE(inrows);
+mc = $SIZE(incols);
+valcols_size = $SIZE(valcols);
+for(i=0;i<mr;i++)
+{
+    for(col=0;col<mc;col++)
+    {
+        $coeffs(inrows => i,incols => col) =
+            $in_coeffs(inrows => i,incols => col);
+    }
+    for(col=0;col<valcols_size;col++)
+    {
+        $values(inrows => i,valcols => col) =
+            $in_values(inrows => i,valcols => col);
+    }
+}
+if (mr != mc)
+{
+    return;
+}
+for(i=0;i<mr;i++)
+{
+    try=i;
+    while($coeffs(incols => i, inrows => try) == 0)
+    {
+        try++;
+        if (try >= mr)
+        {
+            break;
+        }
+    }
+    if (try >= mr)
+    {
+        break;
+    }
+    /* Swap i and try */
+    for(col=0;col<mc;col++)
+    {
+        temp = $coeffs(inrows => i, incols => col);
+        $coeffs(inrows => i, incols => col) = $coeffs(inrows => try, incols => col);
+        $coeffs(inrows => try, incols => col) = temp;
+    }
+    for(col=0;col<valcols_size;col++)
+    {
+        temp = $values(inrows => i, valcols => col);
+        $values(inrows => i, valcols => col) = $values(inrows => try, valcols => col);
+        $values(inrows => try, valcols => col) = temp;        
+    }
+
+    f = 1 / $coeffs(inrows => i, incols => i);
+    for(col=0;col<mc;col++)
+    {
+        $coeffs(inrows => i, incols => col) *= f;
+    }
+    for(col=0;col<valcols_size;col++)
+    {
+        $values(inrows => i, valcols => col) *= f;
+    }
+
+    for(j=0;j<mr;j++)
+    {
+        if (i == j)
+        {
+            continue;
+        }
+        f = $coeffs(incols => i, inrows => j);
+        for(col=0;col<mc;col++)
+        {
+            $coeffs(incols => col, inrows => j) -= 
+                $coeffs(incols => col, inrows => i) * f;
+        }
+        for(col=0;col<valcols_size;col++)
+        {
+            $values(valcols => col, inrows => j) -= 
+                $values(valcols => col, inrows => i) * f;
+        }
+    }
+}
+EOF
+    Doc => <<'EOF',
+=for ref
+
+C<linear_solve> provides a solving algorithm for a set of linear equations.
+
+=for example
+
+    ($result_coeffs, $result_values) = linear_solve($coeffs, $values);
+    
+C<linear_solve> accepts two arguments - a square matrix representing the linear
+equation coefficients, and a tensors that contains their values. As a result
+it produces a matrix with the final coefficients (an identity matrix if 
+the coefficients form a regular matrix), and the final values.
+
+=cut
+
+EOF
+);
+
+my $rejects =  <<'EOF';
+loop(inrows) %{
+    loop(incols) %{
+        $coeffs() = $in_coeffs();
+    %}
+    loop(valcols) %{
+        $values() = $in_values();
+    %}
+%}
+EOF
+
+pp_addpm(<<'EOF');
+
+sub eye
+{
+    my $n = shift;
+    my $pdl = zeroes($n,$n);
+    $pdl->diagonal(0,1)++;
+    return $pdl;
+}
+
+sub myinv
+{
+    my $matrix = shift;
+    my @dims = dims($matrix);    
+    my $I = eye($dims[0]);
+    my ($out_coeffs, $out_values) = linear_solve($matrix, $I);
+    return $out_values;
+}
+
+=head2 myinv
+
+=for ref
+
+Inverts a matrix
+
+=for usage
+
+    $inverted_matrix = inv($matrix);
+
+=for example
+
+    $inverted_matrix = inv($matrix);
+
+=cut
+EOF
+
+pp_add_exported('', 'myinv');
+
+pp_done();
+
+

matrix_inversion.pl

+#!/usr/bin/perl -w
+
+use strict;
+
+use MatInv;
+use PDL;
+
+my ($a, $v);
+my ($new_a, $new_v);
+sub mytest
+{
+    print "\n\n-------------\n\n";
+print "\$a = " , $a;
+print "\$v = ", $v, "\n";    
+($new_a, $new_v) = linear_solver($a,$v);
+print "\$new_a = ", $new_a;
+print "\$new_v = ", $new_v;
+print "\n\n";
+    
+}
+
+sub eye
+{
+    my $n = shift;
+    my $pdl = zeroes($n,$n);
+    $pdl->diagonal(0,1)++;
+    return $pdl;
+}
+
+$a = pdl([[0,2,3], [4, 5 ,6], [7, 8, 9]]);
+$v = pdl([1, 5 ,10]);
+
+mytest();
+
+$a = pdl([[1,0,1],[0,1,0],[0,0,1]]);
+$v = pdl([3,2,1]);
+
+mytest();
+
+$a = pdl([[1,0,1],[0,1,0],[0,0,1]]);
+$v = eye(3);
+mytest();
+print $new_a*$new_v;
+
+$a = pdl([[1,0,1],[0,1,0],[0,0,1]]);
+$v = pdl([
+         [[1,2,3],[4,5,6],[7,8,9]],
+         [[10,20,30],[40,50,60],[70,80,90]],
+         [[100,200,300],[400,500,600],[700,800,900]]
+         ]);
+mytest();
+
+
+#!/usr/bin/perl -w -I./blib/lib/ -I./blib/arch/
+
+use strict;
+
+use PDL;
+use PDL::Lib::Linear::Solve;
+
+sub eye
+{
+    my $n = shift;
+    my $pdl = zeroes($n,$n);
+    $pdl->diagonal(0,1)++;
+    return $pdl;
+}
+
+# my $a = pdl([3,4,5]);
+# print PDL::Lib::Mylib::sumover2($a), "\n";
+
+my ($coeffs, $values);
+# my $coeffs = pdl([0,0,2], [1, 1,3],[0,1,9]);
+# my $values = pdl([1], [5],[7]);
+
+$coeffs = pdl([1, 1],[0, 1]);
+$values = eye(2);
+
+$coeffs = pdl([0,0,2], [1, 1,3],[0,1,9]);
+$values = eye(3);
+
+$coeffs = pdl([[1,1,1],[0,1,1],[1,1,0]]);
+$values = eye(3);
+
+$coeffs = pdl([[1,1,1],[1,1,0],[2,2,0]]);
+$values = eye(3);
+
+$coeffs = pdl([[1,0,1],[0,1,0],[0,0,1]]);
+$values = pdl([
+         [[1,2,3],[4,5,6],[7,8,9]],
+         [[10,20,30],[40,50,60],[70,80,90]],
+         [[100,200,300],[400,500,600],[700,800,900]]
+         ]);
+
+my ($out_coeffs, $out_values) = linear_solve($coeffs, $values);
+print "\$coeffs = ", $coeffs, "\n";
+print "\$out_coeffs = ", $out_coeffs, "\n";
+print "\$out_values = ", $out_values, "\n";
+print "c*inv(c) = ", $coeffs x $out_values, "\n";
+
+
+#!/usr/bin/perl -w
+
+use strict;
+
+use constant NUM_LINEAR_SOLVE_TESTS => 4;
+use constant NUM_INV_TESTS => 5;
+use Test::More tests => (1+NUM_LINEAR_SOLVE_TESTS+NUM_INV_TESTS);
+
+use PDL;
+
+BEGIN
+{
+    use_ok('PDL::Lib::Linear::Solve');
+}
+
+sub eye
+{
+    my $n = shift;
+    my $pdl = zeroes($n,$n);
+    $pdl->diagonal(0,1)++;
+    return $pdl;
+}
+
+sub tapprox {
+    my($a,$b) = @_;
+    my $c = abs($a-$b);
+    my $d = max($c);
+    $d < 0.01;
+}
+
+my @linear_solve_tests =
+(
+    {
+        'in_c' => [[1,1],[0,1]],
+        'in_v' => [[5],[2]],
+        'out_c' => eye(2),
+        'out_v' => [[3],[2]],
+    },
+    {
+        'in_c' => [[1,1,1],[0,1,1],[1,1,0]],
+        'in_v' => eye(3),
+        'out_c' => eye(3),
+        'out_v' => [[1,-1,0],[-1,1,1],[1,0,-1]],
+    },
+    {
+        # This is a singular matrix.
+        'in_c' => [[1,1,1],[1,1,0],[2,2,0]],
+        'in_v' => eye(3),
+        'out_c' => [[1,1,1],[0,0,-1],[0,0,-2]],
+        'out_v' => [[1,0,0],[-1,1,0],[-2,0,1]],
+    },
+    {
+        # Test the identity.
+        'in_c' => eye(2),
+        'in_v' => eye(2),
+        'out_c' => eye(2),
+        'out_v' => eye(2),
+    },
+);
+
+if (scalar(@linear_solve_tests) != NUM_LINEAR_SOLVE_TESTS)
+{
+    die "Incorrect number of linear solve tests!";
+}
+
+foreach my $t (@linear_solve_tests)
+{
+    my @out_array = linear_solve(map { pdl($_) } @$t{'in_c','in_v'});
+    ok(tapprox($out_array[0], pdl($t->{'out_c'})) &&
+       tapprox($out_array[1], pdl($t->{'out_v'})));
+}
+
+my @inv_tests =
+(
+    {
+        'in_mat' => [[1,0],[0,1]],
+    },
+    {
+        'in_mat' => [[0,1],[1,0]],
+    },
+    {
+        'in_mat' => [[1,1],[1,0]],
+    },
+    {
+        'in_mat' => [[1,1,1],[1,0,0],[0,2,0]],
+    },
+    {
+        'in_mat' => [[1,2,3],[4,5,6],[-100,2,100]],
+    },    
+);
+
+if (scalar(@inv_tests) != NUM_INV_TESTS)
+{
+    die "Incorrect number of inv tests!";
+}
+
+foreach my $t (@inv_tests)
+{
+    my $mat = pdl($t->{'in_mat'});
+    my @dims = dims($mat);
+    if ((scalar(@dims) != 2) || ($dims[0] != $dims[1]))
+    {
+        die "Incorrect input for inv:\n$mat\n";
+    }
+    my $result = myinv($mat);
+    ok(tapprox($mat x $result, eye($dims[0])));
+}
+
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.