Devuan fork of gpsd
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 
 
 
 
 

120 lines
2.7 KiB

  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. use Math::MatrixReal;
  13. use Math::Factor::XS qw(prime_factors);
  14. use Getopt::Std;
  15. $opt{n} = 1;
  16. getopts('hsdn:', \%opt);
  17. $usage .= "easyinverse.pl [-hsd] [-n N]\n";
  18. $usage .= " -h this message\n";
  19. $usage .= " -s symmetric\n";
  20. $usage .= " -d diagonally dominant\n";
  21. $usage .= " -n N number to generate (DEF 1)\n";
  22. if ($opt{h}) {
  23. print $usage;
  24. exit;
  25. }
  26. sub randint{ # but not 0
  27. my $min = shift or -5;
  28. my $max = shift or +5;
  29. my $ansa = 0;
  30. while ($ansa == 0) {
  31. my $r = rand;
  32. $r *= ($max - $min + 1);
  33. my $i = int($r);
  34. $ansa = $min + $i;
  35. }
  36. return $ansa;
  37. }
  38. $size = 4;
  39. sub printmateq {
  40. my $lbl = shift;
  41. my $mat = shift;
  42. my $round = shift;
  43. if ($round) {
  44. for my $r (1..$size) {
  45. for my $c (1..$size) {
  46. my $elt = $mat->element($r,$c);
  47. my $pelt = abs($elt);
  48. my $prnd = int($pelt * $round + 0.5) / $round;
  49. my $rnd = $prnd * ($elt < 0 ? -1 : 1);
  50. $mat->assign($r, $c, $rnd);
  51. }}
  52. }
  53. my @rows;
  54. for my $r (1..$size) {
  55. my @row = map {$mat->element($r, $_)} (1..$size);
  56. my $rowstr = '{' . (join ',', @row) . '}';
  57. push @rows, $rowstr;
  58. }
  59. my $matstr = '{' . (join ',', @rows) . '}';
  60. $matstr =~ s!\.?000+\d+!!g;
  61. print "$lbl = $matstr;\n";
  62. }
  63. for (1..$opt{n}) { # requested number of matrices
  64. $M = new Math::MatrixReal($size,$size);
  65. while (1) {
  66. for $r (1..$size) {
  67. if ($opt{d}) { $M->assign($r,$r, randint(5,10)) }
  68. else { $M->assign($r,$r, randint(-3, 3)) }
  69. for $c ($r+1..$size) {
  70. $i = randint(-5,5);
  71. $M->assign($r, $c, $i);
  72. $M->assign($c, $r, ($opt{s} ? $i : randint(-5,5)));
  73. }
  74. }
  75. $det = int($M->det + 0.5);
  76. next if $det == 0;
  77. @factors = prime_factors(abs($det));
  78. $all_good_factors = 1; # cross your fingers!
  79. for $f (@factors) {
  80. if ($f != 2 || $f != 5) {
  81. $all_good_factors = 0;
  82. last;
  83. }
  84. }
  85. next unless $all_good_factors;
  86. # test for positive semidefinite ~ all positive eigenvalues
  87. $evals = $M->sym_eigenvalues;
  88. $all_positive_evals = 1; # keep hope alive!
  89. for $r (1..$size) {
  90. if ($evals->element($r,1) <= 0) {
  91. $all_positive_evals = 0;
  92. last;
  93. }
  94. }
  95. last if $all_positive_evals;
  96. }
  97. printmateq(".mat", $M);
  98. printmateq(".inv", $M->inverse, 1000);
  99. }