[DRE-commits] [ruby-extendmatrix] 01/07: Imported Upstream version 0.3.1
Balint Reczey
rbalint at moszumanska.debian.org
Fri May 22 15:28:03 UTC 2015
This is an automated email from the git hooks/post-receive script.
rbalint pushed a commit to branch master
in repository ruby-extendmatrix.
commit dccbdd5dc5f894001a177711cf9609a7f9a7700e
Author: Balint Reczey <balint at balintreczey.hu>
Date: Fri May 22 16:15:29 2015 +0200
Imported Upstream version 0.3.1
---
History.txt | 27 ++
LICENSE.txt | 175 +++++++
Manifest.txt | 9 +
ORIGINAL_README.txt | 22 +
README.txt | 66 +++
Rakefile | 13 +
data.tar.gz.sig | Bin 0 -> 256 bytes
lib/extendmatrix.rb | 1137 +++++++++++++++++++++++++++++++++++++++++++++
metadata.gz.sig | Bin 0 -> 256 bytes
metadata.yml | 130 ++++++
spec/extendmatrix_spec.rb | 423 +++++++++++++++++
spec/spec.opts | 1 +
12 files changed, 2003 insertions(+)
diff --git a/History.txt b/History.txt
new file mode 100644
index 0000000..d02f5dd
--- /dev/null
+++ b/History.txt
@@ -0,0 +1,27 @@
+=== 0.3.1 / 2010-08-16
+* Bug fix: method mssq returns only the last row squares
+* Added method diagonal
+
+=== 0.3.0 / 2010-08-16
+* Added new methods, inspired on SPSS matrix methods: e_mult, e_quo, mssq, eigen, eigenpairs, ln, sqrt, sscp
+
+=== 0.2.3 / 2010-05-05
+* Added LICENSE.txt
+* Complete README.txt
+* Added Vector#magnitude
+* Bug fix on Vector#normalize (quo doesn't defined)
+
+=== 0.2.2 / 2010-05-05
+* Bug fix on Matrix#dup. Doesn't works on Ruby1.9.2
+
+=== 0.2.1 / 2010-05-05
+
+* Works on Ruby 1.9.2. Fixed conflicts between old and new Matrix implementation
+
+=== 0.2.0 / 2010-05-04
+* More rubyesque implementation of some methods
+* Matrix#new is not altered and previously extended #new is called #new_with_values
+
+=== 0.1.0 / 2010-05-04
+
+* First gem version
diff --git a/LICENSE.txt b/LICENSE.txt
new file mode 100644
index 0000000..e97bf01
--- /dev/null
+++ b/LICENSE.txt
@@ -0,0 +1,175 @@
+ Apache License
+ Version 2.0, January 2004
+ http://www.apache.org/licenses/
+
+ TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION
+
+ 1. Definitions.
+
+ "License" shall mean the terms and conditions for use, reproduction,
+ and distribution as defined by Sections 1 through 9 of this document.
+
+ "Licensor" shall mean the copyright owner or entity authorized by
+ the copyright owner that is granting the License.
+
+ "Legal Entity" shall mean the union of the acting entity and all
+ other entities that control, are controlled by, or are under common
+ control with that entity. For the purposes of this definition,
+ "control" means (i) the power, direct or indirect, to cause the
+ direction or management of such entity, whether by contract or
+ otherwise, or (ii) ownership of fifty percent (50%) or more of the
+ outstanding shares, or (iii) beneficial ownership of such entity.
+
+ "You" (or "Your") shall mean an individual or Legal Entity
+ exercising permissions granted by this License.
+
+ "Source" form shall mean the preferred form for making modifications,
+ including but not limited to software source code, documentation
+ source, and configuration files.
+
+ "Object" form shall mean any form resulting from mechanical
+ transformation or translation of a Source form, including but
+ not limited to compiled object code, generated documentation,
+ and conversions to other media types.
+
+ "Work" shall mean the work of authorship, whether in Source or
+ Object form, made available under the License, as indicated by a
+ copyright notice that is included in or attached to the work
+ (an example is provided in the Appendix below).
+
+ "Derivative Works" shall mean any work, whether in Source or Object
+ form, that is based on (or derived from) the Work and for which the
+ editorial revisions, annotations, elaborations, or other modifications
+ represent, as a whole, an original work of authorship. For the purposes
+ of this License, Derivative Works shall not include works that remain
+ separable from, or merely link (or bind by name) to the interfaces of,
+ the Work and Derivative Works thereof.
+
+ "Contribution" shall mean any work of authorship, including
+ the original version of the Work and any modifications or additions
+ to that Work or Derivative Works thereof, that is intentionally
+ submitted to Licensor for inclusion in the Work by the copyright owner
+ or by an individual or Legal Entity authorized to submit on behalf of
+ the copyright owner. For the purposes of this definition, "submitted"
+ means any form of electronic, verbal, or written communication sent
+ to the Licensor or its representatives, including but not limited to
+ communication on electronic mailing lists, source code control systems,
+ and issue tracking systems that are managed by, or on behalf of, the
+ Licensor for the purpose of discussing and improving the Work, but
+ excluding communication that is conspicuously marked or otherwise
+ designated in writing by the copyright owner as "Not a Contribution."
+
+ "Contributor" shall mean Licensor and any individual or Legal Entity
+ on behalf of whom a Contribution has been received by Licensor and
+ subsequently incorporated within the Work.
+
+ 2. Grant of Copyright License. Subject to the terms and conditions of
+ this License, each Contributor hereby grants to You a perpetual,
+ worldwide, non-exclusive, no-charge, royalty-free, irrevocable
+ copyright license to reproduce, prepare Derivative Works of,
+ publicly display, publicly perform, sublicense, and distribute the
+ Work and such Derivative Works in Source or Object form.
+
+ 3. Grant of Patent License. Subject to the terms and conditions of
+ this License, each Contributor hereby grants to You a perpetual,
+ worldwide, non-exclusive, no-charge, royalty-free, irrevocable
+ (except as stated in this section) patent license to make, have made,
+ use, offer to sell, sell, import, and otherwise transfer the Work,
+ where such license applies only to those patent claims licensable
+ by such Contributor that are necessarily infringed by their
+ Contribution(s) alone or by combination of their Contribution(s)
+ with the Work to which such Contribution(s) was submitted. If You
+ institute patent litigation against any entity (including a
+ cross-claim or counterclaim in a lawsuit) alleging that the Work
+ or a Contribution incorporated within the Work constitutes direct
+ or contributory patent infringement, then any patent licenses
+ granted to You under this License for that Work shall terminate
+ as of the date such litigation is filed.
+
+ 4. Redistribution. You may reproduce and distribute copies of the
+ Work or Derivative Works thereof in any medium, with or without
+ modifications, and in Source or Object form, provided that You
+ meet the following conditions:
+
+ (a) You must give any other recipients of the Work or
+ Derivative Works a copy of this License; and
+
+ (b) You must cause any modified files to carry prominent notices
+ stating that You changed the files; and
+
+ (c) You must retain, in the Source form of any Derivative Works
+ that You distribute, all copyright, patent, trademark, and
+ attribution notices from the Source form of the Work,
+ excluding those notices that do not pertain to any part of
+ the Derivative Works; and
+
+ (d) If the Work includes a "NOTICE" text file as part of its
+ distribution, then any Derivative Works that You distribute must
+ include a readable copy of the attribution notices contained
+ within such NOTICE file, excluding those notices that do not
+ pertain to any part of the Derivative Works, in at least one
+ of the following places: within a NOTICE text file distributed
+ as part of the Derivative Works; within the Source form or
+ documentation, if provided along with the Derivative Works; or,
+ within a display generated by the Derivative Works, if and
+ wherever such third-party notices normally appear. The contents
+ of the NOTICE file are for informational purposes only and
+ do not modify the License. You may add Your own attribution
+ notices within Derivative Works that You distribute, alongside
+ or as an addendum to the NOTICE text from the Work, provided
+ that such additional attribution notices cannot be construed
+ as modifying the License.
+
+ You may add Your own copyright statement to Your modifications and
+ may provide additional or different license terms and conditions
+ for use, reproduction, or distribution of Your modifications, or
+ for any such Derivative Works as a whole, provided Your use,
+ reproduction, and distribution of the Work otherwise complies with
+ the conditions stated in this License.
+
+ 5. Submission of Contributions. Unless You explicitly state otherwise,
+ any Contribution intentionally submitted for inclusion in the Work
+ by You to the Licensor shall be under the terms and conditions of
+ this License, without any additional terms or conditions.
+ Notwithstanding the above, nothing herein shall supersede or modify
+ the terms of any separate license agreement you may have executed
+ with Licensor regarding such Contributions.
+
+ 6. Trademarks. This License does not grant permission to use the trade
+ names, trademarks, service marks, or product names of the Licensor,
+ except as required for reasonable and customary use in describing the
+ origin of the Work and reproducing the content of the NOTICE file.
+
+ 7. Disclaimer of Warranty. Unless required by applicable law or
+ agreed to in writing, Licensor provides the Work (and each
+ Contributor provides its Contributions) on an "AS IS" BASIS,
+ WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or
+ implied, including, without limitation, any warranties or conditions
+ of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A
+ PARTICULAR PURPOSE. You are solely responsible for determining the
+ appropriateness of using or redistributing the Work and assume any
+ risks associated with Your exercise of permissions under this License.
+
+ 8. Limitation of Liability. In no event and under no legal theory,
+ whether in tort (including negligence), contract, or otherwise,
+ unless required by applicable law (such as deliberate and grossly
+ negligent acts) or agreed to in writing, shall any Contributor be
+ liable to You for damages, including any direct, indirect, special,
+ incidental, or consequential damages of any character arising as a
+ result of this License or out of the use or inability to use the
+ Work (including but not limited to damages for loss of goodwill,
+ work stoppage, computer failure or malfunction, or any and all
+ other commercial damages or losses), even if such Contributor
+ has been advised of the possibility of such damages.
+
+ 9. Accepting Warranty or Additional Liability. While redistributing
+ the Work or Derivative Works thereof, You may choose to offer,
+ and charge a fee for, acceptance of support, warranty, indemnity,
+ or other liability obligations and/or rights consistent with this
+ License. However, in accepting such obligations, You may act only
+ on Your own behalf and on Your sole responsibility, not on behalf
+ of any other Contributor, and only if You agree to indemnify,
+ defend, and hold each Contributor harmless for any liability
+ incurred by, or claims asserted against, such Contributor by reason
+ of your accepting any such warranty or additional liability.
+
diff --git a/Manifest.txt b/Manifest.txt
new file mode 100644
index 0000000..992ad27
--- /dev/null
+++ b/Manifest.txt
@@ -0,0 +1,9 @@
+History.txt
+LICENSE.txt
+Manifest.txt
+ORIGINAL_README.txt
+README.txt
+Rakefile
+lib/extendmatrix.rb
+spec/extendmatrix_spec.rb
+spec/spec.opts
diff --git a/ORIGINAL_README.txt b/ORIGINAL_README.txt
new file mode 100755
index 0000000..c96491a
--- /dev/null
+++ b/ORIGINAL_README.txt
@@ -0,0 +1,22 @@
+Extensions to the Ruby Matrix module
+====================================
+
+This README is a small description of the work done by Cosmin Bonchis as a
+Google Summer of Code 2007 project for Ruby Central Inc.
+
+The project consists of some enhancements to the Ruby "Matrix" module and includes: LU and QR (Householder, Givens, Gram Schmidt, Hessenberg) decompositions, bidiagonalization, eigenvalue and eigenvector calculations.
+
+This archive contains in extendmatrix.rb file the source code of the project, an implementation of mapcar used in extending matrix, and all the tests files in the "tests" directory.
+
+The code can also be found on the RubyForge repository at http://matrix.rubyforge.org/svn/trunk/ or the project's SVN repository can be checked out through anonymous access with the following command(s).
+
+svn checkout svn://rubyforge.org/var/svn/matrix
+svn checkout http://matrix.rubyforge.org/svn/trunk/
+
+
+Relevant URLs:
+==============
+
+Project sources:
+http://matrix.rubyforge.org/svn/trunk/
+
diff --git a/README.txt b/README.txt
new file mode 100644
index 0000000..b4f65f9
--- /dev/null
+++ b/README.txt
@@ -0,0 +1,66 @@
+= extendmatrix
+
+* http://github.com/clbustos/extendmatrix
+
+== DESCRIPTION:
+
+The project consists of some enhancements to the Ruby "Matrix" module and includes: LU and QR (Householder, Givens, Gram Schmidt, Hessenberg) decompositions, bidiagonalization, eigenvalue and eigenvector calculations.
+Include some aditional code to obtains marginal for rows and columns.
+
+Original code from http://rubyforge.org/projects/matrix/ , done by Cosmin Bonchis as a Google Summer of Code 2007 project for Ruby Central Inc.
+
+Gem, github repository and current version manteined by Claudio Bustos.
+
+== SYNOPSIS:
+
+ require 'extendmatrix'
+ v = Vector[1, 2, 3, 4]
+
+ v.magnitude # => 5.47722557505166
+
+ v.normalize # => Vector[0.182574185835055, 0.365148371670111, 0.547722557505166, 0.730296743340221]
+
+ # Backport from Ruby 1.9.2+
+ m = Matrix.build(4, 3){|i, j| i * 3 + j}
+ # => Matrix[[0, 1, 2], [3, 4, 5], [6, 7, 8], [9, 10, 11]]
+
+ # You could modify the matrix
+ m = Matrix.I(3)
+ m[0,1]=m[1,0]=0.5
+ m[0,2]=m[2,0]=0.7
+ m[1,2]=m[2,1]=0.9
+ m # => Matrix[[1, 0.5, 0.7], [0.5, 1, 0.9], [0.7, 0.9, 1]]
+
+ # Eigenvalues and EigenVectors. PCA at sight :)
+
+ m.eigenvaluesJacobi => Vector[0.523942339006665, 0.0632833995384682, 2.41277426145487]
+ m.cJacobiV
+ # => Matrix[[0.818814082563014, 0.249617871497675, 0.516947208547894], [-0.550168858227442, 0.598307531004925, 0.58253096551128], [-0.163883268313767, -0.761392813580323, 0.62723461144538]]
+== REQUIREMENTS:
+
+* Only Ruby
+
+== INSTALL:
+
+* sudo gem install matrix-extensions
+
+== LICENSE:
+
+
+
+Copyright [2007] Cosmin Bonchis
+Copyright [2010] Claudio Bustos
+
+Licensed under the Apache License, Version 2.0 (the "License"); you may not use
+this file except in compliance with the License. You may obtain a copy of the
+License at
+
+http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software distributed
+under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR
+CONDITIONS OF ANY KIND, either express or implied. See the License for the
+specific language governing permissions and limitations under the License.
+
+
+See LICENSE.txt for more details
\ No newline at end of file
diff --git a/Rakefile b/Rakefile
new file mode 100644
index 0000000..dc5e941
--- /dev/null
+++ b/Rakefile
@@ -0,0 +1,13 @@
+# -*- ruby -*-
+$:.unshift(File.dirname(__FILE__)+"/lib")
+require 'rubygems'
+require 'hoe'
+require 'extendmatrix.rb'
+Hoe.plugin :git
+Hoe.spec 'extendmatrix' do
+ self.rubyforge_name = 'ruby-statsample'
+ self.version = Matrix::EXTENSION_VERSION
+ self.developer('Cosmin Bonchis', 'cbonchis_info.uvt.ro')
+end
+
+# vim: syntax=ruby
diff --git a/data.tar.gz.sig b/data.tar.gz.sig
new file mode 100644
index 0000000..8097a29
Binary files /dev/null and b/data.tar.gz.sig differ
diff --git a/lib/extendmatrix.rb b/lib/extendmatrix.rb
new file mode 100644
index 0000000..d60c3a0
--- /dev/null
+++ b/lib/extendmatrix.rb
@@ -0,0 +1,1137 @@
+require 'rational'
+require 'matrix'
+class Vector
+ include Enumerable
+ # fix for Vector#coerce on Ruby 1.8.x
+ if RUBY_VERSION<="1.9.0"
+ alias_method :old_coerce, :coerce
+ def coerce(other)
+ case other
+ when Numeric
+ return Matrix::Scalar.new(other), self
+ else
+ raise TypeError, "#{self.class} can't be coerced into #{other.class}"
+ end
+ end
+ end
+
+ module Norm
+ def self.sqnorm(obj, p)
+ obj.inject(0) {|ac,x| ac+(x**p)}
+ end
+ end
+
+ alias :length :size
+ alias :index :[]
+ #
+ # Returns the value of an index vector or
+ # a Vector with the values of a range
+ # v = Vector[1, 2, 3, 4]
+ # v[0] => 1
+ # v[0..2] => Vector[1, 2, 3]
+ #
+ def [](i)
+ case i
+ when Range
+ Vector[*to_a.slice(i)]
+ else
+ index(i)
+ end
+ end
+
+ # Magnitude or length of the vector
+ # Equal to sqrt(sum(x_i^2))
+ def magnitude
+ Math::sqrt(to_a.inject(0) {|ac, v| ac+(v**2)})
+ end
+
+ #
+ # Sets a vector value/(range of values) with a new value/(values from a vector)
+ # v = Vector[1, 2, 3]
+ # v[2] = 9 => Vector[1, 2, 9]
+ # v[1..2] = Vector[9, 9, 9, 9, 9] => v: Vector[1, 9, 9]
+ #
+ def []=(i, v)
+ case i
+ when Range
+ (self.size..i.begin - 1).each{|e| self[e] = 0} # self.size must be in the first place because the size of self can be modified
+ [v.size, i.entries.size].min.times {|e| self[e + i.begin] = v[e]}
+ (v.size + i.begin .. i.end).each {|e| self[e] = 0}
+ else
+ @elements[i]=v
+ end
+ end
+
+ class << self
+ #
+ # Returns a concatenated Vector
+ #
+ # Vector.concat(Vector[1,2,3], Vector[4,5,6])
+ # => Vector[1, 2, 3, 4, 5, 6]
+ #
+ def concat(*args)
+ Vector.elements(args.inject([]){|ac,v| ac+v.to_a}, false)
+ end
+ end
+
+ #
+ # Changes the elements of vector and returns a Vector
+ #
+ def collect!
+ els = @elements.collect! {|v| yield(v)}
+ Vector.elements(els, false)
+ end
+
+ #
+ # Iterates the elements of a vector
+ #
+ def each
+ (0...size).each {|i| yield(self[i])}
+ nil
+ end
+
+ #
+ # Returns the maximum element of a vector
+ #
+ def max
+ to_a.max
+ end
+
+ #
+ # Returns the minimum element of a vector
+ #
+ def min
+ to_a.min
+ end
+ #
+ # Returns the sum of values of the vector
+ #
+ def sum
+ to_a.inject(&:+)
+ end
+ #
+ # Returns the p-norm of a vector
+ #
+ def norm(p = 2)
+ Norm.sqnorm(self, p) ** (1.quo(p))
+ end
+
+ #
+ # Returns the infinite-norm
+ #
+ def norm_inf
+ [min.abs, max.abs].max
+ end
+
+ #
+ # Returns a slice of vector
+ #
+ def slice(*args)
+ Vector[*to_a.slice(*args)]
+ end
+
+ def slice_set(v, b, e)
+ for i in b..e
+ self[i] = v[i-b]
+ end
+ end
+
+ #
+ # Sets a slice of vector
+ #
+ def slice=(args)
+ case args[1]
+ when Range
+ slice_set(args[0], args[1].begin, args[1].last)
+ else
+ slice_set(args[0], args[1], args[2])
+ end
+ end
+
+ #
+ # Return the vector divided by a scalar
+ #
+ def / (c)
+ map {|e| e.quo(c)}
+ end
+ alias :quo :/
+ #
+ # Return the matrix column coresponding to the vector transpose
+ #
+ def transpose
+ Matrix[self.to_a]
+ end
+
+ alias :t :transpose
+
+ #
+ # Computes the Householder vector (MC, Golub, p. 210, algorithm 5.1.1)
+ #
+ def house
+ s = self[1..length-1]
+ sigma = s.inner_product(s)
+ v = clone; v[0] = 1
+ if sigma == 0
+ beta = 0
+ else
+ mu = Math.sqrt(self[0] ** 2 + sigma)
+ if self[0] <= 0
+ v[0] = self[0] - mu
+ else
+ v[0] = - sigma.quo(self[0] + mu)
+ end
+ v2 = v[0] ** 2
+ beta = 2 * v2.quo(sigma + v2)
+ v /= v[0]
+ end
+ return v, beta
+ end
+
+ #
+ # Projection operator
+ # (http://en.wikipedia.org/wiki/Gram-Schmidt_process#The_Gram.E2.80.93Schmidt_process)
+ #
+ def proj(v)
+ vp = v.inner_product(self)
+ vp = Float vp if vp.is_a?(Integer)
+ self * (vp / inner_product(self))
+ end
+
+ #
+ # Return the vector normalized
+ #
+ def normalize
+ self.quo(self.norm)
+ end
+
+ #
+ # Stabilized Gram-Schmidt process
+ # (http://en.wikipedia.org/wiki/Gram-Schmidt_process#Algorithm)
+ #
+ def self.gram_schmidt(*vectors)
+ v = vectors.clone
+ for j in 0...v.size
+ for i in 0..j-1
+ v[j] -= v[i] * v[j].inner_product(v[i])
+ end
+ v[j] /= v[j].norm
+ end
+ v
+ end
+end
+
+class Matrix
+
+ EXTENSION_VERSION="0.3.1"
+ include Enumerable
+
+ attr_reader :rows, :wrap
+ @wrap = nil
+ #
+ # For invoking a method
+ #
+ def initialize_old(init_method, *argv) #:nodoc:
+ self.send(init_method, *argv)
+ end
+ # Return an empty matrix on n=0
+ # else, returns I(n)
+ def self.robust_I(n)
+ n>0 ? self.I(n) : self.[]
+ end
+ alias :ids :[]
+ #
+ # Return a value or a vector/matrix of values depending
+ # if the indexes are ranges or not
+ # m = Matrix.new_with_value(4, 3){|i, j| i * 3 + j}
+ # m: 0 1 2
+ # 3 4 5
+ # 6 7 8
+ # 9 10 11
+ # m[1, 2] => 5
+ # m[3,1..2] => Vector[10, 11]
+ # m[0..1, 0..2] => Matrix[[0, 1, 2], [3, 4, 5]]
+ #
+ def [](i, j)
+ case i
+ when Range
+ case j
+ when Range
+ Matrix[*i.collect{|l| self.row(l)[j].to_a}]
+ else
+ column(j)[i]
+ end
+ else
+ case j
+ when Range
+ row(i)[j]
+ else
+ ids(i, j)
+ end
+ end
+ end
+
+
+
+
+ #
+ # Set the values of a matrix
+ # m = Matrix.build(3, 3){|i, j| i * 3 + j}
+ # m: 0 1 2
+ # 3 4 5
+ # 6 7 8
+ # m[1, 2] = 9 => Matrix[[0, 1, 2], [3, 4, 9], [6, 7, 8]]
+ # m[2,1..2] = Vector[8, 8] => Matrix[[0, 1, 2], [3, 8, 8], [6, 7, 8]]
+ # m[0..1, 0..1] = Matrix[[0, 0, 0],[0, 0, 0]]
+ # => Matrix[[0, 0, 2], [0, 0, 8], [6, 7, 8]]
+ #
+ def []=(i, j, v)
+ case i
+ when Range
+ if i.entries.size == 1
+ self[i.begin, j] = (v.is_a?(Matrix) ? v.row(0) : v)
+ else
+ case j
+ when Range
+ if j.entries.size == 1
+ self[i, j.begin] = (v.is_a?(Matrix) ? v.column(0) : v)
+ else
+ i.each{|l| self.row= l, v.row(l - i.begin), j}
+ end
+ else
+ self.column= j, v, i
+ end
+ end
+ else
+ case j
+ when Range
+ if j.entries.size == 1
+ self[i, j.begin] = (v.is_a?(Vector) ? v[0] : v)
+ else
+ self.row= i, v, j
+ end
+ else
+ @rows[i][j] = (v.is_a?(Vector) ? v[0] : v)
+
+ end
+ end
+ end
+
+ #
+ # Return a duplicate matrix, with all elements copied
+ #
+ def dup
+ (self.class).rows(self.rows.dup)
+ end
+
+ def initialize_copy(orig)
+ init_rows(orig.rows, true)
+ self.wrap=(orig.wrap)
+ end
+
+
+ class << self
+ if !self.respond_to? :build
+ #
+ # Creates a matrix <tt>n</tt> x <tt>m</tt>
+ # If you provide a block, it will be used to set the values.
+ # If not, <tt>val</tt> will be used
+ #
+
+ def build(n,m,val=0)
+ f = (block_given?)? lambda {|i,j| yield(i, j)} : lambda {|i,j| val}
+ Matrix.rows((0...n).collect {|i| (0...m).collect {|j| f.call(i,j)}})
+ end
+ end
+ #
+ # Creates a matrix with the given matrices as diagonal blocks
+ #
+ def diag(*args)
+ dsize = 0
+ sizes = args.collect{|e| x = (e.is_a?(Matrix)) ? e.row_size : 1; dsize += x; x}
+ m = Matrix.zero(dsize)
+ count = 0
+
+ sizes.size.times{|i|
+ range = count..(count+sizes[i]-1)
+ m[range, range] = args[i]
+ count += sizes[i]
+ }
+ m
+ end
+
+ #
+ # Tests if all the elements of two matrix are equal in delta
+ #
+ def equal_in_delta?(m0, m1, delta = 1.0e-10)
+ delta = delta.abs
+ m0.row_size.times {|i|
+ m0.column_size.times {|j|
+ x=m0[i,j]; y=m1[i,j]
+ return false if (x < y - delta or x > y + delta)
+ }
+ }
+ true
+ end
+
+ #
+ # Tests if all the diagonal elements of two matrix are equal in delta
+ #
+ def diag_in_delta?(m1, m0, eps = 1.0e-10)
+ n = m1.row_size
+ return false if n != m0.row_size or m1.column_size != m0.column_size
+ n.times{|i|
+ return false if (m1[i,i]-m0[i,i]).abs > eps
+ }
+ true
+ end
+ end
+
+ #
+ # Returns the matrix divided by a scalar
+ #
+ def quo(v)
+ map {|e| e.quo(v)}
+ end
+
+ alias :old_divition :/
+ #
+ # quo seems always desirable
+ #
+ alias :/ :quo
+
+ #
+ # Set de values of a matrix and the value of wrap property
+ #
+ def set(m)
+ 0.upto(m.row_size - 1) do |i|
+ 0.upto(m.column_size - 1) do |j|
+ self[i, j] = m[i, j]
+ end
+ end
+ self.wrap = m.wrap
+ end
+
+ def wraplate(ijwrap = "")
+ "class << self
+ def [](i, j)
+ #{ijwrap}; @rows[i][j]
+ end
+
+ def []=(i, j, v)
+ #{ijwrap}; @rows[i][j] = v
+ end
+ end"
+ end
+
+ #
+ # Set wrap feature of a matrix
+ #
+ def wrap=(mode = :torus)
+ case mode
+ when :torus then eval(wraplate("i %= row_size; j %= column_size"))
+ when :h_cylinder then eval(wraplate("i %= row_size"))
+ when :v_cylinder then eval(wraplate("j %= column_size"))
+ when :nil then eval(wraplate)
+ end
+ @wrap = mode
+ end
+
+ #
+ # Returns the maximum length of column elements
+ #
+ def max_len_column(j)
+ column_collect(j) {|x| x.to_s.length}.max
+ end
+
+ #
+ # Returns a list with the maximum lengths
+ #
+ def cols_len
+ (0...column_size).collect {|j| max_len_column(j)}
+ end
+
+ alias :to_s_old :to_s
+
+ #
+ # Returns a string for nice printing matrix
+ #
+ def to_s(mode = :pretty, len_col = 3)
+ return "Matrix[]" if empty?
+ if mode == :pretty
+ clen = cols_len
+ to_a.collect {|r|
+ i=0
+ r.map {|x|
+ l=clen[i]
+ i+=1
+ format("%#{l}s ", x.to_s)
+ } << "\n"
+ }.join("")
+ else
+ i = 0; s = ""; cs = column_size
+ each do |e|
+ i = (i + 1) % cs
+ s += format("%#{len_col}s ", e.to_s)
+ s += "\n" if i == 0
+ end
+ s
+ end
+ end
+
+ #
+ # Iterate the elements of a matrix
+ #
+ def each
+ @rows.each {|x| x.each {|e| yield(e)}}
+ nil
+ end
+ #
+ # :section: SPSS like methods
+ #
+
+ # Element wise operation
+ def elementwise_operation(op,other)
+ Matrix.build(row_size,column_size) do |row, column|
+ self[row,column].send(op,other[row,column])
+ end
+ end
+ # Element wise multiplication
+ def e_mult(other)
+ elementwise_operation(:*,other)
+ end
+ # Element wise multiplication
+ def e_quo(other)
+ elementwise_operation(:quo,other)
+ end
+ # Matrix sum of squares
+ def mssq
+ @rows.inject(0){|ac,row| ac+(row.inject(0) {|acr,i| acr+(i**2)})}
+ end
+ def eigenpairs
+ eigval, eigvec= eigenvaluesJacobi, cJacobiV
+ eigenpairs=eigval.size.times.map {|i|
+ [eigval[i], eigvec.column(i)]
+ }
+ eigenpairs=eigenpairs.sort{|a,b| a[0]<=>b[0]}.reverse
+ end
+ # Returns eigenvalues and eigenvectors of a matrix on a Hash
+ # like CALL EIGEN on SPSS.
+ # * _:eigenvectors_: contains the eigenvectors as columns of a new Matrix, ordered in descendent order
+ # * _:eigenvalues_: contains an array with the eigenvalues, ordered in descendent order
+ def eigen
+ ep=eigenpairs
+ { :eigenvalues=>ep.map {|a| a[0]} ,
+ :eigenvectors=>Matrix.columns(ep.map{|a| a[1]})
+ }
+ end
+ # Returns a new matrix with the natural logarithm of initial values
+ def ln
+ map {|v| Math::log(v)}
+ end
+ # Returns a new matrix, with the square root of initial values
+ def sqrt
+ map {|v| Math::sqrt(v)}
+ end
+ # Sum of squares and cross-products. Equal to m.t * m
+ def sscp
+ transpose*self
+ end
+ # Diagonal of matrix.
+ # Returns an array with as many elements as the minimum of the number of rows
+ # and the number of columns in the argument
+ def diagonal
+ min = (row_size<column_size) ? row_size : column_size
+ min.times.map {|i| self[i,i]}
+ end
+ #
+ # :section: Advanced methods
+ #
+
+ #
+ # a hided module of Matrix
+ module MMatrix
+ def self.default_block(block)
+ block ? lambda { |i| block.call(i) } : lambda {|i| i }
+ end
+
+ #
+ # Returns:
+ # 1) the index of row/column and
+ # 2) the values Vector for changing the row/column and
+ # 3) the range of changes
+ #
+ def self.id_vect_range(args, l)
+ i = args[0] # the column(/the row) to be change
+ vect = args[1] # the values vector
+
+ case args.size
+ when 3 then range = args[2] # the range of the elements to be change
+ when 4 then range = args[2]..args[3] #the range by borders
+ else range = 0...l
+ end
+ return i, vect, range
+ end
+
+ end
+
+ #
+ # Returns an array with the elements collected from the row "i".
+ # When a block is given, the elements of that vector are iterated.
+ #
+ def row_collect(i, &block)
+ f = MMatrix.default_block(block)
+ @rows[i].collect {|e| f.call(e)}
+ end
+
+ #
+ # Returns row vector number "i" like Matrix.row as a Vector.
+ # When the block is given, the elements of row "i" are modified
+ #
+ def row!(i)
+ if block_given?
+ @rows[i].collect! {|e| yield e }
+ else
+ Vector.elements(@rows[i], false)
+ end
+ end
+ alias :row_collect! :row!
+
+ #
+ # Returns an array with the elements collected from the column "j".
+ # When a block is given, the elements of that vector are iterated.
+ #
+ def column_collect(j, &block)
+ f = MMatrix.default_block(block)
+ (0...row_size).collect {|r| f.call(self[r, j])}
+ end
+
+#
+# Returns column vector number "j" as a Vector.
+# When the block is given, the elements of column "j" are mmodified
+#
+def column!(j)
+ if block_given?
+ (0...row_size).collect { |i| @rows[i][j] = yield @rows[i][j] }
+ else
+ column(j)
+ end
+end
+alias :column_collect! :column!
+
+#
+# Set a certain column with the values of a Vector
+# m = Matrix.build(3, 3){|i, j| i * 3 + j + 1}
+# m.column= 1, Vector[1, 1, 1], 1..2
+# m => 1 2 3
+# 4 1 6
+# 7 1 9
+#
+def column=(args)
+ m = row_size
+ c, v, r = MMatrix.id_vect_range(args, m)
+ (m..r.begin - 1).each{|i| self[i, c] = 0}
+ [v.size, r.entries.size].min.times{|i| self[i + r.begin, c] = v[i]}
+ ((v.size + r.begin)..r.entries.last).each {|i| self[i, c] = 0}
+end
+
+#
+# Set a certain row with the values of a Vector
+# m = Matrix.new(3, 3){|i, j| i * 3 + j + 1}
+# m.row= 0, Vector[0, 0], 1..2
+# m => 1 0 0
+# 4 5 6
+# 7 8 9
+#
+def row=(args)
+ i, val, range = MMatrix.id_vect_range(args, column_size)
+ row!(i)[range] = val
+end
+
+def norm(p = 2)
+ Vector::Norm.sqnorm(self, p) ** (1.quo(p))
+end
+
+def norm_frobenius
+ norm
+end
+alias :normF :norm_frobenius
+
+#
+# Tests if the matrix is empty or not
+#
+def empty?
+ @rows.empty? if @rows
+end
+
+#
+# Returns row(s) of matrix as a Matrix
+#
+
+def row2matrix(r)
+ if r.is_a? Range
+ a=r.map {|v| v<row_size ? row(v).to_a : nil }.find_all {|v| !v.nil?}
+ else
+ a = row(r).to_a
+ end
+ if r.is_a?(Range) and r.entries.size > 1
+ return Matrix[*a]
+ else
+ return Matrix[a]
+ end
+end
+
+#
+# Returns the column/s of matrix as a Matrix
+#
+def column2matrix(c)
+ if c.is_a?(Range)
+ a=c.map {|v| column(v).to_a}.find_all {|v| v.size>0}
+ else
+ a = column(c).to_a
+ end
+ if c.is_a?(Range)
+ return Matrix.columns(a)
+ else
+ return Matrix[*a.collect{|x| [x]}]
+ end
+end
+
+# Returns the marginal of rows
+def row_sum
+ (0...row_size).collect {|i|
+ row(i).sum
+ }
+end
+# Returns the marginal of columns
+def column_sum
+ (0...column_size).collect {|i|
+ column(i).sum
+ }
+end
+# Calculate sum of cells
+def total_sum
+ row_sum.inject(&:+)
+end
+
+module LU
+ #
+ # Return the Gauss vector, MC, Golub, 3.2.1 Gauss Transformation, p94
+ #
+ def self.gauss_vector(mat, k)
+ t = mat.column2matrix(k)
+ tk = t[k, 0]
+ (0..k).each{|i| t[i, 0] = 0}
+ return t if tk == 0
+ (k+1...mat.row_size).each{|i| t[i, 0] = t[i, 0].to_f / tk}
+ t
+ end
+
+ #
+ # Return the Gauss transformation matrix: M_k = I - tau * e_k^T
+ #
+ def self.gauss(mat, k)
+ i = Matrix.I(mat.column_size)
+ tau = gauss_vector(mat, k)
+ e = i.row2matrix(k)
+ i - tau * e
+ end
+
+ #
+ # LU factorization: A = LU
+ # where L is lower triangular and U is upper triangular
+ #
+ def self.factorization(mat)
+ u = mat.clone
+ n = u.column_size
+ i = Matrix.I(n)
+ l = i.clone
+ (n-1).times {|k|
+ mk = gauss(u, k)
+ u = mk * u # M_{n-1} * ... * M_1 * A = U
+ l += i - mk # L = M_1^{-1} * ... * M_{n-1}^{-1} = I + sum_{k=1}^{n-1} tau * e
+ }
+ return l, u
+ end
+end
+
+#
+# Return the upper triangular matrix of LU factorization
+# M_{n-1} * ... * M_1 * A = U
+#
+def U
+ LU.factorization(self)[1]
+end
+
+#
+# Return the lower triangular matrix of LU factorization
+# L = M_1^{-1} * ... * M_{n-1}^{-1} = I + sum_{k=1}^{n-1} tau * e
+#
+def L
+ LU.factorization(self)[0]
+end
+
+ module Householder
+ #
+ # a QR factorization that uses Householder transformation
+ # Q^T * A = R
+ # MC, Golub & van Loan, pg 224, 5.2.1 Householder QR
+ #
+ def self.QR(mat)
+ h = []
+ a = mat.clone
+ m = a.row_size
+ n = a.column_size
+ n.times do |j|
+ v, beta = a[j..m - 1, j].house
+ h[j] = Matrix.diag(Matrix.robust_I(j), Matrix.I(m-j)- beta * (v * v.t))
+
+ a[j..m-1, j..n-1] = (Matrix.I(m-j) - beta * (v * v.t)) * a[j..m-1, j..n-1]
+ a[(j+1)..m-1,j] = v[2..(m-j)] if j < m - 1
+ end
+ h
+ end
+ #
+ # From the essential part of Householder vector
+ # it returns the coresponding upper(U_j)/lower(V_j) matrix
+ #
+ def self.bidiagUV(essential, dim, beta)
+ v = Vector.concat(Vector[1], essential)
+ dimv = v.size
+ Matrix.diag(Matrix.robust_I(dim - dimv), Matrix.I(dimv) - beta * (v * v.t) )
+ end
+
+ #
+ # Householder Bidiagonalization algorithm. MC, Golub, pg 252, Algorithm 5.4.2
+ # Returns the matrices U_B and V_B such that: U_B^T * A * V_B = B,
+ # where B is upper bidiagonal.
+ #
+ def self.bidiag(mat)
+ a = mat.clone
+ m = a.row_size
+ n = a.column_size
+ ub = Matrix.I(m)
+ vb = Matrix.I(n)
+ n.times do |j|
+ v, beta = a[j..m-1,j].house
+ a[j..m-1, j..n-1] = (Matrix.I(m-j) - beta * (v * v.t)) * a[j..m-1, j..n-1]
+ a[j+1..m-1, j] = v[1..(m-j-1)]
+ ub *= bidiagUV(a[j+1..m-1,j], m, beta) #Ub = U_1 * U_2 * ... * U_n
+ if j < n - 2
+ v, beta = (a[j, j+1..n-1]).house
+ a[j..m-1, j+1..n-1] = a[j..m-1, j+1..n-1] * (Matrix.I(n-j-1) - beta * (v * v.t))
+ a[j, j+2..n-1] = v[1..n-j-2]
+ vb *= bidiagUV(a[j, j+2..n-1], n, beta) #Vb = V_1 * U_2 * ... * V_n-2
+ end
+ end
+ return ub, vb
+ end
+
+ #
+ #Householder Reduction to Hessenberg Form
+ #
+ def self.toHessenberg(mat)
+ h = mat.clone
+ n = h.row_size
+ u0 = Matrix.I(n)
+ for k in (0...n - 2)
+ v, beta = h[k+1..n-1, k].house #the householder matrice part
+ houseV = Matrix.I(n-k-1) - beta * (v * v.t)
+ u0 *= Matrix.diag(Matrix.I(k+1), houseV)
+ h[k+1..n-1, k..n-1] = houseV * h[k+1..n-1, k..n-1]
+ h[0..n-1, k+1..n-1] = h[0..n-1, k+1..n-1] * houseV
+ end
+ return h, u0
+ end
+ end #end of Householder module
+
+ #
+ # Returns the upper bidiagonal matrix obtained with Householder Bidiagonalization algorithm
+ #
+ def bidiagonal
+ ub, vb = Householder.bidiag(self)
+ ub.t * self * vb
+ end
+
+ #
+ # Returns the orthogonal matrix Q of Householder QR factorization
+ # where Q = H_1 * H_2 * H_3 * ... * H_n,
+ #
+ def houseQ
+ h = Householder.QR(self)
+ q = h[0]
+ (1...h.size).each{|i| q *= h[i]}
+ q
+ end
+
+ #
+ # Returns the matrix R of Householder QR factorization
+ # R = H_n * H_n-1 * ... * H_1 * A is an upper triangular matrix
+ #
+ def houseR
+ h = Householder.QR(self)
+ r = self.clone
+ h.size.times{|i| r = h[i] * r}
+ r
+ end
+
+ #
+ # Modified Gram Schmidt QR factorization (MC, Golub, p. 232)
+ # A = Q_1 * R_1
+ #
+ def gram_schmidt
+ a = clone
+ n = column_size
+ m = row_size
+ q = Matrix.build(m, n){0}
+ r = Matrix.zero(n)
+ for k in 0...n
+ r[k,k] = a[0...m, k].norm
+ q[0...m, k] = a[0...m, k] / r[k, k]
+ for j in (k+1)...n
+ r[k, j] = q[0...m, k].t * a[0...m, j]
+ a[0...m, j] -= q[0...m, k] * r[k, j]
+ end
+ end
+ return q, r
+ end
+
+ #
+ # Returns the Q_1 matrix of Modified Gram Schmidt algorithm
+ # Q_1 has orthonormal columns
+ #
+ def gram_schmidtQ
+ gram_schmidt[0]
+ end
+
+ #
+ # Returns the R_1 upper triangular matrix of Modified Gram Schmidt algorithm
+ #
+ def gram_schmidtR
+ gram_schmidt[1]
+ end
+
+
+ module Givens
+ #
+ # Returns the values "c and s" of a Given rotation
+ # MC, Golub, pg 216, Alghorithm 5.1.3
+ #
+ def self.givens(a, b)
+ if b == 0
+ c = 0; s = 0
+ else
+ if b.abs > a.abs
+ tau = Float(-a)/b; s = 1/Math.sqrt(1+tau**2); c = s * tau
+ else
+ tau = Float(-b)/a; c = 1/Math.sqrt(1+tau**2); s = c * tau
+ end
+ end
+ return c, s
+ end
+
+ #
+ # a QR factorization using Givens rotation
+ # Computes the upper triangular matrix R and the orthogonal matrix Q
+ # where Q^t A = R (MC, Golub, p227 algorithm 5.2.2)
+ #
+ def self.QR(mat)
+ r = mat.clone
+ m = r.row_size
+ n = r.column_size
+ q = Matrix.I(m)
+ n.times{|j|
+ m-1.downto(j+1){|i|
+ c, s = givens(r[i - 1, j], r[i, j])
+ qt = Matrix.I(m); qt[i-1..i, i-1..i] = Matrix[[c, s],[-s, c]]
+ q *= qt
+ r[i-1..i, j..n-1] = Matrix[[c, -s],[s, c]] * r[i-1..i, j..n-1]
+ }
+ }
+ return r, q
+ end
+ end
+
+ #
+ # Returns the upper triunghiular matrix R of a Givens QR factorization
+ #
+ def givensR
+ Givens.QR(self)[0]
+ end
+
+ #
+ # Returns the orthogonal matrix Q of Givens QR factorization.
+ # Q = G_1 * ... * G_t where G_j is the j'th Givens rotation
+ # and 't' is the total number of rotations
+ #
+ def givensQ
+ Givens.QR(self)[1]
+ end
+
+ module Hessenberg
+ #
+ # the matrix must be an upper R^(n x n) Hessenberg matrix
+ #
+ def self.QR(mat)
+ r = mat.clone
+ n = r.row_size
+ q = Matrix.I(n)
+ for j in (0...n-1)
+ c, s = Givens.givens(r[j,j], r[j+1, j])
+ cs = Matrix[[c, s], [-s, c]]
+ q *= Matrix.diag(Matrix.robust_I(j), cs, Matrix.robust_I(n - j - 2))
+ r[j..j+1, j..n-1] = cs.t * r[j..j+1, j..n-1]
+ end
+ return q, r
+ end
+ end
+
+ #
+ # Returns the orthogonal matrix Q of Hessenberg QR factorization
+ # Q = G_1 *...* G_(n-1) where G_j is the Givens rotation G_j = G(j, j+1, omega_j)
+ #
+ def hessenbergQ
+ Hessenberg.QR(self)[0]
+ end
+
+ #
+ # Returns the upper triunghiular matrix R of a Hessenberg QR factorization
+ #
+ def hessenbergR
+ Hessenberg.QR(self)[1]
+ end
+
+ #
+ # Return an upper Hessenberg matrix obtained with Householder reduction to Hessenberg Form algorithm
+ #
+ def hessenberg_form_H
+ Householder.toHessenberg(self)[0]
+ end
+
+ #
+ # The real Schur decomposition.
+ # The eigenvalues are aproximated in diagonal elements of the real Schur decomposition matrix
+ #
+ def realSchur(eps = 1.0e-10, steps = 100)
+ h = self.hessenberg_form_H
+ h1 = Matrix[]
+ i = 0
+ loop do
+ h1 = h.hessenbergR * h.hessenbergQ
+ break if Matrix.diag_in_delta?(h1, h, eps) or steps <= 0
+ h = h1.clone
+ steps -= 1
+ i += 1
+ end
+ h1
+ end
+
+
+ module Jacobi
+ #
+ # Returns the nurm of the off-diagonal element
+ #
+ def self.off(a)
+ n = a.row_size
+ sum = 0
+ n.times{|i| n.times{|j| sum += a[i, j]**2 if j != i}}
+ Math.sqrt(sum)
+ end
+
+ #
+ # Returns the index pair (p, q) with 1<= p < q <= n and A[p, q] is the maximum in absolute value
+ #
+ def self.max(a)
+ n = a.row_size
+ max = 0
+ p = 0
+ q = 0
+ n.times{|i|
+ ((i+1)...n).each{|j|
+ val = a[i, j].abs
+ if val > max
+ max = val
+ p = i
+ q = j
+ end
+ }
+ }
+ return p, q
+ end
+
+ #
+ # Compute the cosine-sine pair (c, s) for the element A[p, q]
+ #
+ def self.sym_schur2(a, p, q)
+ if a[p, q] != 0
+ tau = Float(a[q, q] - a[p, p])/(2 * a[p, q])
+ if tau >= 0
+ t = 1./(tau + Math.sqrt(1 + tau ** 2))
+ else
+ t = -1./(-tau + Math.sqrt(1 + tau ** 2))
+ end
+ c = 1./Math.sqrt(1 + t ** 2)
+ s = t * c
+ else
+ c = 1
+ s = 0
+ end
+ return c, s
+ end
+
+ #
+ # Returns the Jacobi rotation matrix
+ #
+ def self.J(p, q, c, s, n)
+ j = Matrix.I(n)
+ j[p,p] = c; j[p, q] = s
+ j[q,p] = -s; j[q, q] = c
+ j
+ end
+ end
+
+ #
+ # Classical Jacobi 8.4.3 Golub & van Loan
+ #
+ def cJacobi(tol = 1.0e-10)
+ a = self.clone
+ n = row_size
+ v = Matrix.I(n)
+ eps = tol * a.normF
+ while Jacobi.off(a) > eps
+ p, q = Jacobi.max(a)
+ c, s = Jacobi.sym_schur2(a, p, q)
+ #print "\np:#{p} q:#{q} c:#{c} s:#{s}\n"
+ j = Jacobi.J(p, q, c, s, n)
+ a = j.t * a * j
+ v = v * j
+ end
+ return a, v
+ end
+
+ #
+ # Returns the aproximation matrix computed with Classical Jacobi algorithm.
+ # The aproximate eigenvalues values are in the diagonal of the matrix A.
+ #
+ def cJacobiA(tol = 1.0e-10)
+ cJacobi(tol)[0]
+ end
+
+ #
+ # Returns a Vector with the eigenvalues aproximated values.
+ # The eigenvalues are computed with the Classic Jacobi Algorithm.
+ #
+ def eigenvaluesJacobi
+ a = cJacobiA
+ Vector[*(0...row_size).collect{|i| a[i, i]}]
+ end
+
+ #
+ # Returns the orthogonal matrix obtained with the Jacobi eigenvalue algorithm.
+ # The columns of V are the eigenvector.
+ #
+ def cJacobiV(tol = 1.0e-10)
+ cJacobi(tol)[1]
+ end
+end
+
+
diff --git a/metadata.gz.sig b/metadata.gz.sig
new file mode 100644
index 0000000..ee44d8e
Binary files /dev/null and b/metadata.gz.sig differ
diff --git a/metadata.yml b/metadata.yml
new file mode 100644
index 0000000..134cb94
--- /dev/null
+++ b/metadata.yml
@@ -0,0 +1,130 @@
+--- !ruby/object:Gem::Specification
+name: extendmatrix
+version: !ruby/object:Gem::Version
+ prerelease: false
+ segments:
+ - 0
+ - 3
+ - 1
+ version: 0.3.1
+platform: ruby
+authors:
+- Cosmin Bonchis
+autorequire:
+bindir: bin
+cert_chain:
+- |
+ -----BEGIN CERTIFICATE-----
+ MIIDMjCCAhqgAwIBAgIBADANBgkqhkiG9w0BAQUFADA/MREwDwYDVQQDDAhjbGJ1
+ c3RvczEVMBMGCgmSJomT8ixkARkWBWdtYWlsMRMwEQYKCZImiZPyLGQBGRYDY29t
+ MB4XDTEwMDMyOTIxMzg1NVoXDTExMDMyOTIxMzg1NVowPzERMA8GA1UEAwwIY2xi
+ dXN0b3MxFTATBgoJkiaJk/IsZAEZFgVnbWFpbDETMBEGCgmSJomT8ixkARkWA2Nv
+ bTCCASIwDQYJKoZIhvcNAQEBBQADggEPADCCAQoCggEBALf8JVMGqE7m5kYb+PNN
+ neZv2pcXV5fQCi6xkyG8bi2/SIFy/LyxuvLzEeOxBeaz1Be93bayIUquOIqw3dyw
+ /KXWa31FxuNuvAm6CN8fyeRYX/ou4cw3OIUUnIvB7RMNIu4wbgeM6htV/QEsNLrv
+ at1/mh9JpqawPrcjIOVMj4BIp67vmzJCaUf+S/H2uYtSO09F+YQE3tv85TPeRmqU
+ yjyXyTc/oJiw1cXskUL8UtMWZmrwNLHXuZWWIMzkjiz3UNdhJr/t5ROk8S2WPznl
+ 0bMy/PMIlAbqWolRn1zl2VFJ3TaXScbqImY8Wf4g62b/1ZSUlGrtnLNsCYXrWiso
+ UPUCAwEAAaM5MDcwCQYDVR0TBAIwADALBgNVHQ8EBAMCBLAwHQYDVR0OBBYEFGu9
+ rrJ1H64qRmNNu3Jj/Qjvh0u5MA0GCSqGSIb3DQEBBQUAA4IBAQCV0Unka5isrhZk
+ GjqSDqY/6hF+G2pbFcbWUpjmC8NWtAxeC+7NGV3ljd0e1SLfoyBj4gnFtFmY8qX4
+ K02tgSZM0eDV8TpgFpWXzK6LzHvoanuahHLZEtk/+Z885lFene+nHadkem1n9iAB
+ cs96JO9/JfFyuXM27wFAwmfHCmJfPF09R4VvGHRAvb8MGzSVgk2i06OJTqkBTwvv
+ JHJdoyw3+8bw9RJ+jLaNoQ+xu+1pQdS2bb3m7xjZpufml/m8zFCtjYM/7qgkKR8z
+ /ZZt8lCiKfFArppRrZayE2FVsps4X6WwBdrKTMZ0CKSXTRctbEj1BAZ67eoTvBBt
+ rpP0jjs0
+ -----END CERTIFICATE-----
+
+date: 2010-08-16 00:00:00 -04:00
+default_executable:
+dependencies:
+- !ruby/object:Gem::Dependency
+ name: rubyforge
+ prerelease: false
+ requirement: &id001 !ruby/object:Gem::Requirement
+ requirements:
+ - - ">="
+ - !ruby/object:Gem::Version
+ segments:
+ - 2
+ - 0
+ - 4
+ version: 2.0.4
+ type: :development
+ version_requirements: *id001
+- !ruby/object:Gem::Dependency
+ name: hoe
+ prerelease: false
+ requirement: &id002 !ruby/object:Gem::Requirement
+ requirements:
+ - - ">="
+ - !ruby/object:Gem::Version
+ segments:
+ - 2
+ - 6
+ - 1
+ version: 2.6.1
+ type: :development
+ version_requirements: *id002
+description: |-
+ The project consists of some enhancements to the Ruby "Matrix" module and includes: LU and QR (Householder, Givens, Gram Schmidt, Hessenberg) decompositions, bidiagonalization, eigenvalue and eigenvector calculations.
+ Include some aditional code to obtains marginal for rows and columns.
+
+ Original code from http://rubyforge.org/projects/matrix/ , done by Cosmin Bonchis as a Google Summer of Code 2007 project for Ruby Central Inc.
+
+ Gem, github repository and current version manteined by Claudio Bustos.
+email:
+- cbonchis_info.uvt.ro
+executables: []
+
+extensions: []
+
+extra_rdoc_files:
+- History.txt
+- LICENSE.txt
+- Manifest.txt
+- ORIGINAL_README.txt
+- README.txt
+files:
+- History.txt
+- LICENSE.txt
+- Manifest.txt
+- ORIGINAL_README.txt
+- README.txt
+- Rakefile
+- lib/extendmatrix.rb
+- spec/extendmatrix_spec.rb
+- spec/spec.opts
+has_rdoc: true
+homepage: http://github.com/clbustos/extendmatrix
+licenses: []
+
+post_install_message:
+rdoc_options:
+- --main
+- README.txt
+require_paths:
+- lib
+required_ruby_version: !ruby/object:Gem::Requirement
+ requirements:
+ - - ">="
+ - !ruby/object:Gem::Version
+ segments:
+ - 0
+ version: "0"
+required_rubygems_version: !ruby/object:Gem::Requirement
+ requirements:
+ - - ">="
+ - !ruby/object:Gem::Version
+ segments:
+ - 0
+ version: "0"
+requirements: []
+
+rubyforge_project: ruby-statsample
+rubygems_version: 1.3.6
+signing_key:
+specification_version: 3
+summary: "The project consists of some enhancements to the Ruby \"Matrix\" module and includes: LU and QR (Householder, Givens, Gram Schmidt, Hessenberg) decompositions, bidiagonalization, eigenvalue and eigenvector calculations"
+test_files: []
+
diff --git a/spec/extendmatrix_spec.rb b/spec/extendmatrix_spec.rb
new file mode 100755
index 0000000..9668be5
--- /dev/null
+++ b/spec/extendmatrix_spec.rb
@@ -0,0 +1,423 @@
+$:.unshift(File.dirname(__FILE__)+"/../lib")
+require 'spec'
+require 'spec/autorun'
+
+require 'extendmatrix'
+
+describe "Vector class extension:" do
+ before do
+ @v = Vector.[](1, 2, 3, 4)
+ end
+ it "[] with range returns correct values" do
+ @v[1..2]==Vector[2,3]
+ end
+ it "[]= should change only specified range" do
+ @v[1..2] = Vector[9, 9, 9, 9, 9]
+ @v.should == Vector[1, 9, 9, 4]
+ end
+
+ it "magnitude methods returns vector length" do
+ @v.magnitude.should==Math::sqrt(30)
+ end
+ it "normalize methods returns the vector normalized" do
+ mag=Math::sqrt(30)
+ @v.normalize.should==@v.quo(mag)
+ end
+ it "Vector.concat method should concat vectors" do
+ Vector.concat(Vector[1], Vector[2, 3]).should == Vector[1, 2, 3]
+ Vector.concat(Vector[1], Vector[2, 3], Vector[4, 5]).should == Vector[1, 2, 3, 4, 5]
+ end
+ it "collect method returns collected vector and doesn't change original" do
+ @v.collect{|i| i*i}.should== Vector.[](1,4,9,16)
+ @v.should==Vector[1,2,3,4]
+ end
+ it "collect! method should change vector" do
+ @v.collect!{|i| i * i}
+ @v.should == Vector.[](1, 4, 9, 16)
+ @v.collect!{|i| 3 * i}
+ @v.should == Vector.[](3, 12, 27, 48)
+ end
+
+ it "each method" do
+ r = []
+ @v.each{|i| r << i+1}
+ r.should == [2, 3, 4, 5]
+ end
+
+ it "max method should return max element" do
+ @v.max.should == 4
+ end
+
+ it "min method should return min element" do
+ @v.min.should == 1
+ end
+
+ it "norm method returns correct norm" do
+ v = Vector.[](3, 4)
+ v.norm.should == 5
+ end
+
+ it "method p_norm(1)" do
+ v = Vector.[](3, 4)
+ v.norm(1).should == 7
+ end
+
+ it "method p_norm(2)" do
+ v = Vector.[](3, 4)
+ v.norm(2).should == 5
+ end
+
+ it "method p_norm(3)" do
+ v = Vector.[](3, 4)
+ v.norm(3).should > 4.49
+ v.norm(3).should < 4.50
+ end
+
+ it "method p_norm(4)" do
+ v = Vector.[](3, 4)
+ v.norm(4).should > 4.28
+ v.norm(4).should < 4.29
+ end
+
+ it "[]= method should change specific value" do
+ @v[3] = 10
+ @v.should == Vector.[](1, 2, 3, 10)
+ end
+
+ it "norm_inf" do
+ @v.norm_inf.should == 4
+ end
+ it "sum method returns the sum of elements" do
+ @v.sum.should==10
+ end
+end
+
+describe "Matrix class extension:" do
+ before do
+ @m = Matrix[[1, 2, 222], [2, 33, 4]]
+ end
+ it "to_s method should return something" do
+ @m.to_s.size.should_not == 0
+ end
+ it "[] method" do
+ m = Matrix.build(4, 3){|i, j| i * 3 + j}
+ m[1, 2].should == 5
+ m[3, 1..2].should == Vector[10, 11]
+ m[0..1, 0..2].should == Matrix[[0, 1, 2], [3, 4, 5]]
+ end
+ it "row_sum method" do
+ @m.row_sum[0].should==225
+ end
+ it "column_sum method" do
+ @m.column_sum[0].should==3
+ end
+ it "total_sum method" do
+ @m.total_sum.should==264
+ end
+
+ it "[]= method" do
+ m = Matrix.build(3, 3){|i, j| i * 3 + j}
+ m[1, 2] = 9
+ m.should == Matrix[[0, 1, 2], [3, 4, 9], [6, 7, 8]]
+ m[1, 1..2] = Vector[8, 8]
+ m.should == Matrix[[0, 1, 2], [3, 8, 8], [6, 7, 8]]
+ m[0..1, 0..1] = Matrix[[0, 0, 0], [0, 0, 0]]
+ m.should == Matrix[[0, 0, 2], [0, 0, 8], [6, 7, 8]]
+ end
+
+ it "set method" do
+ n = Matrix.build(2, 3)
+ n.set(@m)
+ n.should == @m
+ end
+
+ it "set method and wrap value" do
+ @m.wrap = :torus
+ n = Matrix.build(2, 3)
+ n.set(@m)
+ n.wrap.should == :torus
+ end
+
+ it "wrap method" do
+ @m.wrap=:torus
+ @m[2, 3].should == 1
+ @m.wrap=:h_cylinder
+ @m[2, 0].should == 1
+ @m.wrap=:v_cylinder
+ @m[0, 3].should == 1
+ end
+
+ it "maximum length of a column" do
+ @m.max_len_column(1).should == 2
+ end
+
+ it "list of maximum lengths of columns" do
+ @m.cols_len.should == [1, 2, 3]
+ end
+
+ it "matrix each method" do
+ r = []
+ @m.each{|x| r << x + 3}
+ r.should == [4, 5, 225, 5, 36, 7]
+
+ s = 0
+ @m.each{|x| s += x}
+ s.should == 264
+ end
+
+ it "row! method" do
+ @m.row!(0){|x| x+x}.should == [2, 4, 444]
+ @m.should == Matrix[[2, 4, 444], [2, 33, 4]]
+ end
+
+ it "row_collect method" do
+ @m.row_collect(1){|x| x+10}.should == [12, 43, 14]
+ end
+
+ it "column_collect method" do
+ @m.column_collect(0){|x| x*3}.should == [3, 6]
+ end
+
+ it "row_collect! method, identicaly with row!" do
+ @m.row_collect!(0){|x| x+x}.should == [2, 4, 444]
+ @m.should == Matrix[[2, 4, 444], [2, 33, 4]]
+ end
+
+ it "column_collect! method" do
+ @m.column_collect!(2){|x| x+10}.should == [232, 14]
+ @m.should == Matrix[[1, 2, 232], [2, 33, 14]]
+ end
+
+ it "column= " do
+ m = Matrix.build(3, 3){|i, j| i * 3 + j + 1}
+ m.column= 1, Vector[1,1,1,1,1,1]
+ m.should == Matrix[[1, 1, 3],[4, 1, 6],[7, 1, 9]]
+ m.column= 2, Vector[9,9], 0..1
+ m.should == Matrix[[1, 1, 9],[4, 1, 9],[7, 1, 9]]
+ end
+
+ it "row= " do
+ m = Matrix.build(3, 3){|i, j| i * 3 + j + 1}
+ m.row= 1, Vector[1,1,1,1,1,1]
+ m.should == Matrix[[1, 2, 3],[1, 1, 1],[7, 8, 9]]
+ m.row= 2, Vector[9,9], 0..2
+ m.should == Matrix[[1, 2, 3],[1, 1, 1],[9, 9, 0]]
+ end
+
+ it "norm of a matrix" do
+ m = Matrix[[1, 2, 3], [1, 2, 3]]
+ m.norm.should == Math.sqrt(28)
+ end
+
+ it "test empty matrix" do
+ @m.empty?.should == false
+ n = Matrix[]
+ n.empty?.should == true
+ end
+
+ it "row2matrix" do
+ m = Matrix.build(4, 3){|i, j| i * 3 + j + 1}
+
+ m.row2matrix(1..2).should == Matrix[[4, 5, 6],[7, 8, 9]]
+ m.row2matrix(2).should == Matrix[[7, 8, 9]]
+ m.row2matrix(0..4).should == m
+ end
+
+ it "column2matrix" do
+ m = Matrix.build(4, 3){|i, j| i * 3 + j + 1}
+ m.column2matrix(1).should == Matrix[[2], [5], [8], [11]]
+ m.column2matrix(1..1).should == Matrix[[2], [5], [8], [11]]
+
+ m.column2matrix(1..2).should == Matrix[[2,3], [5,6], [8,9], [11,12]]
+
+ end
+
+ it "diag" do
+ m1 = Matrix[[1]]
+ m2 = Matrix[[2, 0], [0, 3]]
+ m3 = Matrix[[4, 0, 0], [0, 5, 0], [0, 0, 6]]
+ a1 = Matrix.build(6, 6){|i, j| i == j ? i + 1: 0}
+ Matrix.diag(m1, m2, m3).should == a1
+ Matrix.diag(m2).should == m2
+ a2 = Matrix[[2, 0, 0, 0],
+ [0, 3, 0, 0],
+ [0, 0, 2, 0],
+ [0, 0, 0, 3]]
+ Matrix.diag(m2, m2).should == a2
+ end
+ it "dup" do
+ m=Matrix.build(4, 3){|i, j| i * 3 + j +1}
+ mm=m.dup
+ mm.object_id.should_not==m.object_id
+ m[0,0]=10
+ mm[0,0].should_not==m[0,0]
+ end
+ it "equal_in_delta" do
+ m = Matrix.build(4, 3){|i, j| i * 3 + j +1}
+ Matrix.equal_in_delta?(m, m).should == true
+ mm = m.clone
+ mm[0,0] += 2
+ Matrix.equal_in_delta?(m, mm, 0.001).should == false
+ Matrix.equal_in_delta?(m, mm, 2).should == true
+ end
+
+ it "diag_in_delta" do
+ Matrix.diag_in_delta?(Matrix.I(5), Matrix.build(4, 4){|i, j| i + j}).should == false
+ m = Matrix.build(5, 5){|i, j| i == j ? 1 + 0.001 * (i+1) : i + j}
+ Matrix.diag_in_delta?(Matrix.I(5), m, 0.01).should == true
+ end
+
+
+ it "e_mult method" do
+ m = Matrix.build(4, 3){|i, j| i * 3 + j +1}
+ n = Matrix.build(4, 3){|i, j| i * 2 + j +2}
+ m.e_mult(n).should==Matrix.build(4,3){|i,j| (i*3+j+1)*(i*2+j+2)}
+ end
+ it "e_mult method" do
+ m = Matrix.build(4, 3){|i, j| i * 3 + j +1}
+ n = Matrix.build(4, 3){|i, j| i * 2 + j +2}
+ m.e_quo(n).should==Matrix.build(4,3){|i,j| (i*3+j+1).quo(i*2+j+2)}
+ end
+ it "mssq method" do
+ m = Matrix[[1,2,3],[4,5,6],[7,8,9]]
+ m.mssq.should==(1..9).each.inject(0) {|ac,v| ac+v**2}
+ end
+ it "eigen" do
+ m=Matrix[[0.95,0.95,0.01,0.01],[0.95,0.95,0.01,0.01],[0.01, 0.01,0.95,0.95], [0.01, 0.01, 0.95, 0.95]]
+ eigenvalues=[1.92,1.88,0.0,0.0]
+ eigen=m.eigen
+ eigen[:eigenvalues].each_with_index do |v,i|
+ v.should be_close(eigenvalues[i],0.01)
+ end
+ eigenvectors=Matrix[[0.5, 0.5, 0.0, 0.707106781186547], [0.5, 0.5, 0.0, -0.707106781186547], [0.5, -0.5, 0.707106781186547, 0.0], [0.5, -0.5, -0.707106781186547, 0.0]]
+ Matrix.equal_in_delta?(eigen[:eigenvectors], eigenvectors).should be_true
+ end
+ it "sqrt" do
+ m=Matrix[[1,4,9],[16,25,36]]
+ m.sqrt.should==Matrix[[1,2,3],[4,5,6]]
+ end
+ it "sscp" do
+ m=Matrix[[1,4,9],[16,25,36]]
+ m.sscp.should== m.t*m
+ end
+ it "diagonal" do
+ m=Matrix.diag(1,4,5,6,7)
+ m.diagonal.should==[1,4,5,6,7]
+ m=Matrix[[1,2,3],[4,5,6],[7,8,9],[10,11,12]]
+ m.diagonal.should==[1,5,9]
+ end
+ it "LU " do
+ m = Matrix[[1, 4, 7],
+ [2, 5, 8],
+ [3, 6, 10]]
+ l = Matrix[[1, 0, 0],[2, 1, 0],[3, 2, 1]]
+ m.L.should == l
+ u = Matrix[[1, 4, 7],[0, -3, -6],[0, 0, 1]]
+ m.U.should == u
+ end
+
+ it "L " do
+ # e.g.: MC, Golub, 3.2 LU factorization, pg 94
+ m = Matrix[[3, 5],
+ [6, 7]]
+ l = Matrix[[1, 0],
+ [2, 1]]
+ m.L.should == l
+ end
+
+ it "U " do
+ # e.g.: MC, Golub, 3.2 LU factorization, pg 94
+ m = Matrix[[3, 5],
+ [6, 7]]
+ u = Matrix[[3, 5],
+ [0, -3]]
+ m.U.should == u
+ end
+
+ it "houseQR " do
+ m = Matrix.build(4, 3){|i, j| i * 3 + j +1}
+ Matrix.equal_in_delta?(m, m.houseQ * m.houseR).should == true
+ q = Matrix[[0.0776, 0.8330, 0.5329, 0.1264],
+ [0.3104, 0.4512, -0.8048, 0.2286],
+ [0.5433, 0.0694, 0.0108, -0.8365],
+ [0.7761, -0.3123, 0.2610, 0.4815]]
+ Matrix.equal_in_delta?(m.houseQ, q, 0.0001).should == true
+ end
+
+ it "houseR " do
+ m = Matrix.build(4, 3){|i, j| i * 3 + j +1}
+ r = Matrix[[12.88409, 14.59162, 16.29916],
+ [ 0, 1.04131, 2.082630],
+ [ 0, 0, 0],
+ [ 0, 0, 0]]
+ Matrix.equal_in_delta?(r, m.houseR, 1.0e-5).should == true
+ end
+
+ it "bidiagonalization" do
+ # MC, Golub, p252, Example 5.4.2
+ m = Matrix.build(4, 3){|i, j| i * 3 + j +1}
+ bidiag = Matrix[[12.884, 21.876, 0 ],
+ [0, 2.246, -0.613],
+ [0, 0, 0 ],
+ [0, 0, 0 ]]
+ Matrix.equal_in_delta?(bidiag, m.bidiagonal, 0.001).should == true
+ end
+
+ it "gram_schmidt" do
+ m = Matrix[[1, 1],
+ [0.001, 0],
+ [0, 0.001]]
+ gsQ = Matrix[[ 1, 0],
+ [0.001, -0.707107],
+ [ 0, 0.707100]]
+ Matrix.equal_in_delta?(gsQ, m.gram_schmidt[0], 0.001).should == true
+ Matrix.equal_in_delta?(m,m.gram_schmidt[0] * m.gram_schmidt[1], 1.0e-5).should == true
+ end
+
+ it "givens " do
+ m = Matrix.build(4, 3){|i, j| i * 3 + j +1}
+ Matrix.equal_in_delta?(m, m.givensQ * m.givensR, 0.001).should == true
+ end
+
+ it "hessenbergQR " do
+ hess = Matrix[[1, 2, 1, 2, 1],
+ [1, 3, 2, 3, 4],
+ [0, 2, 4, 3, 5],
+ [0, 0, 1, 4, 3],
+ [0, 0, 0, 6, 1]]
+ hessR = hess.hessenbergR
+ r = Matrix[[1.41421, 3.53553, 2.12132, 3.53553, 3.53553],
+ [ 0, -2.12132, -4.00693, -3.06412, -5.42115],
+ [ 0, 0, -1.20185, -3.51310, -2.31125],
+ [ 0, 0, 0, -6.30628, -1.54912],
+ [ 0, 0, 0, 0, 1.53929]]
+
+ Matrix.equal_in_delta?(r, hessR, 1.0e-5).should == true
+ Matrix.equal_in_delta?(hessR, hess.hessenbergQ.t * hess, 1.0e-5).should == true
+ end
+
+ it "hessenberg_form " do
+ m = Matrix[[1, 5, 7],[3, 0, 6],[4, 3, 1]]
+ h = Matrix[[1, 8.6, -0.2],[5, 4.96, -0.72],[0, 2.28, -3.96]]
+ Matrix.equal_in_delta?(h, m.hessenberg_form_H, 0.001).should == true
+ end
+
+ it "realSchur" do
+ m = Matrix.build(3, 3){1} + Matrix.diagonal(2, 2, 2)
+ e = Matrix[[5, 0, 0],[0, 2, 0],[0, 0, 2]]
+ Matrix.diag_in_delta?(m.realSchur, e, 1.0e-5).should == true
+ end
+
+ it "Classic Jacobi algorithm" do
+ m = Matrix[[3, 1, 1],[1, 3, 1],[1, 1, 3]]
+ v = Matrix[[2, 0, 0],[0, 5, 0],[0, 0, 2]]
+ Matrix.diag_in_delta?(v, m.cJacobiA, 0.01).should == true
+ a = Matrix[[1, 1, 1, 4],
+ [1, 1, 0, 5],
+ [1, 0, 1, 4],
+ [4, 5, 4, 1]]
+ e = Matrix[[-0.26828, 0, 0, 0], [0, -5.97550, 0, 0], [0, 0, 1.01373, 0], [0, 0, 0, 9.23004]]
+ Matrix.diag_in_delta?(e, a.cJacobiA, 1.0e-5).should == true
+ end
+end
+
diff --git a/spec/spec.opts b/spec/spec.opts
new file mode 100644
index 0000000..4e1e0d2
--- /dev/null
+++ b/spec/spec.opts
@@ -0,0 +1 @@
+--color
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-ruby-extras/ruby-extendmatrix.git
More information about the Pkg-ruby-extras-commits
mailing list