/src/rdkit/Code/Numerics/Vector.h

Source
//
//  Copyright (C) 2004-2025 Greg Landrum and other RDKit contributors
//
//   @@ All Rights Reserved @@
//  This file is part of the RDKit.
//  The contents are covered by the terms of the BSD license
//  which is included in the file license.txt, found at the root
//  of the RDKit source tree.
//
#include <RDGeneral/export.h>
#ifndef RD_VECTOR_H
#define RD_VECTOR_H

#include <RDGeneral/Invariant.h>
#include <RDGeneral/utils.h>
#include <cmath>
#include <iomanip>
#include <cstdlib>
#include <cstring>
#include <ctime>
#include <memory>
#include <boost/random.hpp>

static constexpr double zero_tolerance = 1.e-16;

namespace RDNumeric {

//! A class to represent vectors of numbers.
template <class TYPE>
class Vector {
 public:
  using DATA_SPTR = std::shared_ptr<TYPE[]>;

  //! Initialize with only a size.
  constexpr explicit Vector(unsigned int N) {
    d_size = N;
    TYPE *data = new TYPE[N];
    memset(static_cast<void *>(data), 0, d_size * sizeof(TYPE));
    d_data.reset(data);
  }

  //! Initialize with a size and default value.
  constexpr Vector(unsigned int N, TYPE val) {
    d_size = N;
    d_data.reset(new TYPE[N]);
    for (auto i = 0u; i < N; i++) {
      d_data[i] = val;
    }
  }

  //! Initialize from a smart pointer.
  /*!
    <b>NOTE:</b> the data is not copied in this case
  */
  constexpr Vector(unsigned int N, DATA_SPTR data) {
    d_size = N;
    d_data = data;
  }

  //! copy constructor
  /*! We make a copy of the other vector's data.
   */
  constexpr Vector(const Vector &other) {
    d_size = other.size();
    const TYPE *otherData = other.getData();
    d_data.reset(new TYPE[d_size]);

    memcpy(static_cast<void *>(d_data.get()),
           static_cast<const void *>(otherData), d_size * sizeof(TYPE));
  }

  constexpr ~Vector() = default;

  //! return the size (dimension) of the vector
  constexpr unsigned int size() const { return d_size; }

  //! returns the value at a particular index
  TYPE getVal(unsigned int i) const {
    PRECONDITION(i < d_size, "bad index");
    return d_data[i];
  }

  //! sets the index at a particular value
  void setVal(unsigned int i, TYPE val) {
    PRECONDITION(i < d_size, "bad index");
    d_data[i] = val;
  }

  TYPE operator[](unsigned int i) const {
    PRECONDITION(i < d_size, "bad index");
    return d_data[i];
  }

  TYPE &operator[](unsigned int i) {
    PRECONDITION(i < d_size, "bad index");
    return d_data[i];
  }

  //! returns a pointer to our data array
  constexpr TYPE *getData() { return d_data.get(); }

  //! returns a const pointer to our data array
  constexpr const TYPE *getData() const {
    // return dp_data;
    return d_data.get();
  }

  //! Copy operator.
  /*! We make a copy of the other Vector's data.
   */

  Vector<TYPE> &assign(const Vector<TYPE> &other) {
    PRECONDITION(d_size == other.size(), "Size mismatch in vector copying");
    const TYPE *otherData = other.getData();
    memcpy(static_cast<void *>(d_data.get()),
           static_cast<const void *>(otherData), d_size * sizeof(TYPE));
    return *this;
  }

  //! elementwise addition, vectors must be the same size.
  Vector<TYPE> &operator+=(const Vector<TYPE> &other) {
    PRECONDITION(d_size == other.size(), "Size mismatch in vector addition");
    const TYPE *otherData = other.getData();
    TYPE *data = d_data.get();
    unsigned int i;
    for (i = 0; i < d_size; i++) {
      data[i] += otherData[i];
    }
    return *this;
  }

  //! elementwise subtraction, vectors must be the same size.
  Vector<TYPE> &operator-=(const Vector<TYPE> &other) {
    PRECONDITION(d_size == other.size(), "Size mismatch in vector subtraction");
    const TYPE *otherData = other.getData();
    TYPE *data = d_data.get();
    unsigned int i;
    for (i = 0; i < d_size; i++) {
      data[i] -= otherData[i];
    }
    return *this;
  }

