diff --git a/.travis.yml b/.travis.yml index 537bfb905..41dae76f8 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,16 +1,24 @@ language: rust -sudo: false +# use trusty for newer openblas +sudo: required +dist: trusty matrix: include: - rust: stable + env: + - FEATURES='test' - rust: beta - rust: nightly env: - - FEATURES='assign_ops rustc-serialize' + - FEATURES='test' - BENCH=1 branches: only: - master +addons: + apt: + packages: + - libopenblas-dev script: - | cargo build --verbose && diff --git a/Cargo.toml b/Cargo.toml index 6589d4516..fc04f58e3 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,7 +1,7 @@ [package] name = "ndarray" -version = "0.4.0-alpha.6" +version = "0.4.1" authors = ["bluss"] license = "MIT/Apache-2.0" @@ -16,6 +16,8 @@ build = "build.rs" [lib] name = "ndarray" +bench = false +test = false [build-dependencies] rustc_version = "0.1.3" @@ -33,18 +35,24 @@ version = "0.3.16" optional = true [dependencies] +blas-sys = { version = "0.5", optional = true, default-features = false } +rand = "0.3" +# deprecated! use ndarray-rblas instead rblas = { version = "0.0.13", optional = true } -#[dependencies.serde] -#version = "0.4" -#optional = true - [features] - +# deprecated! Will be default assign_ops = [] +blas = ["blas-sys"] +blas-openblas-sys = [ + "blas", + "blas-sys/openblas", + "blas-sys/openblas-provider/system-openblas", +] + # This feature is used for testing -all = ["assign_ops", "rblas", "rustc-serialize"] +test = ["blas-openblas-sys", "rustc-serialize"] [profile.release] [profile.bench] diff --git a/Documentation/split_at.svg b/Documentation/split_at.svg new file mode 100644 index 000000000..e7a4ec693 --- /dev/null +++ b/Documentation/split_at.svg @@ -0,0 +1,755 @@ + + + + + + + + + + + + + + + + + + + + image/svg+xml + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Axis(0) + Axis(1) + Axis(2) + .split_at(Axis(2), 2) + Input shape: (3, 5, 5)Output shapes: (3, 5, 2) and (3, 5, 3) + + diff --git a/Makefile b/Makefile index b7d83c568..39079dea5 100644 --- a/Makefile +++ b/Makefile @@ -3,7 +3,7 @@ DOCCRATES = ndarray # deps to delete the generated docs RMDOCS = -FEATURES = "assign_ops rustc-serialize rblas" +FEATURES = "assign_ops rustc-serialize rblas blas" VERSIONS = $(patsubst %,target/VERS/%,$(DOCCRATES)) @@ -24,6 +24,7 @@ mkdocs: Cargo.toml cargo doc --no-deps --features=$(FEATURES) rm -rf ./master cp -r ./target/doc ./master + cp ./Documentation/split_at.svg ./master/ndarray/ - cat ./custom.css >> master/main.css $(RMDOCS): mkdocs diff --git a/README.rst b/README.rst index 62bbacb18..5b8cdc864 100644 --- a/README.rst +++ b/README.rst @@ -2,7 +2,7 @@ ndarray ========= The ``ndarray`` crate provides an N-dimensional container similar to numpy’s -ndarray. Requires Rust 1.5. +ndarray. Requires Rust 1.7. Please read the `API documentation here (master)`__, `(0.3)`__, `(0.2)`__ @@ -57,8 +57,7 @@ Status and Lookout + ``.fold()`` and ``.zip_mut_with()`` are the most efficient ways to perform single traversal and lock step traversal respectively. + ``.iter()`` and ``.iter_mut()`` are efficient for contiguous arrays. - -- There is experimental bridging to the linear algebra package ``rblas``. + + Can use BLAS in some operations (``dot`` and ``mat_mul``). Crate Feature Flags ------------------- @@ -68,7 +67,7 @@ your `Cargo.toml`. - ``assign_ops`` - - Optional, requires nightly + - Requires Rust 1.8, will be default soon. - Enables the compound assignment operators - ``rustc-serialize`` @@ -78,17 +77,61 @@ your `Cargo.toml`. - ``rblas`` - - Optional, stable + - **Deprecated:** replaced by separate crate ``ndarray-rblas`` - Enables ``rblas`` integration How to use with cargo:: [dependencies] - ndarray = "0.3" + ndarray = "0.4" Recent Changes -------------- +- 0.4.1 + + - Mark iterators ``Send + Sync`` when possible. + +- **0.4.0** `Release Announcement`__ + + - New array splitting via ``.split_at(Axis, Ix)`` and ``.axis_chunks_iter()`` + - Added traits ``NdFloat``, ``AsArray`` and ``From for ArrayView`` which + improve generic programming. + - Array constructors panic when attempting to create an array whose element + count overflows ``usize``. (Would be a debug assertion for overflow before.) + - Performance improvements for ``.map()``. + - Added ``stack`` and macro ``stack![axis, arrays..]`` to concatenate arrays. + - Added constructor ``OwnedArray::range(start, end, step)``. + - The type alias ``Array`` was renamed to ``RcArray`` (and the old name deprecated). + - Binary operators are not defined when consuming a mutable array view as + the left hand side argument anymore. + - Remove methods and items deprecated since 0.3 or earlier; deprecated methods + have notes about replacements in 0.3 docs. + - See below for full changelog through alphas. + +__ http://bluss.github.io/rust/2016/03/06/ndarray-0.4/ + +- 0.4.0-alpha.8 + + - In debug mode, indexing an array out of bounds now has a detailed + message about index and shape. (In release mode it does not.) + - Enable assign_ops feature automatically when it is supported (Rust 1.8 beta + or later). + - Add trait ``NdFloat`` which makes it easy to be generic over ``f32, f64``. + - Add ``From`` implementations that convert slices or references to arrays + into array views. This replaces ``from_slice`` from a previous alpha. + - Add ``AsArray`` trait, which is simply based on those ``From`` implementations. + - Improve ``.map()`` so that it can autovectorize. + - Use ``Axis`` argument in ``RemoveAxis`` too. + - Require ``DataOwned`` in the raw data methods. + - Merged error types into a single ``ShapeError``, which uses no allocated data. + +- 0.4.0-alpha.7 + + - Fix too strict lifetime bound in arithmetic operations like ``&a @ &b``. + - Rename trait Scalar to ScalarOperand (and improve its docs). + - Implement <<= and >>= for arrays. + - 0.4.0-alpha.6 - All axis arguments must now be wrapped in newtype ``Axis``. diff --git a/benches/bench1.rs b/benches/bench1.rs index 7ccab0395..382ec94fb 100644 --- a/benches/bench1.rs +++ b/benches/bench1.rs @@ -1,6 +1,5 @@ #![feature(test)] #![allow(unused_imports)] -#![cfg_attr(feature = "assign_ops", feature(augmented_assignments))] extern crate test; #[macro_use(s)] @@ -19,36 +18,12 @@ use ndarray::{arr0, arr1, arr2}; use test::black_box; #[bench] -fn small_iter_1d(bench: &mut test::Bencher) +fn map(bench: &mut test::Bencher) { - let a = arr1::(&[1., 2., 2., - 3., 4., 4., - 3., 4., 4., - 3., 4., 4., - 5., 6., 6.]); - bench.iter(|| for &elt in a.iter() { black_box(elt); }) -} - -#[bench] -fn small_iter_1d_raw(bench: &mut test::Bencher) -{ - let a = arr1::(&[1., 2., 2., - 3., 4., 4., - 3., 4., 4., - 3., 4., 4., - 5., 6., 6.]); - bench.iter(|| for &elt in a.raw_data().iter() { black_box(elt); }) -} - -#[bench] -fn small_iter_2d(bench: &mut test::Bencher) -{ - let a = arr2::(&[[1., 2., 2.], - [3., 4., 4.], - [3., 4., 4.], - [3., 4., 4.], - [5., 6., 6.]]); - bench.iter(|| for &elt in a.iter() { black_box(elt); }) + let a = OwnedArray::linspace(0., 127., 128).into_shape((8, 16)).unwrap(); + bench.iter(|| { + a.map(|&x| 2. * x) + }); } #[bench] @@ -270,7 +245,6 @@ fn add_2d_regular(bench: &mut test::Bencher) }); } -#[cfg(feature = "assign_ops")] #[bench] fn add_2d_assign_ops(bench: &mut test::Bencher) { @@ -459,19 +433,52 @@ fn bench_col_iter(bench: &mut test::Bencher) bench.iter(|| for elt in it.clone() { black_box(elt); }) } -#[bench] -fn bench_mat_mul_large(bench: &mut test::Bencher) -{ - let a = OwnedArray::::zeros((64, 64)); - let b = OwnedArray::::zeros((64, 64)); - let a = black_box(a.view()); - let b = black_box(b.view()); - bench.iter(|| a.mat_mul(&b)); +macro_rules! mat_mul { + ($modname:ident, $ty:ident, $(($name:ident, $m:expr, $n:expr, $k:expr))+) => { + mod $modname { + use test::{black_box, Bencher}; + use ndarray::OwnedArray; + $( + #[bench] + fn $name(bench: &mut Bencher) + { + let a = OwnedArray::<$ty, _>::zeros(($m, $n)); + let b = OwnedArray::<$ty, _>::zeros(($n, $k)); + let a = black_box(a.view()); + let b = black_box(b.view()); + bench.iter(|| a.mat_mul(&b)); + } + )+ + } + } +} + +mat_mul!{mat_mul_f32, f32, + (m004, 4, 4, 4) + (m007, 7, 7, 7) + (m008, 8, 8, 8) + (m012, 12, 12, 12) + (m016, 16, 16, 16) + (m032, 32, 32, 32) + (m064, 64, 64, 64) + (m127, 127, 127, 127) + (mix16x4, 32, 4, 32) + (mix32x2, 32, 2, 32) +} + +mat_mul!{mat_mul_i32, i32, + (m004, 4, 4, 4) + (m007, 7, 7, 7) + (m008, 8, 8, 8) + (m012, 12, 12, 12) + (m016, 16, 16, 16) + (m032, 32, 32, 32) + (m064, 64, 64, 64) } #[cfg(feature = "rblas")] #[bench] -fn bench_mat_mul_rblas_large(bench: &mut test::Bencher) +fn bench_mat_mul_rblas_64(bench: &mut test::Bencher) { use rblas::Gemm; use rblas::attribute::Transpose; @@ -541,6 +548,22 @@ fn dot_f32_16(bench: &mut test::Bencher) bench.iter(|| a.dot(&b)); } +#[bench] +fn dot_f32_20(bench: &mut test::Bencher) +{ + let a = OwnedArray::::zeros(20); + let b = OwnedArray::::zeros(20); + bench.iter(|| a.dot(&b)); +} + +#[bench] +fn dot_f32_32(bench: &mut test::Bencher) +{ + let a = OwnedArray::::zeros(32); + let b = OwnedArray::::zeros(32); + bench.iter(|| a.dot(&b)); +} + #[bench] fn dot_f32_256(bench: &mut test::Bencher) { @@ -559,6 +582,26 @@ fn dot_f32_1024(bench: &mut test::Bencher) }); } +#[bench] +fn dot_extended(bench: &mut test::Bencher) { + let m = 10; + let n = 33; + let k = 10; + let av = OwnedArray::::zeros((m, n)); + let bv = OwnedArray::::zeros((n, k)); + let mut res = OwnedArray::::zeros((m, k)); + // make a manual simple matrix multiply to test + bench.iter(|| { + for i in 0..m { + for j in 0..k { + unsafe { + *res.uget_mut((i, j)) = av.row(i).dot(&bv.column(j)); + } + } + } + }) +} + #[bench] fn means(bench: &mut test::Bencher) { let a = OwnedArray::from_iter(0..100_000i64); diff --git a/build.rs b/build.rs index 4f7e61b62..8183ca975 100644 --- a/build.rs +++ b/build.rs @@ -8,6 +8,8 @@ extern crate rustc_version; use rustc_version::Channel; const DEPRECATED_CFG: &'static str = "has_deprecated"; +const ASSIGN_FEATURE: &'static str = r#"feature="assign_ops""#; +const ASSIGN_CFG: &'static str = "has_assign"; fn main() { let version = rustc_version::version_meta(); @@ -23,6 +25,16 @@ fn main() { if ndate >= vec![2015, 12, 18] { println!("cargo:rustc-cfg={}", DEPRECATED_CFG); } + // assign_ops is available from nightly 2016-03-01 + if ndate >= vec![2016, 3, 1] { + println!("cargo:rustc-cfg={}", ASSIGN_CFG); + println!("cargo:rustc-cfg={}", ASSIGN_FEATURE); + } + } + } else { + if rustc_version::version_matches(">= 1.8") { + println!("cargo:rustc-cfg={}", ASSIGN_FEATURE); + println!("cargo:rustc-cfg={}", ASSIGN_CFG); } } } diff --git a/custom.css b/custom.css index 29d60e51c..abe09352c 100644 --- a/custom.css +++ b/custom.css @@ -1,6 +1,6 @@ .docblock pre.rust { background: #eeeeff; } -pre.trait, pre.macro, pre.fn, pre.struct, pre.enum, pre.typedef { background: #fcfefc; } +pre.trait, pre.macro, pre.fn, pre.struct, pre.enum, pre.typedef, pre.const { background: inherit; } /* Small “example” label for doc examples */ .docblock pre.rust::before { diff --git a/examples/convo.rs b/examples/convo.rs index 77b42db6a..64abfb89f 100644 --- a/examples/convo.rs +++ b/examples/convo.rs @@ -59,10 +59,10 @@ fn main() { } } } - println!("{:?}", a); + println!("{:2}", a); let mut res = OwnedArray::zeros(a.dim()); for _ in 0..1000 { conv_3x3(&a.view(), &mut res.view_mut(), &SOBEL_X); } - println!("{:?}", res); + println!("{:2}", res); } diff --git a/examples/life.rs b/examples/life.rs index 66913b178..4b38edea6 100644 --- a/examples/life.rs +++ b/examples/life.rs @@ -36,6 +36,33 @@ fn parse(x: &[u8]) -> Board { // 3 neighbors: birth // otherwise: death +#[cfg(feature = "assign_ops")] +fn iterate(z: &mut Board, scratch: &mut Board) { + // compute number of neighbors + let mut neigh = scratch.view_mut(); + neigh.assign_scalar(&0); + neigh += &z.slice(s![0..-2, 0..-2]); + neigh += &z.slice(s![0..-2, 1..-1]); + neigh += &z.slice(s![0..-2, 2.. ]); + + neigh += &z.slice(s![1..-1, 0..-2]); + neigh += &z.slice(s![1..-1, 2.. ]); + + neigh += &z.slice(s![2.. , 0..-2]); + neigh += &z.slice(s![2.. , 1..-1]); + neigh += &z.slice(s![2.. , 2.. ]); + + // birth where n = 3 and z[i] = 0, + // survive where n = 2 || n = 3 and z[i] = 1 + let mut zv = z.slice_mut(s![1..-1, 1..-1]); + + // this is autovectorized amazingly well! + zv.zip_mut_with(&neigh, |y, &n| { + *y = ((n == 3) || (n == 2 && *y > 0)) as u8 + }); +} + +#[cfg(not(feature = "assign_ops"))] fn iterate(z: &mut Board, scratch: &mut Board) { // compute number of neighbors let mut neigh = scratch.view_mut(); diff --git a/ndarray-rblas/Cargo.toml b/ndarray-rblas/Cargo.toml new file mode 100644 index 000000000..4eb5d1373 --- /dev/null +++ b/ndarray-rblas/Cargo.toml @@ -0,0 +1,18 @@ +[package] +name = "ndarray-rblas" +version = "0.1.0" +authors = ["bluss"] +license = "MIT/Apache-2.0" + +repository = "https://github.com/bluss/rust-ndarray" + +description = "`rblas` bindings for `ndarray`." + +keywords = ["multidimensional", "matrix", "blas"] + +[dependencies] +rblas = "0.0.13" +ndarray = { version = "0.4.1" } + +[dev-dependencies] +num = { version = "0.1", default-features = false } diff --git a/ndarray-rblas/LICENSE-APACHE b/ndarray-rblas/LICENSE-APACHE new file mode 100644 index 000000000..16fe87b06 --- /dev/null +++ b/ndarray-rblas/LICENSE-APACHE @@ -0,0 +1,201 @@ + 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. + +END OF TERMS AND CONDITIONS + +APPENDIX: How to apply the Apache License to your work. + + To apply the Apache License to your work, attach the following + boilerplate notice, with the fields enclosed by brackets "[]" + replaced with your own identifying information. (Don't include + the brackets!) The text should be enclosed in the appropriate + comment syntax for the file format. We also recommend that a + file or class name and description of purpose be included on the + same "printed page" as the copyright notice for easier + identification within third-party archives. + +Copyright [yyyy] [name of copyright owner] + +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. diff --git a/ndarray-rblas/LICENSE-MIT b/ndarray-rblas/LICENSE-MIT new file mode 100644 index 000000000..9203baa05 --- /dev/null +++ b/ndarray-rblas/LICENSE-MIT @@ -0,0 +1,25 @@ +Copyright (c) 2015 + +Permission is hereby granted, free of charge, to any +person obtaining a copy of this software and associated +documentation files (the "Software"), to deal in the +Software without restriction, including without +limitation the rights to use, copy, modify, merge, +publish, distribute, sublicense, and/or sell copies of +the Software, and to permit persons to whom the Software +is furnished to do so, subject to the following +conditions: + +The above copyright notice and this permission notice +shall be included in all copies or substantial portions +of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF +ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED +TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A +PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT +SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY +CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION +OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR +IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +DEALINGS IN THE SOFTWARE. diff --git a/ndarray-rblas/src/lib.rs b/ndarray-rblas/src/lib.rs new file mode 100644 index 000000000..00ffe2016 --- /dev/null +++ b/ndarray-rblas/src/lib.rs @@ -0,0 +1,397 @@ +// Copyright 2014-2016 bluss and ndarray developers. +// +// Licensed under the Apache License, Version 2.0 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. +//! Experimental BLAS (Basic Linear Algebra Subprograms) integration +//! +//! Depends on crate [`rblas`], ([docs]). +//! +//! [`rblas`]: https://crates.io/crates/rblas/ +//! [docs]: http://mikkyang.github.io/rust-blas/doc/rblas/ +//! +//! ``` +//! extern crate ndarray; +//! extern crate ndarray_rblas; +//! extern crate rblas; +//! +//! use rblas::Gemv; +//! use rblas::attribute::Transpose; +//! +//! use ndarray::{arr1, arr2}; +//! use ndarray_rblas::AsBlas; +//! +//! fn main() { +//! // Gemv is the operation y = α a x + β y +//! let alpha = 1.; +//! let mut a = arr2(&[[1., 2., 3.], +//! [4., 5., 6.], +//! [7., 8., 9.]]); +//! let x = [1., 0., 1.]; +//! let beta = 1.; +//! let mut y = arr1(&[0., 0., 0.]); +//! +//! Gemv::gemv(Transpose::NoTrans, &alpha, &a.blas(), &x[..], +//! &beta, &mut y.blas()); +//! +//! assert_eq!(y, arr1(&[4., 10., 16.])); +//! } +//! +//! ``` +//! +//! Use the methods in trait `AsBlas` to convert an array into a view that +//! implements rblas’ `Vector` or `Matrix` traits. +//! +//! Blas supports strided vectors and matrices; Matrices need to be contiguous +//! in their lowest dimension, so they will be copied into c-contiguous layout +//! automatically if needed. You should be able to use blocks sliced out +//! from a larger matrix without copying. Use the transpose flags in blas +//! instead of transposing with `ndarray`. +//! +//! Blas has its own error reporting system and will not panic on errors (that +//! I know), instead output its own error conditions, for example on dimension +//! mismatch in a matrix multiplication. +//! + +extern crate rblas; +extern crate ndarray; + +use std::os::raw::{c_int}; + +use rblas::{ + Matrix, + Vector, +}; +use ndarray::{ + ShapeError, + ErrorKind, + ArrayView, + ArrayViewMut, + Data, + DataOwned, + DataMut, + Dimension, + ArrayBase, + Ix, Ixs, +}; + + +/// ***Requires crate feature `"rblas"`*** +pub struct BlasArrayView<'a, A: 'a, D>(ArrayView<'a, A, D>); +impl<'a, A, D: Copy> Copy for BlasArrayView<'a, A, D> { } +impl<'a, A, D: Clone> Clone for BlasArrayView<'a, A, D> { + fn clone(&self) -> Self { + BlasArrayView(self.0.clone()) + } +} + +/// ***Requires crate feature `"rblas"`*** +pub struct BlasArrayViewMut<'a, A: 'a, D>(ArrayViewMut<'a, A, D>); + +struct Priv(T); + +/// Return `true` if the innermost dimension is contiguous (includes +/// the special cases of 0 or 1 length in that axis). +fn is_inner_contiguous(a: &ArrayBase) -> bool + where S: Data, + D: Dimension, +{ + let ndim = a.ndim(); + if ndim == 0 { + return true; + } + a.shape()[ndim - 1] <= 1 || a.strides()[ndim - 1] == 1 +} + +/// If the array is not in the standard layout, copy all elements +/// into the standard layout so that the array is C-contiguous. +fn ensure_standard_layout(a: &mut ArrayBase) + where S: DataOwned, + D: Dimension, + A: Clone +{ + if !a.is_standard_layout() { + let d = a.dim(); + let v: Vec = a.iter().cloned().collect(); + *a = ArrayBase::from_vec_dim(d, v).unwrap(); + } +} + + +impl Priv> + where S: Data, + D: Dimension +{ + fn size_check(&self) -> Result<(), ShapeError> { + let max = c_int::max_value(); + let self_ = &self.0; + for (&dim, &stride) in self_.shape().iter().zip(self_.strides()) { + if dim > max as Ix || stride > max as Ixs { + return Err(ShapeError::from_kind(ErrorKind::RangeLimited)); + } + } + Ok(()) + } + + fn contiguous_check(&self) -> Result<(), ShapeError> { + // FIXME: handle transposed? + if is_inner_contiguous(&self.0) { + Ok(()) + } else { + Err(ShapeError::from_kind(ErrorKind::IncompatibleLayout)) + } + } +} + +impl<'a, A, D> Priv> + where D: Dimension +{ + pub fn into_blas_view(self) -> Result, ShapeError> { + if self.0.ndim() > 1 { + try!(self.contiguous_check()); + } + try!(self.size_check()); + Ok(BlasArrayView(self.0)) + } +} + +impl<'a, A, D> Priv> + where D: Dimension +{ + fn into_blas_view_mut(self) -> Result, ShapeError> { + if self.0.ndim() > 1 { + try!(self.contiguous_check()); + } + try!(self.size_check()); + Ok(BlasArrayViewMut(self.0)) + } +} +/* +*/ + +/// Convert an array into a blas friendly wrapper. +/// +/// Note that `blas` suppors four different element types: `f32`, `f64`, +/// `Complex`, and `Complex`. +/// +/// ***Requires crate feature `"rblas"`*** +pub trait AsBlas { + /// Return an array view implementing Vector (1D) or Matrix (2D) + /// traits. + /// + /// Elements are copied if needed to produce a contiguous matrix.
+ /// The result is always mutable, due to the requirement of having write + /// access to update the layout either way. Breaks sharing if the array is + /// an `RcArray`. + /// + /// **Errors** if any dimension is larger than `c_int::MAX`. + fn blas_checked(&mut self) -> Result, ShapeError> + where S: DataOwned + DataMut, + A: Clone; + + /// Equivalent to `.blas_checked().unwrap()` + /// + /// **Panics** if there was a an error in `.blas_checked()`. + fn blas(&mut self) -> BlasArrayViewMut + where S: DataOwned + DataMut, + A: Clone + { + self.blas_checked().unwrap() + } + + /// Return a read-only array view implementing Vector (1D) or Matrix (2D) + /// traits. + /// + /// The array must already be in a blas compatible layout: its innermost + /// dimension must be contiguous. + /// + /// **Errors** if any dimension is larger than `c_int::MAX`.
+ /// **Errors** if the inner dimension is not c-contiguous. + /// + /// Layout requirements may be loosened in the future. + fn blas_view_checked(&self) -> Result, ShapeError> + where S: Data; + + /// `bv` stands for **b**las **v**iew. + /// + /// Equivalent to `.blas_view_checked().unwrap()` + /// + /// **Panics** if there was a an error in `.blas_view_checked()`. + fn bv(&self) -> BlasArrayView + where S: Data, + { + self.blas_view_checked().unwrap() + } + + /// Return a read-write array view implementing Vector (1D) or Matrix (2D) + /// traits. + /// + /// The array must already be in a blas compatible layout: its innermost + /// dimension must be contiguous. + /// + /// **Errors** if any dimension is larger than `c_int::MAX`.
+ /// **Errors** if the inner dimension is not c-contiguous. + /// + /// Layout requirements may be loosened in the future. + fn blas_view_mut_checked(&mut self) -> Result, ShapeError> + where S: DataMut; + + /// `bvm` stands for **b**las **v**iew **m**ut. + /// + /// Equivalent to `.blas_view_mut_checked().unwrap()` + /// + /// **Panics** if there was a an error in `.blas_view_mut_checked()`. + fn bvm(&mut self) -> BlasArrayViewMut + where S: DataMut, + { + self.blas_view_mut_checked().unwrap() + } + /* + + /// Equivalent to `.blas_checked().unwrap()`, except elements + /// are not copied to make the array contiguous: instead just + /// dimensions and strides are adjusted, and elements end up in + /// arbitrary location. Useful if the content of the array doesn't matter. + /// + /// **Panics** if there was a an error in `blas_checked`. + fn blas_overwrite(&mut self) -> BlasArrayViewMut + where S: DataMut; + */ +} + +/// ***Requires crate feature `"rblas"`*** +impl AsBlas for ArrayBase + where S: Data, + D: Dimension, +{ + fn blas_checked(&mut self) -> Result, ShapeError> + where S: DataOwned + DataMut, + A: Clone, + { + try!(Priv(self.view()).size_check()); + match self.ndim() { + 0 | 1 => { } + 2 => { + if !is_inner_contiguous(self) { + ensure_standard_layout(self); + } + } + _n => ensure_standard_layout(self), + } + Priv(self.view_mut()).into_blas_view_mut() + } + + fn blas_view_checked(&self) -> Result, ShapeError> + where S: Data + { + Priv(self.view()).into_blas_view() + } + + fn blas_view_mut_checked(&mut self) -> Result, ShapeError> + where S: DataMut, + { + Priv(self.view_mut()).into_blas_view_mut() + } + + /* + fn blas_overwrite(&mut self) -> BlasArrayViewMut + where S: DataMut, + { + self.size_check().unwrap(); + if self.dim.ndim() > 1 { + self.force_standard_layout(); + } + BlasArrayViewMut(self.view_mut()) + } + */ +} + +/// **Panics** if `as_mut_ptr` is called on a read-only view. +impl<'a, A> Vector
for BlasArrayView<'a, A, Ix> { + fn len(&self) -> c_int { + self.0.len() as c_int + } + + fn as_ptr(&self) -> *const A { + self.0.as_ptr() + } + + fn as_mut_ptr(&mut self) -> *mut A { + panic!("ndarray: as_mut_ptr called on BlasArrayView (not mutable)"); + } + + // increment: stride + fn inc(&self) -> c_int { + self.0.strides()[0] as c_int + } +} + +impl<'a, A> Vector for BlasArrayViewMut<'a, A, Ix> { + fn len(&self) -> c_int { + self.0.len() as c_int + } + + fn as_ptr(&self) -> *const A { + self.0.as_ptr() + } + + fn as_mut_ptr(&mut self) -> *mut A { + self.0.as_mut_ptr() + } + + // increment: stride + fn inc(&self) -> c_int { + self.0.strides()[0] as c_int + } +} + +/// **Panics** if `as_mut_ptr` is called on a read-only view. +impl<'a, A> Matrix for BlasArrayView<'a, A, (Ix, Ix)> { + fn rows(&self) -> c_int { + self.0.dim().0 as c_int + } + + fn cols(&self) -> c_int { + self.0.dim().1 as c_int + } + + // leading dimension == stride between each row + fn lead_dim(&self) -> c_int { + debug_assert!(self.cols() <= 1 || self.0.strides()[1] == 1); + self.0.strides()[0] as c_int + } + + fn as_ptr(&self) -> *const A { + self.0.as_ptr() + } + + fn as_mut_ptr(&mut self) -> *mut A { + panic!("ndarray: as_mut_ptr called on BlasArrayView (not mutable)"); + } +} + +impl<'a, A> Matrix for BlasArrayViewMut<'a, A, (Ix, Ix)> { + fn rows(&self) -> c_int { + self.0.dim().0 as c_int + } + + fn cols(&self) -> c_int { + self.0.dim().1 as c_int + } + + // leading dimension == stride between each row + fn lead_dim(&self) -> c_int { + debug_assert!(self.cols() <= 1 || self.0.strides()[1] == 1); + self.0.strides()[0] as c_int + } + + fn as_ptr(&self) -> *const A { + self.0.as_ptr() + } + + fn as_mut_ptr(&mut self) -> *mut A { + self.0.as_mut_ptr() + } +} diff --git a/ndarray-rblas/tests/blas.rs b/ndarray-rblas/tests/blas.rs new file mode 100644 index 000000000..5825e0835 --- /dev/null +++ b/ndarray-rblas/tests/blas.rs @@ -0,0 +1,164 @@ +extern crate rblas; +extern crate num; +#[macro_use] extern crate ndarray; +extern crate ndarray_rblas; + +use rblas::Gemm; +use rblas::attribute::Transpose; + +use num::Float; + +use ndarray::{ + OwnedArray, + ArrayView, + ArrayViewMut, + arr2, + Ix, + ShapeError +}; + +use ndarray_rblas::AsBlas; + +#[test] +fn strided_matrix() { + // smoke test, a matrix multiplication of uneven size + let (n, m) = (45, 33); + let mut a = OwnedArray::linspace(0., ((n * m) - 1) as f32, n as usize * m as usize ).into_shape((n, m)).unwrap(); + let mut b = OwnedArray::eye(m); + let mut res = OwnedArray::zeros(a.dim()); + Gemm::gemm(&1., Transpose::NoTrans, &a.blas(), Transpose::NoTrans, &b.blas(), + &0., &mut res.blas()); + assert_eq!(res, a); + + // matrix multiplication, strided + let mut aprim = a.to_shared(); + aprim.islice(s![0..12, 0..11]); + println!("{:?}", aprim.shape()); + println!("{:?}", aprim.strides()); + let mut b = OwnedArray::eye(aprim.shape()[1]); + let mut res = OwnedArray::zeros(aprim.dim()); + Gemm::gemm(&1., Transpose::NoTrans, &aprim.blas(), Transpose::NoTrans, &b.blas(), + &0., &mut res.blas()); + assert_eq!(res, aprim); + + // Transposed matrix multiply + let (np, mp) = aprim.dim(); + let mut res = OwnedArray::zeros((mp, np)); + let mut b = OwnedArray::eye(np); + Gemm::gemm(&1., Transpose::Trans, &aprim.blas(), Transpose::NoTrans, &b.blas(), + &0., &mut res.blas()); + let mut at = aprim.clone(); + at.swap_axes(0, 1); + assert_eq!(at, res); + + // strided, needs copy + let mut abis = a.to_shared(); + abis.islice(s![0..12, ..;2]); + println!("{:?}", abis.shape()); + println!("{:?}", abis.strides()); + let mut b = OwnedArray::eye(abis.shape()[1]); + let mut res = OwnedArray::zeros(abis.dim()); + Gemm::gemm(&1., Transpose::NoTrans, &abis.blas(), Transpose::NoTrans, &b.blas(), + &0., &mut res.blas()); + assert_eq!(res, abis); +} + +#[test] +fn strided_view() { + // smoke test, a matrix multiplication of uneven size + let (n, m) = (45, 33); + let mut a = OwnedArray::linspace(0., ((n * m) - 1) as f32, n as usize * m as usize ).into_shape((n, m)).unwrap(); + let mut b = OwnedArray::eye(m); + let mut res = OwnedArray::zeros(a.dim()); + Gemm::gemm(&1., + Transpose::NoTrans, &a.blas_view_mut_checked().unwrap(), + Transpose::NoTrans, &b.blas_view_mut_checked().unwrap(), + &0., &mut res.blas_view_mut_checked().unwrap()); + assert_eq!(res, a); + + // matrix multiplication, strided + let aprim = a.slice(s![0..12, 0..11]); + let mut b = OwnedArray::eye(aprim.shape()[1]); + let mut res = OwnedArray::zeros(aprim.dim()); + Gemm::gemm(&1., + Transpose::NoTrans, &aprim.bv(), + Transpose::NoTrans, &b.blas(), + &0., &mut res.blas()); + assert_eq!(res, aprim); + + // test out with matrices where lower axis is strided but has length 1 + let mut a3 = arr2(&[[1., 2., 3.]]); + a3.swap_axes(0, 1); + let mut b = OwnedArray::eye(a3.shape()[1]); + let mut res = arr2(&[[0., 0., 0.]]); + res.swap_axes(0, 1); + Gemm::gemm(&1., + Transpose::NoTrans, &a3.bvm(), + Transpose::NoTrans, &b.blas(), + &0., &mut res.bvm()); + assert_eq!(res, a3); +} + +#[test] +fn as_blas() { + let mut a = OwnedArray::::zeros((4, 4)); + assert!(a.blas_view_mut_checked().is_ok()); + a.swap_axes(0, 1); + assert!(a.blas_view_mut_checked().is_err()); + a.swap_axes(0, 1); + + { + // increased row stride + let mut b = a.slice_mut(s![..;2, ..]); + assert!(b.blas_view_mut_checked().is_ok()); + b.bvm(); // no panic + } + { + // inner dimension is not contig + let mut b = a.slice_mut(s![.., ..;2]); + assert!(b.blas_view_mut_checked().is_err()); + } + { + // inner dimension is length 1, is ok again + let mut b = a.slice_mut(s![.., ..;4]); + assert!(b.blas_view_mut_checked().is_ok()); + b.bvm(); + } +} + +type Ix2 = (Ix, Ix); +fn dot(a: ArrayView, b: ArrayView, + c: &mut ArrayViewMut) -> Result<(), ShapeError> + where F: Gemm + Float, +{ + let at = Transpose::NoTrans; + let bt = Transpose::NoTrans; + + let ba = try!(a.blas_view_checked()); + let bb = try!(b.blas_view_checked()); + let mut bc = try!(c.blas_view_mut_checked()); + F::gemm(&F::one(), + at, &ba, + bt, &bb, + &F::zero(), &mut bc); + Ok(()) +} + +#[test] +fn test_dot() { + let mut a = arr2(&[[1., 2.], + [0., 3.]]); + + let b = a.clone(); + let mut res = a.clone(); + res.assign_scalar(&0.); + + dot(a.view(), b.view(), &mut res.view_mut()).unwrap(); + println!("{:?}", res); + + a.swap_axes(0, 1); + res.assign_scalar(&0.); + + let result = dot(a.view(), b.view(), &mut res.view_mut()); + assert!(result.is_err()); +} diff --git a/src/arrayformat.rs b/src/arrayformat.rs index d0597ef1c..53232d94c 100644 --- a/src/arrayformat.rs +++ b/src/arrayformat.rs @@ -1,3 +1,10 @@ +// Copyright 2014-2016 bluss and ndarray developers. +// +// Licensed under the Apache License, Version 2.0 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. use std::fmt; use super::{ ArrayBase, @@ -8,7 +15,7 @@ use super::{ fn format_array(view: &ArrayBase, f: &mut fmt::Formatter, mut format: F) -> fmt::Result - where F: FnMut(&mut fmt::Formatter, &A) -> fmt::Result, + where F: FnMut(&A, &mut fmt::Formatter) -> fmt::Result, D: Dimension, S: Data, { @@ -63,7 +70,7 @@ fn format_array(view: &ArrayBase, f: &mut fmt::Formatter, try!(write!(f, ", ")); } first = false; - try!(format(f, elt)); + try!(format(elt, f)); if update_index { last_index = index; @@ -85,7 +92,7 @@ impl<'a, A: fmt::Display, S, D: Dimension> fmt::Display for ArrayBase where S: Data, { fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { - format_array(self, f, |f, elt| elt.fmt(f)) + format_array(self, f, <_>::fmt) } } @@ -98,7 +105,10 @@ impl<'a, A: fmt::Debug, S, D: Dimension> fmt::Debug for ArrayBase where S: Data, { fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { - format_array(self, f, |f, elt| elt.fmt(f)) + // Add extra information for Debug + try!(format_array(self, f, <_>::fmt)); + try!(write!(f, " shape={:?}, strides={:?}", self.shape(), self.strides())); + Ok(()) } } @@ -111,7 +121,7 @@ impl<'a, A: fmt::LowerExp, S, D: Dimension> fmt::LowerExp for ArrayBase where S: Data, { fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { - format_array(self, f, |f, elt| elt.fmt(f)) + format_array(self, f, <_>::fmt) } } @@ -124,7 +134,7 @@ impl<'a, A: fmt::UpperExp, S, D: Dimension> fmt::UpperExp for ArrayBase where S: Data, { fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { - format_array(self, f, |f, elt| elt.fmt(f)) + format_array(self, f, <_>::fmt) } } /// Format the array using `LowerHex` and apply the formatting parameters used @@ -136,7 +146,7 @@ impl<'a, A: fmt::LowerHex, S, D: Dimension> fmt::LowerHex for ArrayBase where S: Data, { fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { - format_array(self, f, |f, elt| elt.fmt(f)) + format_array(self, f, <_>::fmt) } } diff --git a/src/arrayserialize.rs b/src/arrayserialize.rs index b443999d3..d4723a131 100644 --- a/src/arrayserialize.rs +++ b/src/arrayserialize.rs @@ -1,3 +1,10 @@ +// Copyright 2014-2016 bluss and ndarray developers. +// +// Licensed under the Apache License, Version 2.0 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. use serde::ser::impls::SeqIteratorVisitor; use serde::{self, Serialize, Deserialize}; diff --git a/src/arraytraits.rs b/src/arraytraits.rs index 8494763be..1c76c2a0f 100644 --- a/src/arraytraits.rs +++ b/src/arraytraits.rs @@ -1,3 +1,10 @@ +// Copyright 2014-2016 bluss and ndarray developers. +// +// Licensed under the Apache License, Version 2.0 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. #[cfg(feature = "rustc-serialize")] use serialize::{Encodable, Encoder, Decodable, Decoder}; @@ -9,16 +16,10 @@ use std::ops::{ IndexMut, }; -use super::{ - Dimension, Ix, +use imp_prelude::*; +use { Elements, ElementsMut, - ArrayBase, - ArrayView, - ArrayViewMut, - Data, - DataMut, - DataOwned, NdIndex, }; @@ -30,6 +31,31 @@ fn array_out_of_bounds() -> ! { panic!("ndarray: index out of bounds"); } +// Macro to insert more informative out of bounds message in debug builds +#[cfg(debug_assertions)] +macro_rules! debug_bounds_check { + ($self_:ident, $index:expr) => { + if let None = $index.index_checked(&$self_.dim, &$self_.strides) { + panic!("ndarray: index {:?} is out of bounds for array of shape {:?}", + $index, $self_.shape()); + } + }; +} + +#[cfg(not(debug_assertions))] +macro_rules! debug_bounds_check { + ($self_:ident, $index:expr) => { }; +} + +#[inline(always)] +pub fn debug_bounds_check(_a: &ArrayBase, _index: &I) + where D: Dimension, + I: NdIndex, + S: Data, +{ + debug_bounds_check!(_a, *_index); +} + /// Access the element at **index**. /// /// **Panics** if index is out of bounds. @@ -41,6 +67,7 @@ impl Index for ArrayBase type Output = S::Elem; #[inline] fn index(&self, index: I) -> &S::Elem { + debug_bounds_check!(self, index); self.get(index).unwrap_or_else(|| array_out_of_bounds()) } } @@ -55,6 +82,7 @@ impl IndexMut for ArrayBase { #[inline] fn index_mut(&mut self, index: I) -> &mut S::Elem { + debug_bounds_check!(self, index); self.get_mut(index).unwrap_or_else(|| array_out_of_bounds()) } } @@ -248,3 +276,84 @@ impl Decodable for ArrayBase }) } } + + +// use "raw" form instead of type aliases here so that they show up in docs +/// Implementation of `ArrayView::from(&S)` where `S` is a slice or slicable. +/// +/// Create a one-dimensional read-only array view of the data in `slice`. +impl<'a, A, Slice: ?Sized> From<&'a Slice> for ArrayBase, Ix> + where Slice: AsRef<[A]> +{ + fn from(slice: &'a Slice) -> Self { + let xs = slice.as_ref(); + unsafe { + Self::new_(xs.as_ptr(), xs.len(), 1) + } + } +} + +/// Implementation of `ArrayView::from(&A)` where `A` is an array. +/// +/// Create a read-only array view of the array. +impl<'a, A, S, D> From<&'a ArrayBase> for ArrayBase, D> + where S: Data, + D: Dimension, +{ + fn from(array: &'a ArrayBase) -> Self { + array.view() + } +} + +/// Implementation of `ArrayViewMut::from(&mut S)` where `S` is a slice or slicable. +/// +/// Create a one-dimensional read-write array view of the data in `slice`. +impl<'a, A, Slice: ?Sized> From<&'a mut Slice> for ArrayBase, Ix> + where Slice: AsMut<[A]> +{ + fn from(slice: &'a mut Slice) -> Self { + let xs = slice.as_mut(); + unsafe { + Self::new_(xs.as_mut_ptr(), xs.len(), 1) + } + } +} + +/// Implementation of `ArrayViewMut::from(&mut A)` where `A` is an array. +/// +/// Create a read-write array view of the array. +impl<'a, A, S, D> From<&'a mut ArrayBase> for ArrayBase, D> + where S: DataMut, + D: Dimension, +{ + fn from(array: &'a mut ArrayBase) -> Self { + array.view_mut() + } +} + +/// Argument conversion into an array view +/// +/// The trait is parameterized over `A`, the element type, and `D`, the +/// dimensionality of the array. `D` defaults to one-dimensional. +/// +/// Use `.into()` to do the conversion. +/// +/// ``` +/// use ndarray::AsArray; +/// +/// fn sum<'a, V: AsArray<'a, f64>>(data: V) -> f64 { +/// let array_view = data.into(); +/// array_view.scalar_sum() +/// } +/// +/// assert_eq!( +/// sum(&[1., 2., 3.]), +/// 6. +/// ); +/// +/// ``` +pub trait AsArray<'a, A: 'a, D = Ix> : Into> where D: Dimension { } +impl<'a, A: 'a, D, T> AsArray<'a, A, D> for T + where T: Into>, + D: Dimension, +{ } diff --git a/src/blas.rs b/src/blas.rs index d7dae530a..4ffabd15f 100644 --- a/src/blas.rs +++ b/src/blas.rs @@ -1,3 +1,10 @@ +// Copyright 2014-2016 bluss and ndarray developers. +// +// Licensed under the Apache License, Version 2.0 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. //! Experimental BLAS (Basic Linear Algebra Subprograms) integration //! //! ***Requires crate feature `"rblas"`*** @@ -48,6 +55,8 @@ //! I know), instead output its own error conditions, for example on dimension //! mismatch in a matrix multiplication. //! +#![cfg_attr(has_deprecated, deprecated(note="`rblas` integration has moved to crate `ndarray-rblas`, use it instead."))] +#![allow(deprecated)] use std::os::raw::{c_int}; @@ -59,6 +68,7 @@ use super::{ ShapeError, zipsl, }; +use error::{from_kind, ErrorKind}; use imp_prelude::*; @@ -82,9 +92,7 @@ impl ArrayBase let max = c_int::max_value(); for (&dim, &stride) in zipsl(self.shape(), self.strides()) { if dim > max as Ix || stride > max as Ixs { - return Err(ShapeError::DimensionTooLarge(self.shape() - .to_vec() - .into_boxed_slice())); + return Err(from_kind(ErrorKind::RangeLimited)); } } Ok(()) @@ -95,7 +103,7 @@ impl ArrayBase if self.is_inner_contiguous() { Ok(()) } else { - Err(ShapeError::IncompatibleLayout) + Err(from_kind(ErrorKind::IncompatibleLayout)) } } } @@ -273,7 +281,7 @@ impl<'a, A> Vector for BlasArrayView<'a, A, Ix> { } fn as_mut_ptr(&mut self) -> *mut A { - panic!("BlasArrayView is not mutable"); + panic!("ndarray: as_mut_ptr called on BlasArrayView (not mutable)"); } // increment: stride @@ -322,7 +330,7 @@ impl<'a, A> Matrix for BlasArrayView<'a, A, (Ix, Ix)> { } fn as_mut_ptr(&mut self) -> *mut A { - panic!("BlasArrayView is not mutable"); + panic!("ndarray: as_mut_ptr called on BlasArrayView (not mutable)"); } } diff --git a/src/data_traits.rs b/src/data_traits.rs index 1fff60a28..530efb12c 100644 --- a/src/data_traits.rs +++ b/src/data_traits.rs @@ -1,3 +1,10 @@ +// Copyright 2014-2016 bluss and ndarray developers. +// +// Licensed under the Apache License, Version 2.0 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. //! The data (inner representation) traits for ndarray @@ -10,33 +17,47 @@ use { ViewRepr, }; -/// Array’s inner representation. +/// Array representation trait. /// /// ***Note:*** `Data` is not an extension interface at this point. /// Traits in Rust can serve many different roles. This trait is public because /// it is used as a bound on public methods. pub unsafe trait Data { + /// The array element type. type Elem; + #[doc(hidden)] fn slice(&self) -> &[Self::Elem]; } -/// Array’s writable inner representation. +/// Array representation trait. +/// +/// For an array with writable elements. +/// +/// ***Internal trait, see `Data`.*** pub unsafe trait DataMut : Data { + #[doc(hidden)] fn slice_mut(&mut self) -> &mut [Self::Elem]; + #[doc(hidden)] #[inline] fn ensure_unique(&mut ArrayBase) where Self: Sized, D: Dimension { } + #[doc(hidden)] #[inline] fn is_unique(&mut self) -> bool { true } } -/// Clone an Array’s storage. +/// Array representation trait. +/// +/// An array representation that can be cloned. +/// +/// ***Internal trait, see `Data`.*** pub unsafe trait DataClone : Data { + #[doc(hidden)] /// Unsafe because, `ptr` must point inside the current storage. unsafe fn clone_with_ptr(&self, ptr: *mut Self::Elem) -> (Self, *mut Self::Elem); } @@ -145,13 +166,23 @@ unsafe impl<'a, A> DataMut for ViewRepr<&'a mut A> { } } -/// Array representation that is a unique or shared owner of its data. +/// Array representation trait. +/// +/// A representation that is a unique or shared owner of its data. +/// +/// ***Internal trait, see `Data`.*** pub unsafe trait DataOwned : Data { + #[doc(hidden)] fn new(elements: Vec) -> Self; + #[doc(hidden)] fn into_shared(self) -> Rc>; } -/// Array representation that is a lightweight view. +/// Array representation trait. +/// +/// A representation that is a lightweight view. +/// +/// ***Internal trait, see `Data`.*** pub unsafe trait DataShared : Clone + DataClone { } unsafe impl DataShared for Rc> {} diff --git a/src/dimension.rs b/src/dimension.rs index 704fafa21..f7c768bc7 100644 --- a/src/dimension.rs +++ b/src/dimension.rs @@ -1,8 +1,17 @@ +// Copyright 2014-2016 bluss and ndarray developers. +// +// Licensed under the Apache License, Version 2.0 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. +use std::cmp::Ordering; +use std::fmt::Debug; use std::slice; use super::{Si, Ix, Ixs}; use super::zipsl; -use stride_error::StrideError; +use error::{from_kind, ErrorKind, ShapeError}; /// Calculate offset from `Ix` stride converting sign properly #[inline] @@ -63,11 +72,11 @@ pub fn dim_stride_overlap(dim: &D, strides: &D) -> bool { /// of the slice. Also, the strides should not allow a same element to be /// referenced by two different index. pub fn can_index_slice(data: &[A], dim: &D, strides: &D) - -> Result<(), StrideError> + -> Result<(), ShapeError> { if strides.slice().iter().cloned().all(stride_is_positive) { if dim.size_checked().is_none() { - return Err(StrideError::OutOfBounds); + return Err(from_kind(ErrorKind::OutOfBounds)); } let mut last_index = dim.clone(); for mut index in last_index.slice_mut().iter_mut() { @@ -80,17 +89,17 @@ pub fn can_index_slice(data: &[A], dim: &D, strides: &D) // offset is guaranteed to be positive so no issue converting // to usize here if (offset as usize) >= data.len() { - return Err(StrideError::OutOfBounds); + return Err(from_kind(ErrorKind::OutOfBounds)); } if dim_stride_overlap(dim, strides) { - return Err(StrideError::Unsupported); + return Err(from_kind(ErrorKind::Unsupported)); } } else { - return Err(StrideError::OutOfBounds); + return Err(from_kind(ErrorKind::OutOfBounds)); } Ok(()) } else { - Err(StrideError::Unsupported) + return Err(from_kind(ErrorKind::Unsupported)); } } @@ -119,13 +128,13 @@ fn stride_offset_checked_arithmetic(dim: &D, strides: &D, index: &D) Some(offset) } -/// Trait for the shape and index types of arrays. +/// Array shape and index trait. /// /// `unsafe` because of the assumptions in the default methods. /// /// ***Don't implement or call methods in this trait, its interface is internal /// to the crate and will evolve at will.*** -pub unsafe trait Dimension : Clone + Eq { +pub unsafe trait Dimension : Clone + Eq + Debug + Send + Sync { /// `SliceArg` is the type which is used to specify slicing for this /// dimension. /// @@ -325,6 +334,51 @@ pub unsafe trait Dimension : Clone + Eq { } offset } + +} + +/// Implementation-specific extensions to `Dimension` +pub trait DimensionExt { +// note: many extensions go in the main trait if they need to be special- +// cased per dimension + /// Get the dimension at `axis`. + /// + /// *Panics* if `axis` is out of bounds. + #[inline] + fn axis(&self, axis: Axis) -> Ix; + + /// Set the dimension at `axis`. + /// + /// *Panics* if `axis` is out of bounds. + #[inline] + fn set_axis(&mut self, axis: Axis, value: Ix); +} + +impl DimensionExt for D + where D: Dimension +{ + #[inline] + fn axis(&self, axis: Axis) -> Ix { + self.slice()[axis.axis()] + } + + #[inline] + fn set_axis(&mut self, axis: Axis, value: Ix) { + self.slice_mut()[axis.axis()] = value; + } +} + +impl<'a> DimensionExt for [Ix] +{ + #[inline] + fn axis(&self, axis: Axis) -> Ix { + self[axis.axis()] + } + + #[inline] + fn set_axis(&mut self, axis: Axis, value: Ix) { + self[axis.axis()] = value; + } } fn abs_index(len: Ixs, index: Ixs) -> Ix { @@ -541,11 +595,13 @@ unsafe impl Dimension for Vec fn slice_mut(&mut self) -> &mut [Ix] { self } } -/// Helper trait to define a larger-than relation for array shapes: +/// Array shape with a next smaller dimension. +/// +/// `RemoveAxis` defines a larger-than relation for array shapes: /// removing one axis from *Self* gives smaller dimension *Smaller*. pub trait RemoveAxis : Dimension { type Smaller: Dimension; - fn remove_axis(&self, axis: usize) -> Self::Smaller; + fn remove_axis(&self, axis: Axis) -> Self::Smaller; } macro_rules! impl_shrink( @@ -555,12 +611,12 @@ impl RemoveAxis for ($from $(,$more)*) type Smaller = ($($more),*); #[allow(unused_parens)] #[inline] - fn remove_axis(&self, axis: usize) -> ($($more),*) { + fn remove_axis(&self, axis: Axis) -> ($($more),*) { let mut tup = ($(0 as $more),*); { let mut it = tup.slice_mut().iter_mut(); for (i, &d) in self.slice().iter().enumerate() { - if i == axis { + if i == axis.axis() { continue; } for rr in it.by_ref() { @@ -588,14 +644,14 @@ impl_shrink_recursive!(Ix, Ix, Ix, Ix, Ix, Ix, Ix, Ix, Ix, Ix, Ix, Ix,); impl RemoveAxis for Vec { type Smaller = Vec; - fn remove_axis(&self, axis: usize) -> Vec { + fn remove_axis(&self, axis: Axis) -> Vec { let mut res = self.clone(); - res.remove(axis); + res.remove(axis.axis()); res } } -/// A tuple or fixed size array that can be used to index an array. +/// Tuple or fixed size arrays that can be used to index an array. /// /// ``` /// use ndarray::arr2; @@ -608,7 +664,7 @@ impl RemoveAxis for Vec { /// /// **Note** the blanket implementation that's not visible in rustdoc: /// `impl NdIndex for D where D: Dimension { ... }` -pub unsafe trait NdIndex { +pub unsafe trait NdIndex : Debug { type Dim: Dimension; #[doc(hidden)] fn index_checked(&self, dim: &Self::Dim, strides: &Self::Dim) -> Option; @@ -680,10 +736,11 @@ unsafe impl<'a> NdIndex for &'a [Ix] { } } +// NOTE: These tests are not compiled & tested #[cfg(test)] mod test { use super::Dimension; - use stride_error::StrideError; + use error::StrideError; #[test] fn fastest_varying_order() { @@ -723,7 +780,7 @@ mod test { /// /// All array axis arguments use this type to make the code easier to write /// correctly and easier to understand. -#[derive(Copy, Clone, PartialEq, Eq, PartialOrd, Ord, Hash, Debug)] +#[derive(Copy, Eq, Ord, Hash, Debug)] pub struct Axis(pub usize); impl Axis { @@ -731,3 +788,27 @@ impl Axis { pub fn axis(&self) -> usize { self.0 } } +macro_rules! clone_from_copy { + ($typename:ident) => { + impl Clone for $typename { + #[inline] + fn clone(&self) -> Self { *self } + } + } +} + +macro_rules! derive_cmp { + ($traitname:ident for $typename:ident, $method:ident -> $ret:ty) => { + impl $traitname for $typename { + #[inline(always)] + fn $method(&self, rhs: &Self) -> $ret { + (self.0).$method(&rhs.0) + } + } + } +} + +derive_cmp!{PartialEq for Axis, eq -> bool} +derive_cmp!{PartialOrd for Axis, partial_cmp -> Option} +clone_from_copy!{Axis} + diff --git a/src/error.rs b/src/error.rs new file mode 100644 index 000000000..e5ff22b92 --- /dev/null +++ b/src/error.rs @@ -0,0 +1,105 @@ +// Copyright 2014-2016 bluss and ndarray developers. +// +// Licensed under the Apache License, Version 2.0 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. +use std::fmt; +use std::error::Error; +use super::{ + Dimension, +}; + +/// An error related to array shape or layout. +#[derive(Clone)] +pub struct ShapeError { + // we want to be able to change this representation later + repr: ErrorKind, +} + +impl ShapeError { + /// Return the `ErrorKind` of this error. + #[inline] + pub fn kind(&self) -> ErrorKind { + self.repr + } + + /// Create a new `ShapeError` + pub fn from_kind(error: ErrorKind) -> Self { + from_kind(error) + } +} + +/// Error code for an error related to array shape or layout. +/// +/// This enumeration is not exhaustive. The representation of the enum +/// is not guaranteed. +#[derive(Copy, Clone, Debug)] +pub enum ErrorKind { + /// incompatible shape + IncompatibleShape, + /// incompatible memory layout + IncompatibleLayout, + /// the shape does not fit inside type limits + RangeLimited, + /// out of bounds indexing + OutOfBounds, + /// aliasing array elements + Unsupported, + #[doc(hidden)] + __Incomplete, +} + +#[inline(always)] +pub fn from_kind(k: ErrorKind) -> ShapeError { + ShapeError { + repr: k + } +} + +impl PartialEq for ErrorKind { + #[inline(always)] + fn eq(&self, rhs: &Self) -> bool { + *self as u8 == *rhs as u8 + } +} + +impl PartialEq for ShapeError { + #[inline(always)] + fn eq(&self, rhs: &Self) -> bool { + self.repr == rhs.repr + } +} + +impl Error for ShapeError { + fn description(&self) -> &str { + match self.kind() { + ErrorKind::IncompatibleShape => "incompatible shapes", + ErrorKind::IncompatibleLayout => "incompatible memory layout", + ErrorKind::RangeLimited => "the shape does not fit in type limits", + ErrorKind::OutOfBounds => "out of bounds indexing", + ErrorKind::Unsupported => "unsupported operation", + ErrorKind::__Incomplete => "this error variant is not in use", + } + } +} + +impl fmt::Display for ShapeError { + fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { + self.description().fmt(f) + } +} + +impl fmt::Debug for ShapeError { + fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { + write!(f, "ShapeError {:?}: {}", self.kind(), self.description()) + } +} + +pub fn incompatible_shapes(_a: &D, _b: &E) -> ShapeError + where D: Dimension, + E: Dimension +{ + from_kind(ErrorKind::IncompatibleShape) +} diff --git a/src/free_functions.rs b/src/free_functions.rs index 85db5df30..2d7bab1a7 100644 --- a/src/free_functions.rs +++ b/src/free_functions.rs @@ -1,3 +1,10 @@ +// Copyright 2014-2016 bluss and ndarray developers. +// +// Licensed under the Apache License, Version 2.0 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. use std::slice; @@ -6,7 +13,7 @@ use imp_prelude::*; /// ***Deprecated: Use `ArrayBase::zeros` instead.*** /// -/// Return an array filled with zeros +/// Create an array filled with zeros #[cfg_attr(has_deprecated, deprecated(note="Use `ArrayBase::zeros` instead."))] pub fn zeros(dim: D) -> OwnedArray where A: Clone + libnum::Zero, D: Dimension, @@ -14,28 +21,28 @@ pub fn zeros(dim: D) -> OwnedArray ArrayBase::zeros(dim) } -/// Return a zero-dimensional array with the element `x`. +/// Create a zero-dimensional array with the element `x`. pub fn arr0(x: A) -> OwnedArray { unsafe { ArrayBase::from_vec_dim_unchecked((), vec![x]) } } -/// Return a one-dimensional array with elements from `xs`. +/// Create a one-dimensional array with elements from `xs`. pub fn arr1(xs: &[A]) -> OwnedArray { ArrayBase::from_vec(xs.to_vec()) } -/// Return a one-dimensional array with elements from `xs`. +/// Create a one-dimensional array with elements from `xs`. pub fn rcarr1(xs: &[A]) -> RcArray { arr1(xs).into_shared() } -/// Return a zero-dimensional array view borrowing `x`. +/// Create a zero-dimensional array view borrowing `x`. pub fn aview0(x: &A) -> ArrayView { unsafe { ArrayView::new_(x, (), ()) } } -/// Return a one-dimensional array view with elements borrowing `xs`. +/// Create a one-dimensional array view with elements borrowing `xs`. /// /// ``` /// use ndarray::aview1; @@ -50,10 +57,10 @@ pub fn aview0(x: &A) -> ArrayView { /// ); /// ``` pub fn aview1(xs: &[A]) -> ArrayView { - ArrayView::from_slice(xs) + ArrayView::from(xs) } -/// Return a two-dimensional array view with elements borrowing `xs`. +/// Create a two-dimensional array view with elements borrowing `xs`. pub fn aview2>(xs: &[V]) -> ArrayView { let cols = V::len(); let rows = xs.len(); @@ -67,7 +74,7 @@ pub fn aview2>(xs: &[V]) -> ArrayView>(xs: &[V]) -> ArrayView(xs: &mut [A]) -> ArrayViewMut { - ArrayViewMut::from_slice(xs) + ArrayViewMut::from(xs) } /// Fixed-size array used for array initialization @@ -114,7 +121,7 @@ macro_rules! impl_arr_init { impl_arr_init!(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,); -/// Return a two-dimensional array with elements from `xs`. +/// Create a two-dimensional array with elements from `xs`. /// /// ``` /// use ndarray::arr2; @@ -126,27 +133,24 @@ impl_arr_init!(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,); /// ); /// ``` pub fn arr2>(xs: &[V]) -> OwnedArray { - // FIXME: Simplify this when V is fix size array - let (m, n) = (xs.len() as Ix, - xs.get(0).map_or(0, |snd| snd.as_init_slice().len() as Ix)); + let (m, n) = (xs.len() as Ix, V::len() as Ix); let dim = (m, n); let mut result = Vec::::with_capacity(dim.size()); for snd in xs { - let snd = snd.as_init_slice(); - result.extend(snd.iter().cloned()); + result.extend_from_slice(snd.as_init_slice()); } unsafe { ArrayBase::from_vec_dim_unchecked(dim, result) } } -/// Return a two-dimensional array with elements from `xs`. +/// Create a two-dimensional array with elements from `xs`. /// pub fn rcarr2>(xs: &[V]) -> RcArray { arr2(xs).into_shared() } -/// Return a three-dimensional array with elements from `xs`. +/// Create a three-dimensional array with elements from `xs`. /// /// **Panics** if the slices are not all of the same length. /// @@ -166,19 +170,11 @@ pub fn rcarr2>(xs: &[V]) -> RcArray, U: FixedInitializer>(xs: &[V]) -> OwnedArray { - // FIXME: Simplify this when U/V are fix size arrays - let m = xs.len() as Ix; - let fst = xs.get(0).map(|snd| snd.as_init_slice()); - let thr = fst.and_then(|elt| elt.get(0).map(|elt2| elt2.as_init_slice())); - let n = fst.map_or(0, |v| v.len() as Ix); - let o = thr.map_or(0, |v| v.len() as Ix); - let dim = (m, n, o); + let dim = (xs.len() as Ix, V::len() as Ix, U::len() as Ix); let mut result = Vec::::with_capacity(dim.size()); for snd in xs { - let snd = snd.as_init_slice(); - for thr in snd.iter() { - let thr = thr.as_init_slice(); - result.extend(thr.iter().cloned()); + for thr in snd.as_init_slice() { + result.extend_from_slice(thr.as_init_slice()); } } unsafe { @@ -186,7 +182,7 @@ pub fn arr3, U: FixedInitializer>( } } -/// Return a three-dimensional array with elements from `xs`. +/// Create a three-dimensional array with elements from `xs`. pub fn rcarr3, U: FixedInitializer>(xs: &[V]) -> RcArray { diff --git a/src/impl_clone.rs b/src/impl_clone.rs index 5d834be8f..b756804cd 100644 --- a/src/impl_clone.rs +++ b/src/impl_clone.rs @@ -1,3 +1,10 @@ +// Copyright 2014-2016 bluss and ndarray developers. +// +// Licensed under the Apache License, Version 2.0 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. use imp_prelude::*; use DataClone; diff --git a/src/impl_constructors.rs b/src/impl_constructors.rs index 8fc899d7a..978563ccd 100644 --- a/src/impl_constructors.rs +++ b/src/impl_constructors.rs @@ -1,35 +1,85 @@ +// Copyright 2014-2016 bluss and ndarray developers. +// +// Licensed under the Apache License, Version 2.0 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. //! Constructor methods for ndarray //! +extern crate rand; +use self::rand::distributions::{Normal, IndependentSample}; +use self::rand::{Rng, XorShiftRng, SeedableRng}; use libnum; use imp_prelude::*; use dimension; use linspace; -use shape_error::{self, ShapeError}; -use stride_error::StrideError; +use error::{self, ShapeError}; /// Constructor methods for one-dimensional arrays. +/// +/// Note that the constructor methods apply to `OwnedArray` and `RcArray`, +/// the two array types that have owned storage. impl ArrayBase where S: DataOwned { - /// Create a one-dimensional array from a vector (no allocation needed). + /// Create a one-dimensional array from a vector (no copying needed). + /// + /// ```rust + /// use ndarray::OwnedArray; + /// + /// let array = OwnedArray::from_vec(vec![1., 2., 3., 4.]); + /// ``` pub fn from_vec(v: Vec) -> ArrayBase { unsafe { Self::from_vec_dim_unchecked(v.len() as Ix, v) } } /// Create a one-dimensional array from an iterable. - pub fn from_iter>(iterable: I) -> ArrayBase { + /// + /// ```rust + /// use ndarray::{OwnedArray, arr1}; + /// + /// let array = OwnedArray::from_iter((0..5).map(|x| x * x)); + /// assert!(array == arr1(&[0, 1, 4, 9, 16])) + /// ``` + pub fn from_iter(iterable: I) -> ArrayBase + where I: IntoIterator + { Self::from_vec(iterable.into_iter().collect()) } - /// Create a one-dimensional array from inclusive interval + /// Create a one-dimensional array from the inclusive interval /// `[start, end]` with `n` elements. `F` must be a floating point type. + /// + /// ```rust + /// use ndarray::{OwnedArray, arr1}; + /// + /// let array = OwnedArray::linspace(0., 1., 5); + /// assert!(array == arr1(&[0.0, 0.25, 0.5, 0.75, 1.0])) + /// ``` pub fn linspace(start: F, end: F, n: usize) -> ArrayBase where S: Data, F: libnum::Float, { - Self::from_iter(linspace::linspace(start, end, n)) + Self::from_vec(::iterators::to_vec(linspace::linspace(start, end, n))) + } + + /// Create a one-dimensional array from the half-open interval + /// `[start, end)` with elements spaced by `step`. `F` must be a floating point type. + /// + /// ```rust + /// use ndarray::{OwnedArray, arr1}; + /// + /// let array = OwnedArray::range(0., 5., 1.); + /// assert!(array == arr1(&[0., 1., 2., 3., 4.])) + /// ``` + pub fn range(start: F, end: F, step: F) -> ArrayBase + where S: Data, + F: libnum::Float, + { + Self::from_vec(::iterators::to_vec(linspace::range(start, end, step))) } } @@ -52,7 +102,7 @@ impl ArrayBase } } -/// Constructor methods for arrays. +/// Constructor methods for n-dimensional arrays. impl ArrayBase where S: DataOwned, D: Dimension, @@ -62,10 +112,10 @@ impl ArrayBase /// **Panics** if the number of elements in `dim` would overflow usize. /// /// ``` - /// use ndarray::RcArray; + /// use ndarray::OwnedArray; /// use ndarray::arr3; /// - /// let a = RcArray::from_elem((2, 2, 2), 1.); + /// let a = OwnedArray::from_elem((2, 2, 2), 1.); /// /// assert!( /// a == arr3(&[[[1., 1.], @@ -73,6 +123,7 @@ impl ArrayBase /// [[1., 1.], /// [1., 1.]]]) /// ); + /// assert!(a.strides() == &[4, 2, 1]); /// ``` pub fn from_elem(dim: D, elem: A) -> ArrayBase where A: Clone @@ -91,17 +142,9 @@ impl ArrayBase /// **Panics** if the number of elements would overflow usize. /// /// ``` - /// use ndarray::RcArray; - /// use ndarray::arr3; + /// use ndarray::OwnedArray; /// - /// let a = RcArray::from_elem_f((2, 2, 2), 1.); - /// - /// assert!( - /// a == arr3(&[[[1., 1.], - /// [1., 1.]], - /// [[1., 1.], - /// [1., 1.]]]) - /// ); + /// let a = OwnedArray::from_elem_f((2, 2, 2), 1.); /// assert!(a.strides() == &[1, 2, 4]); /// ``` pub fn from_elem_f(dim: D, elem: A) -> ArrayBase @@ -140,18 +183,17 @@ impl ArrayBase unsafe { Self::from_vec_dim_unchecked(dim, v) } } - /// Create an array from a vector (with no allocation needed). + /// Create an array from a vector (no copying needed). /// - /// **Errors** if `dim` does not correspond to the number of elements - /// in `v`. + /// **Errors** if `dim` does not correspond to the number of elements in `v`. pub fn from_vec_dim(dim: D, v: Vec) -> Result, ShapeError> { if dim.size_checked() != Some(v.len()) { - return Err(shape_error::incompatible_shapes(&v.len(), &dim)); + return Err(error::incompatible_shapes(&v.len(), &dim)); } unsafe { Ok(Self::from_vec_dim_unchecked(dim, v)) } } - /// Create an array from a vector (with no allocation needed). + /// Create an array from a vector (no copying needed). /// /// Unsafe because dimension is unchecked, and must be correct. pub unsafe fn from_vec_dim_unchecked(dim: D, mut v: Vec) -> ArrayBase { @@ -164,7 +206,7 @@ impl ArrayBase } } - /// Create an array from a vector (with no allocation needed), + /// Create an array from a vector (with no copying needed), /// using fortran memory order to interpret the data. /// /// Unsafe because dimension is unchecked, and must be correct. @@ -187,7 +229,7 @@ impl ArrayBase /// **Errors** if strides and dimensions can point out of bounds of `v`.
/// **Errors** if strides allow multiple indices to point to the same element. pub fn from_vec_dim_stride(dim: D, strides: D, v: Vec
) - -> Result, StrideError> + -> Result, ShapeError> { dimension::can_index_slice(&v, &dim, &strides).map(|_| { unsafe { @@ -214,3 +256,27 @@ impl ArrayBase } +use super::NdFloat; + +impl ArrayBase + where S: DataOwned, + A: NdFloat, + D: Dimension, +{ + pub fn randn_rng(dim: D, rng: &mut R) -> ArrayBase, D> + { + let size = dim.size_checked().expect("Shape too large: overflow in size"); + let normal = Normal::new(1.0, 0.0); + let iter = (0..size).map(|_| A::from(normal.ind_sample(rng)).unwrap()); + let arr = OwnedArray::from_iter(iter); + arr.into_shape(dim).unwrap() + } + + pub fn randn(dim: D) -> ArrayBase, D> + { + // Use the thread_rng once to seed a fast, non-crypto grade rng. + // To use a crypto quality RNG, use randn_rng directly. + let mut rng : XorShiftRng = SeedableRng::from_seed(rand::thread_rng().gen::<[u32;4]>()); + ArrayBase::, D>::randn_rng(dim, &mut rng) + } +} diff --git a/src/impl_linalg.rs b/src/impl_linalg.rs index 036c8bb33..de8298765 100644 --- a/src/impl_linalg.rs +++ b/src/impl_linalg.rs @@ -1,5 +1,12 @@ +// Copyright 2014-2016 bluss and ndarray developers. +// +// Licensed under the Apache License, Version 2.0 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. -use libnum; +use libnum::Zero; use itertools::free::enumerate; use imp_prelude::*; @@ -9,6 +16,28 @@ use { LinalgScalar, }; +#[cfg(feature="blas")] +use std::mem::swap; +#[cfg(feature="blas")] +use std::os::raw::c_int; +#[cfg(feature="blas")] +use std::any::{Any, TypeId}; + +#[cfg(feature="blas")] +use blas_sys::c::{CblasNoTrans, CblasTrans, CblasRowMajor}; +#[cfg(feature="blas")] +use blas_sys; + +/// len of vector before we use blas +#[cfg(feature="blas")] +const DOT_BLAS_CUTOFF: usize = 32; +/// side of matrix before we use blas +#[cfg(feature="blas")] +const GEMM_BLAS_CUTOFF: usize = 7; +#[cfg(feature="blas")] +#[allow(non_camel_case_types)] +type blas_index = c_int; // blas index type + impl ArrayBase where S: Data, { @@ -29,7 +58,8 @@ impl ArrayBase where S2: Data, A: LinalgScalar, { - assert_eq!(self.len(), rhs.len()); + debug_assert_eq!(self.len(), rhs.len()); + assert!(self.len() == rhs.len()); if let Some(self_s) = self.as_slice() { if let Some(rhs_s) = rhs.as_slice() { return numeric_util::unrolled_dot(self_s, rhs_s); @@ -44,7 +74,7 @@ impl ArrayBase sum } - #[cfg(not(feature="rblas"))] + #[cfg(not(feature="blas"))] fn dot_impl(&self, rhs: &ArrayBase) -> A where S2: Data, A: LinalgScalar, @@ -52,39 +82,36 @@ impl ArrayBase self.dot_generic(rhs) } - #[cfg(feature="rblas")] + #[cfg(feature="blas")] fn dot_impl(&self, rhs: &ArrayBase) -> A where S2: Data, A: LinalgScalar, { - use std::any::{Any, TypeId}; - use rblas::vector::ops::Dot; - use linalg::AsBlasAny; - - // Read pointer to type `A` as type `B`. - // - // **Panics** if `A` and `B` are not the same type - fn cast_as(a: &A) -> B { - assert_eq!(TypeId::of::(), TypeId::of::()); - unsafe { - ::std::ptr::read(a as *const _ as *const B) - } - } // Use only if the vector is large enough to be worth it - if self.len() >= 32 { - assert_eq!(self.len(), rhs.len()); - if let Ok(self_v) = self.blas_view_as_type::() { - if let Ok(rhs_v) = rhs.blas_view_as_type::() { - let f_ret = f32::dot(&self_v, &rhs_v); - return cast_as::(&f_ret); - } + if self.len() >= DOT_BLAS_CUTOFF { + debug_assert_eq!(self.len(), rhs.len()); + assert!(self.len() == rhs.len()); + macro_rules! dot { + ($ty:ty, $func:ident) => {{ + if blas_compat_1d::<$ty, _>(self) && blas_compat_1d::<$ty, _>(rhs) { + let n = self.len() as blas_index; + let incx = self.strides()[0] as blas_index; + let incy = rhs.strides()[0] as blas_index; + let ret = unsafe { + blas_sys::c::$func( + n, + self.ptr as *const $ty, + incx, + rhs.ptr as *const $ty, + incy) + }; + return cast_as::<$ty, A>(&ret); } - if let Ok(self_v) = self.blas_view_as_type::() { - if let Ok(rhs_v) = rhs.blas_view_as_type::() { - let f_ret = f64::dot(&self_v, &rhs_v); - return cast_as::(&f_ret); - } + }} } + + dot!{f32, cblas_sdot}; + dot!{f64, cblas_ddot}; } self.dot_generic(rhs) } @@ -154,36 +181,12 @@ impl ArrayBase pub fn mat_mul(&self, rhs: &ArrayBase) -> OwnedArray where A: LinalgScalar, { - // NOTE: Matrix multiplication only defined for Copy types to - // avoid trouble with panicking + and *, and destructors - let ((m, a), (b, n)) = (self.dim, rhs.dim); - let (self_columns, other_rows) = (a, b); - assert!(self_columns == other_rows); + let (lhs_columns, rhs_rows) = (a, b); + assert!(lhs_columns == rhs_rows); + assert!(m.checked_mul(n).is_some()); - // Avoid initializing the memory in vec -- set it during iteration - // Panic safe because A: Copy - let mut res_elems = Vec::::with_capacity(m as usize * n as usize); - unsafe { - res_elems.set_len(m as usize * n as usize); - } - let mut i = 0; - let mut j = 0; - for rr in &mut res_elems { - unsafe { - *rr = (0..a).fold(libnum::zero::(), - move |s, k| s + *self.uget((i, k)) * *rhs.uget((k, j)) - ); - } - j += 1; - if j == n { - j = 0; - i += 1; - } - } - unsafe { - ArrayBase::from_vec_dim_unchecked((m, n), res_elems) - } + mat_mul_impl(self, rhs) } /// Perform the matrix multiplication of the rectangular array `self` and @@ -209,7 +212,7 @@ impl ArrayBase } for (i, rr) in enumerate(&mut res_elems) { unsafe { - *rr = (0..a).fold(libnum::zero::(), + *rr = (0..a).fold(A::zero(), move |s, k| s + *self.uget((i, k)) * *rhs.uget(k) ); } @@ -220,4 +223,210 @@ impl ArrayBase } } +#[cfg(not(feature="blas"))] +use self::mat_mul_general as mat_mul_impl; + +#[cfg(feature="blas")] +fn mat_mul_impl(lhs: &ArrayBase, rhs: &ArrayBase) + -> OwnedArray + where A: LinalgScalar, + S: Data, +{ + // size cutoff for using BLAS + let cut = GEMM_BLAS_CUTOFF; + let ((mut m, a), (_, mut n)) = (lhs.dim, rhs.dim); + if !(m > cut || n > cut || a > cut) || + !(same_type::() || same_type::()) { + return mat_mul_general(lhs, rhs); + } + // Use `c` for c-order and `f` for an f-order matrix + // We can handle c * c, f * f generally and + // c * f and f * c if the `f` matrix is square. + let mut lhs_ = lhs.view(); + let mut rhs_ = rhs.view(); + let lhs_s0 = lhs_.strides()[0]; + let rhs_s0 = rhs_.strides()[0]; + let both_f = lhs_s0 == 1 && rhs_s0 == 1; + let mut lhs_trans = CblasNoTrans; + let mut rhs_trans = CblasNoTrans; + if both_f { + // A^t B^t = C^t => B A = C + lhs_ = lhs_.reversed_axes(); + rhs_ = rhs_.reversed_axes(); + swap(&mut lhs_, &mut rhs_); + swap(&mut m, &mut n); + } else if lhs_s0 == 1 && m == a { + lhs_ = lhs_.reversed_axes(); + lhs_trans = CblasTrans; + } else if rhs_s0 == 1 && a == n { + rhs_ = rhs_.reversed_axes(); + rhs_trans = CblasTrans; + } + + macro_rules! gemm { + ($ty:ty, $gemm:ident) => { + if blas_row_major_2d::<$ty, _>(&lhs_) && blas_row_major_2d::<$ty, _>(&rhs_) { + let mut elems = Vec::::with_capacity(m * n); + let c; + unsafe { + elems.set_len(m * n); + c = OwnedArray::from_vec_dim_unchecked((m, n), elems); + } + { + let (m, k) = match lhs_trans { + CblasNoTrans => lhs_.dim(), + _ => { + let (rows, cols) = lhs_.dim(); + (cols, rows) + } + }; + let n = match rhs_trans { + CblasNoTrans => rhs_.dim().1, + _ => rhs_.dim().0, + }; + // gemm is C ← αA^Op B^Op + βC + // Where Op is notrans/trans/conjtrans + unsafe { + blas_sys::c::$gemm( + CblasRowMajor, + lhs_trans, + rhs_trans, + m as blas_index, // m, rows of Op(a) + n as blas_index, // n, cols of Op(b) + k as blas_index, // k, cols of Op(a) + 1.0, // alpha + lhs_.ptr as *const _, // a + lhs_.strides()[0] as blas_index, // lda + rhs_.ptr as *const _, // b + rhs_.strides()[0] as blas_index, // ldb + 0.0, // beta + c.ptr as *mut _, // c + c.strides()[0] as blas_index, // ldc + ); + } + } + return if both_f { + c.reversed_axes() + } else { + c + }; + } + } + } + gemm!(f32, cblas_sgemm); + gemm!(f64, cblas_dgemm); + return mat_mul_general(lhs, rhs); +} + +fn mat_mul_general(lhs: &ArrayBase, rhs: &ArrayBase) + -> OwnedArray + where A: LinalgScalar, + S: Data, +{ + let ((m, a), (_, n)) = (lhs.dim, rhs.dim); + + let lhs_s0 = lhs.strides()[0]; + let rhs_s0 = rhs.strides()[0]; + let column_major = lhs_s0 == 1 && rhs_s0 == 1; + + // Avoid initializing the memory in vec -- set it during iteration + // Panic safe because A: Copy + let mut res_elems = Vec::::with_capacity(m as usize * n as usize); + unsafe { + res_elems.set_len(m as usize * n as usize); + } + let mut i = 0; + let mut j = 0; + for rr in &mut res_elems { + unsafe { + *rr = (0..a).fold(A::zero(), + move |s, k| s + *lhs.uget((i, k)) * *rhs.uget((k, j))); + } + if !column_major { + j += 1; + if j == n { + j = 0; + i += 1; + } + } else { + i += 1; + if i == m { + i = 0; + j += 1; + } + } + } + unsafe { + if !column_major { + ArrayBase::from_vec_dim_unchecked((m, n), res_elems) + } else { + ArrayBase::from_vec_dim_unchecked_f((m, n), res_elems) + } + } +} +#[cfg(feature="blas")] +#[inline(always)] +/// Return `true` if `A` and `B` are the same type +fn same_type() -> bool { + TypeId::of::() == TypeId::of::() +} + +#[cfg(feature="blas")] +// Read pointer to type `A` as type `B`. +// +// **Panics** if `A` and `B` are not the same type +fn cast_as(a: &A) -> B { + assert!(same_type::()); + unsafe { + ::std::ptr::read(a as *const _ as *const B) + } +} + +#[cfg(feature="blas")] +fn blas_compat_1d(a: &ArrayBase) -> bool + where S: Data, + A: Any, + S::Elem: Any, +{ + if !same_type::() { + return false; + } + if a.len() > blas_index::max_value() as usize { + return false; + } + let stride = a.strides()[0]; + if stride > blas_index::max_value() as isize || + stride < blas_index::min_value() as isize { + return false; + } + true +} + +#[cfg(feature="blas")] +fn blas_row_major_2d(a: &ArrayBase) -> bool + where S: Data, + A: Any, + S::Elem: Any, +{ + if !same_type::() { + return false; + } + let s0 = a.strides()[0]; + let s1 = a.strides()[1]; + if s1 != 1 { + return false; + } + if (s0 > blas_index::max_value() as isize || s0 < blas_index::min_value() as isize) || + (s1 > blas_index::max_value() as isize || s1 < blas_index::min_value() as isize) + { + return false; + } + let (m, n) = a.dim(); + if m > blas_index::max_value() as usize || + n > blas_index::max_value() as usize + { + return false; + } + true +} diff --git a/src/impl_methods.rs b/src/impl_methods.rs index cbc16a8e7..fa9e00227 100644 --- a/src/impl_methods.rs +++ b/src/impl_methods.rs @@ -1,12 +1,20 @@ +// Copyright 2014-2016 bluss and ndarray developers. +// +// Licensed under the Apache License, Version 2.0 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. use std::cmp; use std::slice; use imp_prelude::*; +use arraytraits; use dimension; use iterators; -use shape_error::{self, ShapeError}; +use error::{self, ShapeError}; use super::zipsl; use { NdIndex, @@ -234,9 +242,7 @@ impl ArrayBase where S: Data, D: Dimension /// **Note:** only unchecked for non-debug builds of ndarray. #[inline] pub unsafe fn uget(&self, index: D) -> &A { - debug_assert!(self.dim - .stride_offset_checked(&self.strides, &index) - .is_some()); + arraytraits::debug_bounds_check(self, &index); let off = Dimension::stride_offset(&index, &self.strides); &*self.ptr.offset(off) } @@ -252,9 +258,7 @@ impl ArrayBase where S: Data, D: Dimension where S: DataMut { debug_assert!(self.data.is_unique()); - debug_assert!(self.dim - .stride_offset_checked(&self.strides, &index) - .is_some()); + arraytraits::debug_bounds_check(self, &index); let off = Dimension::stride_offset(&index, &self.strides); &mut *self.ptr.offset(off) } @@ -267,7 +271,7 @@ impl ArrayBase where S: Data, D: Dimension /// **Panics** if `axis` or `index` is out of bounds. /// /// ``` - /// use ndarray::{arr1, arr2, Axis}; + /// use ndarray::{arr2, ArrayView, Axis}; /// /// let a = arr2(&[[1., 2.], // -- axis 0, row 0 /// [3., 4.], // -- axis 0, row 1 @@ -276,8 +280,8 @@ impl ArrayBase where S: Data, D: Dimension /// // \ axis 1, column 1 /// // axis 1, column 0 /// assert!( - /// a.subview(Axis(0), 1) == arr1(&[3., 4.]) && - /// a.subview(Axis(1), 1) == arr1(&[2., 4., 6.]) + /// a.subview(Axis(0), 1) == ArrayView::from(&[3., 4.]) && + /// a.subview(Axis(1), 1) == ArrayView::from(&[2., 4., 6.]) /// ); /// ``` pub fn subview(&self, axis: Axis, index: Ix) @@ -331,7 +335,6 @@ impl ArrayBase where S: Data, D: Dimension where D: RemoveAxis, { self.isubview(axis, index); - let axis = axis.axis(); // don't use reshape -- we always know it will fit the size, // and we can use remove_axis on the strides as well ArrayBase { @@ -558,6 +561,19 @@ impl ArrayBase where S: Data, D: Dimension true } + #[inline(always)] + pub fn as_ptr(&self) -> *const A { + self.ptr + } + + #[inline(always)] + pub fn as_mut_ptr(&mut self) -> *mut A + where S: DataMut + { + self.ensure_unique(); // for RcArray + self.ptr + } + /// Return the array’s data as a slice, if it is contiguous and /// the element order corresponds to the memory order. Return `None` otherwise. pub fn as_slice(&self) -> Option<&[A]> { @@ -608,7 +624,7 @@ impl ArrayBase where S: Data, D: Dimension E: Dimension, { if shape.size_checked() != Some(self.dim.size()) { - panic!("Incompatible shapes in reshape, attempted from: {:?}, to: {:?}", + panic!("ndarray: incompatible shapes in reshape, attempted from: {:?}, to: {:?}", self.dim.slice(), shape.slice()) } @@ -650,7 +666,7 @@ impl ArrayBase where S: Data, D: Dimension where E: Dimension { if shape.size_checked() != Some(self.dim.size()) { - return Err(shape_error::incompatible_shapes(&self.dim, &shape)); + return Err(error::incompatible_shapes(&self.dim, &shape)); } // Check if contiguous, if not => copy all, else just adapt strides if self.is_standard_layout() { @@ -661,7 +677,7 @@ impl ArrayBase where S: Data, D: Dimension dim: shape, }) } else { - Err(ShapeError::IncompatibleLayout) + Err(error::from_kind(error::ErrorKind::IncompatibleLayout)) } } @@ -757,7 +773,7 @@ impl ArrayBase where S: Data, D: Dimension where D: Dimension, E: Dimension, { - panic!("Could not broadcast array from shape: {:?} to: {:?}", + panic!("ndarray: could not broadcast array from shape: {:?} to: {:?}", from.slice(), to.slice()) } @@ -804,8 +820,9 @@ impl ArrayBase where S: Data, D: Dimension /// **Note:** Data memory order may not correspond to the index order /// of the array. Neither is the raw data slice is restricted to just the /// array’s view.
- /// **Note:** the slice may be empty. - pub fn raw_data(&self) -> &[A] { + pub fn raw_data(&self) -> &[A] + where S: DataOwned, + { self.data.slice() } @@ -814,12 +831,11 @@ impl ArrayBase where S: Data, D: Dimension /// **Note:** Data memory order may not correspond to the index order /// of the array. Neither is the raw data slice is restricted to just the /// array’s view.
- /// **Note:** the slice may be empty. /// /// **Note:** The data is uniquely held and nonaliased /// while it is mutably borrowed. pub fn raw_data_mut(&mut self) -> &mut [A] - where S: DataMut, + where S: DataOwned + DataMut, { self.ensure_unique(); self.data.slice_mut() @@ -984,16 +1000,13 @@ impl ArrayBase where S: Data, D: Dimension /// [false, true]]) /// ); /// ``` - pub fn map<'a, B, F>(&'a self, mut f: F) -> OwnedArray + pub fn map<'a, B, F>(&'a self, f: F) -> OwnedArray where F: FnMut(&'a A) -> B, A: 'a, { - let mut res = Vec::with_capacity(self.dim.size()); - for elt in self.iter() { - res.push(f(elt)) - } + let v = ::iterators::to_vec(self.iter().map(f)); unsafe { - ArrayBase::from_vec_dim_unchecked(self.dim.clone(), res) + ArrayBase::from_vec_dim_unchecked(self.dim.clone(), v) } } } diff --git a/src/impl_numeric.rs b/src/impl_numeric.rs index afeee0c54..98486fa92 100644 --- a/src/impl_numeric.rs +++ b/src/impl_numeric.rs @@ -1,3 +1,10 @@ +// Copyright 2014-2016 bluss and ndarray developers. +// +// Licensed under the Apache License, Version 2.0 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. use std::ops::Add; use libnum::{self, Float}; @@ -7,6 +14,7 @@ use numeric_util; use { LinalgScalar, + aview0, }; impl ArrayBase @@ -87,14 +95,12 @@ impl ArrayBase D: RemoveAxis, { let n = self.shape()[axis.axis()]; - let mut sum = self.sum(axis); - let one = libnum::one::
(); - let mut cnt = one; + let sum = self.sum(axis); + let mut cnt = A::one(); for _ in 1..n { - cnt = cnt + one; + cnt = cnt + A::one(); } - sum.idiv_scalar(&cnt); - sum + sum / &aview0(&cnt) } /// Return `true` if the arrays' elementwise differences are all within diff --git a/src/impl_ops.rs b/src/impl_ops.rs index fe43eea90..40a7d591b 100644 --- a/src/impl_ops.rs +++ b/src/impl_ops.rs @@ -1,3 +1,10 @@ +// Copyright 2014-2016 bluss and ndarray developers. +// +// Licensed under the Apache License, Version 2.0 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. use std::any::Any; use libnum::Complex; @@ -12,11 +19,11 @@ use libnum::Complex; /// and let `C` be an array with mutable data. /// /// `ScalarOperand` determines for which scalars `K` operations `&A @ K`, and `B @ K`, -/// and `C @= K` are defined, as **right hand side** operands, for applicable +/// and `C @= K` are defined, as ***right hand side operands***, for applicable /// arithmetic operators (denoted `@`). /// -/// **Left hand side** scalar operands are implemented differently -/// (one `impl` per concrete scalar type); they are +/// ***Left hand side*** scalar operands are not related to this trait. +/// (They need one `impl` per concrete scalar type); they are /// implemented for the default `ScalarOperand` types, allowing /// operations `K @ &A`, and `K @ B`. /// diff --git a/src/impl_ops_inplace.rs b/src/impl_ops_inplace.rs index fab2da1e3..37f2cf3aa 100644 --- a/src/impl_ops_inplace.rs +++ b/src/impl_ops_inplace.rs @@ -1,3 +1,10 @@ +// Copyright 2014-2016 bluss and ndarray developers. +// +// Licensed under the Apache License, Version 2.0 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. use std::ops::{ Add, Sub, Mul, Div, Rem, Neg, Not, Shr, Shl, @@ -41,7 +48,10 @@ macro_rules! impl_binary_op_inplace( ); ); -/// *In-place* arithmetic operations. +/// In-place arithmetic operations. +/// +/// ***Note: These will be deprecated in favour of overloaded `@=` operators +/// when Rust 1.8 is released.*** impl ArrayBase where S: DataMut, D: Dimension, diff --git a/src/impl_views.rs b/src/impl_views.rs index a67ca80f2..b0098c2a4 100644 --- a/src/impl_views.rs +++ b/src/impl_views.rs @@ -1,28 +1,21 @@ +// Copyright 2014-2016 bluss and ndarray developers. +// +// Licensed under the Apache License, Version 2.0 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. use imp_prelude::*; use dimension::{self, stride_offset}; -use stride_error::StrideError; - -use { - ViewRepr, -}; +use error::ShapeError; /// # Methods for Array Views /// /// Methods for read-only array views `ArrayView<'a, A, D>` -impl<'a, A> ArrayBase, Ix> { - /// Create a one-dimensional read-only array view of the data in `xs`. - #[inline] - pub fn from_slice(xs: &'a [A]) -> Self { - ArrayView { - data: ViewRepr::new(), - ptr: xs.as_ptr() as *mut A, - dim: xs.len(), - strides: 1, - } - } -} - +/// +/// Note that array views implement traits like [`From`][f] and `IntoIterator` too. +/// [f]: #method.from impl<'a, A, D> ArrayBase, D> where D: Dimension, { @@ -51,7 +44,7 @@ impl<'a, A, D> ArrayBase, D> /// assert!(a.strides() == &[1, 4, 2]); /// ``` pub fn from_slice_dim_stride(dim: D, strides: D, xs: &'a [A]) - -> Result + -> Result { dimension::can_index_slice(xs, &dim, &strides).map(|_| { unsafe { @@ -64,6 +57,11 @@ impl<'a, A, D> ArrayBase, D> /// split and one view after the split. /// /// **Panics** if `axis` or `index` is out of bounds. + /// + /// Below, an illustration of `.split_at(Axis(2), 2)` on + /// an array with shape 3 × 5 × 5. + /// + /// pub fn split_at(self, axis: Axis, index: Ix) -> (Self, Self) { @@ -99,19 +97,10 @@ impl<'a, A, D> ArrayBase, D> } /// Methods for read-write array views `ArrayViewMut<'a, A, D>` -impl<'a, A> ArrayBase, Ix> { - /// Create a one-dimensional read-write array view of the data in `xs`. - #[inline] - pub fn from_slice(xs: &'a mut [A]) -> Self { - ArrayViewMut { - data: ViewRepr::new(), - ptr: xs.as_mut_ptr(), - dim: xs.len(), - strides: 1, - } - } -} - +/// +/// Note that array views implement traits like [`From`][f] and `IntoIterator` too. +/// +/// [f]: #method.from impl<'a, A, D> ArrayBase, D> where D: Dimension, { @@ -141,7 +130,7 @@ impl<'a, A, D> ArrayBase, D> /// assert!(a.strides() == &[1, 4, 2]); /// ``` pub fn from_slice_dim_stride(dim: D, strides: D, xs: &'a mut [A]) - -> Result + -> Result { dimension::can_index_slice(xs, &dim, &strides).map(|_| { unsafe { diff --git a/src/indexes.rs b/src/indexes.rs index 6fd4a1590..855b082ac 100644 --- a/src/indexes.rs +++ b/src/indexes.rs @@ -1,3 +1,10 @@ +// Copyright 2014-2016 bluss and ndarray developers. +// +// Licensed under the Apache License, Version 2.0 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. use super::Dimension; /// An iterator over the indexes of an array shape. diff --git a/src/iterators.rs b/src/iterators.rs index 4bdc6ea78..3b9e3eb28 100644 --- a/src/iterators.rs +++ b/src/iterators.rs @@ -1,4 +1,12 @@ +// Copyright 2014-2016 bluss and ndarray developers. +// +// Licensed under the Apache License, Version 2.0 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. use std::marker::PhantomData; +use std::ptr; use super::{Dimension, Ix, Ixs}; use super::{Elements, ElementsRepr, ElementsBase, ElementsBaseMut, ElementsMut, Indexed, IndexedMut}; @@ -8,6 +16,7 @@ use super::{ ArrayView, ArrayViewMut, RemoveAxis, + Axis, }; /// Base for array iterators @@ -422,8 +431,8 @@ fn new_outer_core(v: ArrayBase, axis: usize) index: 0, len: shape, stride: stride, - inner_dim: v.dim.remove_axis(axis), - inner_strides: v.strides.remove_axis(axis), + inner_dim: v.dim.remove_axis(Axis(axis)), + inner_strides: v.strides.remove_axis(Axis(axis)), ptr: v.ptr, } } @@ -842,3 +851,71 @@ pub fn new_chunk_iter_mut(v: ArrayViewMut, axis: usize, size: usize) chunk_iter_impl!(AxisChunksIter, ArrayView); chunk_iter_impl!(AxisChunksIterMut, ArrayViewMut); + + +// Send and Sync +// All the iterators are thread safe the same way the slice's iterator are + +// read-only iterators use Sync => Send rules, same as `std::slice::Iter`. +macro_rules! send_sync_read_only { + ($name:ident) => { + unsafe impl<'a, A, D> Send for $name<'a, A, D> where A: Sync, D: Send { } + unsafe impl<'a, A, D> Sync for $name<'a, A, D> where A: Sync, D: Sync { } + } +} + +// read-write iterators use Send => Send rules, same as `std::slice::IterMut`. +macro_rules! send_sync_read_write { + ($name:ident) => { + unsafe impl<'a, A, D> Send for $name<'a, A, D> where A: Send, D: Send { } + unsafe impl<'a, A, D> Sync for $name<'a, A, D> where A: Sync, D: Sync { } + } +} + +send_sync_read_only!(Elements); +send_sync_read_only!(Indexed); +send_sync_read_only!(InnerIter); +send_sync_read_only!(OuterIter); +send_sync_read_only!(AxisChunksIter); + +send_sync_read_write!(ElementsMut); +send_sync_read_write!(IndexedMut); +send_sync_read_write!(InnerIterMut); +send_sync_read_write!(OuterIterMut); +send_sync_read_write!(AxisChunksIterMut); + +/// (Trait used internally) An iterator that we trust +/// to deliver exactly as many items as it said it would. +pub unsafe trait TrustedIterator { } + +use std::iter; +use linspace::Linspace; + +unsafe impl TrustedIterator for Linspace { } +unsafe impl<'a, A, D> TrustedIterator for Elements<'a, A, D> { } +unsafe impl TrustedIterator for iter::Map + where I: TrustedIterator { } + + +/// Like Iterator::collect, but only for trusted length iterators +pub fn to_vec(iter: I) -> Vec + where I: TrustedIterator + ExactSizeIterator +{ + // Use an `unsafe` block to do this efficiently. + // We know that iter will produce exactly .size() elements, + // and the loop can vectorize if it's clean (without branch to grow the vector). + let (size, _) = iter.size_hint(); + let mut result = Vec::with_capacity(size); + let mut out_ptr = result.as_mut_ptr(); + let mut len = 0; + for elt in iter { + unsafe { + ptr::write(out_ptr, elt); + len += 1; + result.set_len(len); + out_ptr = out_ptr.offset(1); + } + } + debug_assert_eq!(size, result.len()); + result +} diff --git a/src/lib.rs b/src/lib.rs index a8dd049b4..33f97b07d 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,3 +1,10 @@ +// Copyright 2014-2016 bluss and ndarray developers. +// +// Licensed under the Apache License, Version 2.0 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. #![crate_name="ndarray"] #![cfg_attr(has_deprecated, feature(deprecated))] #![doc(html_root_url = "http://bluss.github.io/rust-ndarray/master/")] @@ -24,7 +31,7 @@ //! - Iteration and most operations are efficient on arrays with contiguous //! innermost dimension. //! - Array views can be used to slice and mutate any `[T]` data using -//! `ArrayView::from_slice` and `ArrayViewMut::from_slice`. +//! `ArrayView::from` and `ArrayViewMut::from`. //! //! ## Crate Status //! @@ -40,7 +47,7 @@ //! + `.fold()` and `.zip_mut_with()` are the most efficient ways to //! perform single traversal and lock step traversal respectively. //! + `.iter()` and `.iter_mut()` are efficient for contiguous arrays. -//! - There is experimental bridging to the linear algebra package `rblas`. +//! + Can use BLAS in some operations (`dot` and `mat_mul`). //! //! ## Crate Feature Flags //! @@ -48,17 +55,17 @@ //! `Cargo.toml`. //! //! - `assign_ops` -//! - Optional, requires nightly +//! - Requires Rust 1.8, will be default soon. //! - Enables the compound assignment operators //! - `rustc-serialize` //! - Optional, stable //! - Enables serialization support //! - `rblas` -//! - Optional, stable +//! - ***Deprecated:*** replaced by separate crate `ndarray-rblas` //! - Enables `rblas` integration //! -#![cfg_attr(feature = "assign_ops", feature(augmented_assignments, - op_assign_traits))] +#![cfg_attr(all(feature = "assign_ops", not(has_assign)), + feature(augmented_assignments, op_assign_traits))] #[cfg(feature = "serde")] extern crate serde; @@ -67,6 +74,8 @@ extern crate rustc_serialize as serialize; #[cfg(feature = "rblas")] extern crate rblas; +#[cfg(feature="blas")] +extern crate blas_sys; extern crate itertools; extern crate num as libnum; @@ -85,8 +94,7 @@ pub use dimension::{ pub use dimension::NdIndex; pub use indexes::Indexes; -pub use shape_error::ShapeError; -pub use stride_error::StrideError; +pub use error::{ShapeError, ErrorKind}; pub use si::{Si, S}; use iterators::Baseiter; @@ -99,7 +107,9 @@ pub use iterators::{ AxisChunksIterMut, }; -pub use linalg::LinalgScalar; +pub use arraytraits::AsArray; +pub use linalg::{LinalgScalar, NdFloat}; +pub use stacking::stack; mod arraytraits; #[cfg(feature = "serde")] @@ -128,8 +138,8 @@ mod linalg; mod linspace; mod numeric_util; mod si; -mod shape_error; -mod stride_error; +mod error; +mod stacking; /// Implementation's prelude. Common types used everywhere. mod imp_prelude { @@ -147,7 +157,9 @@ mod imp_prelude { DataMut, DataOwned, DataShared, + ViewRepr, }; + pub use dimension::DimensionExt; /// Wrapper type for private methods #[derive(Copy, Clone, Debug)] pub struct Priv(pub T); @@ -164,13 +176,11 @@ pub type Ixs = isize; /// can be sliced into subsets of its data. /// The array supports arithmetic operations by applying them elementwise. /// -/// The `ArrayBase` is parameterized by: - -/// - `S` for the data container -/// - `D` for the number of dimensions +/// The `ArrayBase` is parameterized by `S` for the data container and +/// `D` for the dimensionality. /// /// Type aliases [`OwnedArray`], [`RcArray`], [`ArrayView`], and [`ArrayViewMut`] refer -/// to `ArrayBase` with different types for the data storage. +/// to `ArrayBase` with different types for the data container. /// /// [`OwnedArray`]: type.OwnedArray.html /// [`RcArray`]: type.RcArray.html @@ -358,7 +368,8 @@ pub type Ixs = isize; /// /// The trait [`ScalarOperand`](trait.ScalarOperand.html) marks types that can be used in arithmetic /// with arrays directly. For a scalar `K` the following combinations of operands -/// are supported (scalar can be on either side). +/// are supported (scalar can be on either the left or right side, but +/// `ScalarOperand` docs has the detailed condtions). /// /// - `&A @ K` or `K @ &A` which produces a new `OwnedArray` /// - `B @ K` or `K @ B` which consumes `B`, updates it with the result and returns it @@ -418,11 +429,25 @@ pub type OwnedArray = ArrayBase, D>; /// A lightweight array view. /// -/// `ArrayView` implements `IntoIterator`. +/// An array view represents an array or a part of it, created from +/// an iterator, subview or slice of an array. +/// +/// Array views have all the methods of an array (see [`ArrayBase`][ab]). +/// +/// See also specific [**Methods for Array Views**](struct.ArrayBase.html#methods-for-array-views). +/// +/// [ab]: struct.ArrayBase.html pub type ArrayView<'a, A, D> = ArrayBase, D>; /// A lightweight read-write array view. /// -/// `ArrayViewMut` implements `IntoIterator`. +/// An array view represents an array or a part of it, created from +/// an iterator, subview or slice of an array. +/// +/// Array views have all the methods of an array (see [`ArrayBase`][ab]). +/// +/// See also specific [**Methods for Array Views**](struct.ArrayBase.html#methods-for-array-views). +/// +/// [ab]: struct.ArrayBase.html pub type ArrayViewMut<'a, A, D> = ArrayBase, D>; /// Array view’s representation. diff --git a/src/linalg.rs b/src/linalg.rs index ae979c29b..a4c57e257 100644 --- a/src/linalg.rs +++ b/src/linalg.rs @@ -1,20 +1,25 @@ -use libnum::{Zero, One}; -use std::ops::{Add, Sub, Mul, Div}; +// Copyright 2014-2016 bluss and ndarray developers. +// +// Licensed under the Apache License, Version 2.0 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. +use libnum::{Zero, One, Float}; use std::any::Any; +use std::fmt; +use std::ops::{Add, Sub, Mul, Div}; +#[cfg(feature="assign_ops")] +use std::ops::{ + AddAssign, + SubAssign, + MulAssign, + DivAssign, + RemAssign, +}; +use ScalarOperand; -#[cfg(feature="rblas")] -use std::any::TypeId; - -#[cfg(feature="rblas")] -use ShapeError; - -#[cfg(feature="rblas")] -use blas::{AsBlas, BlasArrayView}; - -#[cfg(feature="rblas")] -use imp_prelude::*; - -/// Trait union for scalars (array elements) that support linear algebra operations. +/// Elements that support linear algebra operations. /// /// `Any` for type-based specialization, `Copy` so that they don't need move /// semantics or destructors, and the rest are numerical traits. @@ -39,30 +44,34 @@ impl LinalgScalar for T Div { } -#[cfg(feature = "rblas")] -pub trait AsBlasAny : AsBlas { - fn blas_view_as_type(&self) -> Result, ShapeError> - where S: Data; -} +/// Floating-point element types `f32` and `f64`. +/// +/// Trait `NdFloat` is only implemented for `f32` and `f64` but encompasses as +/// much float-relevant ndarray functionality as possible, including the traits +/// needed for linear algebra (`Any`) and for *right hand side* scalar +/// operations (`ScalarOperand`). +#[cfg(not(feature="assign_ops"))] +pub trait NdFloat : + Float + + fmt::Display + fmt::Debug + fmt::LowerExp + fmt::UpperExp + + ScalarOperand + LinalgScalar +{ } + +/// Floating-point element types `f32` and `f64`. +/// +/// Trait `NdFloat` is only implemented for `f32` and `f64` but encompasses as +/// much float-relevant ndarray functionality as possible, including the traits +/// needed for linear algebra (`Any`) and for *right hand side* scalar +/// operations (`ScalarOperand`). +/// +/// This trait can only be implemented by `f32` and `f64`. +#[cfg(feature="assign_ops")] +pub trait NdFloat : + Float + + AddAssign + SubAssign + MulAssign + DivAssign + RemAssign + + fmt::Display + fmt::Debug + fmt::LowerExp + fmt::UpperExp + + ScalarOperand + LinalgScalar +{ } -#[cfg(feature = "rblas")] -/// ***Requires `features = "rblas"`*** -impl AsBlasAny for ArrayBase - where S: Data, - D: Dimension, - A: Any, -{ - fn blas_view_as_type(&self) -> Result, ShapeError> - where S: Data - { - if TypeId::of::() == TypeId::of::() { - unsafe { - let v = self.view(); - let u = ArrayView::new_(v.ptr as *const T, v.dim, v.strides); - Priv(u).into_blas_view() - } - } else { - Err(ShapeError::IncompatibleLayout) - } - } -} +impl NdFloat for f32 { } +impl NdFloat for f64 { } diff --git a/src/linspace.rs b/src/linspace.rs index a81e577e5..ba1f481b9 100644 --- a/src/linspace.rs +++ b/src/linspace.rs @@ -1,3 +1,10 @@ +// Copyright 2014-2016 bluss and ndarray developers. +// +// Licensed under the Apache License, Version 2.0 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. use libnum::Float; /// An iterator of a sequence of evenly spaced floats. @@ -78,3 +85,24 @@ pub fn linspace(a: F, b: F, n: usize) -> Linspace len: n, } } + +/// Return an iterator of floats spaced by `step`, from +/// the half-open interval [a, b). +/// Numerical reasons can result in `b` being included +/// in the result. +/// +/// Iterator element type is `F`, where `F` must be +/// either `f32` or `f64`. +#[inline] +pub fn range(a: F, b: F, step: F) -> Linspace + where F: Float +{ + let len = b - a; + let steps = F::ceil(len / step); + Linspace { + start: a, + step: step, + len: steps.to_usize().unwrap(), + index: 0, + } +} diff --git a/src/numeric_util.rs b/src/numeric_util.rs index 7752fa1cb..2e0af2c93 100644 --- a/src/numeric_util.rs +++ b/src/numeric_util.rs @@ -1,3 +1,10 @@ +// Copyright 2014-2016 bluss and ndarray developers. +// +// Licensed under the Apache License, Version 2.0 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. use libnum; use std::cmp; @@ -34,13 +41,16 @@ pub fn unrolled_sum(mut xs: &[A]) -> A sum = sum.clone() + (p1 + p5); sum = sum.clone() + (p2 + p6); sum = sum.clone() + (p3 + p7); - for elt in xs { - sum = sum.clone() + elt.clone(); + + // make it clear to the optimizer that this loop is short + // and can not be autovectorized. + for i in 0..xs.len() { + if i >= 7 { break; } + sum = sum.clone() + xs[i].clone() } sum } - /// Compute the dot product. /// /// `xs` and `ys` must be the same length @@ -75,8 +85,13 @@ pub fn unrolled_dot(xs: &[A], ys: &[A]) -> A sum = sum + (p1 + p5); sum = sum + (p2 + p6); sum = sum + (p3 + p7); + for i in 0..xs.len() { - sum = sum + xs[i] * ys[i]; + if i >= 7 { break; } + unsafe { + // get_unchecked is needed to avoid the bounds check + sum = sum + xs[i] * *ys.get_unchecked(i); + } } sum } diff --git a/src/shape_error.rs b/src/shape_error.rs deleted file mode 100644 index 24fc5d5db..000000000 --- a/src/shape_error.rs +++ /dev/null @@ -1,55 +0,0 @@ -use std::fmt; -use std::error::Error; -use super::{ - Ix, - Dimension, -}; - -/// An error that can be produced by `.into_shape()` -#[derive(Clone, Debug)] -pub enum ShapeError { - /// incompatible shapes in reshape, (from, to) - IncompatibleShapes(Box<[Ix]>, Box<[Ix]>), - /// incompatible layout: not contiguous - IncompatibleLayout, - /// Dimension too large (shape) - DimensionTooLarge(Box<[Ix]>), -} - -#[inline(never)] -#[cold] -pub fn incompatible_shapes(a: &D, b: &E) -> ShapeError - where D: Dimension, - E: Dimension -{ - ShapeError::IncompatibleShapes(a.slice().to_vec().into_boxed_slice(), - b.slice().to_vec().into_boxed_slice()) -} - - -impl Error for ShapeError { - fn description(&self) -> &str { - match *self { - ShapeError::IncompatibleShapes(..) => "incompatible shapes", - ShapeError::IncompatibleLayout => "incompatible layout (not contiguous)", - ShapeError::DimensionTooLarge(..) => "dimension too large", - } - } -} - -impl fmt::Display for ShapeError { - fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { - match *self { - ShapeError::IncompatibleShapes(ref a, ref b) => { - write!(f, "incompatible shapes, attempted from: {:?}, to: {:?}", - a, b) - } - ShapeError::IncompatibleLayout => { - write!(f, "{}", self.description()) - } - ShapeError::DimensionTooLarge(ref a) => { - write!(f, "dimension too large in shape {:?}", a) - } - } - } -} diff --git a/src/si.rs b/src/si.rs index 4b37f0955..12a5d965b 100644 --- a/src/si.rs +++ b/src/si.rs @@ -1,3 +1,10 @@ +// Copyright 2014-2016 bluss and ndarray developers. +// +// Licensed under the Apache License, Version 2.0 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. use std::ops::{Range, RangeFrom, RangeTo, RangeFull}; use super::Ixs; diff --git a/src/stacking.rs b/src/stacking.rs new file mode 100644 index 000000000..9d83093fb --- /dev/null +++ b/src/stacking.rs @@ -0,0 +1,102 @@ +// Copyright 2014-2016 bluss and ndarray developers. +// +// Licensed under the Apache License, Version 2.0 or the MIT license +// , at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. + +use imp_prelude::*; +use error::{ShapeError, ErrorKind, from_kind}; + +/// Stack arrays along the given axis. +/// +/// ``` +/// use ndarray::{arr2, Axis, stack}; +/// +/// let a = arr2(&[[2., 2.], +/// [3., 3.]]); +/// assert!( +/// stack(Axis(0), &[a.view(), a.view()]) +/// == Ok(arr2(&[[2., 2.], +/// [3., 3.], +/// [2., 2.], +/// [3., 3.]])) +/// ); +/// ``` +pub fn stack<'a, A, D>(axis: Axis, arrays: &[ArrayView<'a, A, D>]) + -> Result, ShapeError> + where A: Copy, + D: RemoveAxis +{ + if arrays.len() == 0 { + return Err(from_kind(ErrorKind::Unsupported)); + } + let mut res_dim = arrays[0].dim(); + if axis.axis() >= res_dim.ndim() { + return Err(from_kind(ErrorKind::OutOfBounds)); + } + let common_dim = res_dim.remove_axis(axis); + if arrays.iter().any(|a| a.dim().remove_axis(axis) != common_dim) { + return Err(from_kind(ErrorKind::IncompatibleShape)); + } + + let stacked_dim = arrays.iter() + .fold(0, |acc, a| acc + a.shape().axis(axis)); + res_dim.set_axis(axis, stacked_dim); + + // we can safely use uninitialized values here because they are Copy + // and we will only ever write to them + let size = res_dim.size(); + let mut v = Vec::with_capacity(size); + unsafe { + v.set_len(size); + } + let mut res = try!(OwnedArray::from_vec_dim(res_dim, v)); + + { + let mut assign_view = res.view_mut(); + for array in arrays { + let len = array.shape().axis(axis); + let (mut front, rest) = assign_view.split_at(axis, len); + front.assign(array); + assign_view = rest; + } + } + Ok(res) +} + +/// Stack arrays along the given axis. +/// +/// Uses the [`stack`][1] function, calling `ArrayView::from(&a)` on each +/// argument `a`. +/// +/// [1]: fn.stack.html +/// +/// ***Panics*** if the `stack` function would return an error. +/// +/// ``` +/// #[macro_use(stack)] +/// extern crate ndarray; +/// +/// use ndarray::{arr2, Axis, stack}; +/// +/// # fn main() { +/// +/// let a = arr2(&[[2., 2.], +/// [3., 3.]]); +/// assert!( +/// stack![Axis(0), a, a] +/// == arr2(&[[2., 2.], +/// [3., 3.], +/// [2., 2.], +/// [3., 3.]]) +/// ); +/// # } +/// ``` +#[macro_export] +macro_rules! stack { + ($axis:expr, $( $array:expr ),+ ) => { + ndarray::stack($axis, &[ $(ndarray::ArrayView::from(&$array) ),* ]).unwrap() + } +} diff --git a/src/stride_error.rs b/src/stride_error.rs deleted file mode 100644 index e42d85b4e..000000000 --- a/src/stride_error.rs +++ /dev/null @@ -1,26 +0,0 @@ -use std::fmt; -use std::error::Error; - -/// An error to describe invalid stride states -#[derive(Clone, Debug, PartialEq)] -pub enum StrideError { - /// stride leads to out of bounds indexing - OutOfBounds, - /// stride leads to aliasing array elements - Unsupported, -} - -impl Error for StrideError { - fn description(&self) -> &str { - match *self { - StrideError::OutOfBounds => "stride leads to out of bounds indexing", - StrideError::Unsupported => "stride leads to aliasing array elements", - } - } -} - -impl fmt::Display for StrideError { - fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { - self.description().fmt(f) - } -} diff --git a/tests/array.rs b/tests/array.rs index 654b66726..9bdefd236 100644 --- a/tests/array.rs +++ b/tests/array.rs @@ -1,8 +1,8 @@ #![allow(non_snake_case)] -#![cfg_attr(feature = "assign_ops", feature(augmented_assignments))] #[macro_use] extern crate ndarray; +extern crate itertools; use ndarray::{RcArray, S, Si, OwnedArray, @@ -68,6 +68,13 @@ fn test_slice() assert!(vi.iter().zip(A.iter()).all(|(a, b)| a == b)); } +#[should_panic] +#[test] +fn index_out_of_bounds() { + let mut a = OwnedArray::::zeros((3, 4)); + a[[3, 2]] = 1; +} + #[should_panic] #[test] fn slice_oob() @@ -564,7 +571,7 @@ fn reshape() { } #[test] -#[should_panic(expected = "IncompatibleShapes")] +#[should_panic(expected = "IncompatibleShape")] fn reshape_error1() { let data = [1, 2, 3, 4, 5, 6, 7, 8]; let v = aview1(&data); @@ -706,3 +713,51 @@ fn deny_split_at_index_out_of_bounds() { let a = arr2(&[[1., 2.], [3., 4.]]); a.view().split_at(Axis(1), 3); } + +#[test] +fn test_range() { + let a = OwnedArray::range(0., 5., 1.); + assert_eq!(a.len(), 5); + assert_eq!(a[0], 0.); + assert_eq!(a[4], 4.); + + let b = OwnedArray::range(0., 2.2, 1.); + assert_eq!(b.len(), 3); + assert_eq!(b[0], 0.); + assert_eq!(b[2], 2.); + + let c = OwnedArray::range(0., 5., 2.); + assert_eq!(c.len(), 3); + assert_eq!(c[0], 0.); + assert_eq!(c[1], 2.); + assert_eq!(c[2], 4.); + + let d = OwnedArray::range(1.0, 2.2, 0.1); + assert_eq!(d.len(), 13); + assert_eq!(d[0], 1.); + assert_eq!(d[10], 2.); + assert_eq!(d[12], 2.2); +} + +#[test] +fn test_f_order() { + // Test that arrays are logically equal in every way, + // even if the underlying memory order is different + let c = arr2(&[[1, 2, 3], + [4, 5, 6]]); + let mut f = OwnedArray::zeros_f(c.dim()); + f.assign(&c); + assert_eq!(f, c); + assert_eq!(f.shape(), c.shape()); + assert_eq!(c.strides(), &[3, 1]); + assert_eq!(f.strides(), &[1, 2]); + itertools::assert_equal(f.iter(), c.iter()); + itertools::assert_equal(f.inner_iter(), c.inner_iter()); + itertools::assert_equal(f.outer_iter(), c.outer_iter()); + itertools::assert_equal(f.axis_iter(Axis(0)), c.axis_iter(Axis(0))); + itertools::assert_equal(f.axis_iter(Axis(1)), c.axis_iter(Axis(1))); + + let dupc = &c + &c; + let dupf = &f + &f; + assert_eq!(dupc, dupf); +} diff --git a/tests/blas.rs b/tests/blas.rs index 7779e7190..7336e61b0 100644 --- a/tests/blas.rs +++ b/tests/blas.rs @@ -1,4 +1,5 @@ #![cfg(feature = "rblas")] +#![allow(deprecated)] extern crate rblas; extern crate num; diff --git a/tests/dimension.rs b/tests/dimension.rs index 84f72a421..0f423eb6f 100644 --- a/tests/dimension.rs +++ b/tests/dimension.rs @@ -11,12 +11,12 @@ use ndarray::{ #[test] fn remove_axis() { - assert_eq!(3.remove_axis(0), ()); - assert_eq!((1, 2).remove_axis(0), 2); - assert_eq!((4, 5, 6).remove_axis(1), (4, 6)); + assert_eq!(3.remove_axis(Axis(0)), ()); + assert_eq!((1, 2).remove_axis(Axis(0)), 2); + assert_eq!((4, 5, 6).remove_axis(Axis(1)), (4, 6)); - assert_eq!(vec![1,2].remove_axis(0), vec![2]); - assert_eq!(vec![4, 5, 6].remove_axis(1), vec![4, 6]); + assert_eq!(vec![1,2].remove_axis(Axis(0)), vec![2]); + assert_eq!(vec![4, 5, 6].remove_axis(Axis(1)), vec![4, 6]); let a = RcArray::::zeros((4,5)); a.subview(Axis(1), 0); @@ -40,4 +40,3 @@ fn dyn_dimension() let z = OwnedArray::::zeros(dim.clone()); assert_eq!(z.shape(), &dim[..]); } - diff --git a/tests/format.rs b/tests/format.rs index 9bfc293fa..c3adb2a2c 100644 --- a/tests/format.rs +++ b/tests/format.rs @@ -7,25 +7,25 @@ use ndarray::{arr0, rcarr1, aview1}; fn formatting() { let a = rcarr1::(&[1., 2., 3., 4.]); - assert_eq!(format!("{:?}", a), + assert_eq!(format!("{}", a), //"[ 1, 2, 3, 4]"); "[1, 2, 3, 4]"); - assert_eq!(format!("{:4?}", a), + assert_eq!(format!("{:4}", a), "[ 1, 2, 3, 4]"); let a = a.reshape((4, 1, 1)); - assert_eq!(format!("{:4?}", a), + assert_eq!(format!("{:4}", a), "[[[ 1]],\n [[ 2]],\n [[ 3]],\n [[ 4]]]"); let a = a.reshape((2, 2)); assert_eq!(format!("{}", a), "[[1, 2],\n [3, 4]]"); - assert_eq!(format!("{:?}", a), + assert_eq!(format!("{}", a), "[[1, 2],\n [3, 4]]"); - assert_eq!(format!("{:#4?}", a), + assert_eq!(format!("{:#4}", a), "[[ 1, 2], [ 3, 4]]"); let b = arr0::(3.5); - assert_eq!(format!("{:?}", b), + assert_eq!(format!("{}", b), "3.5"); let s = format!("{:.3e}", aview1::(&[1.1, 2.2, 33., 440.])); diff --git a/tests/iterators.rs b/tests/iterators.rs index 06c606b3b..d50bf9944 100644 --- a/tests/iterators.rs +++ b/tests/iterators.rs @@ -12,6 +12,7 @@ use ndarray::{ arr2, arr3, Axis, + Indexes, }; use itertools::assert_equal; @@ -454,3 +455,28 @@ fn outer_iter_mut_split_at() { assert_eq!(a[[1, 2, 1]], 12); assert_eq!(a[[4, 2, 1]], 28); } + +#[test] +fn iterators_are_send_sync() { + // When the element type is Send + Sync, then the iterators and views + // are too. + fn _send_sync(_: &T) { } + + let mut a = RcArray::from_iter(0..30).into_shape((5, 3, 2)).unwrap(); + + _send_sync(&a.view()); + _send_sync(&a.view_mut()); + _send_sync(&a.iter()); + _send_sync(&a.iter_mut()); + _send_sync(&a.indexed_iter()); + _send_sync(&a.indexed_iter_mut()); + _send_sync(&a.inner_iter()); + _send_sync(&a.inner_iter_mut()); + _send_sync(&a.outer_iter()); + _send_sync(&a.outer_iter_mut()); + _send_sync(&a.axis_iter(Axis(1))); + _send_sync(&a.axis_iter_mut(Axis(1))); + _send_sync(&a.axis_chunks_iter(Axis(1), 1)); + _send_sync(&a.axis_chunks_iter_mut(Axis(1), 1)); + _send_sync(&Indexes::new(a.dim())); +} diff --git a/tests/oper.rs b/tests/oper.rs index 76fc3e648..69930ec82 100644 --- a/tests/oper.rs +++ b/tests/oper.rs @@ -5,6 +5,7 @@ use ndarray::RcArray; use ndarray::{arr0, rcarr1, rcarr2}; use ndarray::{ OwnedArray, + Ix, }; use std::fmt; @@ -120,14 +121,79 @@ fn assert_approx_eq(f: F, g: F, tol: F) -> bool { #[test] fn dot_product() { - let a = OwnedArray::linspace(0., 63., 64); - let b = OwnedArray::linspace(0., 63., 64); - let dot = 85344.; + let a = OwnedArray::range(0., 69., 1.); + let b = &a * 2. - 7.; + let dot = 197846.; assert_approx_eq(a.dot(&b), dot, 1e-5); let a = a.map(|f| *f as f32); - let b = a.map(|f| *f as f32); + let b = b.map(|f| *f as f32); assert_approx_eq(a.dot(&b), dot as f32, 1e-5); let a = a.map(|f| *f as i32); - let b = a.map(|f| *f as i32); + let b = b.map(|f| *f as i32); assert_eq!(a.dot(&b), dot as i32); } + +fn range_mat(m: Ix, n: Ix) -> OwnedArray { + OwnedArray::linspace(0., (m * n - 1) as f32, m * n).into_shape((m, n)).unwrap() +} + +#[cfg(has_assign)] +#[test] +fn mat_mul() { + let (m, n, k) = (8, 8, 8); + let a = range_mat(m, n); + let b = range_mat(n, k); + let mut b = b / 4.; + { + let mut c = b.column_mut(0); + c += 1.0; + } + let ab = a.mat_mul(&b); + + let mut af = OwnedArray::zeros_f(a.dim()); + let mut bf = OwnedArray::zeros_f(b.dim()); + af.assign(&a); + bf.assign(&b); + + assert_eq!(ab, a.mat_mul(&bf)); + assert_eq!(ab, af.mat_mul(&b)); + assert_eq!(ab, af.mat_mul(&bf)); + + let (m, n, k) = (10, 5, 11); + let a = range_mat(m, n); + let b = range_mat(n, k); + let mut b = b / 4.; + { + let mut c = b.column_mut(0); + c += 1.0; + } + let ab = a.mat_mul(&b); + + let mut af = OwnedArray::zeros_f(a.dim()); + let mut bf = OwnedArray::zeros_f(b.dim()); + af.assign(&a); + bf.assign(&b); + + assert_eq!(ab, a.mat_mul(&bf)); + assert_eq!(ab, af.mat_mul(&b)); + assert_eq!(ab, af.mat_mul(&bf)); +} + +// Check that matrix multiplication of contiguous matrices returns a +// matrix with the same order +#[test] +fn mat_mul_order() { + let (m, n, k) = (8, 8, 8); + let a = range_mat(m, n); + let b = range_mat(n, k); + let mut af = OwnedArray::zeros_f(a.dim()); + let mut bf = OwnedArray::zeros_f(b.dim()); + af.assign(&a); + bf.assign(&b); + + let cc = a.mat_mul(&b); + let ff = af.mat_mul(&bf); + + assert_eq!(cc.strides()[1], 1); + assert_eq!(ff.strides()[0], 1); +} diff --git a/tests/stacking.rs b/tests/stacking.rs new file mode 100644 index 000000000..65be6993f --- /dev/null +++ b/tests/stacking.rs @@ -0,0 +1,44 @@ + +#[macro_use(stack)] +extern crate ndarray; + + +use ndarray::{ + aview1, + arr2, + Axis, + Ix, + OwnedArray, + ErrorKind, +}; + +#[test] +fn stacking() { + let a = arr2(&[[2., 2.], + [3., 3.]]); + let b = ndarray::stack(Axis(0), &[a.view(), a.view()]).unwrap(); + assert_eq!(b, arr2(&[[2., 2.], + [3., 3.], + [2., 2.], + [3., 3.]])); + + let c = stack![Axis(0), a, b]; + assert_eq!(c, arr2(&[[2., 2.], + [3., 3.], + [2., 2.], + [3., 3.], + [2., 2.], + [3., 3.]])); + + let d = stack![Axis(0), a.row(0), &[9., 9.]]; + assert_eq!(d, aview1(&[2., 2., 9., 9.])); + + let res = ndarray::stack(Axis(1), &[a.view(), c.view()]); + assert_eq!(res.unwrap_err().kind(), ErrorKind::IncompatibleShape); + + let res = ndarray::stack(Axis(2), &[a.view(), c.view()]); + assert_eq!(res.unwrap_err().kind(), ErrorKind::OutOfBounds); + + let res: Result, _> = ndarray::stack(Axis(0), &[]); + assert_eq!(res.unwrap_err().kind(), ErrorKind::Unsupported); +}