openkernel/0000750021576700116100000000000011600656073012400 5ustar allauzenengopenkernel/AUTHORS0000640021576700116100000000017511600655010013442 0ustar allauzenengCyril Allauzen Mehryar Mohri Afshin Rostamizadeh openkernel/COPYING0000640021576700116100000000101611600655014013424 0ustar allauzenengLicensed under the Apache License, Version 2.0 (the "License"); you may not use these files 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. openkernel/Makedefs0000640021576700116100000000101111600655023014026 0ustar allauzeneng# Assumes that the PREFIX variable is set by the makefile including # this file. # # Include dependencies definitions and default compiler options include $(PREFIX)/Makedefs.custom # Architecture specific definitions ARCH= $(shell uname -s) ifeq ($(ARCH),Darwin) LIBTYPE= dylib SOFLAGS= -dynamiclib -flat_namespace -undefined suppress else LIBTYPE= so SOFLAGS= -shared endif # Internal variables FARDIR=$(PREFIX)/far KERNELLIBDIR=$(PREFIX)/kernel/lib LIBLINEARDIR=$(PREFIX)/liblinear-1.33 LIBSVMDIR=$(PREFIX)/libsvm-2.82 openkernel/Makedefs.custom0000640021576700116100000000110011600655023015336 0ustar allauzeneng# Defines path to Fst install directories FSTPREFIX= /usr/local FSTINCDIR= $(FSTPREFIX)/include FSTLIBDIR= $(FSTPREFIX)/lib # Flags for linking the Fst Library LIBFST= -lfst LIBFSTMAIN= -fstmain # Flags for linking the ICU Unicode library LIBICU= -licuuc # "-licucore" on Mac OS X with Fink # Flags for linking BLAS and LAPACK LIBLAPACK= -llapack -lblas /usr/lib/libgfortran.so.3 # Flags for linking GLPK LIBGLPK= -lglpk # Compiler options CC= g++ OPT= -O2 -fPIC EXTRA_CCFLAGS= # "-I/sw/include" on Mac OS X with Fink EXTRA_LDFLAGS= # "-Wl,-L/sw/lib" on Mac OS X with Fink openkernel/Makefile0000640021576700116100000000040411600655023014031 0ustar allauzeneng all: ( cd liblinear-1.33 ; $(MAKE) all ) ( cd libsvm-2.82 ; $(MAKE) all ) ( cd kernel ; $(MAKE) all ) clean: rm -f *~ ( cd far ; $(MAKE) clean ) ( cd liblinear-1.33 ; $(MAKE) clean ) ( cd libsvm-2.82 ; $(MAKE) clean ) ( cd kernel ; $(MAKE) clean ) openkernel/README0000640021576700116100000000210611600660661013256 0ustar allauzenengOpenKernel - Beta release PLATFORMS: This version is known to work under Linux and MacOS X using g++. The makefiles are unlikely to work on other platforms unmodified. REQUIREMENT: The version 1.2.7 release of the OpenFst is required and can be downloaded at http://www.openfst.org/FstDownload. The library needs to be compiled with --enable-far. Other requirements are: * The ICU library: http://site.icu-project.org/ Tested with versions 3.8 and 4.0, might work with older versions. Only needed if OpenFst was compiled with --with-icu. * The GLPK (GNU Linear Programming Kit) library: http://www.gnu.org/software/glpk/ Tested with version 4.25, might work with older versions. * An implementation of BLAS and LAPACK. The package also include a modified version of LIBSVM 2.82 and an unmodified version of LIBLINEAR 1.33. BUILDING: To build (rudimentary Makefiles for now): * Edit the Makedefs.custom file at the top level if necessary. * Run at the top level: make all DOCUMENTATION: At http://www.openkernel.org. openkernel/kernel/0000750021576700116100000000000011600656322013655 5ustar allauzenengopenkernel/kernel/lib/0000750021576700116100000000000011600656322014423 5ustar allauzenengopenkernel/kernel/lib/sparse-vector.h0000640021576700116100000002643411600655760017410 0ustar allauzeneng// sparse-vector.h // // 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. // // Author: allauzen@google.com (Cyril Allauzen) // // \file // A simple reference-counted sparse vector representation #ifndef KERNEL_LIB_SPARSE_VECTOR_H__ #define KERNEL_LIB_SPARSE_VECTOR_H__ #include #include #include "kernel/lib/kernel.h" // change to min req namespace kernel { // Forward Declarations template class DenseMatrix; // // SPARSE VECTOR CLASS // // Forward declarations template class SparseVectorData; template class DenseVector; // SparseVector: a simple reference-counted dense vector representation template > class SparseVector { public: typedef T ElementType; typedef I IndexType; typedef D DataType; SparseVector() : data_(new DataType) {} SparseVector(const SparseVector &vec) : data_(vec.data_) { if (data_) data_->IncrRefCount(); } ~SparseVector() { if (data_ && !data_->DecrRefCount()) delete data_; } SparseVector &operator=(const SparseVector &v) { if (this != &v) { if (!data_->DecrRefCount()) delete data_; v.data_->IncrRefCount(); data_ = v.data_; } return *this; } size_t Size() const { return data_->Size(); } size_t MaxIndex() const { if(data_ && data_->Size() > 0) return data_->Pair(data_->Size() - 1).first; else return 0; } void Clear() { MutateCheck(); data_->Clear(); } void Reserve(size_t n) { data_->Reserve(n); } const pair &Pair(size_t i) const { return data_->Pair(i); } T Element(I i) const { return data_->Element(i); } void AddElement(I i, const T &x) { if (x == 0.0F) return; MutateCheck(); data_->AddElement(i, x); } void SetElement(I i, const T &x) { MutateCheck(); data_->SetElement(i, x); } void SetPair(size_t i, const pair &p) { MutateCheck(); data_->SetPair(i, p); } static SparseVector *Read(istream &strm) { return new SparseVector(SparseVectorData::Read(strm)); } bool Write(ostream &strm) const { return data_->Write(strm); } private: SparseVector(DataType *data) : data_(data) {} void MutateCheck() { // Copy on write if (data_->RefCount() > 1) { data_->DecrRefCount(); data_ = new DataType(*data_); } } DataType *data_; }; template class SparseVectorData { public: SparseVectorData() {} explicit SparseVectorData(const SparseVectorData &data) : pairs_(data.pairs_) {} ~SparseVectorData() {} size_t Size() const { return pairs_.size(); } void Reserve(size_t n) { pairs_.reserve(n); } void Clear() { pairs_.clear(); } const pair &Pair(size_t i) const { return pairs_[i]; } T Element(I i) const { typename vector >::const_iterator iter = lower_bound( pairs_.begin(), pairs_.end(), make_pair(i, 0.0F), Compare()); return (iter != pairs_.end()) && ((*iter).first == i) ? (*iter).second : 0.0F; } void AddElement(I i, const T &x) { if (!pairs_.empty() && (pairs_.back().first >= i)) LOG(FATAL) << "SparseVector::AddElement: non-zero elements must" << " be added in order"; pairs_.push_back(make_pair(i, x)); } void SetElement(I i, const T &x) { typename vector >::iterator iter = lower_bound( pairs_.begin(), pairs_.end(), make_pair(i, 0.0F), Compare()); if ((*iter).first == i) { if (x == 0.0F) pairs_.erase(iter); else (*iter).second = x; } else { pairs_.insert(iter, make_pair(i, x)); } } void SetPair(size_t i, const pair &p) { pairs_[i] = p; } int RefCount() const { return ref_count_.count(); } int IncrRefCount() { return ref_count_.Incr(); } int DecrRefCount() { return ref_count_.Decr(); } static SparseVectorData *Read(istream &strm) { SparseVectorData *data = new SparseVectorData; fst::ReadType(strm, &(data->pairs_)); if (!strm) { LOG(ERROR) << "DenseVector::Read: Read failed!"; return 0; } return data; } bool Write(ostream &strm) const { fst::WriteType(strm, pairs_); if (!strm) { LOG(ERROR) << "DenseVector::Write: write failed!"; return false; } return true; } private: struct Compare { bool operator()(const pair &x, const pair &y) { return x.first < y.first; } }; vector > pairs_; RefCounter ref_count_; void operator=(const SparseVectorData &); // disallow }; // // VECTOR OPERATIONS // // Vector equality template inline bool Equal(const V1 &a, const V2 &b) { if (a.Size() != b.Size()) return false; for (size_t i = 0; i < a.Size(); ++i) if (a.Pair(i) != b.Pair(i)) return false; return true; } // Vector approximate equality template inline bool ApproxEqual(const V &a, const V &b, typename V::ElementType delta) { if (a.Size() != b.Size()) return false; for (size_t i = 0; i < a.Size(); ++i) { if (a.Pair(i).first != b.Pair(i).first) return false; if ((a.Pair(i).second > b.Pair(i).second + delta) || (a.Pair(i).second < b.Pair(i).second - delta)) return false; } return true; } // // Multiplication by a scalar // // y = alpha * x template inline void ScalarProduct(const T &a, const V &x, SparseVector *y) { y->Clear(); for (size_t i = 0; i < x.Size(); ++i) y->AddElement(x.Pair(i).first, x.Pair(i).second * a); } // x = alpha * x template inline void ScalarProduct(const T &alpha, SparseVector *x) { for(size_t i = 0; i < x->Size(); ++i) x->SetPair(i, make_pair(x.Pair(i).first, alpha * x.Pair(i).second)); } // // Addition // // z = alpha * x + y template inline void Sum(const T &alpha, const V1 &x, const V2 &y, SparseVector *z) { z->Clear(); for (size_t i = 0, j = 0; i < x.Size(); ++i) { while ((j < y.Size()) && (y.Pair(j).first < x.Pair(i).first)) { z->AddElement(y.Pair(j).first, y.Pair(j).second); ++j; } if (j < y.Size() && (x.Pair(i).first == y.Pair(j).first)) { z->AddElement(x.Pair(i).first, alpha *x.Pair(i).second + y.Pair(j).second); ++j; } else { z->AddElement(x.Pair(i).first, alpha * x.Pair(i).second); } } } // y = alpha * x + y template inline void Sum(const T &alpha, const V &x, SparseVector *y) { SparseVector z; Sum(alpha, x, *y, &z); *y = z; } // z = x + y template inline void Sum(const V1 &x, const V2 &y, SparseVector *z) { Sum(1.0F, x, y, z); } // y = x + y template inline void Sum(const V &x, SparseVector *y) { Sum(1.0F, x, y); } // // Pointwise and dot products // // Point-wise product between two sparse vectors: // z = x .* y template inline void PointwiseProduct(const V1 &x, const V2 &y, SparseVector *z) { z->Clear(); for (size_t i = 0, j = 0; (i < x.Size()) && (j < y.Size()); ++i) { while ((j < y.Size()) && (y.Pair(j).first < x.Pair(i).first)) ++j; if (j == y.Size()) break; if (x.Pair(i).first == y.Pair(j).first) { z->AddElement(x.Pair(i).first, x.Pair(i).second * y.Pair(j).second); ++j; } } } // Point-wise product between one sparse and one dense vector: // z = x .* y template inline void PointwiseProduct(const V &x, const DenseVector &y, SparseVector *z) { z->Clear(); for (size_t i = 0; i < x.Size(); ++i) { if (x.Pair(i).first >= y.Size()) // append zeros to y continue; z->AddElement(x.Pair(i).first, x.Pair(i).second * y.Element(x.Pair(i).first)); } } // Dot product between two sparse vectors: // dot = x^T * y template inline typename V1::ElementType DotProduct(const V1 &x, const V2 &y) { typename V1::ElementType dot = 0.0F; for (size_t i = 0, j = 0; (i < x.Size()) && (j < y.Size()); ++i) { while ((j < y.Size()) && (y.Pair(j).first < x.Pair(i).first)) ++j; if (j == y.Size()) break; if (x.Pair(i).first == y.Pair(j).first) { dot += x.Pair(i).second * y.Pair(j).second; ++j; } } return dot; } // Dot product between one sparse and one dense vector: // dot = x^T * y template inline T DotProduct(const V &x, const DenseVector &y) { T dot = 0.0F; for (size_t i = 0; (i < x.Size()) && (x.Pair(i).first < y.Size()); ++i) dot += x.Pair(i).second * y.Element(x.Pair(i).first); return dot; } // Dot product between two sparse vector: // dot = x^T * y // Uses binary search to find element in larger vector. template inline T DotProductSkew(const SparseVector &big, const SparseVector &small) { T dot = 0.0F; T val; for (size_t i = 0; i < small.Size(); ++i) if ((val = big.Element(small.Pair(i).first)) != 0) dot += val * small.Pair(i).second; return dot; } // Weighted dot product: // dot = x^T * Diag(y) * z template inline T DotProduct(const V1 &x, const DenseVector &y, const V2 &z) { T dot = 0.0F; for (size_t i = 0, j = 0; (i < x.Size()) && (j < z.Size()); ++i) { while ((j < z.Size()) && (z.Pair(j).first < x.Pair(i).first)) ++j; if (j == z.Size()) break; if ((x.Pair(i).first == z.Pair(j).first) && (x.Pair(i).first < y.Size())) { dot += x.Pair(i).second * y.Element(x.Pair(i).first) * z.Pair(j).second; ++j; } } return dot; } // Weighted dot product: // dot = x^T * Diag(y) * z template inline typename V::ElementType DotProduct(const V &x, const V &y, const V &z) { SparseVector t; PointwiseProduct(x, y, &t); return DotProduct(t, z); } template inline T OuterProduct(const SparseVector &x, size_t size_x, const SparseVector &y, size_t size_y, DenseMatrix *mat) { mat->Reset(size_x, size_y, 0); for (size_t i = 0; i < x.Size(); ++i) for (size_t j = 0; j < y.Size(); ++j) mat->SetElement(x.Pair(i).first, y.Pair(j).first, x.Pair(i).second * y.Pair(j).second); } } // namespace kernel #endif // KERNEL_LIB_SPARSE_VECTOR_H__ openkernel/kernel/lib/properties.h0000640021576700116100000000367311600655760017007 0ustar allauzeneng// properties.h // // 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. // // Author: allauzen@cs.nyu.edu (Cyril Allauzen) // // \file // Definitions of property bits. #ifndef KERNEL_LIB_PROPERTIES_H__ #define KERNEL_LIB_PROPERTIES_H__ #include "kernel/lib/compat.h" namespace kernel { const uint64 kFeatured = 0x0000000000000001ULL; const uint64 kNormalized = 0x0000000000010000ULL; const uint64 kUnnormalized = 0x0000000000020000ULL; const uint64 kSymmetric = 0x0000000000040000ULL; const uint64 kAsymmetric = 0x0000000000080000ULL; const uint64 kPosDefinite = 0x0000000000100000ULL; const uint64 kNegDefinite = 0x0000000000200000ULL; const uint64 kNotDefinite = 0x0000000000400000ULL; // // Composite Properties // // PDS: Positive Definite Symmetric const uint64 kPDS = kSymmetric | kPosDefinite; // NDS: Negative Definite Symmetric const uint64 kNDS = kSymmetric | kNegDefinite; const uint64 kKernelProperties = kFeatured | kSymmetric | kAsymmetric | kPosDefinite | kNegDefinite | kNotDefinite | kNormalized | kUnnormalized; // // Properties for Kernel operations // uint64 ExternalProductProperties(uint64 inprops, double a); uint64 ProductProperties(uint64 inprops1, uint64 inprops2); uint64 SumProperties(uint64 inprops1, uint64 inprops2, double a1, double a2); } // namespace kernel #endif // KERNEL_LIB_PROPERTIES_H__ openkernel/kernel/lib/sparse-matrix.h0000640021576700116100000002725311600655760017412 0ustar allauzeneng// sparse-matrix.h // // 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. // // Author: allauzen@google.com (Cyril Allauzen) // // \file // A simple reference-counted sparse matrix representation. #ifndef KERNEL_LIB_SPARSE_MATRIX_H__ #define KERNEL_LIB_SPARSE_MATRIX_H__ #include #include "kernel/lib/dense-matrix.h" #include "kernel/lib/sparse-vector.h" #include "kernel/lib/kernel.h" // change to min req namespace kernel { enum MajorType { ROW_MAJOR = 0, COLUMN_MAJOR = 1 }; #define TRANSPOSE_MAJOR_TYPE(A) \ ((A) == ROW_MAJOR ? COLUMN_MAJOR : ROW_MAJOR) // // SPARSE MATRIX CLASS // // Forward declaration template class SparseMatrix; template class SparseMatrixData; template class SparseMatrixSplice; // SparseMatrixBase: template class SparseMatrixBase { public: typedef T ElementType; typedef I IndexType; typedef SparseMatrixBase TransposeMatrix; typedef SparseMatrixSplice VectorType; friend class SparseMatrixBase; SparseMatrixBase() : data_(new SparseMatrixData) {} SparseMatrixBase(const SparseMatrixBase &base) : data_(base.data_) { data_->IncrRefCount(); } virtual ~SparseMatrixBase() { if (!data_->DecrRefCount()) delete data_; } SparseMatrixBase &operator=(const SparseMatrixBase &m) { if (this != &m) { if (!data_->DecrRefCount()) delete data_; m.data_->IncrRefCount(); data_ = m.data_; } return *this; } void Transpose(TransposeMatrix *m) const { *m = TransposeMatrix(data_); } VectorType Vector(size_t i) const { return SparseMatrixSplice(data_, data_->Index(i), data_->Index(i + 1)); } //VectorType Vector(size_t i) const { return data_->Vector(i); } size_t Size() const { return data_->Size(); } size_t NumNonZeros() const { return data_->NumNonZeros(); } size_t MaxSparseIndex() const { return data_->MaxSparseIndex(); } template void AddVector(const V &vec) { data_->AddVector(vec); } void Clear() { data_->Clear(); } void Reserve(size_t num_nonzeros, size_t num_vectors) { data_->Reserve(num_nonzeros, num_vectors); } protected: SparseMatrixBase(SparseMatrixData *data) : data_(data) { data_->IncrRefCount(); } void MutateCheck() { // Copy on write if (data_->RefCount() > 1) { data_->DecrRefCount(); data_ = new SparseMatrixData(*data_); } } private: SparseMatrixData *data_; }; template class SparseMatrix : public SparseMatrixBase { public: typedef T ElementType; typedef I IndexType; typedef SparseMatrix TransposeMatrix; typedef SparseMatrixSplice RowType; using SparseMatrixBase::Transpose; using SparseMatrixBase::operator=; using SparseMatrixBase::Size; using SparseMatrixBase::MaxSparseIndex; using SparseMatrixBase::Vector; SparseMatrix() {} SparseMatrix(const SparseMatrixBase &base) : SparseMatrixBase(base) {} size_t NumRow() const { return Size(); } size_t NumColumn() const { return MaxSparseIndex() + 1; } RowType Row(size_t i) const { return Vector(i); } template void AddRow(const V &row) { AddVector(row); } }; template class SparseMatrix : public SparseMatrixBase { public: typedef T ElementType; typedef I IndexType; typedef SparseMatrix TransposeMatrix; typedef SparseMatrixSplice ColumnType; using SparseMatrixBase::Transpose; using SparseMatrixBase::operator=; using SparseMatrixBase::Size; using SparseMatrixBase::MaxSparseIndex; using SparseMatrixBase::Vector; SparseMatrix() {}; SparseMatrix(const SparseMatrixBase &base) : SparseMatrixBase(base) {} size_t NumColumn() const { return Size(); } size_t NumRow() const { return MaxSparseIndex() + 1; } ColumnType Column(size_t i) const { return Vector(i); } template void AddColumn(const V &col) { AddVector(col); } }; template class SparseMatrixData { public: typedef SparseMatrixSplice VectorType; SparseMatrixData() : max_sparse_index_(0) { indexes_.push_back(0); } explicit SparseMatrixData(const SparseMatrixData &data) : pairs_(data.pairs_), indexes_(data.indexes_), max_sparse_index_(data.max_sparse_index_) {} ~SparseMatrixData() {} size_t Size() const { return indexes_.size() - 1; } size_t MaxSparseIndex() const { return max_sparse_index_; } size_t NumNonZeros() const { return pairs_.size(); } void Clear() { pairs_.clear(); indexes_.resize(1); } void Reserve(size_t num_nonzeros, size_t num_vectors) { pairs_.reserve(num_nonzeros); indexes_.reserve(num_vectors + 1); } I Index(size_t i) const { return indexes_[i]; } const pair &Pair(size_t i) const { return pairs_[i]; } template void AddVector(const V &v) { for (size_t i = 0; i < v.Size(); ++i) pairs_.push_back(v.Pair(i)); indexes_.push_back(pairs_.size()); if (v.MaxIndex() > max_sparse_index_) max_sparse_index_ = v.MaxIndex(); } // TODO: Reserve() int RefCount() const { return ref_count_.count(); } int IncrRefCount() { return ref_count_.Incr(); } int DecrRefCount() { return ref_count_.Decr(); } private: vector > pairs_; vector indexes_; size_t max_sparse_index_; RefCounter ref_count_; void operator=(const SparseMatrixData &); // disallow }; template class SparseMatrixSplice { public: typedef T ElementType; typedef I IndexType; SparseMatrixSplice(SparseMatrixData *data, size_t begin, size_t end) : data_(data), begin_(begin), end_(end) { data_->IncrRefCount(); } SparseMatrixSplice(const SparseMatrixSplice &splice) : data_(splice.data_), begin_(splice.begin_), end_(splice.end_) { data_->IncrRefCount(); } ~SparseMatrixSplice() { if (!data_->DecrRefCount()) delete data_; } size_t Size() const { return end_ - begin_; } const pair &Pair(size_t i) const { return data_->Pair(begin_ + i); } private: SparseMatrixData *data_; size_t begin_; size_t end_; }; // // MATRIX OPERATIONS // // Equality tests template inline bool Equal(const SparseMatrix &a, const SparseMatrix &b) { if (a.Size() != b.Size()) return false; bool flag = true; for (size_t i = 0; i < a.Size(); ++i) if (!Equal(a.Vector(i), b.Vector(i))) flag = false; return flag; } template inline bool ApproxEqual(const SparseMatrix &a, const SparseMatrix &b, T delta) { if (a.Size() != b.Size()) return false; bool flag = true; for (size_t i = 0; i < a.Size(); ++i) if (!ApproxEqual(a.Vector(i), b.Vector(i), delta)) flag = false; return flag; } // Construct a sparse diagonal matrix, with diagonal equal to // diag_vector; template inline void Diagonal(const DenseVector &diag_vector, SparseMatrixBase *m) { for (size_t i = 0; i < diag_vector.Size(); ++i) { SparseVector temp; temp.AddElement(i, diag_vector.Element(i)); m->AddVector(temp); } } // C = A * B // with A dense by row and B dense by column template inline void Product(const SparseMatrix &a, const SparseMatrix &b, DenseMatrix *c) { c->Reset(a.NumRow(), b.NumColumn(), 0.0F); for (size_t i = 0; i < a.NumRow(); ++i) for (size_t j = 0; j < b.NumColumn(); ++j ) c->SetElement(i, j, DotProduct(a.Row(i), b.Column(j))); } // C = A * B // with A dense by row and B dense by column // return a semidense matrix template inline void Product(const SparseMatrix &a, const SparseMatrix &b, SparseMatrix *c) { c->Clear(); for (size_t i = 0; i < a.NumRow(); ++i) { SparseVector temp; for (size_t j = 0; j < b.NumColumn(); ++j ) { SparseMatrixSplice row(a.Row(i)); SparseMatrixSplice column(b.Column(j)); temp.AddElement(j, DotProduct(row, column)); // temp.AddElement(j, DotProduct, // SparseMatrixSplice >(a.Row(i), b.Column(j))); } c->AddRow(temp); } } // C = A * B // with A dense by column and B dense by row template inline void Product(const SparseMatrix &a, const SparseMatrix &b, DenseMatrix *c) { CHECK(a.NumColumn() == b.NumRow()); c->Reset(a.NumRow(), b.NumColumn(), 0); for (size_t k = 0; k < a.NumColumn() ; ++k) { const SparseMatrixSplice &column = a.Column(k); const SparseMatrixSplice &row = b.Row(k); for (size_t i = 0; i < column.Size(); ++i) for (size_t j = 0; j < row.Size(); ++j) { pair pair_i = column.Pair(i), pair_j = row.Pair(j); c->SetElement(pair_i.first, pair_j.first, c->Element(pair_i.first, pair_j.first) + pair_i.second * pair_j.second); } } } // // MATRIX-VECTOR OPERATIONS // // B = A * diag(v) template inline void DiagonalProduct(const SparseMatrix &A, const DenseVector &v, SparseMatrix *B) { B->Clear(); for (size_t i = 0; i < A.NumRow(); ++i) { SparseVector temp; PointwiseProduct(A.Row(i), v, &temp); B->AddRow(temp); } } // u = B * v // B is dense by column. template inline void Product(const SparseMatrix &a, const DenseVector &v, DenseVector *u) { CHECK(v.Size() == a.NumColumn()); u->Clear(); if (u->Size() < a.NumRow()) u->Resize(a.NumRow(), 0); pair p; for (size_t j = 0; j < a.NumColumn(); ++j) { SparseMatrixSplice col(a.Column(j)); for (size_t i = 0; i < col.Size(); ++i) { p = col.Pair(i); u->SetElement(p.first, u->Element(p.first) + v.Element(j) * p.second); } } } // u = B * v // B is dense by row. template inline void Product(const SparseMatrix &a, const DenseVector &v, DenseVector *u) { if (u->Size() != a.NumRow()) u->Resize(a.NumRow(), 0); SparseVector row; pair p; for (size_t i = 0; i < a.NumRow(); ++i) u->SetElement(i, DotProduct(a.Row(i), v)); } } // namespace kernel #endif // KERNEL_LIB_SPARSE_MATRIX_H__ openkernel/kernel/lib/product-kernel.h0000640021576700116100000001342111600655760017541 0ustar allauzeneng// product-kernel.h // // 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. // // Author: allauzen@cs.nyu.edu (Cyril Allauzen) // // \file // Product kernel representation. #ifndef KERNEL_LIB_PRODUCT_KERNEL_H__ #define KERNEL_LIB_PRODUCT_KERNEL_H__ #include "kernel/lib/kernel.h" namespace kernel { template class ProductKernelImpl : public KernelImpl { public: typedef O Object; ProductKernelImpl(const Kernel &k1, const Kernel &k2) : kernel1_(k1.Copy()), kernel2_(k2.Copy()) { KernelImpl::SetType("product"); Init(); KernelImpl::SetProperties( ProductProperties(kernel1_->Properties(), kernel2_->Properties())); } ~ProductKernelImpl() { delete kernel1_; delete kernel2_; } double Eval(size_t i, size_t j) { return kernel1_->Eval(i, j) * kernel2_->Eval(i,j); } void SetData(DataSet *ds) { kernel1_->SetData(ds); kernel2_->SetData(ds); KernelImpl::SetData(ds); } static ProductKernelImpl *Read(istream &strm, const KernelReadOptions &opts) { ProductKernelImpl *impl = new ProductKernelImpl(); KernelHeader hdr; if (!impl->ReadHeader(strm, opts, kMinFileVersion, &hdr)) return 0; Kernel *kernel1 = down_cast*>( Kernel::Read(strm, KernelReadOptions(opts.source))); if (!kernel1) { LOG(ERROR) << "ProductKernel::Read: read failed: " << opts.source; return 0; } Kernel *kernel2 = down_cast*>( Kernel::Read(strm, KernelReadOptions(opts.source))); if (!kernel2) { LOG(ERROR) << "ProductKernel::Read: read failed: " << opts.source; return 0; } impl->kernel1_ = kernel1; impl->kernel2_ = kernel2; impl->Init(); return impl; } bool Write(ostream &strm, const KernelWriteOptions &opts) const { KernelHeader hdr; KernelImpl::WriteHeader(strm, opts, kFileVersion, &hdr); if (!strm) return false; kernel1_->Write(strm, KernelWriteOptions(opts.source)); kernel2_->Write(strm, KernelWriteOptions(opts.source)); if (!strm) LOG(ERROR) << "ProductKernel::Write: write failed: " << opts.source; return strm; } private: void Init() { if (kernel1_->Data() == kernel2_->Data()) KernelImpl::SetData(kernel1_->Data()); else LOG(FATAL) << "ProductKernel: incompatible datasets"; } ProductKernelImpl() : kernel1_(0), kernel2_(0) { KernelImpl::SetType("product"); } private: Kernel *kernel1_; Kernel *kernel2_; // Identifies current file format version static const int32 kFileVersion = 1; // Identifies minimum compatible file format version static const int32 kMinFileVersion = 1; }; template class ProductKernel :public Kernel { public: typedef O Object; ProductKernel(const Kernel &k1, const Kernel &k2) : impl_(new ProductKernelImpl(k1, k2)) { } ProductKernel(const ProductKernel &kernel) : impl_(kernel.impl_) { impl_->IncrRefCount(); } virtual ~ProductKernel() { if(!impl_->DecrRefCount()) delete impl_;} ProductKernel &operator=(const ProductKernel &kernel) { if (this != &kernel) { if(!impl_->DecrRefCount()) delete impl_; kernel.impl_->IncrRefCount(); impl_ = kernel.impl_; } return *this; } virtual double Eval(size_t i, size_t j) { return impl_->Eval(i,j); } virtual const string& Type() const { return impl_->Type(); } virtual ProductKernel *Copy() const { impl_->IncrRefCount(); return new ProductKernel(impl_); } // Read and Write static ProductKernel *Read(istream &strm, const KernelReadOptions &opts) { ProductKernelImpl *impl = ProductKernelImpl::Read(strm, opts); return impl ? new ProductKernel(impl) : 0; } static ProductKernel *Read(const string &filename) { if (!filename.empty()) { ifstream strm(filename.c_str()); if (!strm) { LOG(ERROR) << "ProductKernel::Read: Can't open file: " << filename; return 0; } return Read(strm, KernelReadOptions(filename)); } else { return Read(std::cin, KernelReadOptions("standard input")); } } virtual bool Write(ostream &strm, const KernelWriteOptions &opts) const { impl_->Write(strm, opts); return strm; } virtual bool Write(const string &filename) const { if (!filename.empty()) { ofstream strm(filename.c_str()); if (!strm) { LOG(ERROR) << "ProductKernel::Write: Can't open file: " << filename; return false; } return Write(strm, KernelWriteOptions(filename)); } else { return Write(std::cout, KernelWriteOptions("standard output")); } } virtual void SetData(DataSet* ds) { impl_->SetData(ds); } virtual const DataSet* Data() const { return impl_->Data(); } virtual DataSet* Data() { return impl_->Data(); } virtual uint64 Properties(uint64 mask) const { return impl_->Properties(mask); } virtual void SetProperties(uint64 props, uint64 mask) { impl_->SetProperties(props, mask); } protected: explicit ProductKernel(ProductKernelImpl *impl) : impl_(impl) {} private: ProductKernelImpl *impl_; }; } // namespace kernel #endif // KERNEL_LIB_PRODUCT_KERNEL_H__ openkernel/kernel/lib/dense-vector.h0000640021576700116100000002523011600655760017202 0ustar allauzeneng// dense-vector.h // // 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. // // Author: allauzen@google.com (Cyril Allauzen) // // \file // A simple reference-counted dense vector representation #ifndef KERNEL_LIB_DENSE_VECTOR_H_ #define KERNEL_LIB_DENSE_VECTOR_H_ #include "kernel/lib/blas.h" #include "kernel/lib/kernel.h" // change to min req namespace kernel { // // DENSE VECTOR CLASS // // Forward declaration template class DenseVectorData; // DenseVector: a simple reference-counted dense vector representation template class DenseVector { public: typedef T ElementType; typedef size_t IndexType; DenseVector() : data_(new DenseVectorData()) {} DenseVector(size_t n, const T& x = T()) : data_(new DenseVectorData()) { data_->Resize(n, x); } DenseVector(const DenseVector &v) : data_(v.data_) { if (data_) data_->IncrRefCount(); } ~DenseVector() { if (!data_->DecrRefCount()) delete data_; } DenseVector &operator=(const DenseVector &v) { if (this != &v) { if (!data_->DecrRefCount()) delete data_; v.data_->IncrRefCount(); data_ = v.data_; } return *this; } size_t Size() const { return data_->Size(); } void Clear() { MutateCheck(); data_->Clear(); } void Resize(size_t n, const T &x = 0.0F) { MutateCheck(); data_->Resize(n, x); } const T &Element(size_t i) const { CHECK(i < data_->Size()); return data_->Element(i); } void SetElement(size_t i, const T &x) { CHECK(i < data_->Size()); MutateCheck(); data_->SetElement(i, x); } void AddElement(size_t i, const T &x) { MutateCheck(); if (i >= Size()) data_->Resize(i + 1, static_cast(0.0F)); data_->SetElement(i, x); } void SetBlasInfo(BlasVectorInfo *info) const { info->dimension = Size(); info->increment = 1; info->array = data_->CArray(); } static DenseVector *Read(istream &strm) { return new DenseVector(DenseVectorData::Read(strm)); } bool Write(ostream &strm) const { return data_->Write(strm); } private: DenseVector(DenseVectorData *data) : data_(data) {} void MutateCheck() { // Copy on write if (data_->RefCount() > 1) { data_->DecrRefCount(); data_ = new DenseVectorData(*data_); } } DenseVectorData *data_; }; // Provides the data representation of a dense vector as a STL vector template class DenseVectorData { public: DenseVectorData() {} explicit DenseVectorData(const DenseVectorData &data) : elements_(data.elements_) {} ~DenseVectorData() {} size_t Size() const { return elements_.size(); } void Reserve(size_t n) { elements_.reserve(n); } void Resize(size_t n, const T &x) { elements_.resize(n, x); } void Clear() { elements_.clear(); } const T &Element(size_t i) const { return elements_[i]; } void SetElement(size_t i, const T &x) { elements_[i] = x; } int RefCount() const { return ref_count_.count(); } int IncrRefCount() { return ref_count_.Incr(); } int DecrRefCount() { return ref_count_.Decr(); } T *CArray() const { return const_cast(&(elements_[0])); } static DenseVectorData *Read(istream &strm) { DenseVectorData *data = new DenseVectorData; fst::ReadType(strm, &(data->elements_)); if (!strm) { LOG(ERROR) << "DenseVector::Read: Read failed!"; return 0; } return data; } bool Write(ostream &strm) const { fst::WriteType(strm, elements_); if (!strm) { LOG(ERROR) << "DenseVector::Write: write failed!"; return false; } return true; } private: vector elements_; RefCounter ref_count_; void operator=(const DenseVectorData &); // disallow }; // // VECTOR OPERATIONS // // Fill vector with integers (i, i + 1, ..., j). template inline void MakeRange(size_t i, size_t j, DenseVector *v) { CHECK(j > i); v->Resize(j - i, 0); for (size_t k = 0; k < j - i; ++k) v->SetElement(k, i + k); } // If j > size, then fill vector with integers // (0, ..., j - size, i, ..., size) template inline void MakeRangeWrap(size_t i, size_t j, size_t size, DenseVector *v) { CHECK(j > i); CHECK(i <= size); v->Clear(); size_t vi = 0; if (j >= size) { CHECK(j - size < size); //check does not wrap around twice for (size_t k = 0; k < j - size; ++k) { v->AddElement(vi, k); ++vi; } } for (size_t k = i; k < min(size, j); ++k) { v->AddElement(vi, k); ++vi; } } // Copy a subset of a vector. template inline void SubsetCopy(size_t i, size_t j, const DenseVector &source, DenseVector *subset) { CHECK(j > i); CHECK(i >= 0); CHECK(j < source.Size()); if (subset->Size() != j - i) subset->Resize(j - i, 0); for (size_t k = i; k < j; ++k) subset->SetElement(k - i, source.Element(k)); } // Copy a subset of a vector. template inline void SubsetCopy(const DenseVector &source, const DenseVector index, DenseVector *subset) { subset->Resize(index.Size(), 0); for (size_t i = 0; i < index.Size(); ++i) subset->SetElement(i, source.Element(index.Element(i))); } // Subtract the mean of v from element of v. template inline void Center(DenseVector *v) { T mean = 0; for (size_t i = 0; i < v->Size(); ++i) mean += v->Element(i); mean = mean / v->Size(); for (size_t i = 0; i < v->Size(); ++i) v->SetElement(i, v->Element(i) - mean); } // Vector equality template inline bool Equal(const DenseVector &a, const DenseVector &b) { if (a.Size() != b.Size()) return false; for (size_t i = 0; i < a.Size(); ++i) if (a.Element(i) != b.Element(i)) return false; return true; } // Vector approximate equality template inline bool ApproxEqual(const DenseVector &a, const DenseVector &b, T delta) { if (a.Size() != b.Size()) return false; for (size_t i = 0; i < a.Size(); ++i) if ((a.Element(i) > b.Element(i) + delta) || (a.Element(i) < b.Element(i) - delta)) { VLOG(1) << "Element " << i << " differs: (" << a.Element(i) << "," << b.Element(i) << ")" << endl; return false; } return true; } // Return largest element in vector, linear time template inline T Max(const DenseVector &v) { T max = 0; for (size_t i = 0; i < v.Size(); i++) if (max < v.Element(i)) max = v.Size(); return max; } // Return the l2 norm template inline T Norm(const DenseVector &v) { T norm = 0; for (size_t i = 0; i < v.Size(); ++i) norm += pow(v.Element(i), 2); return sqrt(norm); } // // Multiplication by a scalar // // y = alpha * x template inline void ScalarProduct(const T &alpha, const DenseVector &x, DenseVector *y) { if(y->Size() != x.Size()) y->Resize(x.Size(), 0); for (size_t i = 0; i < x.Size(); ++i) y->SetElement(i, alpha * x.Element(i)); } // x = alpha * x template inline void ScalarProduct(const T &alpha, DenseVector *x) { ScalarProduct(alpha, *x, x); } // // Addition // // z = alpha * x + y template inline void Sum(const T &alpha, const DenseVector &x, const DenseVector &y, DenseVector *z){ size_t m = x.Size() < y.Size() ? x.Size() : y.Size(); size_t n = x.Size() > y.Size() ? x.Size() : y.Size(); if (z->Size() < n) z->Resize(n, 0.0F); for (size_t i = 0; i < m; ++i) z->SetElement(i, alpha * x.Element(i) + y.Element(i)); for (size_t i = m; i < x.Size(); ++i) z->SetElement(i, alpha * x.Element(i)); for (size_t i = m; i < y.Size(); ++i) z->SetElement(i, y.Element(i)); } // y = alpha * x + y template inline void Sum(const T &alpha, const DenseVector &x, DenseVector *y) { Sum(alpha, x, *y, y); } // z = x + y template inline void Sum(const DenseVector &x, const DenseVector &y, DenseVector *z) { Sum(1.0F, x, y, z); } // y = x + y template inline void Sum(const DenseVector &x, DenseVector *y) { Sum(1.0F, x, *y, y); } // // Dot product // // dot = x^T * y template inline T DotProduct(const DenseVector &x, const DenseVector &y) { T dot = 0.0F; size_t n = x.Size() <= y.Size() ? x.Size() : y.Size(); for (size_t i = 0; i < n; ++i) dot += x.Element(i) * y.Element(i); return dot; } // dot = x^T * Diag(y) * z template inline T DotProduct(const DenseVector &x, const DenseVector &y, const DenseVector &z) { T dot = 0.0F; size_t n = x.Size() <= y.Size() ? (z.Size() <= x.Size() ? z.Size() :x.Size()) : (z.Size() <= y.Size() ? z.Size() : y.Size()); for (size_t i = 0; i < n; ++i) dot += x.Element(i) * y.Element(i) * z.Element(i); return dot; } // z = x .* y, pointwise product of x and y. template inline void PointwiseProduct(const DenseVector &x, const DenseVector &y, DenseVector *z) { CHECK(x.Size() == y.Size()); z->Resize(x.Size(), 0); for (size_t i = 0; i < z->Size(); ++i) z->SetElement(i, x.Element(i) * y.Element(i)); } // Sum of entries template inline T Sum(const DenseVector &v) { T sum = 0; for (size_t i = 0; i < v.Size(); ++i) sum += v.Element(i); return sum; } // Return +/-1 vector w, indicating the signs of entries of v. template inline void SignVector(const DenseVector &v, DenseVector *w) { w->Resize(v.Size(), 0); for (size_t i = 0; i < v.Size(); ++i) { if (v.Element(i) > 0) w->SetElement(i, 1); else w->SetElement(i, -1); } } // Mean of entries. template inline T Mean(const DenseVector &v) { T mean = Sum(v); mean = mean / v.Size(); return mean; } // Standard deviation of entries. template inline T StdDev(const DenseVector &v) { T std; T mean = Mean(v); T ss = 0; for (size_t i = 0; i < v.Size(); ++i) ss += pow(v.Element(i), 2); ss = ss / v.Size(); std = sqrt(ss - pow(mean, 2)); return std; } } // namespace kernel #endif // KERNEL_LIB_DENSE_VECTOR_H_ openkernel/kernel/lib/liblinear.cc0000640021576700116100000001152211600655760016702 0ustar allauzeneng // 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. // // Author: rostami@google.com (Afshin Rostamizadeh) #include #include "kernel/lib/liblinear.h" // This file contains functions that are copied directly from LIBLINEAR, // and modified only to also return the dual optimal variables (i.e. the // double *alpha parameter in solve_linear_c_svc). namespace kernel { #define INF HUGE_VAL static void info(const char *fmt,...) { va_list ap; va_start(ap,fmt); vprintf(fmt,ap); va_end(ap); } static void info_flush() { fflush(stdout); } void solve_linear_c_svc(const problem *prob, double *w, double *alpha, double eps, double Cp, double Cn, int solver_type) { int l = prob->l; int n = prob->n; int i, s, iter = 0; double C, d, G; double *QD = new double[l]; int max_iter = 20000; int *index = new int[l]; //double *alpha = new double[l]; //rostami: allocated outside of function schar *y = new schar[l]; int active_size = l; // PG: projected gradient, for shrinking and stopping double PG; double PGmax_old = INF; double PGmin_old = -INF; double PGmax_new, PGmin_new; // default solver_type: L2LOSS_SVM_DUAL double diag_p = 0.5/Cp, diag_n = 0.5/Cn; double upper_bound_p = INF, upper_bound_n = INF; if(solver_type == L1LOSS_SVM_DUAL) { diag_p = 0; diag_n = 0; upper_bound_p = Cp; upper_bound_n = Cn; } for(i=0; iy[i] > 0) { y[i] = +1; QD[i] = diag_p; } else { y[i] = -1; QD[i] = diag_n; } feature_node *xi = prob->x[i]; while (xi->index != -1) { QD[i] += (xi->value)*(xi->value); xi++; } index[i] = i; } while (iter < max_iter) { PGmax_new = -INF; PGmin_new = INF; for (i=0; ix[i]; while(xi->index!= -1) { G += w[xi->index-1]*(xi->value); xi++; } G = G*yi-1; if(yi == 1) { C = upper_bound_p; G += alpha[i]*diag_p; } else { C = upper_bound_n; G += alpha[i]*diag_n; } PG = 0; if (alpha[i] == 0) { if (G > PGmax_old) { active_size--; swap(index[s], index[active_size]); s--; continue; } else if (G < 0) PG = G; } else if (alpha[i] == C) { if (G < PGmin_old) { active_size--; swap(index[s], index[active_size]); s--; continue; } else if (G > 0) PG = G; } else PG = G; PGmax_new = max(PGmax_new, PG); PGmin_new = min(PGmin_new, PG); if(fabs(PG) > 1.0e-12) { double alpha_old = alpha[i]; alpha[i] = min(max(alpha[i] - G/QD[i], 0.0), C); d = (alpha[i] - alpha_old)*yi; xi = prob->x[i]; while (xi->index != -1) { w[xi->index-1] += d*xi->value; xi++; } } } iter++; if(iter % 10 == 0) { info("."); info_flush(); } if(PGmax_new - PGmin_new <= eps) { if(active_size == l) break; else { active_size = l; info("*"); info_flush(); PGmax_old = INF; PGmin_old = -INF; continue; } } PGmax_old = PGmax_new; PGmin_old = PGmin_new; if (PGmax_old <= 0) PGmax_old = INF; if (PGmin_old >= 0) PGmin_old = -INF; } info("\noptimization finished, #iter = %d\n",iter); if (iter >= max_iter) info("Warning: reaching max number of iterations\n"); // calculate objective value double v = 0; int nSV = 0; for(i=0; i 0) ++nSV; } info("Objective value = %lf\n",v/2); info("nSV = %d\n",nSV); delete [] QD; //delete [] alpha; //rostami deleted outside of function. delete [] y; delete [] index; } } // namespace kernel openkernel/kernel/lib/basic-kernel.h0000640021576700116100000001726511600655760017154 0ustar allauzeneng// basic-kernel.h // // 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. // // Author: allauzen@cs.nyu.edu (Cyril Allauzen) // // \file // Simple and modular cached kernel representation. #ifndef KERNEL_LIB_BASIC_KERNEL_H__ #define KERNEL_LIB_BASIC_KERNEL_H__ #include "kernel/lib/kernel.h" #include "kernel/lib/cache.h" namespace kernel { template class BasicKernelImpl : public KernelImpl { public: typedef F Function; typedef typename Function::Object Object; typedef C Cache; using KernelImpl::Data; using KernelImpl::SetData; BasicKernelImpl(const F *function, C* cache, bool normalize) : normalize_(normalize) { KernelImpl::SetType("basic"); cache_ = cache; function_ = function; uint64 props = normalize ? function->Properties() | kNormalized : function->Properties(); KernelImpl::SetProperties(props); } ~BasicKernelImpl() { VLOG(1) << "~BasicKernelImpl[" << KernelImpl::RefCount() << "]"; delete cache_; delete function_; } double Eval(size_t i, size_t j) { return Eval(i, j, normalize_); } double Eval(size_t i, size_t j, bool normalize) { if (normalize) { if (i == j) return 1.0F; double nu = sqrt(Eval(i, i, false) * Eval(j, j, false)); if (nu == 0.0) LOG(ERROR) << "BasicKernel::Eval(" << i << "," << j << "): cannot normalize kernel when" << "k(x,x) = 0 for some x"; return Eval(i, j, false) / nu; } uint32 type = cache_->CacheType(); if ((type & kCacheValue) && cache_->HasValue(i, j)) return cache_->GetValue(i, j); if ((i == j) && (type & kCacheNorm) && cache_->HasNorm(i)) return cache_->GetNorm(i); double value; DataSet *dataset = KernelImpl::Data(); if ((type & kCacheFeature)) { typename F::Feature *fi, *fj; if (cache_->HasFeature(i)) { fi = cache_->GetFeature(i); } else { fi = function_->GetFeature((*dataset)[i]); cache_->SetFeature(i, fi); } if (cache_->HasFeature(j)) { fj = cache_->GetFeature(j); } else { fj = function_->GetFeature((*dataset)[j]); cache_->SetFeature(j, fj); } value = function_->Product(*fi, *fj); } else { value = (*function_)((*dataset)[i], (*dataset)[j]); } if (type & kCacheValue) cache_->SetValue(i, j, value); else if ((i == j) && (type & kCacheNorm)) cache_->SetNorm(i, value); return value; } static BasicKernelImpl *Read(istream &strm, const KernelReadOptions &opts); bool Write(ostream &strm, const KernelWriteOptions &opts) const; void ClearCache() { if (cache_) cache_->Clear(); } const F* GetFunction() { return function_; } private: BasicKernelImpl() : cache_(0), function_(0), normalize_(false) {} private: C *cache_; const F *function_; bool normalize_; // Identifies current file format version static const int32 kFileVersion = 2; // Identifies minimum compatible file format version static const int32 kMinFileVersion = 2; }; template BasicKernelImpl *BasicKernelImpl::Read( istream &strm, const KernelReadOptions &opts) { BasicKernelImpl *impl = new BasicKernelImpl(); impl->SetType(opts.type); KernelHeader hdr; if (!impl->ReadHeader(strm, opts, kMinFileVersion, &hdr)) return 0; ReadType(strm, &impl->normalize_); impl->function_ = F::Read(strm, opts); impl->cache_ = C::Read(strm, opts); if (!impl->function_ || !impl->cache_) { LOG(ERROR) << "BasicKernel::Read: read failed: " << opts.source; return 0; } bool has_ds; ReadType(strm, &has_ds); if (has_ds) { DataSet *ds = DataSet::Read(strm, opts); impl->SetData(ds); } return impl; } template bool BasicKernelImpl::Write(ostream &strm, const KernelWriteOptions &opts) const { KernelHeader hdr; KernelImpl::WriteHeader(strm, opts, kFileVersion, &hdr); if (!strm) return false; WriteType(strm, normalize_); function_->Write(strm, opts); cache_->Write(strm, opts); const DataSet *ds = this->Data(); bool has_ds = (ds != 0); WriteType(strm, has_ds); if (has_ds) ds->Write(strm, opts); strm.flush(); if (!strm) LOG(ERROR) << "BasicKernel::Write: write failed: " << opts.source; return strm; } template class BasicKernel : public Kernel { public: typedef F Function; typedef C Cache; typedef typename Function::Object Object; BasicKernel(const F* f, C* c, bool n) : impl_(new BasicKernelImpl(f,c,n)) {} BasicKernel(const BasicKernel &kernel) : impl_(kernel.impl_) { impl_->IncrRefCount(); } virtual ~BasicKernel() { if(!impl_->DecrRefCount()) delete impl_;} BasicKernel &operator=(const BasicKernel &kernel) { if (this != &kernel) { if(!impl_->DecrRefCount()) delete impl_; kernel.impl_->IncrRefCount(); impl_ = kernel.impl_; } return *this; } virtual double Eval(size_t i, size_t j) { return impl_->Eval(i,j); } virtual const string& Type() const { return impl_->Type(); } virtual BasicKernel *Copy() const { impl_->IncrRefCount(); return new BasicKernel(impl_); } static BasicKernel *Read(istream &strm, const KernelReadOptions &opts) { BasicKernelImpl* impl = BasicKernelImpl::Read(strm, opts); return impl ? new BasicKernel(impl) : 0; } static BasicKernel *Read(const string &filename) { if (!filename.empty()) { ifstream strm(filename.c_str()); if (!strm) { LOG(ERROR) << "BasicKernel::Read: Can't open file: " << filename; return 0; } return Read(strm, KernelReadOptions(filename)); } else { return Read(std::cin, KernelReadOptions("standard input")); } } virtual bool Write(ostream &strm, const KernelWriteOptions &opts) const { return impl_->Write(strm, opts); } virtual bool Write(const string &filename) const { if (!filename.empty()) { ofstream strm(filename.c_str()); if (!strm) { LOG(ERROR) << "BasicKernel::Write: Can't open file: " << filename; return false; } return Write(strm, KernelWriteOptions(filename)); } else { return Write(std::cout, KernelWriteOptions("standard output")); } } virtual void SetData(DataSet* ds) { impl_->SetData(ds); impl_->ClearCache(); } virtual const DataSet* Data() const { return impl_->Data(); } virtual DataSet* Data() { return impl_->Data(); } virtual uint64 Properties(uint64 mask) const { return impl_->Properties(mask); } virtual void SetProperties(uint64 props, uint64 mask) { impl_->SetProperties(props, mask); } const F* GetFunction() { return impl_->GetFunction(); } protected: BasicKernelImpl *Impl() {return impl_;} explicit BasicKernel(BasicKernelImpl *impl) : impl_(impl) {} private: BasicKernelImpl *impl_; }; } // namespace kernel #endif // KERNEL_LIB_BASIC_KERNEL_H__ openkernel/kernel/lib/dense-matrix.h0000640021576700116100000005476511600655760017223 0ustar allauzeneng// dense-matrix.h // // 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. // // Author: allauzen@google.com (Cyril Allauzen) // // \file // A simple reference-counted dense matrix representation #ifndef KERNEL_LIB_DENSE_MATRIX_H_ #define KERNEL_LIB_DENSE_MATRIX_H_ #include "kernel/lib/blas.h" #include "kernel/lib/sparse-vector.h" #include "kernel/lib/kernel.h" // change to min req namespace kernel { // // DENSE MATRIX CLASS // // Forward declarations template class DenseMatrixData; template class DenseVector; // DenseMatrix: a simple reference-counted dense matrix representation template class DenseMatrix { public: typedef T ElementType; typedef size_t IndexType; typedef DenseMatrix TransposeMatrix; DenseMatrix() : data_(new DenseMatrixData), transpose_(false) {} DenseMatrix(size_t num_row, size_t num_col, const T &x = T()) : data_(new DenseMatrixData(num_row, num_col, x)), transpose_(false) {} DenseMatrix(const DenseMatrix &m) : data_(m.data_), transpose_(m.transpose_) { if (data_) data_->IncrRefCount(); } explicit DenseMatrix(const DenseVector &v) : data_(new DenseMatrixData(v.Size(), 1, T())), transpose_(false) { for (size_t i = 0; i < v.Size(); ++i) data_->SetElement(i, 0, v.Element(i)); } ~DenseMatrix() { if (!data_->DecrRefCount()) delete data_; } DenseMatrix &operator=(const DenseMatrix &m) { if (this != &m) { if (!data_->DecrRefCount()) delete data_; m.data_->IncrRefCount(); data_ = m.data_; transpose_ = m.transpose_; } return *this; } size_t NumRow() const { return transpose_ ? data_->NumColumn() : data_->NumRow(); } size_t NumColumn() const { return transpose_ ? data_->NumRow() : data_->NumColumn(); } void Reset(size_t num_row, size_t num_col, const T &x = T()) { MutateCheck(); data_->Reset(num_row, num_col, x); transpose_ = false; } void Clear() { MutateCheck(); data_->Reset(0, 0, 0); transpose_ = false; } const T& Element(size_t i, size_t j) const { return transpose_ ? data_->Element(j, i) : data_->Element(i, j); } void SetElement(size_t i, size_t j, const T &x) { MutateCheck(); CHECK(i < this->NumRow()); CHECK(j < this->NumColumn()); if (transpose_) data_->SetElement(j, i, x); else data_->SetElement(i, j, x); } void Transpose(TransposeMatrix *m) const { *m = TransposeMatrix(data_, !transpose_); } bool IsTranspose() { return transpose_; } void SetBlasInfo(BlasMatrixInfo *info) const { info->transpose = transpose_ ? 't' : 'n'; info->leading_dimension = data_->NumRow(); info->other_dimension = data_->NumColumn(); info->array = data_->CArray(); } void SetBlasInfo(MutableBlasMatrixInfo *info) { MutateCheck(); info->transpose = transpose_ ? 't' : 'n'; info->leading_dimension = data_->NumRow(); info->other_dimension = data_->NumColumn(); info->array = data_->CArray(); } private: DenseMatrix(DenseMatrixData *data, bool transpose) : data_(data), transpose_(transpose) { data_->IncrRefCount(); } void MutateCheck() { // Copy on write if (data_->RefCount() > 1) { data_->DecrRefCount(); data_ = new DenseMatrixData(*data_); } } DenseMatrixData *data_; bool transpose_; }; // Provides the data representation of a dense matrix as // a STL vector in column-major ordering. template class DenseMatrixData { public: DenseMatrixData() : num_row_(0), num_col_(0) {} DenseMatrixData(size_t num_row, size_t num_col, const T &x) : elements_(num_row * num_col, x), num_row_(num_row), num_col_(num_col) {} explicit DenseMatrixData(const DenseMatrixData &data) : elements_(data.elements_), num_row_(data.num_row_), num_col_(data.num_col_) {} ~DenseMatrixData() {} size_t Size() const { return elements_.size(); } size_t NumRow() const { return num_row_; } size_t NumColumn() const { return num_col_; } void Reset(size_t num_row, size_t num_col, const T &x) { elements_.clear(); elements_.resize(num_row * num_col, x); num_row_ = num_row; num_col_ = num_col; } const T &Element(size_t i, size_t j) const { return elements_[i + num_row_ * j]; } void SetElement(size_t i, size_t j, const T &x) { elements_[i + num_row_ * j] = x; } int RefCount() const { return ref_count_.count(); } int IncrRefCount() { return ref_count_.Incr(); } int DecrRefCount() { return ref_count_.Decr(); } const T *CArray() const { return &(elements_[0]); } T *CArray() { return const_cast(&(elements_[0])); } private: vector elements_; size_t num_row_; size_t num_col_; RefCounter ref_count_; void operator=(const DenseMatrix &); // disallow }; // // MATRIX OPERATIONS // // Copy the submatrix defined by the indices row_index, column_index. // Result matrix: // submatrix(i,j) = mat(row_index(i), column_index(j)) template void SubMatrixCopy(const DenseMatrix &mat, const DenseVector &row_index, const DenseVector &column_index, DenseMatrix *submatrix) { if (submatrix->NumRow() != row_index.Size() || submatrix->NumColumn() != column_index.Size()) submatrix->Reset(row_index.Size(), column_index.Size(), 0); int row, col; for (size_t i = 0; i < row_index.Size(); ++i) { row = row_index.Element(i); CHECK(row < mat.NumRow()); for (size_t j = 0; j < column_index.Size(); ++j) { col = column_index.Element(j); CHECK(col < mat.NumColumn()); submatrix->SetElement(i, j, mat.Element(row, col)); } } } template void SubMatrixCopy(const DenseMatrix &mat, const DenseVector &row_index, DenseMatrix *submatrix) { DenseVector all_col; MakeRange(0, mat.NumColumn(), &all_col); SubMatrixCopy(mat, row_index, all_col, submatrix); } // Return a vector of sparse vectors representing the columns of mat template void MatrixToColumnVectors(const DenseMatrix &mat, vector< SparseVector > *columns) { columns->clear(); for (size_t i = 0; i < mat.NumColumn(); ++i) { SparseVector column; for (size_t j = 0; j push_back(column); } } // Matrix Trace: sum of diagonal entries template T Trace(const DenseMatrix &m) { T sum = 0; size_t min_dim = (m.NumRow() < m.NumColumn()) ? m.NumRow() : m.NumColumn(); for (size_t i = 0; i < min_dim; ++i) sum += m.Element(i, i); return sum; } // Frobenius Norm (squared): sum of squares of all entries template T FrobeniusNorm(const DenseMatrix &m) { T sum = 0; for (size_t i = 0; i < m.NumRow(); ++i) for (size_t j = 0; j < m.NumColumn(); ++j) sum += pow(m.Element(i,j), 2); return sum; } // Frobeinus Product template T FrobeniusProductSym(const DenseMatrix &a, const DenseMatrix &b) { CHECK(a.NumRow() == a.NumColumn()); CHECK(b.NumRow() == b.NumColumn()); CHECK(a.NumRow() == b.NumRow()); T sum = 0; for (size_t i = 0; i < a.NumRow(); ++i) for (size_t j = 0; j < i; ++j) sum += a.Element(i,j) * b.Element(i,j) * 2; for (size_t i = 0; i < a.NumRow(); ++i) sum += a.Element(i,i) * b.Element(i,i); return sum; } // Given a DenseVector, output a diagonal DenseMatrix with diagonal equal // to the input vector. template inline void Diagonal(const DenseVector &diag_vector, DenseMatrix *I) { I->Reset(diag_vector.Size(), diag_vector.Size(), 0); for (size_t i = 0; i < diag_vector.Size(); ++i) I->SetElement(i, i, diag_vector.Element(i)); } // Matrix equality: A == B ? template bool Equal(const DenseMatrix &a, const DenseMatrix &b) { if (a.NumRow() != b.NumRow()) return false; if (a.NumColumn() != b.NumColumn()) return false; for (size_t i = 0; i < a.NumRow(); ++i) for (size_t j = 0; j < a.NumColumn(); ++j) if (a.Element(i, j) != b.Element(i, j)) return false; return true; } // Matrix approximate equality template bool ApproxEqual(const DenseMatrix &a, const DenseMatrix &b, T tol) { if (a.NumRow() != b.NumRow()) return false; if (a.NumColumn() != b.NumColumn()) return false; for (size_t i = 0; i < a.NumRow(); ++i) for (size_t j = 0; j < a.NumColumn(); ++j) if ( a.Element(i, j) > b.Element(i, j) + tol || a.Element(i, j) < b.Element(i, j) - tol) { VLOG(2) << "element (" << i << "," << j << ") differ: " << a.Element(i, j) << ", " << b.Element(i, j) << endl; return false; } return true; } // // Matrix addition // // C = alpha * A + beta * B template void Sum(const T &alpha, const DenseMatrix &a, const T &beta, const DenseMatrix &b, DenseMatrix *c) { CHECK(a.NumRow() == b.NumRow()); CHECK(a.NumColumn() == b.NumColumn()); if((c->NumRow() != a.NumRow()) || (c->NumColumn() != a.NumColumn())) c->Reset(a.NumRow(), a.NumColumn(), 0); for (size_t i = 0; i < a.NumRow(); ++i) for (size_t j = 0; j < a.NumColumn(); ++j) c->SetElement(i, j, alpha * a.Element(i, j) + beta * b.Element(i, j)); } // C = alpha * A + beta * C template inline void Sum(const T &alpha, const DenseMatrix &a, const T &beta, DenseMatrix *c) { Sum(alpha, a, beta, *c, c); } // C = A + B template inline void Sum(const DenseMatrix &a, const DenseMatrix &b, DenseMatrix *c) { return Sum(1.0F, a, 1.0F, b, c); } // C = A + C template inline void Sum(const DenseMatrix &a, DenseMatrix *c) { return Sum(1.0F, a, 1.0F, c); } // Sum of all entries of the matrix template inline T Sum(const DenseMatrix &a) { DenseVector ones_row(a.NumRow(), 1.0); DenseVector ones_col(a.NumColumn(), 1.0); DenseVector temp(a.NumRow(), 0.0); Product(a, ones_col, &temp); return DotProduct(ones_row, temp); } // Sum of the indicated row of the matrix template inline T Sum(const DenseMatrix &a, size_t row) { double sum = 0; for (size_t i = 0; i < a.NumColumn(); ++i) sum += a.Element(row, i); return sum; } // Sum the indicated row and columns of the matrix template inline T Sum(const DenseMatrix &a, const DenseVector &rows, const DenseVector &columns) { double sum = 0; for (size_t i = 0; i < rows.Size(); ++i) for (size_t j = 0; j < columns.Size(); ++j) sum += a.Element(rows.Element(i), columns.Element(j)); return sum; } // // Matrix multiplication: // // C = alpha * A * B + beta * C template void Product(T alpha, const DenseMatrix &a, const DenseMatrix &b, T beta, DenseMatrix *c) { CHECK(a.NumColumn() == b.NumRow()); if ((a.NumRow() != c->NumRow() || b.NumColumn() == c->NumColumn())) c->Reset(a.NumRow(), b.NumColumn(), static_cast(0.0F)); BlasMatrixInfo info_a, info_b; MutableBlasMatrixInfo info_c; a.SetBlasInfo(&info_a); b.SetBlasInfo(&info_b); c->SetBlasInfo(&info_c); CHECK(info_c.transpose == 'n'); int32 m = a.NumRow(); int32 n = b.NumColumn(); int32 k = a.NumColumn(); Gemm()(&(info_a.transpose), &(info_b.transpose), &m, &n, &k, &alpha, info_a.array, &(info_a.leading_dimension), info_b.array, &(info_b.leading_dimension), &beta, info_c.array, &(info_c.leading_dimension)); } // C = A * B template inline void Product(const DenseMatrix &a, const DenseMatrix &b, DenseMatrix *c) { if ((a.NumRow() != c->NumRow() || b.NumColumn() == c->NumColumn())) c->Reset(a.NumRow(), b.NumColumn(), static_cast(0.0F)); Product(static_cast(1.0F), a, b, static_cast(0.0F), c); } // B = A * B // In place product, slower but useful if B is large and A is small. // A must be square. // TODO(rostami): use blas to compute chunks of the product, this could // be the best time/space tradeoff. template inline void Product(const DenseMatrix &a, DenseMatrix *b) { CHECK(a.NumColumn() == b->NumRow()); CHECK(a.NumColumn() == a.NumRow()); DenseVector scratch(b->NumRow(), 0); T temp; for (size_t i = 0; i < b->NumColumn(); ++i) { for (size_t j = 0; j < b->NumRow(); ++j) scratch.SetElement(j, b->Element(j, i)); for (size_t j = 0; j < a.NumRow(); ++j) { temp = 0; for (size_t k = 0; k < a.NumColumn(); ++k) temp += a.Element(j, k) * scratch.Element(k); b->SetElement(j, i, temp); } } } // Point-wise matrix product. // C = A .* B template inline void PointwiseProduct(const DenseMatrix &a, const DenseMatrix &b, DenseMatrix *c) { CHECK(a.NumRow() == b.NumRow()); CHECK(a.NumColumn() == b.NumColumn()); if ((c->NumRow() != a.NumRow()) || (c->NumColumn() != a.NumColumn())) c->Reset(a.NumRow(), a.NumColumn(), 0); for (size_t i = 0; i < a.NumRow(); ++i) for (size_t j = 0; j < a.NumColumn(); ++j) c->SetElement(i, j, a.Element(i, j) * b.Element(i, j)); } // Point-wise (in-place) matrix product. // B = A .* B template inline void PointwiseProduct(const DenseMatrix &a, DenseMatrix *b) { PointwiseProduct(a, *b, b); } // // Multiplication by a scalar // // C = alpha * A template void ScalarProduct(const T &alpha, const DenseMatrix &a, DenseMatrix *c) { if ((c->NumRow() != a.NumRow()) || (c->NumColumn() != a.NumColumn())) c->Reset(a.NumRow(), a.NumColumn(), 0); for (size_t i = 0; i < a.NumRow(); ++i) for (size_t j = 0; j < a.NumColumn(); ++j) c->SetElement(i, j, alpha * a.Element(i,j)); } // C = alpha * C template inline void ScalarProduct(const T &alpha, DenseMatrix *c) { ScalarProduct(alpha, *c, c); } // // MATRIX-VECTOR OPERATIONS // // // Multiplication by a vector // // B = A * diag(v) template void DiagonalProduct(const DenseMatrix &a, const DenseVector &v, DenseMatrix *b) { CHECK(a.NumColumn() == v.Size()); if ((b->NumColumn() != a.NumColumn()) || (b->NumRow() != a.NumRow())) b->Reset(a.NumRow(), a.NumColumn(), 0); //TODO(rostami): See if lapack banded matrix multiply is faster. for (size_t i = 0; i < v.Size(); ++i) for (size_t j = 0; j < a.NumRow(); ++j) b->SetElement(j, i, a.Element(j, i) * v.Element(i)); } // y = alpha * A * x + beta * y template void Product(T alpha, const DenseMatrix &a, const DenseVector &x, T beta, DenseVector *y) { CHECK(a.NumColumn() == x.Size()); y->Resize(a.NumRow(), 0); // force MutateCheck BlasMatrixInfo info_a; BlasVectorInfo info_x, info_y; a.SetBlasInfo(&info_a); x.SetBlasInfo(&info_x); y->SetBlasInfo(&info_y); Gemv()(&(info_a.transpose), &(info_a.leading_dimension), &(info_a.other_dimension), &alpha, info_a.array, &(info_a.leading_dimension), info_x.array, &(info_x.increment), &beta, info_y.array, &(info_y.increment)); } // y = A * x template void Product(const DenseMatrix &a, const DenseVector &x, DenseVector *y) { Product(static_cast(1.0), a, x, static_cast(0.0), y); } // // LAPACK OPERATIONS // // // Solve linear equations, symmetric A // // A * X = B // template void SolveLinearEq(const DenseMatrix &a, const DenseMatrix &b, DenseMatrix *x) { CHECK(a.NumRow() == b.NumRow()); DenseMatrix mut_a(a); DenseMatrix mut_b(b); MutableBlasMatrixInfo info_a, info_b; mut_a.SetBlasInfo(&info_a); mut_b.SetBlasInfo(&info_b); int32 n = mut_a.NumRow(); int32 nrhs = mut_b.NumColumn(); int32 info; char uplo = 'U'; Posv()(&uplo, &n, &nrhs, info_a.array, &(info_a.leading_dimension), info_b.array, &(info_b.leading_dimension), &info); if(info < 0) LOG(ERROR) << "Error in solving linear equations:" << " BAD_FORMAT" << endl; if(info > 0) LOG(ERROR) << "Error in solving linear equations:" << " Non-PDS, increase ridge parameter" << endl; *x = mut_b; } // // Compute eigenvalue/eigenvectors of real symmetric matrix A // // Evec * Eval = A * Evec // template void SymEigen(const DenseMatrix &a, DenseMatrix *evectors, DenseVector *evalues) { CHECK(a.NumRow() == a.NumColumn()); DenseMatrix mut_a(a); MutableBlasMatrixInfo info_a; mut_a.SetBlasInfo(&info_a); int32 n = mut_a.NumRow(); int32 info; char uplo = 'U'; char jobz = 'V'; T *w = new T[n]; T *work_query = new T[1]; int32 lwork = -1; //first call queries for optimal work array size Syev()(&jobz, &uplo, &n, info_a.array, &(info_a.leading_dimension), w, work_query, &lwork, &info); if(info != 0) LOG(ERROR) << "Error (Syev, query phase): code -- " << info << endl; //second call does the actual work T *work = new T[static_cast(work_query[0])]; lwork = static_cast(work_query[0]); Syev()(&jobz, &uplo, &n, info_a.array, &(info_a.leading_dimension), w, work, &lwork, &info); if(info != 0) LOG(ERROR) << "Error (Syev, work phase): code -- " << info << endl; *evectors = mut_a; evalues->Resize(n, 0); for (size_t i = 0; i < n; i++) evalues->SetElement(i, w[i]); delete[] w; delete[] work_query; delete[] work; } // // Compute the low-rank approximation to a real symmetric matrix, given // the sampled indices and data. // // Let K = [W, G12; G12^T, G22] and C = [W, G12], // where K is n-by-n and W is l-by-l. // // Nystrom approximation: K_approx = C * W^+ * C^T // // Approximate SVD: Evec_approx = sqrt(l / n) * C * Evec(W) * Eval(W)^+ // Eval_approx = n / l * Eval(W) // // So that: K_approx = Evec_approx * Eval_approx * Evec_approx^T // // On input, C contains columns sampled from the kernel matrix, and the // vector index indicates which columns were sampled. // On output C contains the approximate eigenvectors, and the DenseVector // evalues holds the approximate eigenvalues. template void ApproxSymEigen(const DenseVector &index, DenseMatrix *C, DenseVector *evalues) { CHECK(C->NumColumn() == index.Size()); DenseMatrix W(C->NumColumn(), C->NumColumn(), 0); DenseMatrix Empty(0, 0, 0); //acts as NULL DenseVector range; MakeRange(0, C->NumColumn(), &range); SubMatrixCopy(*C, index, range, &W); DenseMatrix V; SymEigen(W, &V, evalues); W = Empty; DenseVector evalues_inv(*evalues); //eigenvalues of pseudoinverse //TODO(rostami): decide if this is a good threshold //double thresh = pow(2, -23) * Max(*evalues) * W.NumRow(); double thresh = 0; for (size_t i = 0; i < evalues->Size(); i++) if (evalues->Element(i) > thresh) evalues_inv.SetElement(i, 1 / evalues->Element(i)); else evalues_inv.SetElement(i, 0); ScalarProduct(static_cast(C->NumRow()) / static_cast(C->NumColumn()), evalues); DenseMatrix temp, evalues_inv_mat; Diagonal(evalues_inv, &evalues_inv_mat); Product( static_cast( sqrt(static_cast(C->NumColumn()) / static_cast(C->NumRow()))), V, evalues_inv_mat, static_cast(0), &temp); V = Empty; evalues_inv_mat = Empty; DenseMatrix C_t, temp_t; C->Transpose(&C_t); temp.Transpose(&temp_t); *C = Empty; temp = Empty; //In place product Product(temp_t, &C_t); C_t.Transpose(&temp); *C = temp; } // Compute the inverse of a general matrix, via the LU factorization. template void Inverse(const DenseMatrix &a, DenseMatrix *a_inv) { CHECK(a.NumRow() == a.NumColumn()); DenseMatrix mut_a(a); MutableBlasMatrixInfo info_a; mut_a.SetBlasInfo(&info_a); int32 n = mut_a.NumRow(); int32 info; int32 *ipiv; ipiv = new int32[n]; Getrf()(&n, &n, info_a.array, &(info_a.leading_dimension), ipiv, &info); if(info != 0) LOG(ERROR) << "Error (Inverse, LU): code -- " << info << endl; int32 lwork = -1; T *work; work = new T[1]; Getri()(&n, info_a.array, &(info_a.leading_dimension), ipiv, work, &lwork, &info); if(info != 0) LOG(ERROR) << "Error (Inverse, Blocksize Query): code -- " << info << endl; lwork = static_cast(work[0]); delete[] work; work = new T[lwork]; Getri()(&n, info_a.array, &(info_a.leading_dimension), ipiv, work, &lwork, &info); if(info != 0) LOG(ERROR)<< "Error (Inverse, Inv): code -- " << info << endl; delete[] ipiv; delete[] work; *a_inv = mut_a; } // Compute the inverse of a positive definite symmetric matrix, via // the Cholesky factorization. template void PDSInverse(const DenseMatrix &a, DenseMatrix *a_inv) { CHECK(a.NumRow() == a.NumColumn()); DenseMatrix mut_a(a); MutableBlasMatrixInfo info_a; mut_a.SetBlasInfo(&info_a); int32 n = mut_a.NumRow(); int32 info; char uplo = 'U'; //Cholesky Potrf()(&uplo, &n, info_a.array, &(info_a.leading_dimension), &info); if(info != 0) LOG(ERROR) << "Error (PDSInverse, Chol): code -- " << info << endl; //Inversion Potri()(&uplo, &n, info_a.array, &(info_a.leading_dimension), &info); if(info != 0) LOG(ERROR) << "Error (PDSInverse, Inv): code -- " << info << endl; *a_inv = mut_a; //fill in lower triangular region for (size_t j = 0; j < a_inv->NumColumn() - 1; ++j) for (size_t i = j + 1; i < a_inv->NumRow(); ++i) a_inv->SetElement(i, j, a_inv->Element(j, i)); } } // namespace kernel #endif // KERNEL_LIB_DENSE_MATRIX_H_ openkernel/kernel/lib/liblinear.h0000640021576700116100000000764611600655760016560 0ustar allauzeneng // 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. // // Author: rostami@google.com (Afshin Rostamizadeh) #ifndef KERNEL_LIB_LIBLINEAR_H_ #define KERNEL_LIB_LIBLINEAR_H_ #include "kernel/lib/dense-matrix.h" #include "kernel/lib/semidense-matrix.h" #include "kernel/lib/dense-vector.h" #include "kernel/lib/sparse-vector.h" // TODO(allauzen): remove this include and add relevant forward // declarations instead. #include "linear.h" // A wrapper to call the liblinear train command, using OpenKernel // data structures. namespace kernel { // Modified versions of functions and structures taken from liblinear. void solve_linear_c_svc(const problem *prob, double *w, double *alpha, double eps, double Cp, double Cn, int solver_type); // OpenKernel wrapper function. // Solves the l1-loss svm linear svm problem with no additional offset // feature. // INPUT: // data - a Dense or SemiDense matrix with instances as rows. // labels - a DenseVector of labels for each instance, // reg_param - the SVM regularization parameter "C", which upper bounds // the size of each dual variable. // OUTPUT: // dual_vars - hold the optimal dual variables on output. template void LibLinearTrain(const M &data, const DenseVector &labels, const double ®_param, DenseVector *dual_vars) { CHECK(labels.Size() == data.NumRow()); struct problem *prob = new struct problem; prob->l = data.NumRow(); prob->n = data.NumColumn(); prob->bias = -1; prob->y = new int[labels.Size()]; for (size_t i = 0; i < labels.Size(); ++i) if (labels.Element(i) > 0) prob->y[i] = 1; else prob->y[i] = -1; prob->x = new struct feature_node*[data.NumRow()]; DataToLibLinear(data, prob->x); double *w; w = new double[prob->n]; double *alpha; alpha = new double[prob->l]; //solve_linear_c_svc(prob, w, alpha, eps, Cp, Cn, L1LOSS_SVM_DUAL); solve_linear_c_svc(prob, w, alpha, 0.1, reg_param, reg_param, L1LOSS_SVM_DUAL); dual_vars->Resize(data.NumRow(), 0.0); for (size_t i = 0; i < dual_vars->Size(); ++i) dual_vars->SetElement(i, alpha[i]); for (size_t i = 0; i < prob->l; ++i) delete[] prob->x[i]; delete[] prob->x; delete[] prob->y; delete prob; delete[] w; delete[] alpha; } // Convert an OpenKernel DenseMatrix to an array of liblinear feature_nodes. template inline void DataToLibLinear(const DenseMatrix &data, feature_node **x) { for (size_t i = 0; i < data.NumRow(); ++i) { x[i] = new struct feature_node[data.NumColumn() + 1]; for (size_t j = 0; j < data.NumColumn(); ++j) { x[i][j].index = j + 1; x[i][j].value = data.Element(i,j); } x[i][data.NumColumn()].index = -1; } } // Convert an OpenKernel SparseMatrix (i.e. SemiDenseMatrix) // to an array of liblinear feature_nodes. template inline void DataToLibLinear( const SemiDenseMatrix, DENSE_BY_ROW> &data, feature_node **x) { SparseVector row; for (size_t i = 0; i < data.NumRow(); ++i) { row = data.Row(i); x[i] = new struct feature_node[row.Size() + 1]; for (size_t j = 0; j < row.Size(); ++j) { x[i][j].index = row.Pair(j).first + 1; x[i][j].value = row.Pair(j).second; } x[i][row.Size()].index = -1; } } } // namespace kernel #endif // KERNEL_LIB_LIBLINEAR_H_ openkernel/kernel/lib/object.h0000640021576700116100000001146611600655760016060 0ustar allauzeneng// object.h // // 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. // // Author: allauzen@cs.nyu.edu (Cyril Allauzen) // // \file // Object representation. #ifndef KERNEL_LIB_OBJECT_H__ #define KERNEL_LIB_OBJECT_H__ #include #include #include namespace kernel { struct KernelReadOptions; struct KernelWriteOptions; using fst::Fst; using fst::FstReadOptions; using fst::FstWriteOptions; using fst::SplitToVector; using fst::StrToWeight; using fst::istream; using fst::ostream; using fst::istringstream; ////////////////////////////////////////////////////////////////////// // // Functors handling the manipulations of objects // ////////////////////////////////////////////////////////////////////// template struct ObjectWrite { bool operator()(ostream &strm, const O &object, const KernelWriteOptions &opts = KernelWriteOptions()) { return object.Write(strm); // TODO(allauzen): add opts or source? } }; template struct ObjectRead { O *operator()(istream &strm, const KernelReadOptions &opts = KernelReadOptions()) { return O::Read(strm); // TODO(allauzen): add opts or source? } }; template struct ObjectReadText { O *operator()(const string &txt, const string &src, size_t nline) { LOG(ERROR) << "Unsupported"; // TODO(allauzen): clean this mess! return 0; } }; template struct ObjectType { static const string type; }; template struct ObjectCopy { O *operator()(const O&object) { return new O(object); } }; ////////////////////////////////////////////////////////////////////// // // Specializations for the "Fst" object type // ////////////////////////////////////////////////////////////////////// template <> template struct ObjectWrite< Fst > { bool operator()(ostream &strm, const Fst &fst, const KernelWriteOptions &opts = KernelWriteOptions()) { return fst.Write(strm, FstWriteOptions()); // TODO(allauzen): opts.source } }; template <> template struct ObjectRead< Fst > { Fst* operator()(istream &strm, const KernelReadOptions &opts = KernelReadOptions()) { // TODO(allauzen): opts.source return Fst::Read(strm, FstReadOptions()); } }; template <> template struct ObjectReadText< Fst > { Fst* operator()(const string &text, const string &src, size_t nline) { ifstream is(text.c_str()); Fst *fst = Fst::Read(is, FstReadOptions(text)); return fst; } }; template <> template struct ObjectCopy< Fst > { Fst* operator()(const Fst &fst) { return fst.Copy(); } }; ////////////////////////////////////////////////////////////////////// // // Specializations for the "SparseVector" object type // ////////////////////////////////////////////////////////////////////// // This should be added to io-utilities.h and the functions there // should be be rewritten to re-use to use this template void StringToVector(const string &str, V *vec, const string &src, size_t nline) { char *c_str = new char[str.size() + 1]; str.copy(c_str, str.size()); c_str[str.size()] = 0; vec->Clear(); vector entries; SplitToVector(c_str, " \t", &entries, true); vec->Reserve(entries.size()); vector pair; for (size_t i = 0; i < entries.size(); ++i) { pair.clear(); SplitToVector(entries[i], ":", &pair, false); if ((i == 0) && (pair.size() == 1)) continue; CHECK_EQ(pair.size(), 2); vec->AddElement( StrToWeight(string(pair[0]), src, nline), StrToWeight(string(pair[1]), src, nline)); } // TODO: all this mess need to be fixed delete[] c_str; } // Forward declaration of sparse vector template class SparseVector; template <> template struct ObjectReadText > { SparseVector *operator()(const string &text, const string &src, size_t nline) { SparseVector *sv = new SparseVector; StringToVector(text, sv, src, nline); return sv; } }; } // namespace kernel #endif // KERNEL_LIB_OBJECT_H__ openkernel/kernel/lib/kernel.h0000640021576700116100000002236011600655760016065 0ustar allauzeneng// kernel.h // // 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. // // Author: allauzen@cs.nyu.edu (Cyril Allauzen) // // \file // Abstract base class for kernel representation. #ifndef KERNEL_LIB_KERNEL_H__ #define KERNEL_LIB_KERNEL_H__ #include #include #include "kernel/lib/compat.h" #include "kernel/lib/dataset.h" #include "kernel/lib/register.h" #include "kernel/lib/properties.h" namespace kernel { class KernelHeader; struct KernelReadOptions { string source; const KernelHeader *header; string type; explicit KernelReadOptions(const string& src = "", const KernelHeader *hdr = 0) : source(src), header(hdr), type() {} }; struct KernelWriteOptions { string source; bool write_header; bool write_cache; explicit KernelWriteOptions(const string& src = "", bool cache = true, bool hdr = true): source(src), write_header(hdr), write_cache(cache) {} }; // // Kernel HEADER CLASS // // This is the recommended Kernel file header representation. // class KernelHeader { public: KernelHeader() : version_(0), flags_(0), properties_(0) {} const string &KernelType() const { return kerneltype_; } const string &ObjectType() const { return objecttype_; } int32 Version() const { return version_; } int32 GetFlags() const { return flags_; } uint64 Properties() const { return properties_; } void SetKernelType(const string& type) { kerneltype_ = type; } void SetObjectType(const string& type) { objecttype_ = type; } void SetVersion(int32 version) { version_ = version; } void SetFlags(int32 flags) { flags_ = flags; } void SetProperties(uint64 properties) { properties_ = properties; } bool Read(istream &strm, const string &source); bool Write(ostream &strm, const string &source) const; private: string kerneltype_; // E.g. "rational" string objecttype_; // E.g. "log_fst" int32 version_; // File format version int32 flags_; uint64 properties_; // Kernel property bits }; // // Kernel INTERFACE CLASS DEFINITION // // A generic abstract Kernel with common-denominator methods. class KernelBase { public: virtual ~KernelBase() {} virtual const string &Type() const = 0; // Kernel type name virtual const string &ObjectType() const = 0; // Object type name // Compute the Kernel between objects of index i and j in the // dataset the Kernel is defined on. virtual double Eval(size_t i, size_t j) = 0; // Return stored properties fr mask. virtual uint64 Properties(uint64 mask = kKernelProperties) const = 0; // Set properties with respect to mask. virtual void SetProperties(uint64 props, uint64 mask = kKernelProperties) = 0; // Get a copy of this Kernel virtual KernelBase *Copy() const = 0; // Read a Kernel from an input stream; return NULL on error static KernelBase *Read(istream &strm, const KernelReadOptions &opts) { KernelReadOptions ropts(opts); KernelHeader hdr; if (ropts.header) hdr = *opts.header; else { if (!hdr.Read(strm, opts.source)) return 0; ropts.header = &hdr; } KernelRegister *registr = KernelRegister::GetRegister(); const KernelRegister::Reader reader = registr->GetReader(hdr.ObjectType(), hdr.KernelType()); if (!reader) { LOG(ERROR) << "KernelBase::Read: Unknown kernel type \"" << hdr.KernelType() << "\" (object type = \"" << hdr.ObjectType() << "\"): " << ropts.source; return 0; } return reader(strm, ropts); }; // Read a Kernel from file; return NULL on error // Empty filename reads from standard input static KernelBase *Read(const string &filename) { if (!filename.empty()) { ifstream strm(filename.c_str()); if (!strm) { LOG(ERROR) << "Kernel::Read: Can't open file: " << filename; return 0; } return Read(strm, KernelReadOptions(filename)); } else { return Read(std::cin, KernelReadOptions("standard input")); } } // Write an Kernel to an output stream; return false on error virtual bool Write(ostream &strm, const KernelWriteOptions &opts) const { LOG(ERROR) << "Fst::Kernel: No write method for " << Type() << " kernel type"; return false; } // Write an Kernel to a file; return false on error virtual bool Write(const string &filename) const { LOG(ERROR) << "Fst::Write: No write method for " << Type() << " Fst type: " << (filename.empty() ? "standard output" : filename); return false; } // Return the dataset over which the Kernel is defined virtual const DataSetBase *Data() const = 0; }; // A generic Kernel, templated on the object definition, with // common-denominator methods. template class Kernel : public KernelBase { public: typedef O Object; virtual ~Kernel() {} const string &ObjectType() const { return kernel::ObjectType::type;} virtual Kernel *Copy() const = 0; virtual const DataSet *Data() const = 0; virtual DataSet *Data() = 0; virtual void SetData(DataSet *) = 0; static Kernel *Read(const string &filename) { return down_cast*>(KernelBase::Read(filename)); } static Kernel *Read(istream &strm, const KernelReadOptions &opts) { return down_cast*>(KernelBase::Read(strm, opts)); } }; // // Kernel IMPLEMENTATION BASE CLASS // // This is the recommended Kernel implementation base class. It will // handle reference counts, property bits, type information and // dataset. // template class KernelImpl { public: typedef O Object; KernelImpl() : properties_(0), type_("null"), dataset_(0), ref_count_(1) {} KernelImpl(const KernelImpl &impl) : properties_(impl.properties_), type_(impl.type_), dataset_(impl.dataset_), ref_count_(1) {} ~KernelImpl() { VLOG(1) << "~KernelImpl[" << (dataset_ ? dataset_->RefCount() : -1) << "]"; if (dataset_ && (dataset_->DecrRefCount() == 0)) delete dataset_; } const string& Type() const { return type_; } void SetType(const string &type) { type_ = type; } uint64 Properties() const { return properties_; } uint64 Properties(uint64 mask) const { return properties_ & mask; } void SetProperties(uint64 props) { properties_ = props; } void SetProperties(uint64 props, uint64 mask) { properties_ &= ~mask; properties_ |= props & mask; } const DataSet* Data() const { return dataset_; } DataSet* Data() { return dataset_; } void SetData(DataSet* ds) { if (ds == dataset_) return; if (dataset_ && (dataset_->DecrRefCount() == 0)) delete dataset_; if(ds) ds->IncrRefCount(); dataset_ = ds; } int RefCount() const { return ref_count_; } int IncrRefCount() { return ++ref_count_; } int DecrRefCount() { return --ref_count_; } // Read-in header and symbols, initialize kernel, and return the header. // If opts.header is non-null, skip read-in and use the option value. bool ReadHeader(istream &strm, const KernelReadOptions& opts, int min_version, KernelHeader *hdr) { if (opts.header) *hdr = *opts.header; else if (!hdr->Read(strm, opts.source)) return false; if (hdr->KernelType() != type_) { LOG(ERROR) << "KernelImpl::ReadHeader: Kernel not of type \"" << type_ << "\": " << opts.source; return false; } if (hdr->ObjectType() != ObjectType::type) { LOG(ERROR) << "KernelImpl::ReadHeader: Object not of type \"" << ObjectType::type << "\": " << opts.source; return false; } if (hdr->Version() < min_version) { LOG(ERROR) << "KernelImpl::ReadHeader: Obsolete " << type_ << " kernel version: " << opts.source; return false; } properties_ = hdr->Properties(); return true; } // Write-out header and symbols. // If a opts.header is false, skip writing header. void WriteHeader(ostream &strm, const KernelWriteOptions& opts, int version, KernelHeader *hdr) const { if (opts.write_header) { hdr->SetKernelType(type_); hdr->SetObjectType(ObjectType::type); hdr->SetVersion(version); hdr->SetProperties(properties_); int32 file_flags = 0; hdr->SetFlags(file_flags); hdr->Write(strm, opts.source); } } protected: uint64 properties_; // Property bits private: string type_; // Unique name of Kernel class DataSet *dataset_; // Data set int ref_count_; // Reference count void operator=(const KernelImpl &impl); // disallow }; } // namespace kernel #endif // KERNEL_LIB_KERNEL_H__ openkernel/kernel/lib/io-utilities.h0000640021576700116100000004257211600655760017234 0ustar allauzeneng // 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. // // Author: rostami@google.com (Afshin Rostamizadeh) #ifndef KERNEL_LIB____LIB_IO_UTILITIES_H_ #define KERNEL_LIB____LIB_IO_UTILITIES_H_ #include #include #include "kernel/lib/dense-matrix.h" #include "kernel/lib/semidense-matrix.h" #include "kernel/lib/dense-vector.h" #include "kernel/lib/sparse-vector.h" #include "kernel/lib/kernel.h" namespace kernel { #define ERR_BADFORMAT 1 #define ERR_OUTOFBOUNDS 2 #define ERR_FILENOTFOUND 3 typedef vector > Matrix; typedef vector Vector; void AsciiToVector(const string &filename, Vector *vec); void AsciiToMatrix(const string &filename, Matrix *mat); // Expects each line to be of the form index_1 index_2 val // (r,c) is the dimension of the data. // Outputs a DENSE matrix. int SparseAsciiToMatrix(const string &filename, int r, int c, Matrix *mat); void KernelPrintLibSVM(const Matrix &mat, const Vector &labels); template void AsciiToDenseVector(const string &filename, DenseVector *vec); template void MatrixToDenseMatrix(const Matrix &m, DenseMatrix *dm); template void AsciiToDenseMatrix(const string &filename, DenseMatrix *mat); template int SparseAsciiToDenseMatrix(const string &filename, int r, int c, DenseMatrix *mat); template int LibSVMToDenseMatrix(const string &filename, DenseMatrix *mat, DenseVector *labels); template int LibSVMToSparseMatrix(const string &filename, S *mat, DenseVector *labels); template void KernelPrintLibSVM(const DenseMatrix &mat, const DenseVector &labels, ostream &out_stream); template void KernelPrintLibSVM(const DenseMatrix &mat, const DenseVector &labels, const DenseVector &ids, ostream &out_stream); template //deprecated int LibSVMKernelToDense(const string &filename, int size, DenseMatrix *mat, DenseVector *labels); template int LibSVMKernelToDense(const string &filename, DenseMatrix *mat, DenseVector *labels, DenseVector *point_ids); template int KarToDense(const string &filename, DenseMatrix *mat); template void PrintVector(const V &vec, ostream &out); template void PrintDenseMatrix(const DenseMatrix &mat, ostream &out); template void AsciiToDenseVector(const string &filename, DenseVector *vec) { string line; ifstream file(filename.c_str()); CHECK(file.is_open()); uint64 i = 0; while(1) { getline(file, line); if(file.eof()) break; vec->AddElement(i, strtod(line.c_str(), NULL)); i++; } file.close(); } template void MatrixToDenseMatrix(const Matrix &m, DenseMatrix *dm) { size_t nr = m.size(); size_t nc = m[0].size(); dm->Reset(nr, nc, 0.0); for(size_t i = 0; i < nr; ++i) for(size_t j = 0; j < nc; ++j) dm->SetElement(i, j, m[i][j]); } template void AsciiToDenseMatrix(const string &filename, DenseMatrix *dm) { Matrix m = Matrix(); AsciiToMatrix(filename, &m); MatrixToDenseMatrix(m, dm); } template int SparseAsciiToDenseMatrix(const string &filename, int r, int c, DenseMatrix *mat) { string line; char *index, *oindex; ifstream file(filename.c_str()); CHECK(file.is_open()); double temp[3]; if(mat->NumRow() != r || mat->NumColumn() != c) mat->Reset(r, c, 0); while(1) { getline(file, line); if(file.eof()) break; temp[0] = strtod(line.c_str(), &index); oindex = index; temp[1] = strtod(oindex, &index); oindex = index; temp[2] = strtod(oindex, &index); if(oindex == index) return ERR_BADFORMAT; if(temp[0] >= r || temp[1] >= c) return ERR_OUTOFBOUNDS; mat->SetElement(static_cast(temp[0]), static_cast(temp[1]), temp[2]); } file.close(); return 0; } template int SparseAsciiToDenseMatrix(const string &filename, DenseMatrix *mat) { string line; char *index, *oindex; int max_r = 0; int max_c = 0; ifstream file(filename.c_str()); CHECK(file.is_open()); int row, col; //count maximum number of points and features while(1) { getline(file, line); if(file.eof()) break; row = static_cast(strtod(line.c_str(), &index)); oindex = index; col = static_cast(strtod(oindex, &index)); if(row > max_r) max_r = row; if(col > max_c) max_c = col; } file.close(); //do a second pass, now that size of matrix is determined. return SparseAsciiToDenseMatrix(filename, max_r + 1, max_c + 1, mat); } template int LibSVMGetLabels(const string &filename, DenseVector *labels, uint64 *max_feature) { string line; char *cstr; vector tokens, subtokens; uint64 i, id; *max_feature = 0; ifstream file(filename.c_str()); CHECK(file.is_open()); labels->Resize(1, 0); i = 0; while (1) { getline(file, line); if (file.eof()) break; cstr = new char[line.size()+1]; strcpy(cstr, line.c_str()); tokens.clear(); SplitToVector(cstr, " \t", &tokens, true); if (tokens.size() < 1) { delete[] cstr; LOG(ERROR) << "Error -- line: " << i << " missing label" << endl; return ERR_BADFORMAT; } labels->AddElement(i, strtod(tokens[0], NULL)); ++i; if (tokens.size() < 2) { // no data delete[] cstr; continue; } subtokens.clear(); SplitToVector(tokens[tokens.size() - 1], ":", &subtokens, true); if ((subtokens.size() != 2) || (atol(subtokens[0]) <= 0)) { delete[] cstr; LOG(ERROR) << "Error -- line: " << i << " subtokens.size(): " << subtokens.size() << " id_val: " << atol(subtokens[0]) << endl; return ERR_BADFORMAT; } id = atol(subtokens[0]); if (id == 0) { delete[] cstr; return ERR_BADFORMAT; } if (*max_feature < id) *max_feature = id; delete[] cstr; } file.close(); return 0; } // Read a libSVM formatted data file and populate a dense data matrix // and label vector. // // Each line is of the form: //