  //! multiplication by a scalar
  constexpr Vector<TYPE> &operator*=(TYPE scale) {
    unsigned int i;
    for (i = 0; i < d_size; i++) {
      d_data[i] *= scale;
    }
    return *this;
  }

  //! division by a scalar
  constexpr Vector<TYPE> &operator/=(TYPE scale) {
    unsigned int i;
    for (i = 0; i < d_size; i++) {
      d_data[i] /= scale;
    }
    return *this;
  }

  //! L2 norm squared
  constexpr inline TYPE normL2Sq() const {
    TYPE res = (TYPE)0.0;
    unsigned int i;
    TYPE *data = d_data.get();
    for (i = 0; i < d_size; i++) {
      res += data[i] * data[i];
    }
    return res;
  }

  //! L2 norm
  constexpr TYPE normL2() const { return sqrt(this->normL2Sq()); }

  //! L1 norm
  constexpr TYPE normL1() const {
    TYPE res = (TYPE)0.0;
    unsigned int i;
    TYPE *data = d_data.get();
    for (i = 0; i < d_size; i++) {
      res += fabs(data[i]);
    }
    return res;
  }

  //! L-infinity norm
  constexpr TYPE normLinfinity() const {
    TYPE res = (TYPE)(-1.0);
    unsigned int i;
    TYPE *data = d_data.get();
    for (i = 0; i < d_size; i++) {
      if (fabs(data[i]) > res) {
        res = fabs(data[i]);
      }
    }
    return res;
  }

  //! \brief Gets the ID of the entry that has the largest absolute value
  //! i.e. the entry being used for the L-infinity norm
  constexpr unsigned int largestAbsValId() const {
    TYPE res = (TYPE)(-1.0);
    unsigned int i, id = d_size;
    TYPE *data = d_data.get();
    for (i = 0; i < d_size; i++) {
      if (fabs(data[i]) > res) {
        res = fabs(data[i]);
        id = i;
      }
    }
    return id;
  }

  //! \brief Gets the ID of the entry that has the largest value
  constexpr unsigned int largestValId() const {
    TYPE res = (TYPE)(-1.e8);
    unsigned int i, id = d_size;
    TYPE *data = d_data.get();
    for (i = 0; i < d_size; i++) {
      if (data[i] > res) {
        res = data[i];
        id = i;
      }
    }
    return id;
  }

  //! \brief Gets the ID of the entry that has the smallest value
  constexpr unsigned int smallestValId() const {
    TYPE res = (TYPE)(1.e8);
    unsigned int i, id = d_size;
    TYPE *data = d_data.get();
    for (i = 0; i < d_size; i++) {
      if (data[i] < res) {
        res = data[i];
        id = i;
      }
    }
    return id;
  }

  //! returns the dot product between two Vectors
  TYPE dotProduct(const Vector<TYPE> other) const {
    PRECONDITION(d_size == other.size(),
                 "Size mismatch in vector doct product");
    const TYPE *oData = other.getData();
    unsigned int i;
    TYPE res = (TYPE)(0.0);
    TYPE *data = d_data.get();
    for (i = 0; i < d_size; i++) {
      res += (data[i] * oData[i]);
    }
    return res;
  }

  //! Normalize the vector using the L2 norm
  constexpr void normalize() {
    TYPE val = this->normL2();
    if (val < zero_tolerance) {
      throw std::runtime_error("Cannot normalize a zero length vector");
    }
    (*this) /= val;
  }

  //! Set to a random unit vector
  void setToRandom(unsigned int seed = 0) {
    // we want to get our own RNG here instead of using the global
    // one.  This is related to Issue285.
    RDKit::rng_type generator(42u);
    RDKit::uniform_double dist(0, 1.0);
    RDKit::double_source_type randSource(generator, dist);
    if (seed > 0) {
      generator.seed(seed);
    } else {
      // we can't initialize using only clock(), because it's possible
      // that we'll get here fast enough that clock() will return 0
      // and generator.seed(0) is an error:
      generator.seed(clock() + 1);
    }

    unsigned int i;
    TYPE *data = d_data.get();
    for (i = 0; i < d_size; i++) {
      data[i] = randSource();
    }
    this->normalize();
  }

