1#! /bin/perl
2#
3# From: Reuben Settergren <rjs@jhu.edu>
4#
5# I wrote a perl script to randomly generate matrices with smallish integers,
6# and test for a suitable determinant (typically +/-1), and ouput the matrix
7# and its inverse.
8#
9# There are switches to make it -s[ymmetric] and/or -d[iagonally strong],
10# which are necessary and good properties for covariance matrices. You may
11# need to grab Math::MatrixReal and Math::Factor::XS from cpan.
12
13use Math::MatrixReal;
14use Math::Factor::XS qw(prime_factors);
15use Getopt::Std;
16$opt{n} = 1;
17getopts('hsdn:', \%opt);
18
19$usage .= "easyinverse.pl [-hsd] [-n N]\n";
20$usage .= "  -h     this message\n";
21$usage .= "  -s     symmetric\n";
22$usage .= "  -d     diagonally dominant\n";
23$usage .= "  -n N   number to generate (DEF 1)\n";
24if ($opt{h}) {
25  print $usage;
26  exit;
27}
28
29sub randint{ # but not 0
30  my $min = shift or -5;
31  my $max = shift or +5;
32
33  my $ansa = 0;
34  while ($ansa == 0) {
35    my $r = rand;
36    $r *= ($max - $min + 1);
37    my $i = int($r);
38    $ansa = $min + $i;
39  }
40  return $ansa;
41}
42
43$size = 4;
44
45sub printmateq {
46  my $lbl = shift;
47  my $mat = shift;
48  my $round = shift;
49
50  if ($round) {
51    for my $r (1..$size) {
52    for my $c (1..$size) {
53      my $elt = $mat->element($r,$c);
54      my $pelt = abs($elt);
55      my $prnd = int($pelt * $round + 0.5) / $round;
56      my $rnd = $prnd * ($elt < 0 ? -1 : 1);
57      $mat->assign($r, $c, $rnd);
58    }}
59  }
60
61  my @rows;
62  for my $r (1..$size) {
63    my @row = map {$mat->element($r, $_)} (1..$size);
64    my $rowstr = '{' . (join ',', @row) . '}';
65    push @rows, $rowstr;
66  }
67  my $matstr = '{' . (join ',', @rows) . '}';
68  $matstr =~ s!\.?000+\d+!!g;
69  print "$lbl = $matstr;\n";
70}
71
72
73
74for (1..$opt{n}) { # requested number of matrices
75
76  $M = new Math::MatrixReal($size,$size);
77
78  while (1) {
79    for $r (1..$size) {
80      if ($opt{d}) { $M->assign($r,$r, randint(5,10)) }
81      else         { $M->assign($r,$r, randint(-3, 3)) }
82
83      for $c ($r+1..$size) {
84	$i = randint(-5,5);
85	$M->assign($r, $c, $i);
86	$M->assign($c, $r, ($opt{s} ? $i : randint(-5,5)));
87      }
88    }
89
90    $det = int($M->det + 0.5);
91    next if $det == 0;
92    @factors = prime_factors(abs($det));
93    $all_good_factors = 1; # cross your fingers!
94    for $f (@factors) {
95      if ($f != 2 || $f != 5) {
96	$all_good_factors = 0;
97	last;
98      }
99    }
100    next unless $all_good_factors;
101
102    # test for positive semidefinite ~ all positive eigenvalues
103    $evals = $M->sym_eigenvalues;
104    $all_positive_evals = 1; # keep hope alive!
105    for $r (1..$size) {
106      if ($evals->element($r,1) <= 0) {
107	$all_positive_evals = 0;
108	last;
109      }
110    }
111    last if $all_positive_evals;
112  }
113
114  printmateq(".mat", $M);
115  printmateq(".inv", $M->inverse, 1000);
116
117}
118
119
120