diff options
-rw-r--r-- | devtools/README | 5 | ||||
-rwxr-xr-x | devtools/easyinverse.pl | 119 |
2 files changed, 124 insertions, 0 deletions
diff --git a/devtools/README b/devtools/README index 41886497..390d29ea 100644 --- a/devtools/README +++ b/devtools/README @@ -34,6 +34,11 @@ which will be executed then. Very simple wrapper around scons clean, build, and check, saving output in a file named with "git describe". +== easyinverse.pl == + +Generate random pairirs of mutually inverse covariance matrices +suitable for testing a matrix-inversion algorithm. + == editcomment == Give this a commit-ID specification. It will edit the associated comment. diff --git a/devtools/easyinverse.pl b/devtools/easyinverse.pl new file mode 100755 index 00000000..29d42ebd --- /dev/null +++ b/devtools/easyinverse.pl @@ -0,0 +1,119 @@ +#! /bin/perl +# +# From: Reuben Settergren <rjs@jhu.edu> +# +# I wrote a perl script to randomly generate matrices with smallish integers, +# and test for a suitable determinant (typically +/-1), and ouput the matrix +# and its inverse. +# +# There are switches to make it -s[ymmetric] and/or -d[iagonally strong], +# which are necessary and good properties for covariance matrices. You may +# need to grab Math::MatrixReal and Math::Factor::XS from cpan. + +use Math::MatrixReal; +use Math::Factor::XS qw(prime_factors); +use Getopt::Std; +$opt{n} = 1; +getopts('hsdn:', \%opt); + +$usage .= "easyinverse.pl [-hsd] [-n N]\n"; +$usage .= " -h this message\n"; +$usage .= " -s symmetric\n"; +$usage .= " -d diagonally dominant\n"; +$usage .= " -n N number to generate (DEF 1)\n"; +if ($opt{h}) { + print $usage; + exit; +} + +sub randint{ # but not 0 + my $min = shift or -5; + my $max = shift or +5; + + my $ansa = 0; + while ($ansa == 0) { + my $r = rand; + $r *= ($max - $min + 1); + my $i = int($r); + $ansa = $min + $i; + } + return $ansa; +} + +$size = 4; + +sub printmateq { + my $lbl = shift; + my $mat = shift; + my $round = shift; + + if ($round) { + for my $r (1..$size) { + for my $c (1..$size) { + my $elt = $mat->element($r,$c); + my $pelt = abs($elt); + my $prnd = int($pelt * $round + 0.5) / $round; + my $rnd = $prnd * ($elt < 0 ? -1 : 1); + $mat->assign($r, $c, $rnd); + }} + } + + my @rows; + for my $r (1..$size) { + my @row = map {$mat->element($r, $_)} (1..$size); + my $rowstr = '{' . (join ',', @row) . '}'; + push @rows, $rowstr; + } + my $matstr = '{' . (join ',', @rows) . '}'; + $matstr =~ s!\.?000+\d+!!g; + print "$lbl = $matstr;\n"; +} + + + +for (1..$opt{n}) { # requested number of matrices + + $M = new Math::MatrixReal($size,$size); + + while (1) { + for $r (1..$size) { + if ($opt{d}) { $M->assign($r,$r, randint(5,10)) } + else { $M->assign($r,$r, randint(-3, 3)) } + + for $c ($r+1..$size) { + $i = randint(-5,5); + $M->assign($r, $c, $i); + $M->assign($c, $r, ($opt{s} ? $i : randint(-5,5))); + } + } + + $det = int($M->det + 0.5); + next if $det == 0; + @factors = prime_factors(abs($det)); + $all_good_factors = 1; # cross your fingers! + for $f (@factors) { + if ($f != 2 || $f != 5) { + $all_good_factors = 0; + last; + } + } + next unless $all_good_factors; + + # test for positive semidefinite ~ all positive eigenvalues + $evals = $M->sym_eigenvalues; + $all_positive_evals = 1; # keep hope alive! + for $r (1..$size) { + if ($evals->element($r,1) <= 0) { + $all_positive_evals = 0; + last; + } + } + last if $all_positive_evals; + } + + printmateq(".mat", $M); + printmateq(".inv", $M->inverse, 1000); + +} + + |