diff --git a/src/Data/Matlab/Formatter.cs b/src/Data/Matlab/Formatter.cs index 28b5cfad..e1cfd0ac 100644 --- a/src/Data/Matlab/Formatter.cs +++ b/src/Data/Matlab/Formatter.cs @@ -28,7 +28,6 @@ // OTHER DEALINGS IN THE SOFTWARE. // - using System; using System.Collections.Generic; using System.IO; @@ -51,11 +50,6 @@ namespace MathNet.Numerics.Data.Matlab /// const string HeaderText = "MATLAB 5.0 MAT-file, Platform: .NET 4 - Math.NET Numerics, Created on: "; - /// - /// The length of the header text. - /// - const int HeaderTextLength = 116; - /// /// Format a matrix block byte array /// @@ -77,490 +71,104 @@ namespace MathNet.Numerics.Data.Matlab throw new ArgumentException(string.Format(Resources.NameCannotContainASpace, name), "name"); } - if (typeof(T) == typeof(double)) - { - var sparse = matrix as LinearAlgebra.Double.SparseMatrix; - return sparse != null - ? GetSparseDataArray(sparse, name) - : GetDenseDataArray((LinearAlgebra.Double.Matrix)(object)matrix, name); - } - - if (typeof(T) == typeof(float)) - { - var sparse = matrix as LinearAlgebra.Single.SparseMatrix; - return sparse != null - ? GetSparseDataArray(sparse, name) - : GetDenseDataArray((LinearAlgebra.Single.Matrix)(object)matrix, name); - } - - if (typeof(T) == typeof(Complex)) - { - var sparse = matrix as LinearAlgebra.Complex.SparseMatrix; - return sparse != null - ? GetSparseDataArray(sparse, name) - : GetDenseDataArray((LinearAlgebra.Complex.Matrix)(object)matrix, name); - } - - if (typeof(T) == typeof(Complex32)) - { - var sparse = matrix as LinearAlgebra.Complex32.SparseMatrix; - return sparse != null - ? GetSparseDataArray(sparse, name) - : GetDenseDataArray((LinearAlgebra.Complex32.Matrix)(object)matrix, name); - } - - throw new NotSupportedException(); - } - - /// - /// Writes all matrix blocks to a stream. - /// - internal static void FormatFile(Stream stream, IEnumerable matrices) - { - using (var buffer = new BufferedStream(stream)) - using (var writer = new BinaryWriter(buffer)) - { - WriteHeader(writer); - - foreach (var matrix in matrices) - { - // write data type - writer.Write((int)DataType.Compressed); - - WriteCompressedData(writer, matrix.Data); - } - - writer.Flush(); - writer.Close(); - } - } - - /// - /// Writes the matrix tag and name. - /// - /// The writer we are using. - /// The array class we are writing. - /// if set to true if this a complex matrix. - /// The name of the matrix. - /// The number of rows. - /// The columns of columns. - /// The maximum number of non-zero elements. - static void WriteMatrixTagAndName(BinaryWriter writer, ArrayClass arrayClass, bool isComplex, - string name, int rows, int columns, int nzmax) - { - writer.Write((int)DataType.Matrix); - - // add place holder for data size - writer.Write(0); - - // write flag, data type and size - writer.Write((int)DataType.UInt32); - writer.Write(8); - - // write array class and flags - writer.Write((byte)arrayClass); - if (isComplex) - { - writer.Write((byte)ArrayFlags.Complex); - } - else - { - writer.Write((byte)0); - } - - writer.Write((short)0); - writer.Write(nzmax); - - // write dimensions - writer.Write((int)DataType.Int32); - writer.Write(8); - writer.Write(rows); - writer.Write(columns); - - var nameBytes = Encoding.ASCII.GetBytes(name); - - // write name - if (nameBytes.Length > 4) - { - writer.Write((int)DataType.Int8); - writer.Write(nameBytes.Length); - writer.Write(nameBytes); - var pad = 8 - (nameBytes.Length%8); - PadData(writer, pad); - } - else - { - writer.Write((short)DataType.Int8); - writer.Write((short)nameBytes.Length); - writer.Write(nameBytes); - PadData(writer, 4 - nameBytes.Length); - } - } - - /// - /// Gets the dense data array. - /// - /// The matrix to get the data from. - /// The name of the matrix. - /// The matrix data as an array. - static MatlabMatrix GetDenseDataArray(Matrix matrix, string name) - { - using (var stream = new MemoryStream()) - using (var writer = new BinaryWriter(stream)) - { - WriteMatrixTagAndName(writer, ArrayClass.Double, false, name, matrix.RowCount, matrix.ColumnCount, 0); - - // write data - writer.Write((int)DataType.Double); - writer.Write(matrix.RowCount*matrix.ColumnCount*8); - - for (var j = 0; j < matrix.ColumnCount; j++) - { - var column = matrix.Column(j); - foreach (var value in column) - { - writer.Write(value); - } - } - - writer.Flush(); - return new MatlabMatrix(name, stream.ToArray()); - } - - } + var typeT = typeof (T); + bool sparse = matrix.Storage.GetType().GetGenericTypeDefinition() == typeof (SparseCompressedRowMatrixStorage<>); + bool doublePrecision = typeT == typeof (double) || typeT == typeof (Complex); + bool complex = typeT == typeof (Complex) || typeT == typeof (Complex32); - /// - /// Gets the dense data array. - /// - /// The matrix to get the data from. - /// The name of the matrix. - /// The matrix data as an array. - static MatlabMatrix GetDenseDataArray(Matrix matrix, string name) - { - using (var stream = new MemoryStream()) - using (var writer = new BinaryWriter(stream)) + int sparseNonZeroValues = 0; + if (sparse) { - WriteMatrixTagAndName(writer, ArrayClass.Single, false, name, matrix.RowCount, matrix.ColumnCount, 0); - - // write data - int size = matrix.RowCount*matrix.ColumnCount*4; - writer.Write((int)DataType.Single); - writer.Write(size); - - for (var j = 0; j < matrix.ColumnCount; j++) - { - var column = matrix.Column(j); - foreach (var value in column) - { - writer.Write(value); - } - } - - PadData(writer, size%8); - - writer.Flush(); - return new MatlabMatrix(name, stream.ToArray()); + var sparseStorage = (SparseCompressedRowMatrixStorage)matrix.Storage; + sparseNonZeroValues = sparseStorage.ValueCount; } - } - /// - /// Gets the dense data array. - /// - /// The matrix to get the data from. - /// The name of the matrix. - /// The matrix data as an array. - static MatlabMatrix GetDenseDataArray(Matrix matrix, string name) - { using (var stream = new MemoryStream()) using (var writer = new BinaryWriter(stream)) { - WriteMatrixTagAndName(writer, ArrayClass.Double, true, name, matrix.RowCount, matrix.ColumnCount, 0); + // Array Flags tag: data type + size (8 bytes) + writer.Write((int)DataType.UInt32); + writer.Write(8); - // write data - int size = matrix.RowCount*matrix.ColumnCount*8; - writer.Write((int)DataType.Double); - writer.Write(size); + // Array Flags data: flags (byte 3), class (byte 4) (8 bytes) + writer.Write((byte)(sparse ? ArrayClass.Sparse : doublePrecision ? ArrayClass.Double : ArrayClass.Single)); + writer.Write((byte)(complex ? ArrayFlags.Complex : 0)); + writer.Write((short)0); + writer.Write((int)sparseNonZeroValues); - for (var j = 0; j < matrix.ColumnCount; j++) - { - var column = matrix.Column(j); - foreach (var value in column) - { - writer.Write(value.Real); - } - } + // Dimensions Array tag: data type + size (8 bytes) + writer.Write((int)DataType.Int32); + writer.Write(8); - writer.Write((int)DataType.Double); - writer.Write(size); + // Dimensions Array data: row and column count (8 bytes) + writer.Write(matrix.RowCount); + writer.Write(matrix.ColumnCount); - for (var j = 0; j < matrix.ColumnCount; j++) + // Array Name: + var nameBytes = Encoding.ASCII.GetBytes(name); + if (nameBytes.Length > 4) { - var column = matrix.Column(j); - foreach (var value in column) - { - writer.Write(value.Imaginary); - } + // long format + writer.Write((int)DataType.Int8); + writer.Write(nameBytes.Length); + writer.Write(nameBytes); + PadData(writer, 8 - (nameBytes.Length%8)); } - - writer.Flush(); - return new MatlabMatrix(name, stream.ToArray()); - } - } - - /// - /// Gets the dense data array. - /// - /// The matrix to get the data from. - /// The name of the matrix. - /// The matrix data as an array. - static MatlabMatrix GetDenseDataArray(Matrix matrix, string name) - { - using (var stream = new MemoryStream()) - using (var writer = new BinaryWriter(stream)) - { - WriteMatrixTagAndName(writer, ArrayClass.Single, true, name, matrix.RowCount, matrix.ColumnCount, 0); - - // write data - int size = matrix.RowCount*matrix.ColumnCount*4; - writer.Write((int)DataType.Single); - writer.Write(size); - - for (var j = 0; j < matrix.ColumnCount; j++) + else { - var column = matrix.Column(j); - foreach (var value in column) - { - writer.Write(value.Real); - } + // small format + writer.Write((short)DataType.Int8); + writer.Write((short)nameBytes.Length); + writer.Write(nameBytes); + PadData(writer, 4 - nameBytes.Length); } - PadData(writer, size%8); - - writer.Write((int)DataType.Single); - writer.Write(size); - - for (var j = 0; j < matrix.ColumnCount; j++) + if (doublePrecision && !complex) { - var column = matrix.Column(j); - foreach (var value in column) + var sparseMatrix = matrix as LinearAlgebra.Double.SparseMatrix; + if (sparseMatrix != null) { - writer.Write(value.Real); + SparseArrayFormatter.Write(writer, sparseMatrix); } - } - - PadData(writer, size%8); - - writer.Flush(); - return new MatlabMatrix(name, stream.ToArray()); - } - } - - /// - /// Gets the sparse data array. - /// - /// The matrix to get the data from. - /// The name of the matrix. - /// The matrix data as an array. - static MatlabMatrix GetSparseDataArray(LinearAlgebra.Double.SparseMatrix matrix, string name) - { - using (var stream = new MemoryStream()) - using (var writer = new BinaryWriter(stream)) - { - var nzmax = matrix.NonZerosCount; - WriteMatrixTagAndName(writer, ArrayClass.Sparse, false, name, matrix.RowCount, matrix.ColumnCount, nzmax); - - // write ir - writer.Write((int)DataType.Int32); - writer.Write(nzmax*4); - - foreach (var column in matrix.EnumerateColumns()) - { - foreach (var row in column.EnumerateNonZeroIndexed()) + else { - writer.Write(row.Item1); + NumericArrayFormatter.Write(writer, (LinearAlgebra.Double.Matrix)(object)matrix); } } - - // add pad if needed - if (nzmax%2 == 1) + else if (!doublePrecision && !complex) { - writer.Write(0); - } - - // write jc - writer.Write((int)DataType.Int32); - writer.Write((matrix.ColumnCount + 1)*4); - writer.Write(0); - var count = 0; - foreach (var column in matrix.EnumerateColumns()) - { - count += ((SparseVectorStorage)column.Storage).ValueCount; - writer.Write(count); - } - - // add pad if needed - if (matrix.ColumnCount%2 == 0) - { - writer.Write(0); - } - - // write data - writer.Write((int)DataType.Double); - writer.Write(nzmax*8); - - foreach (var column in matrix.EnumerateColumns()) - { - foreach (var row in column.EnumerateNonZeroIndexed()) + var sparseMatrix = matrix as LinearAlgebra.Single.SparseMatrix; + if (sparseMatrix != null) { - writer.Write(row.Item2); + SparseArrayFormatter.Write(writer, sparseMatrix); } - } - - writer.Flush(); - return new MatlabMatrix(name, stream.ToArray()); - } - } - - /// - /// Gets the sparse data array. - /// - /// The matrix to get the data from. - /// The name of the matrix. - /// The matrix data as an array. - static MatlabMatrix GetSparseDataArray(LinearAlgebra.Single.SparseMatrix matrix, string name) - { - using (var stream = new MemoryStream()) - using (var writer = new BinaryWriter(stream)) - { - var nzmax = matrix.NonZerosCount; - WriteMatrixTagAndName(writer, ArrayClass.Sparse, false, name, matrix.RowCount, matrix.ColumnCount, - nzmax); - - // write ir - writer.Write((int)DataType.Int32); - writer.Write(nzmax*4); - - foreach (var column in matrix.EnumerateColumns()) - { - foreach (var row in column.EnumerateNonZeroIndexed()) + else { - writer.Write(row.Item1); + NumericArrayFormatter.Write(writer, (LinearAlgebra.Single.Matrix)(object)matrix); } } - - // add pad if needed - if (nzmax%2 == 1) + else if (doublePrecision) { - writer.Write(0); - } - - // write jc - writer.Write((int)DataType.Int32); - writer.Write((matrix.ColumnCount + 1)*4); - writer.Write(0); - var count = 0; - foreach (var column in matrix.EnumerateColumns()) - { - count += ((SparseVectorStorage)column.Storage).ValueCount; - writer.Write(count); - } - - // add pad if needed - if (matrix.ColumnCount%2 == 0) - { - writer.Write(0); - } - - // write data - writer.Write((int)DataType.Single); - writer.Write(nzmax*4); - - foreach (var column in matrix.EnumerateColumns()) - { - foreach (var row in column.EnumerateNonZeroIndexed()) + var sparseMatrix = matrix as LinearAlgebra.Complex.SparseMatrix; + if (sparseMatrix != null) { - writer.Write(row.Item2); + SparseArrayFormatter.Write(writer, sparseMatrix); } - } - - var pad = nzmax*4%8; - PadData(writer, pad); - - writer.Flush(); - return new MatlabMatrix(name, stream.ToArray()); - } - } - - /// - /// Gets the sparse data array. - /// - /// The matrix to get the data from. - /// The name of the matrix. - /// The matrix data as an array. - static MatlabMatrix GetSparseDataArray(LinearAlgebra.Complex.SparseMatrix matrix, string name) - { - using (var stream = new MemoryStream()) - using (var writer = new BinaryWriter(stream)) - { - var nzmax = matrix.NonZerosCount; - WriteMatrixTagAndName(writer, ArrayClass.Sparse, true, name, matrix.RowCount, matrix.ColumnCount, - nzmax); - - // write ir - writer.Write((int)DataType.Int32); - writer.Write(nzmax*4); - - foreach (var column in matrix.EnumerateColumns()) - { - foreach (var row in column.EnumerateNonZeroIndexed()) + else { - writer.Write(row.Item1); + NumericArrayFormatter.Write(writer, (LinearAlgebra.Complex.Matrix)(object)matrix); } } - - // add pad if needed - if (nzmax%2 == 1) - { - writer.Write(0); - } - - // write jc - writer.Write((int)DataType.Int32); - writer.Write((matrix.ColumnCount + 1)*4); - writer.Write(0); - var count = 0; - foreach (var column in matrix.EnumerateColumns()) + else { - count += ((SparseVectorStorage)column.Storage).ValueCount; - writer.Write(count); - } - - // add pad if needed - if (matrix.ColumnCount%2 == 0) - { - writer.Write(0); - } - - // write data - writer.Write((int)DataType.Double); - writer.Write(nzmax*8); - - foreach (var column in matrix.EnumerateColumns()) - { - foreach (var row in column.EnumerateNonZeroIndexed()) + var sparseMatrix = matrix as LinearAlgebra.Complex32.SparseMatrix; + if (sparseMatrix != null) { - writer.Write(row.Item2.Real); + SparseArrayFormatter.Write(writer, sparseMatrix); } - } - - writer.Write((int)DataType.Double); - writer.Write(nzmax*8); - - foreach (var column in matrix.EnumerateColumns()) - { - foreach (var row in column.EnumerateNonZeroIndexed()) + else { - writer.Write(row.Item2.Imaginary); + NumericArrayFormatter.Write(writer, (LinearAlgebra.Complex32.Matrix)(object)matrix); } } @@ -570,105 +178,43 @@ namespace MathNet.Numerics.Data.Matlab } /// - /// Gets the sparse data array. + /// Writes all matrix blocks to a stream. /// - /// The matrix to get the data from. - /// The name of the matrix. - /// The matrix data as an array. - static MatlabMatrix GetSparseDataArray(LinearAlgebra.Complex32.SparseMatrix matrix, string name) + internal static void FormatFile(Stream stream, IEnumerable matrices) { - using (var stream = new MemoryStream()) - using (var writer = new BinaryWriter(stream)) + using (var buffer = new BufferedStream(stream)) + using (var writer = new BinaryWriter(buffer)) { - var nzmax = matrix.NonZerosCount; - WriteMatrixTagAndName(writer, ArrayClass.Sparse, true, name, matrix.RowCount, matrix.ColumnCount, - nzmax); - - // write ir - writer.Write((int)DataType.Int32); - writer.Write(nzmax*4); - - foreach (var column in matrix.EnumerateColumns()) - { - foreach (var row in column.EnumerateNonZeroIndexed()) - { - writer.Write(row.Item1); - } - } - - // add pad if needed - if (nzmax%2 == 1) - { - writer.Write(0); - } - - // write jc - writer.Write((int)DataType.Int32); - writer.Write((matrix.ColumnCount + 1)*4); - writer.Write(0); - var count = 0; - foreach (var column in matrix.EnumerateColumns()) - { - count += ((SparseVectorStorage)column.Storage).ValueCount; - writer.Write(count); - } + // write header and subsystem data offset (all space) + var header = Encoding.ASCII.GetBytes(HeaderText + DateTime.Now.ToString(Resources.MatlabDateHeaderFormat)); + writer.Write(header); + PadData(writer, 116 - header.Length + 8, 32); - // add pad if needed - if (matrix.ColumnCount%2 == 0) - { - writer.Write(0); - } + // write version + writer.Write((short)0x100); - // write data - writer.Write((int)DataType.Single); - writer.Write(nzmax*4); + // write little endian indicator + writer.Write((byte)0x49); + writer.Write((byte)0x4D); - foreach (var column in matrix.EnumerateColumns()) + foreach (var matrix in matrices) { - foreach (var row in column.EnumerateNonZeroIndexed()) - { - writer.Write(row.Item2.Real); - } - } - - var pad = nzmax*4%8; - PadData(writer, pad); + // write data type + writer.Write((int)DataType.Compressed); - writer.Write((int)DataType.Single); - writer.Write(nzmax*4); + // compress data + var compressedData = PackCompressedBlock(matrix.Data, DataType.Matrix); - foreach (var column in matrix.EnumerateColumns()) - { - foreach (var row in column.EnumerateNonZeroIndexed()) - { - writer.Write(row.Item2.Imaginary); - } + // write compressed data to file + writer.Write(compressedData.Length); + writer.Write(compressedData); } - PadData(writer, pad); - writer.Flush(); - return new MatlabMatrix(name, stream.ToArray()); + writer.Close(); } } - /// - /// Writes the file header. - /// - static void WriteHeader(BinaryWriter writer) - { - var header = Encoding.ASCII.GetBytes(HeaderText + DateTime.Now.ToString(Resources.MatlabDateHeaderFormat)); - writer.Write(header); - PadData(writer, HeaderTextLength - header.Length + 8, 32); - - // write version - writer.Write((short)0x100); - - // write little endian indicator - writer.Write((byte)0x49); - writer.Write((byte)0x4D); - } - /// /// Pads the data with the given byte. /// @@ -684,41 +230,22 @@ namespace MathNet.Numerics.Data.Matlab } /// - /// Writes the compressed data. + /// Packs a compressed block /// - /// The data to write. - static void WriteCompressedData(BinaryWriter writer, byte[] data) - { - // fill in data size - var size = BitConverter.GetBytes(data.Length); - data[4] = size[0]; - data[5] = size[1]; - data[6] = size[2]; - data[7] = size[3]; - - // compress data - var compressedData = CompressData(data); - - // write compressed data to file - writer.Write(compressedData.Length); - writer.Write(compressedData); - } - - /// - /// Compresses the data array. - /// - /// The data to compress. - /// The compressed data. - static byte[] CompressData(byte[] data) + static byte[] PackCompressedBlock(byte[] data, DataType dataType) { var adler = BitConverter.GetBytes(Adler32.Compute(data)); using (var compressedStream = new MemoryStream()) { compressedStream.WriteByte(0x58); compressedStream.WriteByte(0x85); + using (var outputStream = new DeflateStream(compressedStream, CompressionMode.Compress, true)) { + outputStream.Write(BitConverter.GetBytes((int)dataType), 0, 4); + outputStream.Write(BitConverter.GetBytes(data.Length), 0, 4); outputStream.Write(data, 0, data.Length); + outputStream.Flush(); } compressedStream.WriteByte(adler[3]); diff --git a/src/Data/Matlab/Matlab.csproj b/src/Data/Matlab/Matlab.csproj index 76631e7a..df786a5c 100644 --- a/src/Data/Matlab/Matlab.csproj +++ b/src/Data/Matlab/Matlab.csproj @@ -53,13 +53,13 @@ - - + + diff --git a/src/Data/Matlab/MatlabMatrix.cs b/src/Data/Matlab/MatlabMatrix.cs index 0d25634d..6ed6ed49 100644 --- a/src/Data/Matlab/MatlabMatrix.cs +++ b/src/Data/Matlab/MatlabMatrix.cs @@ -30,14 +30,18 @@ namespace MathNet.Numerics.Data.Matlab { + /// + /// MATLAB Matrix Data Element + /// public class MatlabMatrix { + /// Sub-elements of the matrix data element (not including the data element tag) internal byte[] Data { get; private set; } /// Name of the matrix public string Name { get; private set; } - /// Size of the packed matrix in bytes + /// Size of the matrix in bytes public int ByteSize { get { return Data.Length; } diff --git a/src/Data/Matlab/MatlabReader.cs b/src/Data/Matlab/MatlabReader.cs index 5ad36534..eb4a6cf0 100644 --- a/src/Data/Matlab/MatlabReader.cs +++ b/src/Data/Matlab/MatlabReader.cs @@ -37,7 +37,7 @@ using MathNet.Numerics.LinearAlgebra; namespace MathNet.Numerics.Data.Matlab { /// - /// Creates matrices from MATLAB 5 files. + /// Creates matrices from MATLAB Level-5 Mat files. /// public static class MatlabReader { diff --git a/src/Data/Matlab/MatlabWriter.cs b/src/Data/Matlab/MatlabWriter.cs index a12359ee..ece9edbd 100644 --- a/src/Data/Matlab/MatlabWriter.cs +++ b/src/Data/Matlab/MatlabWriter.cs @@ -37,7 +37,7 @@ using MathNet.Numerics.LinearAlgebra; namespace MathNet.Numerics.Data.Matlab { /// - /// Writes matrices to a MATLAB 5 file. + /// Writes matrices to a MATLAB Level-5 Mat file. /// public static class MatlabWriter { diff --git a/src/Data/Matlab/NumericArrayFormatter.cs b/src/Data/Matlab/NumericArrayFormatter.cs new file mode 100644 index 00000000..143e8f21 --- /dev/null +++ b/src/Data/Matlab/NumericArrayFormatter.cs @@ -0,0 +1,150 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2014 Math.NET +// +// 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. +// + +using System.IO; +using System.Numerics; +using MathNet.Numerics.LinearAlgebra; + +namespace MathNet.Numerics.Data.Matlab +{ + internal static class NumericArrayFormatter + { + internal static void Write(BinaryWriter writer, Matrix matrix) + { + // write data + writer.Write((int)DataType.Double); + writer.Write(matrix.RowCount*matrix.ColumnCount*8); + + for (var j = 0; j < matrix.ColumnCount; j++) + { + var column = matrix.Column(j); + foreach (var value in column) + { + writer.Write(value); + } + } + } + + internal static void Write(BinaryWriter writer, Matrix matrix) + { + // write data + int size = matrix.RowCount*matrix.ColumnCount*4; + writer.Write((int)DataType.Single); + writer.Write(size); + + for (var j = 0; j < matrix.ColumnCount; j++) + { + var column = matrix.Column(j); + foreach (var value in column) + { + writer.Write(value); + } + } + + PadData(writer, size%8); + } + + internal static void Write(BinaryWriter writer, Matrix matrix) + { + // write data + int size = matrix.RowCount*matrix.ColumnCount*8; + writer.Write((int)DataType.Double); + writer.Write(size); + + for (var j = 0; j < matrix.ColumnCount; j++) + { + var column = matrix.Column(j); + foreach (var value in column) + { + writer.Write(value.Real); + } + } + + writer.Write((int)DataType.Double); + writer.Write(size); + + for (var j = 0; j < matrix.ColumnCount; j++) + { + var column = matrix.Column(j); + foreach (var value in column) + { + writer.Write(value.Imaginary); + } + } + } + + internal static void Write(BinaryWriter writer, Matrix matrix) + { + // write data + int size = matrix.RowCount*matrix.ColumnCount*4; + writer.Write((int)DataType.Single); + writer.Write(size); + + for (var j = 0; j < matrix.ColumnCount; j++) + { + var column = matrix.Column(j); + foreach (var value in column) + { + writer.Write(value.Real); + } + } + + PadData(writer, size%8); + + writer.Write((int)DataType.Single); + writer.Write(size); + + for (var j = 0; j < matrix.ColumnCount; j++) + { + var column = matrix.Column(j); + foreach (var value in column) + { + writer.Write(value.Real); + } + } + + PadData(writer, size%8); + } + + /// + /// Pads the data with the given byte. + /// + /// Where to write the pad values. + /// The number of bytes to pad. + /// What value to pad with. + static void PadData(BinaryWriter writer, int bytes, byte pad = (byte)0) + { + for (var i = 0; i < bytes; i++) + { + writer.Write(pad); + } + } + } +} diff --git a/src/Data/Matlab/NumericArrayReader.cs b/src/Data/Matlab/NumericArrayReader.cs deleted file mode 100644 index 799ce4c2..00000000 --- a/src/Data/Matlab/NumericArrayReader.cs +++ /dev/null @@ -1,257 +0,0 @@ -// -// Math.NET Numerics, part of the Math.NET Project -// http://numerics.mathdotnet.com -// http://github.com/mathnet/mathnet-numerics -// http://mathnetnumerics.codeplex.com -// -// Copyright (c) 2009-2014 Math.NET -// -// 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. -// - -using System; -using System.IO; -using System.Numerics; -using MathNet.Numerics.LinearAlgebra; - -namespace MathNet.Numerics.Data.Matlab -{ - internal static class NumericArrayReader - where TDataType : struct, IEquatable, IFormattable - { - /// - /// Populates a dense matrix. - /// - /// The MATLAB data type. - /// The reader to read from. - /// if set to true if the MATLAB complex flag is set. - /// The number of rows. - /// The number of columns. - /// The length of the stored data. - /// Returns a populated dense matrix. - public static Matrix PopulateDenseMatrix(DataType type, BinaryReader reader, bool isComplex, int rows, int columns, int size) - { - var dataType = typeof (TDataType); - Matrix matrix; - if (type == DataType.Double && dataType == typeof (double)) - { - var count = rows*columns; - var data = new double[count]; - Buffer.BlockCopy(reader.ReadBytes(count*Constants.SizeOfDouble), 0, data, 0, count*Constants.SizeOfDouble); - matrix = (Matrix)(object)new LinearAlgebra.Double.DenseMatrix(rows, columns, data); - } - else if (type == DataType.Single && dataType == typeof (float)) - { - var count = rows*columns; - var data = new float[count]; - Buffer.BlockCopy(reader.ReadBytes(count*Constants.SizeOfFloat), 0, data, 0, count*Constants.SizeOfFloat); - matrix = (Matrix)(object)new LinearAlgebra.Single.DenseMatrix(rows, columns, data); - } - else - { - matrix = Matrix.Build.Dense(rows, columns); - - if (dataType == typeof (double)) - { - if (isComplex) - { - throw new ArgumentException("Invalid TDataType. Matrix is stored as a complex matrix, but a real data type was given."); - } - - PopulateDoubleDenseMatrix((Matrix)(object)matrix, type, reader, rows, columns); - } - else if (dataType == typeof (float)) - { - if (isComplex) - { - throw new ArgumentException("Invalid TDataType. Matrix is stored as a complex matrix, but a real data type was given."); - } - - PopulateSingleDenseMatrix((Matrix)(object)matrix, type, reader, rows, columns); - } - else if (dataType == typeof (Complex)) - { - PopulateComplexDenseMatrix((Matrix)(object)matrix, type, isComplex, reader, rows, columns, size); - } - else if (dataType == typeof (Complex32)) - { - PopulateComplex32DenseMatrix((Matrix)(object)matrix, type, isComplex, reader, rows, columns, size); - } - else - { - throw new NotSupportedException(); - } - } - - return matrix; - } - - /// - /// Populates the double dense matrix. - /// - /// The matrix to populate. - /// The MATLAB data type. - /// The reader to read from. - /// The number of rows. - /// The number of columns. - static void PopulateDoubleDenseMatrix(Matrix matrix, DataType type, BinaryReader reader, int rows, int columns) - { - for (var j = 0; j < columns; j++) - { - for (var i = 0; i < rows; i++) - { - matrix.At(i, j, ReadDoubleValue(type, reader)); - } - } - } - - /// - /// Populates the complex dense matrix. - /// - /// The matrix to populate. - /// The MATLAB data type. - /// if set to true if the MATLAB complex flag is set. - /// The reader to read from. - /// The number of rows. - /// The number of columns. - /// The length of the stored data. - static void PopulateComplexDenseMatrix(Matrix matrix, DataType type, bool isComplex, BinaryReader reader, int rows, int columns, int dataSize) - { - for (var j = 0; j < columns; j++) - { - for (var i = 0; i < rows; i++) - { - matrix.At(i, j, ReadDoubleValue(type, reader)); - } - } - - if (isComplex) - { - var skip = dataSize%8; - - // skip pad - reader.ReadBytes(skip); - - // skip header - type = (DataType)reader.ReadInt32(); - reader.ReadInt32(); - - for (var j = 0; j < columns; j++) - { - for (var i = 0; i < rows; i++) - { - matrix.At(i, j, new Complex(matrix.At(i, j).Real, ReadDoubleValue(type, reader))); - } - } - } - } - - /// - /// Populates the complex32 dense matrix. - /// - /// The matrix to populate. - /// The MATLAB data type. - /// if set to true if the MATLAB complex flag is set. - /// The reader to read from. - /// The number of rows. - /// The number of columns. - /// The length of the stored data. - static void PopulateComplex32DenseMatrix(Matrix matrix, DataType type, bool isComplex, BinaryReader reader, int rows, int columns, int dataSize) - { - for (var j = 0; j < columns; j++) - { - for (var i = 0; i < rows; i++) - { - matrix.At(i, j, (float)ReadDoubleValue(type, reader)); - } - } - - if (isComplex) - { - var skip = dataSize%8; - - // skip pad - reader.ReadBytes(skip); - - // skip header - type = (DataType)reader.ReadInt32(); - reader.ReadInt32(); - - for (var j = 0; j < columns; j++) - { - for (var i = 0; i < rows; i++) - { - matrix.At(i, j, new Complex32(matrix.At(i, j).Real, (float)ReadDoubleValue(type, reader))); - } - } - } - } - - /// - /// Populates the float dense matrix. - /// - /// The matrix to populate. - /// The MATLAB data type. - /// The reader to read from. - /// The number of rows. - /// The number of columns. - static void PopulateSingleDenseMatrix(Matrix matrix, DataType type, BinaryReader reader, int rows, int columns) - { - for (var j = 0; j < columns; j++) - { - for (var i = 0; i < rows; i++) - { - matrix.At(i, j, (float)ReadDoubleValue(type, reader)); - } - } - } - - static double ReadDoubleValue(DataType type, BinaryReader reader) - { - switch (type) - { - case DataType.Double: - return reader.ReadDouble(); - case DataType.Int8: - return reader.ReadSByte(); - case DataType.UInt8: - return reader.ReadByte(); - case DataType.Int16: - return reader.ReadInt16(); - case DataType.UInt16: - return reader.ReadUInt16(); - case DataType.Int32: - return reader.ReadInt32(); - case DataType.UInt32: - return reader.ReadUInt32(); - case DataType.Single: - return reader.ReadSingle(); - case DataType.Int64: - return reader.ReadInt64(); - case DataType.UInt64: - return reader.ReadUInt64(); - default: - throw new NotSupportedException(); - } - } - } -} diff --git a/src/Data/Matlab/Parser.cs b/src/Data/Matlab/Parser.cs index b4533cc5..8bdf9d31 100644 --- a/src/Data/Matlab/Parser.cs +++ b/src/Data/Matlab/Parser.cs @@ -32,8 +32,10 @@ using System; using System.Collections.Generic; using System.IO; using System.IO.Compression; +using System.Numerics; using System.Text; using MathNet.Numerics.LinearAlgebra; +using MathNet.Numerics.LinearAlgebra.Storage; using MathNet.Numerics.Properties; namespace MathNet.Numerics.Data.Matlab @@ -43,11 +45,6 @@ namespace MathNet.Numerics.Data.Matlab /// internal static class Parser { - /// - /// Large Block Size - /// - const int LargeBlockSize = 8; - /// /// Little Endian Indicator /// @@ -59,74 +56,9 @@ namespace MathNet.Numerics.Data.Matlab const int SmallBlockSize = 4; /// - /// Parse a matrix block byte array + /// Large Block Size /// - internal static Matrix ParseMatrix(byte[] data) - where T : struct, IEquatable, IFormattable - { - using (var stream = new MemoryStream(data)) - using (var reader = new BinaryReader(stream)) - { - // skip tag - doesn't tell us anything we don't already know - reader.BaseStream.Seek(8, SeekOrigin.Current); - - var arrayClass = (ArrayClass)reader.ReadByte(); - var flags = reader.ReadByte(); - var isComplex = (flags & (byte)ArrayFlags.Complex) == (byte)ArrayFlags.Complex; - - // skip unneeded bytes - reader.BaseStream.Seek(10, SeekOrigin.Current); - - var numDimensions = reader.ReadInt32()/8; - if (numDimensions > 2) - { - throw new NotSupportedException(Resources.MoreThan2D); - } - - var rows = reader.ReadInt32(); - var columns = reader.ReadInt32(); - - // skip name and unneeded bytes - reader.BaseStream.Seek(2, SeekOrigin.Current); - int size = reader.ReadInt16(); - var smallBlock = true; - if (size == 0) - { - size = reader.ReadInt32(); - smallBlock = false; - } - - reader.BaseStream.Seek(size, SeekOrigin.Current); - AlignData(reader.BaseStream, size, smallBlock); - - var type = (DataType)reader.ReadInt16(); - size = reader.ReadInt16(); - if (size == 0) - { - size = reader.ReadInt32(); - } - - Matrix matrix; - switch (arrayClass) - { - case ArrayClass.Sparse: - matrix = SparseArrayReader.PopulateSparseMatrix(reader, isComplex, rows, columns, size); - break; - case ArrayClass.Function: - case ArrayClass.Character: - case ArrayClass.Object: - case ArrayClass.Structure: - case ArrayClass.Cell: - case ArrayClass.Unknown: - throw new NotSupportedException(); - default: - matrix = NumericArrayReader.PopulateDenseMatrix(type, reader, isComplex, rows, columns, size); - break; - } - - return matrix; - } - } + const int LargeBlockSize = 8; /// /// Extracts all matrix blocks in a format we support from a stream. @@ -137,38 +69,44 @@ namespace MathNet.Numerics.Data.Matlab using (var reader = new BinaryReader(stream)) { + // skip header (116 bytes) + // skip subsystem data offset (8 bytes) + // skip version (2 bytes) reader.BaseStream.Position = 126; + + // endian indicator (2 bytes) if (reader.ReadByte() != LittleEndianIndicator) { throw new NotSupportedException(Resources.BigEndianNotSupported); } - // skip version since it is always 0x0100. + // set position to first data element, right after full file header (128 bytes) reader.BaseStream.Position = 128; var length = stream.Length; - // for each data block add a MATLAB object to the file. + // for each data element add a MATLAB object to the file. while (reader.BaseStream.Position < length) { - var type = (DataType)reader.ReadInt16(); - int size = reader.ReadInt16(); - var smallBlock = true; - if (size == 0) - { - size = reader.ReadInt32(); - smallBlock = false; - } + // small format: size (2 bytes), type (2 bytes), data (4 bytes) + // long format: type (4 bytes), size (4 bytes), data (size, aligned to 8 bytes) + DataType type; + int size; + bool smallBlock; + ReadElementTag(reader, out type, out size, out smallBlock); + + // read element data of the size provided in the element header + // uncompress if compressed byte[] data; if (type == DataType.Compressed) { - data = DecompressBlock(reader.ReadBytes(size), out type); + data = UnpackCompressedBlock(reader.ReadBytes(size), out type); } else { data = new byte[size]; reader.Read(data, 0, size); - AlignData(reader.BaseStream, size, smallBlock); + SkipElementPadding(reader.BaseStream, size, smallBlock); } if (type == DataType.Matrix) @@ -202,12 +140,330 @@ namespace MathNet.Numerics.Data.Matlab } /// - /// Aligns the data. + /// Parse a matrix block byte array + /// + internal static Matrix ParseMatrix(byte[] data) + where T : struct, IEquatable, IFormattable + { + using (var stream = new MemoryStream(data)) + using (var reader = new BinaryReader(stream)) + { + // Array Flags tag (8 bytes) + reader.BaseStream.Seek(8, SeekOrigin.Current); + + // Array Flags data: flags (byte 3), class (byte 4) (8 bytes) + var arrayClass = (ArrayClass)reader.ReadByte(); + var flags = reader.ReadByte(); + var complex = (flags & (byte)ArrayFlags.Complex) == (byte)ArrayFlags.Complex; + reader.BaseStream.Seek(6, SeekOrigin.Current); + + // Dimensions Array tag (8 bytes) + reader.BaseStream.Seek(4, SeekOrigin.Current); + var numDimensions = reader.ReadInt32()/8; + if (numDimensions > 2) + { + throw new NotSupportedException(Resources.MoreThan2D); + } + + // Dimensions Array data: row and column count (8 bytes) + var rows = reader.ReadInt32(); + var columns = reader.ReadInt32(); + + // Array name + DataType type; + int size; + bool smallBlock; + ReadElementTag(reader, out type, out size, out smallBlock); + reader.BaseStream.Seek(size, SeekOrigin.Current); + SkipElementPadding(reader.BaseStream, size, smallBlock); + + Matrix matrix; + switch (arrayClass) + { + case ArrayClass.Sparse: + matrix = PopulateSparseMatrix(reader, complex, rows, columns); + break; + case ArrayClass.Function: + case ArrayClass.Character: + case ArrayClass.Object: + case ArrayClass.Structure: + case ArrayClass.Cell: + case ArrayClass.Unknown: + throw new NotSupportedException(); + default: + matrix = PopulateDenseMatrix(reader, complex, rows, columns); + break; + } + + return matrix; + } + } + + /// + /// Populates a dense matrix. + /// + /// The reader to read from. + /// if set to true if the MATLAB complex flag is set. + /// The number of rows. + /// The number of columns. + /// Returns a populated dense matrix. + static Matrix PopulateDenseMatrix(BinaryReader reader, bool complex, int rows, int columns) + where T : struct, IEquatable, IFormattable + { + var dataType = typeof(T); + var count = rows*columns; + var data = new T[count]; + + DataType type; + int size; + bool smallBlock; + + // read real part array + ReadElementTag(reader, out type, out size, out smallBlock); + + // direct copy if possible + if (type == DataType.Double && dataType == typeof(double) || type == DataType.Single && dataType == typeof(float)) + { + Buffer.BlockCopy(reader.ReadBytes(size), 0, data, 0, size); + } + else if (dataType == typeof(double)) + { + if (complex) + { + throw new ArgumentException("Invalid TDataType. Matrix is stored as a complex matrix, but a real data type was given."); + } + + PopulateDoubleArray(reader, (double[])(object)data, type); + } + else if (dataType == typeof(float)) + { + if (complex) + { + throw new ArgumentException("Invalid TDataType. Matrix is stored as a complex matrix, but a real data type was given."); + } + + PopulateSingleArray(reader, (float[])(object)data, type); + } + else if (dataType == typeof(Complex)) + { + PopulateComplexArray(reader, (Complex[])(object)data, complex, type, ref size, ref smallBlock); + } + else if (dataType == typeof(Complex32)) + { + PopulateComplex32Array(reader, (Complex32[])(object)data, complex, type, ref size, ref smallBlock); + } + else + { + throw new NotSupportedException(); + } + + SkipElementPadding(reader.BaseStream, size, smallBlock); + return Matrix.Build.Dense(rows, columns, data); + } + + /// + /// Populates a sparse matrix. /// - /// The stream. - /// The size of the array. - /// if set to true if reading from a small block. - internal static void AlignData(Stream stream, int size, bool smallBlock) + /// The reader. + /// if set to true if the MATLAB complex flag is set. + /// The number of rows. + /// The number of columns. + /// A populated sparse matrix. + static Matrix PopulateSparseMatrix(BinaryReader reader, bool complex, int rows, int columns) + where T : struct, IEquatable, IFormattable + { + // Create matrix with CSR storage. + var matrix = Matrix.Build.Sparse(columns, rows); + + // MATLAB sparse matrices are actually stored as CSC, so just read the data and then transpose. + var storage = matrix.Storage as SparseCompressedRowMatrixStorage; + + DataType type; + int size; + bool smallBlock; + + // populate the row data array + ReadElementTag(reader, out type, out size, out smallBlock); + var ir = storage.ColumnIndices = new int[size/4]; + for (var i = 0; i < ir.Length; i++) + { + ir[i] = reader.ReadInt32(); + } + + SkipElementPadding(reader.BaseStream, size, smallBlock); + + // populate the column data array + ReadElementTag(reader, out type, out size, out smallBlock); + var jc = storage.RowPointers; + if (jc.Length != size/4) + { + throw new Exception("invalid jcsize"); + } + + for (var j = 0; j < jc.Length; j++) + { + jc[j] = reader.ReadInt32(); + } + + SkipElementPadding(reader.BaseStream, size, smallBlock); + + // populate the values + ReadElementTag(reader, out type, out size, out smallBlock); + var dataType = typeof(T); + var data = storage.Values = new T[jc[columns]]; + + if (dataType == typeof(double)) + { + if (complex) + { + throw new ArgumentException("Invalid TDataType. Matrix is stored as a complex matrix, but a real data type was given."); + } + + PopulateDoubleArray(reader, (double[])(object)data, type); + } + else if (dataType == typeof(float)) + { + if (complex) + { + throw new ArgumentException("Invalid TDataType. Matrix is stored as a complex matrix, but a real data type was given."); + } + + PopulateSingleArray(reader, (float[])(object)data, type); + } + else if (dataType == typeof(Complex)) + { + PopulateComplexArray(reader, (Complex[])(object)data, complex, type, ref size, ref smallBlock); + } + else if (dataType == typeof(Complex32)) + { + PopulateComplex32Array(reader, (Complex32[])(object)data, complex, type, ref size, ref smallBlock); + } + else + { + throw new NotSupportedException(); + } + + SkipElementPadding(reader.BaseStream, size, smallBlock); + return matrix.Transpose(); + } + + /// + /// Populates the double dense matrix. + /// + static void PopulateDoubleArray(BinaryReader reader, double[] data, DataType type) + { + for (int i = 0; i < data.Length; i++) + { + data[i] = ReadDoubleValue(reader, type); + } + } + + /// + /// Populates the float dense matrix. + /// + static void PopulateSingleArray(BinaryReader reader, float[] data, DataType type) + { + for (int i = 0; i < data.Length; i++) + { + data[i] = (float)ReadDoubleValue(reader, type); + } + } + + /// + /// Populates the complex dense matrix. + /// + static void PopulateComplexArray(BinaryReader reader, Complex[] data, bool complex, DataType type, ref int size, ref bool smallBlock) + { + for (int i = 0; i < data.Length; i++) + { + data[i] = ReadDoubleValue(reader, type); + } + + if (complex) + { + SkipElementPadding(reader.BaseStream, size, smallBlock); + ReadElementTag(reader, out type, out size, out smallBlock); + + for (int i = 0; i < data.Length; i++) + { + data[i] = new Complex(data[i].Real, ReadDoubleValue(reader, type)); + } + } + } + + /// + /// Populates the complex32 dense matrix. + /// + static void PopulateComplex32Array(BinaryReader reader, Complex32[] data, bool complex, DataType type, ref int size, ref bool smallBlock) + { + for (int i = 0; i < data.Length; i++) + { + data[i] = (float)ReadDoubleValue(reader, type); + } + + if (complex) + { + SkipElementPadding(reader.BaseStream, size, smallBlock); + ReadElementTag(reader, out type, out size, out smallBlock); + + for (int i = 0; i < data.Length; i++) + { + data[i] = new Complex32(data[i].Real, (float)ReadDoubleValue(reader, type)); + } + } + } + + static double ReadDoubleValue(BinaryReader reader, DataType type) + { + switch (type) + { + case DataType.Double: + return reader.ReadDouble(); + case DataType.Int8: + return reader.ReadSByte(); + case DataType.UInt8: + return reader.ReadByte(); + case DataType.Int16: + return reader.ReadInt16(); + case DataType.UInt16: + return reader.ReadUInt16(); + case DataType.Int32: + return reader.ReadInt32(); + case DataType.UInt32: + return reader.ReadUInt32(); + case DataType.Single: + return reader.ReadSingle(); + case DataType.Int64: + return reader.ReadInt64(); + case DataType.UInt64: + return reader.ReadUInt64(); + default: + throw new NotSupportedException(); + } + } + + static void ReadElementTag(BinaryReader reader, out DataType dataType, out int size, out bool smallBlock) + { + // assume small format + smallBlock = true; + + // small type (2 bytes) + dataType = (DataType)reader.ReadInt16(); + + // small size (2 bytes) + size = reader.ReadInt16(); + + if (size == 0) + { + // long format detected + smallBlock = false; + + // long size (4 bytes) + size = reader.ReadInt32(); + } + } + + static void SkipElementPadding(Stream stream, int size, bool smallBlock) { var blockSize = smallBlock ? SmallBlockSize : LargeBlockSize; var offset = 0; @@ -221,19 +477,22 @@ namespace MathNet.Numerics.Data.Matlab } /// - /// Decompresses the block. + /// Unpacks a compressed block. /// /// The compressed data. /// The type data type contained in the block. /// The decompressed block. - static byte[] DecompressBlock(byte[] compressed, out DataType type) + static byte[] UnpackCompressedBlock(byte[] compressed, out DataType type) { byte[] data; - using (var compressedStream = new MemoryStream(compressed, 2, compressed.Length - 6)) - using (var decompressor = new DeflateStream(compressedStream, CompressionMode.Decompress)) using (var decompressed = new MemoryStream()) { - decompressor.CopyTo(decompressed); + using (var compressedStream = new MemoryStream(compressed, 2, compressed.Length - 6)) + using (var decompressor = new DeflateStream(compressedStream, CompressionMode.Decompress)) + { + decompressor.CopyTo(decompressed); + } + decompressed.Position = 0; var buf = new byte[4]; decompressed.Read(buf, 0, 4); diff --git a/src/Data/Matlab/SparseArrayFormatter.cs b/src/Data/Matlab/SparseArrayFormatter.cs new file mode 100644 index 00000000..b836081c --- /dev/null +++ b/src/Data/Matlab/SparseArrayFormatter.cs @@ -0,0 +1,301 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2014 Math.NET +// +// 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. +// + +using System.IO; +using System.Numerics; +using MathNet.Numerics.LinearAlgebra.Storage; + +namespace MathNet.Numerics.Data.Matlab +{ + internal static class SparseArrayFormatter + { + //public static void FormatSparseMatrix(BinaryWriter writer, Matrix matrix, string name, bool isComplex, int rows, int columns, int size) + // where T : struct, IEquatable, IFormattable + //{ + // var transposed = matrix.Transpose(); + // var storage = (SparseCompressedRowMatrixStorage)transposed.Storage; + + // WriteMatrixTagAndName(writer, ArrayClass.Sparse, isComplex, name, storage); + //} + + + internal static void Write(BinaryWriter writer, LinearAlgebra.Double.SparseMatrix matrix) + { + var nzmax = matrix.NonZerosCount; + + // write ir + writer.Write((int)DataType.Int32); + writer.Write(nzmax*4); + + foreach (var column in matrix.EnumerateColumns()) + { + foreach (var row in column.EnumerateNonZeroIndexed()) + { + writer.Write(row.Item1); + } + } + + // add pad if needed + if (nzmax%2 == 1) + { + writer.Write(0); + } + + // write jc + writer.Write((int)DataType.Int32); + writer.Write((matrix.ColumnCount + 1)*4); + writer.Write(0); + var count = 0; + foreach (var column in matrix.EnumerateColumns()) + { + count += ((SparseVectorStorage)column.Storage).ValueCount; + writer.Write(count); + } + + // add pad if needed + if (matrix.ColumnCount%2 == 0) + { + writer.Write(0); + } + + // write data + writer.Write((int)DataType.Double); + writer.Write(nzmax*8); + + foreach (var column in matrix.EnumerateColumns()) + { + foreach (var row in column.EnumerateNonZeroIndexed()) + { + writer.Write(row.Item2); + } + } + } + + internal static void Write(BinaryWriter writer, LinearAlgebra.Single.SparseMatrix matrix) + { + var nzmax = matrix.NonZerosCount; + + // write ir + writer.Write((int)DataType.Int32); + writer.Write(nzmax*4); + + foreach (var column in matrix.EnumerateColumns()) + { + foreach (var row in column.EnumerateNonZeroIndexed()) + { + writer.Write(row.Item1); + } + } + + // add pad if needed + if (nzmax%2 == 1) + { + writer.Write(0); + } + + // write jc + writer.Write((int)DataType.Int32); + writer.Write((matrix.ColumnCount + 1)*4); + writer.Write(0); + var count = 0; + foreach (var column in matrix.EnumerateColumns()) + { + count += ((SparseVectorStorage)column.Storage).ValueCount; + writer.Write(count); + } + + // add pad if needed + if (matrix.ColumnCount%2 == 0) + { + writer.Write(0); + } + + // write data + writer.Write((int)DataType.Single); + writer.Write(nzmax*4); + + foreach (var column in matrix.EnumerateColumns()) + { + foreach (var row in column.EnumerateNonZeroIndexed()) + { + writer.Write(row.Item2); + } + } + + var pad = nzmax*4%8; + PadData(writer, pad); + } + + internal static void Write(BinaryWriter writer, LinearAlgebra.Complex.SparseMatrix matrix) + { + var nzmax = matrix.NonZerosCount; + + // write ir + writer.Write((int)DataType.Int32); + writer.Write(nzmax*4); + + foreach (var column in matrix.EnumerateColumns()) + { + foreach (var row in column.EnumerateNonZeroIndexed()) + { + writer.Write(row.Item1); + } + } + + // add pad if needed + if (nzmax%2 == 1) + { + writer.Write(0); + } + + // write jc + writer.Write((int)DataType.Int32); + writer.Write((matrix.ColumnCount + 1)*4); + writer.Write(0); + var count = 0; + foreach (var column in matrix.EnumerateColumns()) + { + count += ((SparseVectorStorage)column.Storage).ValueCount; + writer.Write(count); + } + + // add pad if needed + if (matrix.ColumnCount%2 == 0) + { + writer.Write(0); + } + + // write data + writer.Write((int)DataType.Double); + writer.Write(nzmax*8); + + foreach (var column in matrix.EnumerateColumns()) + { + foreach (var row in column.EnumerateNonZeroIndexed()) + { + writer.Write(row.Item2.Real); + } + } + + writer.Write((int)DataType.Double); + writer.Write(nzmax*8); + + foreach (var column in matrix.EnumerateColumns()) + { + foreach (var row in column.EnumerateNonZeroIndexed()) + { + writer.Write(row.Item2.Imaginary); + } + } + } + + internal static void Write(BinaryWriter writer, LinearAlgebra.Complex32.SparseMatrix matrix) + { + var nzmax = matrix.NonZerosCount; + + // write ir + writer.Write((int)DataType.Int32); + writer.Write(nzmax*4); + + foreach (var column in matrix.EnumerateColumns()) + { + foreach (var row in column.EnumerateNonZeroIndexed()) + { + writer.Write(row.Item1); + } + } + + // add pad if needed + if (nzmax%2 == 1) + { + writer.Write(0); + } + + // write jc + writer.Write((int)DataType.Int32); + writer.Write((matrix.ColumnCount + 1)*4); + writer.Write(0); + var count = 0; + foreach (var column in matrix.EnumerateColumns()) + { + count += ((SparseVectorStorage)column.Storage).ValueCount; + writer.Write(count); + } + + // add pad if needed + if (matrix.ColumnCount%2 == 0) + { + writer.Write(0); + } + + // write data + writer.Write((int)DataType.Single); + writer.Write(nzmax*4); + + foreach (var column in matrix.EnumerateColumns()) + { + foreach (var row in column.EnumerateNonZeroIndexed()) + { + writer.Write(row.Item2.Real); + } + } + + var pad = nzmax*4%8; + PadData(writer, pad); + + writer.Write((int)DataType.Single); + writer.Write(nzmax*4); + + foreach (var column in matrix.EnumerateColumns()) + { + foreach (var row in column.EnumerateNonZeroIndexed()) + { + writer.Write(row.Item2.Imaginary); + } + } + + PadData(writer, pad); + } + + /// + /// Pads the data with the given byte. + /// + /// Where to write the pad values. + /// The number of bytes to pad. + /// What value to pad with. + static void PadData(BinaryWriter writer, int bytes, byte pad = (byte)0) + { + for (var i = 0; i < bytes; i++) + { + writer.Write(pad); + } + } + } +} diff --git a/src/Data/Matlab/SparseArrayReader.cs b/src/Data/Matlab/SparseArrayReader.cs deleted file mode 100644 index 33e07067..00000000 --- a/src/Data/Matlab/SparseArrayReader.cs +++ /dev/null @@ -1,251 +0,0 @@ -// -// Math.NET Numerics, part of the Math.NET Project -// http://numerics.mathdotnet.com -// http://github.com/mathnet/mathnet-numerics -// http://mathnetnumerics.codeplex.com -// -// Copyright (c) 2009-2014 Math.NET -// -// 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. -// - -using System; -using System.IO; -using System.Numerics; -using MathNet.Numerics.LinearAlgebra; -using MathNet.Numerics.LinearAlgebra.Storage; - -namespace MathNet.Numerics.Data.Matlab -{ - internal static class SparseArrayReader - where TDataType : struct, IEquatable, IFormattable - { - /// - /// Populates a sparse matrix. - /// - /// The reader. - /// if set to true if the MATLAB complex flag is set. - /// The number of rows. - /// The number of columns. - /// The size of the block. - /// A populated sparse matrix. - public static Matrix PopulateSparseMatrix(BinaryReader reader, bool isComplex, int rows, int columns, int size) - { - // Create matrix with CSR storage. - var matrix = Matrix.Build.Sparse(columns, rows); - - // MATLAB sparse matrices are actually stored as CSC, so just read the data and then transpose. - var storage = matrix.Storage as SparseCompressedRowMatrixStorage; - - // populate the row data array - var ir = storage.ColumnIndices = new int[size/4]; - for (var i = 0; i < ir.Length; i++) - { - ir[i] = reader.ReadInt32(); - } - - Parser.AlignData(reader.BaseStream, size, false); - - // skip data type since it will always be int32 - reader.BaseStream.Seek(4, SeekOrigin.Current); - - // populate the column data array - var jcsize = reader.ReadInt32(); - var jc = storage.RowPointers; - if (jc.Length != jcsize/4) - { - throw new Exception("invalid jcsize"); - } - - for (var j = 0; j < jc.Length; j++) - { - jc[j] = reader.ReadInt32(); - } - - Parser.AlignData(reader.BaseStream, jcsize, false); - - var type = (DataType)reader.ReadInt32(); - var dataSize = reader.ReadInt32(); - - var dataType = typeof (TDataType); - - // Allocate memory for matrix values - var data = storage.Values = new TDataType[jc[columns]]; - - if (dataType == typeof (double)) - { - if (isComplex) - { - throw new ArgumentException("Invalid TDataType. Matrix is stored as a complex matrix, but a real data type was given."); - } - - PopulateDoubleSparseMatrix(type, (double[])(object)data, reader); - } - else if (dataType == typeof (float)) - { - if (isComplex) - { - throw new ArgumentException("Invalid TDataType. Matrix is stored as a complex matrix, but a real data type was given."); - } - - PopulateSingleSparseMatrix(type, (float[])(object)data, reader); - } - else if (dataType == typeof (Complex)) - { - PopulateComplexSparseMatrix(type, isComplex, (Complex[])(object)data, reader, dataSize); - } - else if (dataType == typeof (Complex32)) - { - PopulateComplex32SparseMatrix(type, isComplex, (Complex32[])(object)data, reader, dataSize); - } - else - { - throw new NotSupportedException(); - } - - return matrix.Transpose(); - } - - /// - /// Populates the double sparse matrix. - /// - /// The MATLAB data type. - /// The matrix values array. - /// The reader to read from. - static void PopulateDoubleSparseMatrix(DataType type, double[] data, BinaryReader reader) - { - for (var i = 0; i < data.Length; i++) - { - data[i] = ReadDoubleValue(type, reader); - } - } - - /// - /// Populates the float sparse matrix. - /// - /// The MATLAB data type. - /// The matrix values array. - /// The reader to read from. - static void PopulateSingleSparseMatrix(DataType type, float[] data, BinaryReader reader) - { - for (var i = 0; i < data.Length; i++) - { - data[i] = (float)ReadDoubleValue(type, reader); - } - } - - /// - /// Populates the complex sparse matrix. - /// - /// The MATLAB data type. - /// if set to true if the MATLAB complex flag is set. - /// The matrix values array. - /// The reader to read from. - /// The length of the stored data. - static void PopulateComplexSparseMatrix(DataType type, bool isComplex, Complex[] data, BinaryReader reader, int dataSize) - { - for (var i = 0; i < data.Length; i++) - { - data[i] = ReadDoubleValue(type, reader); - } - - if (isComplex) - { - var skip = dataSize%8; - - // skip pad - reader.ReadBytes(skip); - - // skip header - type = (DataType)reader.ReadInt32(); - reader.ReadInt32(); - - for (var i = 0; i < data.Length; i++) - { - data[i] += new Complex(0.0, ReadDoubleValue(type, reader)); - } - } - } - - /// - /// Populates the complex32 sparse matrix. - /// - /// The MATLAB data type. - /// if set to true if the MATLAB complex flag is set. - /// The matrix values array. - /// The reader to read from. - /// The length of the stored data. - static void PopulateComplex32SparseMatrix(DataType type, bool isComplex, Complex32[] data, BinaryReader reader, int dataSize) - { - for (var i = 0; i < data.Length; i++) - { - data[i] = (float)ReadDoubleValue(type, reader); - } - - if (isComplex) - { - var skip = dataSize%8; - - // skip pad - reader.ReadBytes(skip); - - // skip header - type = (DataType)reader.ReadInt32(); - reader.ReadInt32(); - - for (var i = 0; i < data.Length; i++) - { - data[i] += new Complex32(0.0f, (float)ReadDoubleValue(type, reader)); - } - } - } - - static double ReadDoubleValue(DataType type, BinaryReader reader) - { - switch (type) - { - case DataType.Double: - return reader.ReadDouble(); - case DataType.Int8: - return reader.ReadSByte(); - case DataType.UInt8: - return reader.ReadByte(); - case DataType.Int16: - return reader.ReadInt16(); - case DataType.UInt16: - return reader.ReadUInt16(); - case DataType.Int32: - return reader.ReadInt32(); - case DataType.UInt32: - return reader.ReadUInt32(); - case DataType.Single: - return reader.ReadSingle(); - case DataType.Int64: - return reader.ReadInt64(); - case DataType.UInt64: - return reader.ReadUInt64(); - default: - throw new NotSupportedException(); - } - } - } -}