 private:
  unsigned int d_size;  //!< our length
  DATA_SPTR d_data;
  Vector<TYPE> &operator=(const Vector<TYPE> &other);
};

using DoubleVector = Vector<double>;

//! returns the algebraic tanimoto similarity [defn' from JCIM 46:587-96 (2006)]
template <typename T>
constexpr double TanimotoSimilarity(const Vector<T> &v1, const Vector<T> &v2) {
  double numer = v1.dotProduct(v2);
  if (numer == 0.0) {
    return 0.0;
  }
  double denom = v1.normL2Sq() + v2.normL2Sq() - numer;
  if (denom == 0.0) {
    return 0.0;
  }
  return numer / denom;
}
}  // end of namespace RDNumeric

//! ostream operator for Vectors
template <typename TYPE>
constexpr std::ostream &operator<<(std::ostream &target,
                                   const RDNumeric::Vector<TYPE> &vec) {
  unsigned int siz = vec.size();
  target << "Size: " << siz << " [";
  unsigned int i;
  for (i = 0; i < siz; i++) {
    target << std::setw(7) << std::setprecision(3) << vec.getVal(i) << ", ";
  }
  target << "]\n";
  return target;
}

#endif

Coverage Report

Created: 2026-11-12 07:33

Line	Count	Source
1		//
2		// Copyright (C) 2004-2025 Greg Landrum and other RDKit contributors
3		//
4		// @@ All Rights Reserved @@
5		// This file is part of the RDKit.
6		// The contents are covered by the terms of the BSD license
7		// which is included in the file license.txt, found at the root
8		// of the RDKit source tree.
9		//
10		#include <RDGeneral/export.h>
11		#ifndef RD_VECTOR_H
12		#define RD_VECTOR_H
13
14		#include <RDGeneral/Invariant.h>
15		#include <RDGeneral/utils.h>
16		#include <cmath>
17		#include <iomanip>
18		#include <cstdlib>
19		#include <cstring>
20		#include <ctime>
21		#include <memory>
22		#include <boost/random.hpp>
23
24		static constexpr double zero_tolerance = 1.e-16;
25
26		namespace RDNumeric {
27
28		//! A class to represent vectors of numbers.
29		template <class TYPE>
30		class Vector {
31		public:
32		using DATA_SPTR = std::shared_ptr<TYPE[]>;
33
34		//! Initialize with only a size.
35		constexpr explicit Vector(unsigned int N) {
36		d_size = N;
37		TYPE *data = new TYPE[N];
38		memset(static_cast<void >(data), 0, d_size sizeof(TYPE));
39		d_data.reset(data);
40		}
41
42		//! Initialize with a size and default value.
43	0	constexpr Vector(unsigned int N, TYPE val) {
44	0	d_size = N;
45	0	d_data.reset(new TYPE[N]);
46	0	for (auto i = 0u; i < N; i++) {
47	0	d_data[i] = val;
48	0	}
49	0	}
50
51		//! Initialize from a smart pointer.
52		/*!
53		<b>NOTE:</b> the data is not copied in this case
54		*/
55		constexpr Vector(unsigned int N, DATA_SPTR data) {
56		d_size = N;
57		d_data = data;
58		}
59
60		//! copy constructor
61		/*! We make a copy of the other vector's data.
62		*/
63	0	constexpr Vector(const Vector &other) {
64	0	d_size = other.size();
65	0	const TYPE *otherData = other.getData();
66	0	d_data.reset(new TYPE[d_size]);
67
68	0	memcpy(static_cast<void *>(d_data.get()),
69	0	static_cast<const void >(otherData), d_size sizeof(TYPE));
70	0	}
71
72	0	constexpr ~Vector() = default;
73
74		//! return the size (dimension) of the vector
75	0	constexpr unsigned int size() const { return d_size; }
76
77		//! returns the value at a particular index
78	0	TYPE getVal(unsigned int i) const {
79	0	PRECONDITION(i < d_size, "bad index");
80	0	return d_data[i];
81	0	}
82
83		//! sets the index at a particular value
84		void setVal(unsigned int i, TYPE val) {
85		PRECONDITION(i < d_size, "bad index");
86		d_data[i] = val;
87		}
88
89		TYPE operator[](unsigned int i) const {
90		PRECONDITION(i < d_size, "bad index");
91		return d_data[i];
92		}
93
94	0	TYPE &operator[](unsigned int i) {
95	0	PRECONDITION(i < d_size, "bad index");
96	0	return d_data[i];
97	0	}
98
99		//! returns a pointer to our data array
100	0	constexpr TYPE *getData() { return d_data.get(); }
101
102		//! returns a const pointer to our data array
103	0	constexpr const TYPE *getData() const {
104		// return dp_data;
105	0	return d_data.get();
106	0	}
107
108		//! Copy operator.
109		/*! We make a copy of the other Vector's data.
110		*/
111
112		Vector<TYPE> &assign(const Vector<TYPE> &other) {
113		PRECONDITION(d_size == other.size(), "Size mismatch in vector copying");
114		const TYPE *otherData = other.getData();
115		memcpy(static_cast<void *>(d_data.get()),
116		static_cast<const void >(otherData), d_size sizeof(TYPE));
117		return *this;
118		}
119
120		//! elementwise addition, vectors must be the same size.
121	0	Vector<TYPE> &operator+=(const Vector<TYPE> &other) {
122	0	PRECONDITION(d_size == other.size(), "Size mismatch in vector addition");
123	0	const TYPE *otherData = other.getData();
124	0	TYPE *data = d_data.get();
125	0	unsigned int i;
126	0	for (i = 0; i < d_size; i++) {
127	0	data[i] += otherData[i];
128	0	}
129	0	return *this;
130	0	}
131
132		//! elementwise subtraction, vectors must be the same size.
133	0	Vector<TYPE> &operator-=(const Vector<TYPE> &other) {
134	0	PRECONDITION(d_size == other.size(), "Size mismatch in vector subtraction");
135	0	const TYPE *otherData = other.getData();
136	0	TYPE *data = d_data.get();
137	0	unsigned int i;
138	0	for (i = 0; i < d_size; i++) {
139	0	data[i] -= otherData[i];
140	0	}
141	0	return *this;
142	0	}
143
144		//! multiplication by a scalar
145	0	constexpr Vector<TYPE> &operator*=(TYPE scale) {
146	0	unsigned int i;
147	0	for (i = 0; i < d_size; i++) {
148	0	d_data[i] *= scale;
149	0	}
150	0	return *this;
151	0	}
152
153		//! division by a scalar
154	0	constexpr Vector<TYPE> &operator/=(TYPE scale) {
155	0	unsigned int i;
156	0	for (i = 0; i < d_size; i++) {
157	0	d_data[i] /= scale;
158	0	}
159	0	return *this;
160	0	}
161
162		//! L2 norm squared
163	0	constexpr inline TYPE normL2Sq() const {
164	0	TYPE res = (TYPE)0.0;
165	0	unsigned int i;
166	0	TYPE *data = d_data.get();
167	0	for (i = 0; i < d_size; i++) {
168	0	res += data[i] * data[i];
169	0	}
170	0	return res;
171	0	}
172
173		//! L2 norm
174	0	constexpr TYPE normL2() const { return sqrt(this->normL2Sq()); }
175
176		//! L1 norm
177		constexpr TYPE normL1() const {
178		TYPE res = (TYPE)0.0;
179		unsigned int i;
180		TYPE *data = d_data.get();
181		for (i = 0; i < d_size; i++) {
182		res += fabs(data[i]);
183		}
184		return res;
185		}
186
187		//! L-infinity norm
188		constexpr TYPE normLinfinity() const {
189		TYPE res = (TYPE)(-1.0);
190		unsigned int i;
191		TYPE *data = d_data.get();
192		for (i = 0; i < d_size; i++) {
193		if (fabs(data[i]) > res) {
194		res = fabs(data[i]);
195		}
196		}
197		return res;
198		}
199
200		//! \brief Gets the ID of the entry that has the largest absolute value
201		//! i.e. the entry being used for the L-infinity norm
202		constexpr unsigned int largestAbsValId() const {
203		TYPE res = (TYPE)(-1.0);
204		unsigned int i, id = d_size;
205		TYPE *data = d_data.get();
206		for (i = 0; i < d_size; i++) {
207		if (fabs(data[i]) > res) {
208		res = fabs(data[i]);
209		id = i;
210		}
211		}
212		return id;
213		}
214
215		//! \brief Gets the ID of the entry that has the largest value
216		constexpr unsigned int largestValId() const {
217		TYPE res = (TYPE)(-1.e8);
218		unsigned int i, id = d_size;
219		TYPE *data = d_data.get();
220		for (i = 0; i < d_size; i++) {
221		if (data[i] > res) {
222		res = data[i];
223		id = i;
224		}
225		}
226		return id;
227		}
228
229		//! \brief Gets the ID of the entry that has the smallest value
230		constexpr unsigned int smallestValId() const {
231		TYPE res = (TYPE)(1.e8);
232		unsigned int i, id = d_size;
233		TYPE *data = d_data.get();
234		for (i = 0; i < d_size; i++) {
235		if (data[i] < res) {
236		res = data[i];
237		id = i;
238		}
239		}
240		return id;
241		}
242
243		//! returns the dot product between two Vectors
244	0	TYPE dotProduct(const Vector<TYPE> other) const {
245	0	PRECONDITION(d_size == other.size(),
246	0	"Size mismatch in vector doct product");
247	0	const TYPE *oData = other.getData();
248	0	unsigned int i;
249	0	TYPE res = (TYPE)(0.0);
250	0	TYPE *data = d_data.get();
251	0	for (i = 0; i < d_size; i++) {
252	0	res += (data[i] * oData[i]);
253	0	}
254	0	return res;
255	0	}
256
257		//! Normalize the vector using the L2 norm
258	0	constexpr void normalize() {
259	0	TYPE val = this->normL2();
260	0	if (val < zero_tolerance) {
261	0	throw std::runtime_error("Cannot normalize a zero length vector");
262	0	}
263	0	(*this) /= val;
264	0	}
265
266		//! Set to a random unit vector
267		void setToRandom(unsigned int seed = 0) {
268		// we want to get our own RNG here instead of using the global
269		// one. This is related to Issue285.
270		RDKit::rng_type generator(42u);
271		RDKit::uniform_double dist(0, 1.0);
272		RDKit::double_source_type randSource(generator, dist);
273		if (seed > 0) {
274		generator.seed(seed);
275		} else {
276		// we can't initialize using only clock(), because it's possible
277		// that we'll get here fast enough that clock() will return 0
278		// and generator.seed(0) is an error:
279		generator.seed(clock() + 1);
280		}
281
282		unsigned int i;
283		TYPE *data = d_data.get();
284		for (i = 0; i < d_size; i++) {
285		data[i] = randSource();
286		}
287		this->normalize();
288		}
289
290		private:
291		unsigned int d_size; //!< our length
292		DATA_SPTR d_data;
293		Vector<TYPE> &operator=(const Vector<TYPE> &other);
294		};
295
296		using DoubleVector = Vector<double>;
297
298		//! returns the algebraic tanimoto similarity [defn' from JCIM 46:587-96 (2006)]
299		template <typename T>
300		constexpr double TanimotoSimilarity(const Vector<T> &v1, const Vector<T> &v2) {
301		double numer = v1.dotProduct(v2);
302		if (numer == 0.0) {
303		return 0.0;
304		}
305		double denom = v1.normL2Sq() + v2.normL2Sq() - numer;
306		if (denom == 0.0) {
307		return 0.0;
308		}
309		return numer / denom;
310		}
311		} // end of namespace RDNumeric
312
313		//! ostream operator for Vectors
314		template <typename TYPE>
315		constexpr std::ostream &operator<<(std::ostream &target,
316		const RDNumeric::Vector<TYPE> &vec) {
317		unsigned int siz = vec.size();
318		target << "Size: " << siz << " [";
319		unsigned int i;
320		for (i = 0; i < siz; i++) {
321		target << std::setw(7) << std::setprecision(3) << vec.getVal(i) << ", ";
322		}
323		target << "]\n";
324		return target;
325		}
326
327		#endif