watzon
/
apatite
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases 1 Wiki Activity

Got vector class pretty much complete

This commit is contained in:
Chris Watson 2019-06-11 16:13:05 -07:00
parent 2cb01cd404
commit da8fe56901
No known key found for this signature in database
GPG Key ID: 37DAEF5F446370A4
9 changed files with 1326 additions and 294 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

31
README.md
View File

@ -1,6 +1,6 @@
# Apatite
Apatite is a fundimental package for scientific computing in Crystal. If that sounds like a modified version of the first line from the NumPy homepage, that's because it is. Apatite has (ok, will have) all of the goodness of NumPy sitting atop the blazing speed and beautiful syntax of Crystal.
Apatite is meant to be a collecion of mathmatical and scientific computing algorithms for the Crystal programming language. I don't expect it to ever reach the level of completeness as numpy, but hopefully it can save some people the trouble of implementing these methods on their own.
## Installation
@ -25,25 +25,16 @@ TODO: Write usage instructions here, but first write the library...
## Roadmap
- [ ] Apetite
- [ ] Array Objects
- [ ] NArray (see [numpy.ndarray](https://docs.scipy.org/doc/numpy/reference/arrays.ndarray.html))
- [ ] DType (see [numpy.dtype](https://docs.scipy.org/doc/numpy/reference/arrays.dtypes.html))
- [ ] Scalars
- [ ] Indexing
- [ ] Routines
- [ ] Binary Operations
- [ ] String Operations
- [ ] FFT (see [numpy.fft](https://docs.scipy.org/doc/numpy/reference/routines.fft.html))
- [ ] Financial Functions
- [ ] LinAlg (see [numpy.linalg](https://docs.scipy.org/doc/numpy/reference/routines.linalg.html))
- [ ] Logic Functions
- [ ] Mathematical Functions
- [ ] Matlib (see [numpy.matlib](https://docs.scipy.org/doc/numpy/reference/routines.matlib.html))
- [ ] Padding Arrays
- [ ] Polynomials
- [ ] Random (see [numpy.random](https://docs.scipy.org/doc/numpy/reference/routines.random.html))
- [ ] Sorting, Searching, and Counting
- [ ] Statistics
- [ ] Vector Math
- [ ] Matrix
- [ ] NDArray
- [x] Vector
- [ ] Line
- [ ] Plane
- [ ] Polygon
- [ ] LinkedList
## Contributing

161
spec/apatite/vector_spec.cr Normal file
View File

@ -0,0 +1,161 @@
require "../spec_helper"
describe "Apatite::Vector" do
x = Apatite::Vector.new([3, 4])
describe "#==" do
it "is true if elements are the same" do
(x == Apatite::Vector.new([3, 4])).should be_true
end
it "is false if elements are different" do
(x == Apatite::Vector.new([7, 1])).should be_false
end
end
describe "#e" do
it "returns the `ith` element in the vector" do
x.e(2).should eq(4)
end
it "returns `nil` if `n` is out of range" do
x.e(7).should be_nil
end
end
describe "#to_unit_vector" do
it "successfully converts a vector to a unit vector" do
res = [0.6, 0.8]
x.to_unit_vector.each_with_index do |e, i|
e.should be_close(res[i], 0.1)
end
end
end
describe "#dimensions" do
it "gets vector's dimensions" do
x.dimensions.should eq({1, 2})
end
end
describe "#rows" do
it "gets vector's rows" do
x.rows.should eq(1)
end
end
describe "#cols" do
it "gets vector's columns" do
x.cols.should eq(2)
end
end
describe "#product" do
it "computes the product of a vector" do
x.product.should eq(12)
end
end
describe "#angle_from" do
it "should compute the angle between `y` and `z`" do
y = Apatite::Vector.create([1, 1])
z = Apatite::Vector.create([1, 0])
y.angle_from(z).should be_close(Math::PI / 4, 0.1)
end
end
describe "#parallel_to?" do
it "correctly determines if a vector is parallel to another" do
x.parallel_to?(Apatite::Vector.create([6, 8])).should be_true
x.parallel_to?(Apatite::Vector.create([1, 1])).should be_false
end
end
describe "#antiparallel_to?" do
it "correctly determines if a vector is antiparallel to another" do
x.antiparallel_to?(Apatite::Vector.create([-3, -4])).should be_true
x.antiparallel_to?(x).should be_false
end
end
describe "#perpendicular_to?" do
it "correctly determines if a vector is antiparallel to another" do
x.perpendicular_to?(Apatite::Vector.create([-4, 3])).should be_true
x.perpendicular_to?(x).should be_false
end
end
describe "#dot" do
it "calculates the dot product of the vector and another vector" do
x.dot(Apatite::Vector.create([2, 3])).should eq(18)
end
end
describe "#add" do
it "adds a number to every item in a vector" do
x.add(2).should eq([5, 6])
end
it "adds an enumerable to a vector" do
x.add([3, 2]).should eq([6, 6])
end
end
describe "#subtract" do
it "subtracts a number from every item in a vector" do
x.subtract(2).should eq([1, 2])
end
it "subtracts an enumerable from a vector" do
x.subtract([3, 2]).should eq([0, 2])
end
end
describe "#multiply" do
it "multiplies a number with every item in a vector" do
x.multiply(2).should eq([6, 8])
end
it "multiplies an enumerable with a vector" do
x.multiply([3, 2]).should eq([9, 8])
end
end
describe "#sum" do
it "sums all items in a vector" do
x.sum.should eq(7)
end
end
describe "#chomp" do
it "returns a new vector with the first `n` items of the old vector" do
x.chomp(1).should eq([4])
end
end
describe "#top" do
it "returns a new vector with the last `n` items of the old vector" do
x.top(1).should eq([3])
end
end
describe "#augment" do
it "creates a new vector with the elements fro vector b appended to those from vector a" do
y = x.clone
y.augment(Apatite::Vector.create([5])).should eq([3, 4, 5])
end
end
describe ".log" do
it "should calculate the log of the vector" do
pp x
x.log.should eq([1.0986122886681098, 1.3862943611198906])
end
end
it "should allow for scalar addition" do
a = Apatite::Vector.create([2, 3, 4])
b = a.add(1)
b.should eq([3, 4, 5])
end
end

5
spec/apatite_spec.cr
View File

@ -1,9 +1,4 @@
require "./spec_helper"
describe Apatite do
# TODO: Write tests
it "works" do
false.should eq(true)
end
end

41
src/apatite.cr
View File

@ -1,5 +1,8 @@
require "./apatite/core_ext/*"
require "./apatite/narray"
require "./apatite/core_ext/array"
require "./apatite/matrix/ndarray"
require "./apatite/matrix/vector"
require "./apatite/matrix/matrix"
# Apatite is a fundimental package for scientific computing in Crystal. If that
# sounds like a modified version of the first line from the NumPy homepage,
@ -9,16 +12,8 @@ require "./apatite/narray"
module Apatite
VERSION = "0.1.0"
# def self.zeros(shape : Array(Int32))
# curr = shape.shift
# Array.new(curr) do
# if shape.empty?
# Array.new(curr, 0)
# else
# self.zeros(shape)
# end
# end
# end
class_property precision = 1e-6
class_property approx_precision = 1e-5
## ## ## ## ## ## ## ## ## ## ## ## ##
# # Array Creation
@ -133,6 +128,18 @@ module Apatite
# # Array Manipulation
## ## ## ## ## ## ## ## ## ## ## ## ##
def self.shape(input : Iterable(U) | U) forall U
if input.is_a?(Iterable(U))
if input.size > 0
[input.size].concat(self.shape(input[0]))
else
[0]
end
else
[] of U
end
end
def self.copyto
end
@ -531,10 +538,16 @@ module Apatite
# # Logic Functions
## ## ## ## ## ## ## ## ## ## ## ## ##
def self.all?(arr : Enumerable)
def self.all?(arr : Enumerable, axis : Int32? | Array(Int32)? = nil, keepdims = false)
raise "axis of #{axis} out of bounds for array of size #{arr.size}" unless axis.nil? || axis < arr.size
if axis.nil? && !keepdims
arr.flatten.each { |i| return false unless i == true }
return true
end
end
def self.any?
def self.any?(arr : Enumerable, axis : Int32? | Array(Int32) = nil, keepdims = false)
end
def self.finite?

4
src/apatite/core_ext/array.cr
View File

@ -20,4 +20,8 @@ class Array(T)
max = max_by { |i| i.is_a?(Array) ? i.size : i }
max.is_a?(Array) ? [size, max.size] : [max]
end
def to_vec
Apatite::Vector.create(self)
end
end

269
src/apatite/matrix/matrix.cr Normal file
View File

@ -0,0 +1,269 @@
require "./vector"
module Apatite
class Matrix(T) < Array(T)
include Indexable(T)
@rows : Int32
@cols : Int32
# Creates a Matrix with each element initialized as *value*.
def self.of(value : T, rows : Int32, cols : Int32)
Matrix(T).new(shape) { value }
end
# Creates a Matrix, invoking *initializer* with each pair of indices.
def self.build(shape, &initializer : UInt32, UInt32 -> T)
Matrix(T).new(shape) do |idx|
i = (idx / N).to_u32
j = (idx % N).to_u32
initializer.call(i, j)
end
end
# Creates a single-column `Matrix` from a `Vector`
def self.col_vector(vector)
col = vector.to_a
Matrix(T).from(col, {col.size, 1})
end
# Creates a Matrix from elements contained within a StaticArray.
#
# The matrix will be filled rows first, such that an array of
#
# [1, 2, 3, 4]
#
# becomes
#
# | 1 2 |
# | 3 4 |
#
def self.from(list : Array(T), shape)
raise("Not enough elements to fill matrix") if list.size < shape[0] * shape[1]
Matrix(T).new(shape) do |idx|
list[idx]
end
end
# Creates a single-row `Matrix` from a `Vector`
def self.row_vector(vector)
row = vector.to_a
Matrix(T).from(row, {1, row.size})
end
# Build a zero matrix (all elements populated with zero value of the type
# isntance).
def self.zero(shape)
Matrix(T).new(shape) { T.zero }
end
# Build the idenity matrix for the instanced type and dimensions.
#
# `id` may be used to specify an identity element for the type. If unspecifed
# a numeric identity will be assumed.
# def self.identity(id = T.zero + 1, shape)
# {{ raise("Matrix dimensions must be square") unless M == N }}
# Matrix(T).(shape) do |i, j|
# i == j ? id : T.zero
# end
# end
# Creates Matrix, yielding the linear index for each element to provide an
# initial value.
def initialize(shape, &block : Int32 -> T)
raise("Matrix dimensions must be positive") if shape[0] < 0 || shape[1] < 0
@rows = shape[0]
@cols = shape[1]
@buffer = Pointer(T).malloc(size, &block)
end
# Equality. Returns `true` if each element in `self` is equal to each
# corresponding element in *other*.
def ==(other : Matrix(U)) forall U
{% if other.rows == rows && other.cols == cols %}
each_with_index do |e, i|
return false unless e == other[i]
end
true
{% else %}
false
{% end %}
end
# Equality with another object, or differently sized matrix. Always `false`.
def ==(other)
false
end
# Returns a new Matrix that is the result of performing a matrix addition with
# *other*
def +(other : Matrix)
merge(other) { |a, b| a + b }
end
# Returns a new Matrix that is the result of performing a matrix subtraction
# with *other*
def -(other : Matrix)
merge(other) { |a, b| a - b }
end
# Performs a matrix multiplication with *other*.
def *(other : Matrix)
raise("Dimension mismatch, cannot multiply a #{rows}x#{cols} by a #{other.rows}x#{other.cols}") unless cols == other.cols
Matrix(typeof(self[0] * other[0])).build(other.rows, other.cols) do |i, j|
pairs = row(i).zip other.col(j)
pairs.map(&.product).sum
end
end
# Performs a scalar multiplication with *other*.
def *(other)
map { |x| x * other }
end
# Retrieves the value of the element at *i*,*j*.
#
# Indicies are zero-based. Negative values may be passed for *i* and *j* to
# enable reverse indexing such that `self[-1, -1] == self[M - 1, N - 1]`
# (same behaviour as arrays).
def [](i : Int, j : Int) : T
idx = index i, j
to_unsafe[idx]
end
# Sets the value of the element at *i*,*j*.
def []=(i : Int, j : Int, value : T)
idx = index i, j
to_unsafe[idx] = value
end
# Gets the contents of row *i*.
def row(i : Int)
Vector.new(cols) { |j| self[i, j] }
end
# Gets the contents of column *j*.
def col(j : Int)
Vector.new(rows) { |i| self[i, j] }
end
# Yields the current element at *i*,*j* and updates the value with the
# block's return value.
def update(i, j, &block : T -> T)
idx = index i, j
to_unsafe[idx] = yield to_unsafe[idx]
end
# Apply a morphism to all elements, returning a new Matrix with the result.
def map(&block : T -> U) forall U
Matrix(U).new({rows, cols}) do |idx|
block.call to_unsafe[idx]
end
end
# ditto
def map_with_indices(&block : T, UInt32, UInt32 -> U) forall U
Matrix(U).from({rows, cols}) do |i, j|
block.call self[i, j], i, j
end
end
# ditto
def map_with_index(&block : T, Int32 -> U) forall U
Matrix(U).new({rows, cols}) do |idx|
block.call to_unsafe[idx], idx
end
end
# Apply an endomorphism to `self`, mutating all elements in place.
def map!(&block : T -> T)
each_with_index do |e, idx|
to_unsafe[idx] = yield e
end
self
end
# Merge with another similarly dimensions matrix, apply the passed block to
# each elemenet pair.
def merge(other : Matrix(U), &block : T, U -> _) forall U
raise("Dimension mismatch") unless other.rows == rows && other.cols == cols
map_with_index do |e, i|
block.call e, other[i]
end
end
# Gets an `Vector` of `Vector` representing rows
def row_vectors
Vector.new(rows) { |i| row(i) }
end
# Gets an `Array` of `Vector` representing columns
def col_vectors
Vector.new(cols) { |i| col(i) }
end
# Creates a new matrix that is the result of inverting the rows and columns
# of `self`.
def transpose
Matrix(T).build({cols, rows}) do |i, j|
self[j, i]
end
end
# Returns the dimensions of `self` as a tuple of `{rows, cols}`.
def dimensions
{rows, cols}
end
# Gets the capacity (total number of elements) of `self`.
def size
@rows * @cols
end
# Count of rows.
def rows
@rows
end
# Count of columns.
def cols
@cols
end
# Returns the element at the given linear index, without doing any bounds
# check.
#
# Used by `Indexable`
@[AlwaysInline]
def unsafe_fetch(index : Int)
to_unsafe[index]
end
# Returns the pointer to the underlying element data.
def to_unsafe : Pointer(T)
@buffer
end
# Map *i*,*j* coords to an index within the buffer.
def index(i : Int, j : Int) : T
i += rows if i < 0
j += cols if j < 0
raise IndexError.new if i < 0 || j < 0
raise IndexError.new unless i < rows && j < cols
i * cols + j
end
def inspect
::String.build do |s|
s << "Matrix{" << self.join(", ") << "}"
end
end
end
end

32
src/apatite/matrix/ndarray.cr Normal file
View File

@ -0,0 +1,32 @@
module Apatite
class NDArray
include Enumerable(Float64)
include Indexable(Float64)
include Comparable(NDArray)
getter data : Array(Float64)
getter shape : Array(Int32)
delegate :[], to: @data
delegate :[]?, to: @data
delegate :[]=, to: @data
delegate :unsafe_fetch, to: @data
delegate :to_unsafe, to: @data
delegate :size, to: @data
def initialize(data : Array(Number), shape : Array(Int32)? = nil)
@data = data.is_a?(Array(Float64)) ? data.flatten : data.flatten.map(&.to_f64)
@shape = shape || [@data.size]
end
# Returns the absolute value of every item in the array
def abs
map { |e| e.abs }
end
# Returns the arccosine of each element in the current array.
def acos
end
end
end

822
src/apatite/matrix/vector.cr Normal file
View File

@ -0,0 +1,822 @@
require "../../apatite"
module Apatite
# Represents a mathematical vector, and also constitutes a row or column
# of a `Matrix`
class Vector
include Enumerable(Float64)
include Indexable(Float64)
include Comparable(Vector)
I = Vector.create([1.0, 0.0, 0.0])
J = Vector.create([0.0, 1.0, 0.0])
K = Vector.create([0.0, 0.0, 1.0])
@buffer : Pointer(Float64)
@capacity : Int32
# Returns the number of elements in the vector
getter size : Int32
# Create a new empty `Vector`.
def initialize
@size = 0
@capacity = 0
@buffer = Pointer(Float64).null
end
# Creates a new empty `Vector` backed by a buffer that is initially
# `initial_capacity` big.
#
# The *initial_capacity* is useful to avoid unnecessary reallocations
# of the internal buffer in case of growth. If you have an estimate
# of the maximum number of elements an vector will hold, the vector should
# be initialized with that capacity for improved performance.
#
# ```
# vec = Vector.new(5)
# vec.size # => 0
# ```
def initialize(initial_capacity : Int)
if initial_capacity < 0
raise ArgumentError.new("Negative array size: #{initial_capacity}")
end
@size = 0
@capacity = initial_capacity.to_i
if initial_capacity == 0
@buffer = Pointer(Float64).null
else
@buffer = Pointer(Float64).malloc(initial_capacity)
end
end
# Creates a new `Vector` of the given *size* filled with the same *value* in each position.
#
# ```
# Vector.new(3, 1.0) # => Vector{1.0, 1.0, 1.0}
# ```
def initialize(size : Int, value : Float64)
if size < 0
raise ArgumentError.new("Negative vector size: #{size}")
end
@size = size.to_i
@capacity = size.to_i
if size == 0
@buffer = Pointer(Float64).null
else
@buffer = Pointer(Float64).malloc(size, value)
end
end
# Creates a new `Vector` of the given *size* and invokes the given block once
# for each index of `self`, assigning the block's value in that index.
#
# ```
# Vector.new(3) { |i| (i + 1.0) ** 2.0 } # => Vector{1.0, 4.0, 9.0}
#
# vec = Vector.new(3) { 5.0 }
# vec # => Vector{5.0, 5.0, 5.0}
# ```
def self.new(size : Int, &block : Int32 -> Float64)
Vector.build(size) do |buffer|
size.to_i.times do |i|
buffer[i] = yield i
end
size
end
end
# Creates a new `Vector`, allocating an internal buffer with the given *capacity*,
# and yielding that buffer. The given block must return the desired size of the vector.
#
# This method is **unsafe**.
#
# ```
# Vector.build(3) do |buffer|
# LibSome.fill_buffer_and_return_number_of_elements_filled(buffer)
# end
# ```
def self.build(capacity : Int) : self
vec = Vector.new(capacity)
vec.size = (yield vec.to_unsafe).to_i
vec
end
# Creates a new `Vector` from the elements of another `Indexable`
# collection.
#
# ```
# Vector.create([1, 2, 3]) # Vector{1.0, 2.0, 3.0}
# ```
def self.create(elements : Indexable(Number))
vec = Vector.new
elements.each do |el|
vec.push(el.to_f64)
end
vec
end
# Generates a random vector of size `n` with elements
# in `range`.
def self.random(n, range = Float64::MIN..Float64::MAX)
random = Random.new
vec = Vector.new(n)
n.times do
vec.push random.rand(range)
end
vec
end
# Returns a standard basis n-vector.
def self.basis(size, index)
raise ArgumentError.new("invalid size (#{size} for 1..)") if size < 1
raise ArgumentError.new("invalid index (#{index} for 0...#{size})") unless 0 <= index && index < size
vec = Vector.new(size, 0.0)
vec[index] = 1.0
vec
end
# Create a new vector of size `n` filled with zeros.
def self.zeros(n)
Vector.new(n, 0.0)
end
# Create a new vector of size `n` filled with ones.
def self.ones(n)
Vector.new(n, 1.0)
end
def self.[](*array)
Vector.create(array)
end
# Returns all elements that are within the given range.
#
# Negative indices count backward from the end of the vector (-1 is the last
# element). Additionally, an empty vector is returned when the starting index
# for an element range is at the end of the vector.
#
# Raises `IndexError` if the starting index is out of range.
#
# ```
# v = Vector{1.0, 2.0, 3.0, 4.0, 5.0}
# v[1..3] # => Vector{2.0, 3.0, 4.0}
# v[6..10] # raise IndexError
# ```
def [](range : Range(Int, Int))
self[*Indexable.range_to_index_and_count(range, size)]
end
# Returns count or less (if there aren't enough) elements starting at the
# given start index.
#
# Negative indices count backward from the end of the vector (-1 is the last
# element). Additionally, an empty vector is returned when the starting index
# for an element range is at the end of the vector.
#
# Raises `IndexError` if the starting index is out of range.
#
# ```
# v = Vector{1.0, 2.0, 3.0, 4.0, 5.0}
# v[1, 3] # => Vector{2.0, 3.0, 4.0}
# v[6, 10] # raise IndexError
# ```
def [](start : Int, count : Int)
raise ArgumentError.new "Negative count: #{count}" if count < 0
if start == size
return Vector.new
end
start += size if start < 0
raise IndexError.new unless 0 <= start <= size
if count == 0
return Vector.new
end
count = Math.min(count, size - start)
Vector.build(count) do |buffer|
buffer.copy_from(@buffer + start, count)
count
end
end
# Sets the given value at the given index.
#
# Negative indices can be used to start counting from the end of the array.
# Raises `IndexError` if trying to set an element outside the array's range.
#
# ```
# vec = Vector{1.0, 2.0, 3.0}
# vec[0] = 5.0
# p vec # => Vec{5.0, 2.0, 3.0}
#
# vec[3] = 5.0 # raises IndexError
# ```
@[AlwaysInline]
def []=(index : Int, value)
index = check_index_out_of_bounds index
@buffer[index] = value.to_f64
end
# Replaces a subrange with a single value. All elements in the range
# `index...index+count` are removed and replaced by a single element
# *value*.
#
# If *count* is zero, *value* is inserted at *index*.
#
# Negative values of *index* count from the end of the vector.
#
# ```
# vec = Vector{1.0, 2.0, 3.0, 4.0, 5.0}
# vec[1, 3] = 6.0
# vec # => Vector{1.0, 6.0, 5.0}
#
# vec = Vector{1.0, 2.0, 3.0, 4.0, 5.0}
# vec[1, 0] = 6.0
# vec # => Vector{1.0, 6.0, 2.0, 3.0, 4.0, 5.0}
# ```
def []=(index : Int, count : Int, value)
raise ArgumentError.new "Negative count: #{count}" if count < 0
value = value.to_f64
index = check_index_out_of_bounds index
count = (index + count <= size) ? count : size - index
case count
when 0
insert index, value
when 1
@buffer[index] = value
else
diff = count - 1
(@buffer + index + 1).move_from(@buffer + index + count, size - index - count)
(@buffer + @size - diff).clear(diff)
@buffer[index] = value
@size -= diff
end
value
end
# Replaces a subrange with a single value.
#
# ```
# vec = Vector{1.0, 2.0, 3.0, 4.0, 5.0}
# vec[1..3] = 6.0
# vec # = Vector{1.0, 6.0, 5.0}
#
# vec = Vector{1.0, 2.0, 3.0, 4.0, 5.0}
# vec[1...1] = 6
# vec # = Vector{1.0, 6.0, 2.0, 3.0, 4.0, 5.0}
# ```
def []=(range : Range(Int, Int), value)
self[*Indexable.range_to_index_and_count(range, size)] = value.to_f64
end
# Replaces a subrange with the elements of the given vector.
#
# ```
# vec = Vector{1.0, 2.0, 3.0, 4.0, 5.0}
# vec[1, 3] = Vector{6.0, 7.0, 8.0}
# vec # = Vector{1.0, 6.0, 7.0, 8.0, 5.0}
#
# vec = Vector{1.0, 2.0, 3.0, 4.0, 5.0}
# vec[1, 3] = Vector{6.0, 7.0}
# vec # = Vector{1.0, 6.0, 7.0, 5.0}
# ```
def []=(index : Int, count : Int, values : Vector)
raise ArgumentError.new "Negative count: #{count}" if count < 0
index = check_index_out_of_bounds index
count = (index + count <= size) ? count : size - index
diff = values.size - count
if diff == 0
# Replace values directly
(@buffer + index).copy_from(values.to_unsafe, values.size)
elsif diff < 0
# Need to shrink
diff = -diff
(@buffer + index).copy_from(values.to_unsafe, values.size)
(@buffer + index + values.size).move_from(@buffer + index + count, size - index - count)
(@buffer + @size - diff).clear(diff)
@size -= diff
else
# Need to grow
resize_to_capacity(Math.pw2ceil(@size + diff))
(@buffer + index + values.size).move_from(@buffer + index + count, size - index - count)
(@buffer + index).copy_from(values.to_unsafe, values.size)
@size += diff
end
values
end
# Combined comparison operator. Returns *0* if `self` equals *other*, *1* if
# `self` is greater than *other* and *-1* if `self` is smaller than *other*.
#
# It compares the elements of both vectors in the same position using the
# `<=>` operator. As soon as one of such comparisons returns a non-zero
# value, that result is the return value of the comparison.
#
# If all elements are equal, the comparison is based on the size of the vectors.
#
# ```
# Vector{8.0} <=> Vector(1.0, 2.0, 3.0} # => 1
# Vector{2.0} <=> Vector{4.0, 2.0, 3.0} # => -1
# Vector{1.0, 2.0} <=> Vector{1.0, 2.0} # => 0
# ```
def <=>(other)
min_size = Math.min(size, other.size)
0.upto(min_size - 1) do |i|
n = @buffer[i] <=> other.to_unsafe[i]
return n if n != 0
end
size <=> other.size
end
# Alias for `#push`
def <<(value : Float64)
push(value)
end
def ==(other : Vector)
equals?(other) { |x, y| x == y }
end
# :nodoc:
def ==(other)
false
end
def +(other)
add(other)
end
def -(other)
subtract(other)
end
def *(other)
multiply(other)
end
def clone
Vector.create(@elements.clone)
end
# Invokes the given block for each element of `self`.
def map(&block)
Vector.new(size) { |i| yield @buffer[i] }
end
# Returns true if all vectors are linearly independent.
def independent?(*vs)
vs.each do |v|
raise "Dimension mismatch. Vectors not all the same size." unless v.size == vs.first.size
end
return false if vs.size > sv.first.size
Matrix[*vs].rank.equal?(vs.count)
end
# Invokes the given block for each element of `self`, replacing the element
# with the value returned by the block. Returns `self`.
#
# ```
# vec = Vector{1.0, 2.0, 3.0}
# vec.map! { |x| x * x }
# a # => Vector{1.0, 4.0, 9.0}
# ```
def map!
@buffer.map!(size) { |e| yield e }
self
end
# Optimized version of `Enumerable#map_with_index`.
def map_with_index(&block)
Vector.new(size) { |i| yield @buffer[i], i }
end
# Like `map_with_index`, but mutates `self` instead of allocating a new object.
def map_with_index!(&block)
to_unsafe.map_with_index!(size) { |e, i| yield e, i }
self
end
# Returns the magnitude/euclidian norm of this vector.
#
# [https://en.wikipedia.org/wiki/Euclidean_distance](https://en.wikipedia.org/wiki/Euclidean_distance)
def magnitude
sum = reduce(0.0) { |acc, e| acc += e * e }
Math.sqrt(sum)
end
# Returns the `ith` element of the vector. Returns `nil` if `i`
# is out of bounds. Indexing starts from 1.
def e(i)
(i < 1 || i > @size) ? nil : self[i - 1]
end
# Returns a new vector created by normalizing this one
# to have a magnitude of `1`. If the vector is a zero
# vector, it will not be modified.
def to_unit_vector
r = magnitude
r == 0 ? dup : map { |x| x.to_f64 / r }
end
# Returns the angle between this vector and another in radians.
# If the vectors are mirrored across their axes this will return `nil`.
def angle_from(vector)
v = vector.is_a?(Vector) ? vector : Vector.create(vector)
unless size == v.size
raise "Cannot compute the angle between vectors with different dimensionality"
end
dot = 0_f64
mod1 = 0_f64
mod2 = 0_f64
zip(vector).each do |x, v|
dot += x * v
mod1 += x * x
mod2 += v * v
end
mod1 = Math.sqrt(mod1)
mod2 = Math.sqrt(mod2)
if mod2 * mod2 == 0
return 0.0
end
theta = (dot / (mod1 * mod2)).clamp(-1, 1)
Math.acos(theta)
end
# Returns whether the vectors are parallel to each other.
def parallel_to?(vector)
angle = angle_from(vector)
angle <= Apatite.precision
end
# Returns whether the vectors are antiparallel to each other.
def antiparallel_to?(vector)
angle = angle_from(vector)
(angle - Math::PI).abs <= Apatite.precision
end
# Returns whether the vectors are perpendicular to each other.
def perpendicular_to?(vector)
(dot(vector)).abs <= Apatite.precision
end
# When the input is a number, this returns the result of adding
# it to all vector elements. When it's a vector, the vectors
# will be added together.
def add(value)
run_binary_op(value) { |a, b| a + b }
end
# When the input is a number, this returns the result of subtracting
# it to all vector elements. When it's a vector, the vectors
# will be subtracted.
def subtract(value)
run_binary_op(value) { |a, b| a - b }
end
# When the input is a number, this returns the result of multiplying
# it to all vector elements. When it's a vector, the vectors
# will be element-wise multiplied.
def multiply(value)
run_binary_op(value) { |a, b| a * b }
end
# Returns the sum of all elements in the vector.
def sum
reduce(0) { |acc, i| acc + i }
end
# Returns the cross product of this vector with the others.
#
# ```
# v1 = Vector{1.0, 0.0, 0.0}
# v2 = Vector{0.0, 1.0, 0.0}
# v1.cross(v2) => Vector{0.0, 0.0, 1.0}
# ```
def cross(*vs)
raise "cross product is not defined on vectors of dimension #{size}" unless size >= 2
raise ArgumentError.new("wrong number of arguments (#{vs.size} for #{size - 2})") unless vs.size == size - 2
vs.each do |v|
raise "Dimension mismatch. Vectors not all the same size." unless v.size == size
end
case size
when 2
Vector[-@buffer[1], @buffer[0]]
when 3
v = vs[0]
Vector[ v[2] * @buffer[1] - v[1] * @buffer[2],
v[0] * @buffer[2] - v[2] * @buffer[0],
v[1] * @buffer[0] - v[0] * @buffer[1] ]
else
# TODO
# rows = [self, *vs, Vector.new(size) {|i| Vector.basis(size: size, index: i) }]
# Matrix.rows(rows).laplace_expansion(row: size - 1)
end
end
# Returns a new vector with the first `n` elements removed from
# the beginning.
def chomp(n)
elements = [] of Float64
each_with_index { |e, i| elements << e if i >= n }
Vector.create(elements)
end
# Returns a vector containing only the first `n` elements.
def top(n)
elements = [] of Float64
each_with_index { |e, i| elements << e if i < n }
Vector.create(elements)
end
# Returns a new vector with the provided `elements` concatenated
# on the end.
def augment(elements)
elements = elements.is_a?(Vector) ? elements : Vector.create(elements)
concat(elements)
end
# Return a new `Vector` with the log of every item in `self`.
def log
map { |x| Math.log(x) }
end
# Get the product of all elements in this vector.
def product
reduce { |acc, v| acc *= v }
end
# Get the scalar (dot) product of this vector with `vector`.
#
# [https://en.wikipedia.org/wiki/Scalar_product](https://en.wikipedia.org/wiki/Scalar_product)
def dot(other)
other = other.is_a?(Vector) ? other : Vector.create(other)
unless size == other.size
raise "Cannot compute the dot product of vectors with different dimensionality"
end
product = 0
(0...size).each do |i|
product += self[i] * other[i]
end
product
end
# Returns the (absolute) largest element in this vector.
def max
reduce { |acc, i| i.abs > acc.abs ? i : acc }
end
# Gets the index of the largest element in this vector.
def max_index
idx = 0
each_with_index { |e, i| idx = e.abs > self[idx].abs ? i : idx }
idx
end
# Creates a single-row matrix from this vector.
def covector
Matrix.row_vector(self)
end
# Returns a diagonal `Matrix` with the vectors elements as its
# diagonal elements.
def to_diagonal_matrix
Matrix::Diagonal.new(elements)
end
# Gets the result of rounding the elements of the vector.
def round
map { |x| x.round }
end
# Transpose this vector into a 1xn `Matrix`
def transpose
Matrix.column_vector(self)
end
# Returns a copy of the vector with elements set to `value` if
# they differ from it by less than `Apetite.precision`
def snap_to(value)
map { |y| (y - x).abs <= Apetite.precision ? value : y }
end
# Gets this vector's distance from the argument, when considered
# a point in space.
def distance_from(obj)
if object.is_a?(Plane) || object.is_a?(Line)
return object.distance_from(self)
end
v = elements.is_a?(Vector) ? elements.elements : elements
unless v.size == @elements.size
return nil
end
sum = 0
part = 0
each_with_index do |x, i|
part = x - v[i - 1]
sum += part * part
end
Math.sqrt(sum)
end
# Returns true if the vector is a point on the given line
def lies_on(line)
line.contains(self)
end
# Returns true if the vector is a point on the given plane.
def lies_in(plane)
plane.contains(self)
end
def push(value : Float64)
check_needs_resize
@buffer[@size] = value
@size += 1
self
end
def push(*values : Float64)
new_size = @size + values.size
resize_to_capacity(Math.pw2ceil(new_size)) if new_size > @capacity
values.each_with_index do |value, i|
@buffer[@size + i] = value
end
@size = new_size
self
end
# Rotates the vector about the given `object`. The object should
# be a point if the vector is 2D, and a line if it is 3D. Be
# careful with line directions!
def rotate(t, object)
# TODO
end
# Returns the result of reflecting the point in the given `object`
# (point, line, or plane).
def reflection_in(object)
if object.is_a?(Plane) || object.is_a?(Line)
# object is a plane or line
p = @elements.dup
c = object.point_closest_to(p).elements
return Vector.create([
C[0] + (C[0] - P[0]),
C[1] + (C[1] - P[1]),
C[2] + (C[2] - (P[2] || 0)),
])
end
# object is a point
q = object.is_a?(Vector) ? object.elements : object
unless @elements.size == q.size
return nil
end
map_with_index { |x, i| q[i - 1] + (q[i - 1] - x) }
end
# Utility to make sure vectors are 3D. If they are 2D, a zero
# z-component is added.
def to_3d
v = dup
case v.elements.size
when 3
break
when 2
v.elements.push(0)
else
return nil
end
return v
end
def set_elements(elements)
@elements = elements.is_a?(Vector) ? elements.elements : elements
self
end
def concat(other : Vector)
other_size = other.size
new_size = size + other_size
if new_size > @capacity
resize_to_capacity(Math.pw2ceil(new_size))
end
(@buffer + @size).copy_from(other.to_unsafe, other_size)
@size = new_size
self
end
def to_a
Array(Float64).new(size) { |i| self[i] }
end
def to_s(io)
io << "Vector{"
join ", ", io, &.inspect(io)
io << "}"
end
def pretty_print(pp) : Nil
pp.list("Vector{", self, "}")
end
@[AlwaysInline]
def unsafe_fetch(index : Int)
@buffer[index]
end
# Returns a pointer to the internal buffer where `self`'s elements are stored.
#
# This method is **unsafe** because it returns a pointer, and the pointed might eventually
# not be that of `self` if the array grows and its internal buffer is reallocated.
#
# ```
# vec = Vector{1.0, 2.0, 3.0}
# vec.to_unsafe[0] # => 1.0
# ```
def to_unsafe : Pointer(Float64)
@buffer
end
# Removes all elements from self.
#
# ```
# vec = Vector{1.0, 2.0, 3.0}
# vec.clear # => Vector{}
# ```
def clear
@buffer.clear(@size)
@size = 0
self
end
def zero?
all?(&.zero?)
end
# :nodoc:
protected def size=(size : Int)
@size = size.to_i
end
private def check_needs_resize
double_capacity if @size == @capacity
end
private def double_capacity
resize_to_capacity(@capacity == 0 ? 3 : (@capacity * 2))
end
private def resize_to_capacity(capacity)
@capacity = capacity
if @buffer
@buffer = @buffer.realloc(@capacity)
else
@buffer = Pointer(Float64).malloc(@capacity)
end
end
# Call a block on the value
private def run_binary_op(value, &block : (T, T) -> T)
if value.is_a?(Number)
return map { |v| yield(v, value) }
end
values = value.is_a?(Vector) ? value.elements : value
unless @elements.size == values.size
raise "Cannot perform operations on vectors with different dimensions."
end
map_with_index { |x, i| yield(x, values[i]) }
end
end
end

255
src/apatite/narray.cr
View File

@ -1,255 +0,0 @@
module Apatite
class NArray(T) < Array(T)
include Indexable(T)
include Comparable(NArray)
getter size : Int32
getter shape : Array(Int32)
getter ndim : Int32
def initialize(@selection : T, @dtype = Int32)
@size = @selection.size
@shape = @selection.shape
@ndim = @shape.size
super(1, @selection)
end
# def self.new(shape : Array(Int32), &block : Int32 -> T)
# end
# def ==(other : NArray)
# end
# # :nodoc:
# def ==(other)
# false
# end
# def <=>(other : NArray)
# end
# def &(other : NArray(U)) forall U
# end
# def |(other : NArray(U)) forall U
# end
# def +(other : NArray(U)) forall U
# end
# def -(other : NArray(U)) forall U
# end
# def *(times : Int)
# end
# def <<(value : T)
# end
# @[AlwaysInline]
# def []=(index : Int, value : T)
# end
# def []=(index : Int, count : Int, value : T)
# end
# def []=(range : Range(Int, Int), value : T)
# end
# def []=(index : Int, count : Int, values : NArray(T))
# end
# def []=(range : Range(Int, Int), values : NArray(T))
# end
# def [](range : Range(Int, Int))
# end
# def [](start : Int, count : Int)
# end
def all?
flatten.all?
end
def any?
flatten.any?
end
# def arg_max
# end
# def arg_min
# end
# def arg_partition
# end
# def arg_sort
# end
# def as_type
# end
# def byte_swap(inplace = false)
# end
# def choose
# end
# def clip
# end
# def compress
# end
# def conj
# end
# def conjugate
# end
# def copy
# end
# def cum_prod
# end
# def cum_sum
# end
# def diagonal
# end
# def dot
# end
# def dump
# end
# def dumps
# end
# def fill
# end
# def flatten
# end
# def get_field
# end
# def item
# end
# def item_set
# end
# def max
# end
# def mean
# end
# def min
# end
# def new_byte_order
# end
# def non_zero
# end
# def partition
# end
# def prod
# end
# def ptp
# end
# def put
# end
# def ravel
# end
# def repeat
# end
# def repeat
# end
# def reshape
# end
# def resize
# end
# def round
# end
# def search_sorted
# end
# def set_field
# end
# def set_flags
# end
# def sort
# end
# def squeeze
# end
# def std
# end
# def sum
# end
# def swap_axes
# end
# def take
# end
# def to_bytes
# end
# def to_a
# end
# def to_json
# end
# def to_s
# end
# def to_unsafe
# end
# def to_yaml
# end
# def trace
# end
# def transpose
# end
# @[AlwaysInline]
# def unsafe_fetch(index : Int)
# @buffer[index]
# end
# def var
# end
# def view
# end
end
end