Skip to content

L
libcifpp
Commit be19e4a9 by Maarten L. Hekkelman
Options

Merge branch 'develop' of github.com:PDB-REDO/libcifpp into develop

parents f83850e3 61ce91a9
Showing with 337 additions and 427 deletions
CMakeLists.txt
......@@ -250,7 +250,6 @@ set(project_headers
${PROJECT_SOURCE_DIR}/include/cif++/CifUtils.hpp
${PROJECT_SOURCE_DIR}/include/cif++/CifValidator.hpp
${PROJECT_SOURCE_DIR}/include/cif++/Compound.hpp
${PROJECT_SOURCE_DIR}/include/cif++/Matrix.hpp
${PROJECT_SOURCE_DIR}/include/cif++/PDB2Cif.hpp
${PROJECT_SOURCE_DIR}/include/cif++/PDB2CifRemark3.hpp
${PROJECT_SOURCE_DIR}/include/cif++/Point.hpp
......
include/cif++/Matrix.hpp deleted 100644 → 0
/*-
* SPDX-License-Identifier: BSD-2-Clause
*
* Copyright Maarten L. Hekkelman, Radboud University 2008-2011.
* Copyright (c) 2021 NKI/AVL, Netherlands Cancer Institute
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* 1. Redistributions of source code must retain the above copyright notice, this
* list of conditions and the following disclaimer
* 2. Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
* ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
* (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
* ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
* SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
// --------------------------------------------------------------------
// uBlas compatible matrix types
#pragma once
#include <iostream>
#include <vector>
// matrix is m x n, addressing i,j is 0 <= i < m and 0 <= j < n
// element m i,j is mapped to [i * n + j] and thus storage is row major
template <typename T>
class MatrixBase
{
public:
using value_type = T;
virtual ~MatrixBase() {}
virtual uint32_t dim_m() const = 0;
virtual uint32_t dim_n() const = 0;
virtual value_type &operator()(uint32_t i, uint32_t j) { throw std::runtime_error("unimplemented method"); }
virtual value_type operator()(uint32_t i, uint32_t j) const = 0;
MatrixBase &operator*=(const value_type &rhs);
MatrixBase &operator-=(const value_type &rhs);
};
template <typename T>
MatrixBase<T> &MatrixBase<T>::operator*=(const T &rhs)
{
for (uint32_t i = 0; i < dim_m(); ++i)
{
for (uint32_t j = 0; j < dim_n(); ++j)
{
operator()(i, j) *= rhs;
}
}
return *this;
}
template <typename T>
MatrixBase<T> &MatrixBase<T>::operator-=(const T &rhs)
{
for (uint32_t i = 0; i < dim_m(); ++i)
{
for (uint32_t j = 0; j < dim_n(); ++j)
{
operator()(i, j) -= rhs;
}
}
return *this;
}
template <typename T>
std::ostream &operator<<(std::ostream &lhs, const MatrixBase<T> &rhs)
{
lhs << '[' << rhs.dim_m() << ',' << rhs.dim_n() << ']' << '(';
for (uint32_t i = 0; i < rhs.dim_m(); ++i)
{
lhs << '(';
for (uint32_t j = 0; j < rhs.dim_n(); ++j)
{
if (j > 0)
lhs << ',';
lhs << rhs(i, j);
}
lhs << ')';
}
lhs << ')';
return lhs;
}
template <typename T>
class Matrix : public MatrixBase<T>
{
public:
using value_type = T;
template <typename T2>
Matrix(const MatrixBase<T2> &m)
: m_m(m.dim_m())
, m_n(m.dim_n())
{
m_data = new value_type[m_m * m_n];
for (uint32_t i = 0; i < m_m; ++i)
{
for (uint32_t j = 0; j < m_n; ++j)
operator()(i, j) = m(i, j);
}
}
Matrix()
: m_data(nullptr)
, m_m(0)
, m_n(0)
{
}
Matrix(const Matrix &m)
: m_m(m.m_m)
, m_n(m.m_n)
{
m_data = new value_type[m_m * m_n];
std::copy(m.m_data, m.m_data + (m_m * m_n), m_data);
}
Matrix &operator=(const Matrix &m)
{
value_type *t = new value_type[m.m_m * m.m_n];
std::copy(m.m_data, m.m_data + (m.m_m * m.m_n), t);
delete[] m_data;
m_data = t;
m_m = m.m_m;
m_n = m.m_n;
return *this;
}
Matrix(uint32_t m, uint32_t n, T v = T())
: m_m(m)
, m_n(n)
{
m_data = new value_type[m_m * m_n];
std::fill(m_data, m_data + (m_m * m_n), v);
}
virtual ~Matrix()
{
delete[] m_data;
}
virtual uint32_t dim_m() const { return m_m; }
virtual uint32_t dim_n() const { return m_n; }
virtual value_type operator()(uint32_t i, uint32_t j) const
{
assert(i < m_m);
assert(j < m_n);
return m_data[i * m_n + j];
}
virtual value_type &operator()(uint32_t i, uint32_t j)
{
assert(i < m_m);
assert(j < m_n);
return m_data[i * m_n + j];
}
template <typename Func>
void each(Func f)
{
for (uint32_t i = 0; i < m_m * m_n; ++i)
f(m_data[i]);
}
template <typename U>
Matrix &operator/=(U v)
{
for (uint32_t i = 0; i < m_m * m_n; ++i)
m_data[i] /= v;
return *this;
}
private:
value_type *m_data;
uint32_t m_m, m_n;
};
// --------------------------------------------------------------------
template <typename T>
class SymmetricMatrix : public MatrixBase<T>
{
public:
typedef typename MatrixBase<T>::value_type value_type;
SymmetricMatrix(uint32_t n, T v = T())
: m_owner(true)
, m_n(n)
{
uint32_t N = (m_n * (m_n + 1)) / 2;
m_data = new value_type[N];
std::fill(m_data, m_data + N, v);
}
SymmetricMatrix(const T *data, uint32_t n)
: m_owner(false)
, m_data(const_cast<T *>(data))
, m_n(n)
{
}
virtual ~SymmetricMatrix()
{
if (m_owner)
delete[] m_data;
}
virtual uint32_t dim_m() const { return m_n; }
virtual uint32_t dim_n() const { return m_n; }
T operator()(uint32_t i, uint32_t j) const;
virtual T &operator()(uint32_t i, uint32_t j);
// erase two rows, add one at the end (for neighbour joining)
void erase_2(uint32_t i, uint32_t j);
template <typename Func>
void each(Func f)
{
uint32_t N = (m_n * (m_n + 1)) / 2;
for (uint32_t i = 0; i < N; ++i)
f(m_data[i]);
}
template <typename U>
SymmetricMatrix &operator/=(U v)
{
uint32_t N = (m_n * (m_n + 1)) / 2;
for (uint32_t i = 0; i < N; ++i)
m_data[i] /= v;
return *this;
}
private:
bool m_owner;
value_type *m_data;
uint32_t m_n;
};
template <typename T>
inline T SymmetricMatrix<T>::operator()(uint32_t i, uint32_t j) const
{
return i < j
? m_data[(j * (j + 1)) / 2 + i]
: m_data[(i * (i + 1)) / 2 + j];
}
template <typename T>
inline T &SymmetricMatrix<T>::operator()(uint32_t i, uint32_t j)
{
if (i > j)
std::swap(i, j);
assert(j < m_n);
return m_data[(j * (j + 1)) / 2 + i];
}
template <typename T>
void SymmetricMatrix<T>::erase_2(uint32_t di, uint32_t dj)
{
uint32_t s = 0, d = 0;
for (uint32_t i = 0; i < m_n; ++i)
{
for (uint32_t j = 0; j < i; ++j)
{
if (i != di and j != dj and i != dj and j != di)
{
if (s != d)
m_data[d] = m_data[s];
++d;
}
++s;
}
}
--m_n;
}
template <typename T>
class IdentityMatrix : public MatrixBase<T>
{
public:
typedef typename MatrixBase<T>::value_type value_type;
IdentityMatrix(uint32_t n)
: m_n(n)
{
}
virtual uint32_t dim_m() const { return m_n; }
virtual uint32_t dim_n() const { return m_n; }
virtual value_type operator()(uint32_t i, uint32_t j) const
{
value_type result = 0;
if (i == j)
result = 1;
return result;
}
private:
uint32_t m_n;
};
// --------------------------------------------------------------------
// matrix functions
template <typename T>
Matrix<T> operator*(const MatrixBase<T> &lhs, const MatrixBase<T> &rhs)
{
Matrix<T> result(std::min(lhs.dim_m(), rhs.dim_m()), std::min(lhs.dim_n(), rhs.dim_n()));
for (uint32_t i = 0; i < result.dim_m(); ++i)
{
for (uint32_t j = 0; j < result.dim_n(); ++j)
{
for (uint32_t li = 0, rj = 0; li < lhs.dim_m() and rj < rhs.dim_n(); ++li, ++rj)
result(i, j) += lhs(li, j) * rhs(i, rj);
}
}
return result;
}
template <typename T>
Matrix<T> operator*(const MatrixBase<T> &lhs, T rhs)
{
Matrix<T> result(lhs);
result *= rhs;
return result;
}
template <typename T>
Matrix<T> operator-(const MatrixBase<T> &lhs, const MatrixBase<T> &rhs)
{
Matrix<T> result(std::min(lhs.dim_m(), rhs.dim_m()), std::min(lhs.dim_n(), rhs.dim_n()));
for (uint32_t i = 0; i < result.dim_m(); ++i)
{
for (uint32_t j = 0; j < result.dim_n(); ++j)
{
result(i, j) = lhs(i, j) - rhs(i, j);
}
}
return result;
}
template <typename T>
Matrix<T> operator-(const MatrixBase<T> &lhs, T rhs)
{
Matrix<T> result(lhs.dim_m(), lhs.dim_n());
result -= rhs;
return result;
}
// template <typename T>
// symmetric_matrix<T> hammingDistance(const MatrixBase<T> &lhs, T rhs);
// template <typename T>
// std::vector<T> sum(const MatrixBase<T> &m);
include/cif++/Point.hpp
......@@ -382,7 +382,16 @@ Quaternion Normalize(Quaternion q);
std::tuple<double, Point> QuaternionToAngleAxis(Quaternion q);
Point Centroid(const std::vector<Point> &Points);
Point CenterPoints(std::vector<Point> &Points);
/// \brief Returns how the two sets of points \a a and \b b can be aligned
///
/// \param a The first set of points
/// \param b The second set of points
/// \result The quaternion which should be applied to the points in \a a to
/// obtain the best superposition.
Quaternion AlignPoints(const std::vector<Point> &a, const std::vector<Point> &b);
/// \brief The RMSd for the points in \a a and \a b
double RMSd(const std::vector<Point> &a, const std::vector<Point> &b);
// --------------------------------------------------------------------
......
src/Point.cpp
......@@ -28,12 +28,271 @@
#include <valarray>
#include "cif++/Point.hpp"
#include "cif++/Matrix.hpp"
namespace mmcif
{
// --------------------------------------------------------------------
// uBlas compatible matrix types
// matrix is m x n, addressing i,j is 0 <= i < m and 0 <= j < n
// element m i,j is mapped to [i * n + j] and thus storage is row major
template <typename M>
class MatrixExpression
{
public:
uint32_t dim_m() const { return static_cast<const M&>(*this).dim_m(); }
uint32_t dim_n() const { return static_cast<const M&>(*this).dim_n(); }
double &operator()(uint32_t i, uint32_t j)
{
return static_cast<M&>(*this).operator()(i, j);
}
double operator()(uint32_t i, uint32_t j) const
{
return static_cast<const M&>(*this).operator()(i, j);
}
};
class Matrix : public MatrixExpression<Matrix>
{
public:
template <typename M2>
Matrix(const MatrixExpression<M2> &m)
: m_m(m.dim_m())
, m_n(m.dim_n())
{
m_data = new double[m_m * m_n];
for (uint32_t i = 0; i < m_m; ++i)
{
for (uint32_t j = 0; j < m_n; ++j)
operator()(i, j) = m(i, j);
}
}
Matrix()
: m_data(nullptr)
, m_m(0)
, m_n(0)
{
}
Matrix(const Matrix &m)
: m_m(m.m_m)
, m_n(m.m_n)
{
m_data = new double[m_m * m_n];
std::copy(m.m_data, m.m_data + (m_m * m_n), m_data);
}
Matrix &operator=(const Matrix &m)
{
double *t = new double[m.m_m * m.m_n];
std::copy(m.m_data, m.m_data + (m.m_m * m.m_n), t);
delete[] m_data;
m_data = t;
m_m = m.m_m;
m_n = m.m_n;
return *this;
}
Matrix(uint32_t m, uint32_t n, double v = 0)
: m_m(m)
, m_n(n)
{
m_data = new double[m_m * m_n];
std::fill(m_data, m_data + (m_m * m_n), v);
}
~Matrix()
{
delete[] m_data;
}
uint32_t dim_m() const { return m_m; }
uint32_t dim_n() const { return m_n; }
double operator()(uint32_t i, uint32_t j) const
{
assert(i < m_m);
assert(j < m_n);
return m_data[i * m_n + j];
}
double &operator()(uint32_t i, uint32_t j)
{
assert(i < m_m);
assert(j < m_n);
return m_data[i * m_n + j];
}
private:
double *m_data;
uint32_t m_m, m_n;
};
// --------------------------------------------------------------------
class SymmetricMatrix : public MatrixExpression<SymmetricMatrix>
{
public:
SymmetricMatrix(uint32_t n, double v = 0)
: m_n(n)
{
uint32_t N = (m_n * (m_n + 1)) / 2;
m_data = new double[N];
std::fill(m_data, m_data + N, v);
}
~SymmetricMatrix()
{
delete[] m_data;
}
uint32_t dim_m() const { return m_n; }
uint32_t dim_n() const { return m_n; }
double operator()(uint32_t i, uint32_t j) const
{
return i < j
? m_data[(j * (j + 1)) / 2 + i]
: m_data[(i * (i + 1)) / 2 + j];
}
double &operator()(uint32_t i, uint32_t j)
{
if (i > j)
std::swap(i, j);
assert(j < m_n);
return m_data[(j * (j + 1)) / 2 + i];
}
private:
double *m_data;
uint32_t m_n;
};
class IdentityMatrix : public MatrixExpression<IdentityMatrix>
{
public:
IdentityMatrix(uint32_t n)
: m_n(n)
{
}
uint32_t dim_m() const { return m_n; }
uint32_t dim_n() const { return m_n; }
double operator()(uint32_t i, uint32_t j) const
{
return i == j ? 1 : 0;
}
private:
uint32_t m_n;
};
// --------------------------------------------------------------------
// matrix functions, implemented as expression templates
template<typename M1, typename M2>
class MatrixSubtraction : public MatrixExpression<MatrixSubtraction<M1, M2>>
{
public:
MatrixSubtraction(const M1 &m1, const M2 &m2)
: m_m1(m1), m_m2(m2)
{
assert(m_m1.dim_m() == m_m2.dim_m());
assert(m_m1.dim_n() == m_m2.dim_n());
}
uint32_t dim_m() const { return m_m1.dim_m(); }
uint32_t dim_n() const { return m_m1.dim_n(); }
double operator()(uint32_t i, uint32_t j) const
{
return m_m1(i, j) - m_m2(i, j);
}
private:
const M1 &m_m1;
const M2 &m_m2;
};
template<typename M1, typename M2>
MatrixSubtraction<M1, M2> operator-(const MatrixExpression<M1> &m1, const MatrixExpression<M2> &m2)
{
return MatrixSubtraction(*static_cast<const M1*>(&m1), *static_cast<const M2*>(&m2));
}
template<typename M>
class MatrixMultiplication : public MatrixExpression<MatrixMultiplication<M>>
{
public:
MatrixMultiplication(const M &m, double v)
: m_m(m), m_v(v)
{
}
uint32_t dim_m() const { return m_m.dim_m(); }
uint32_t dim_n() const { return m_m.dim_n(); }
double operator()(uint32_t i, uint32_t j) const
{
return m_m(i, j) * m_v;
}
private:
const M &m_m;
double m_v;
};
template<typename M>
MatrixMultiplication<M> operator*(const MatrixExpression<M> &m, double v)
{
return MatrixMultiplication(*static_cast<const M*>(&m), v);
}
// --------------------------------------------------------------------
template<class M1>
void cofactors(const M1& m, SymmetricMatrix& cf)
{
const size_t ixs[4][3] =
{
{ 1, 2, 3 },
{ 0, 2, 3 },
{ 0, 1, 3 },
{ 0, 1, 2 }
};
for (size_t x = 0; x < 4; ++x)
{
const size_t* ix = ixs[x];
for (size_t y = x; y < 4; ++y)
{
const size_t* iy = ixs[y];
cf(x, y) =
m(ix[0], iy[0]) * m(ix[1], iy[1]) * m(ix[2], iy[2]) +
m(ix[0], iy[1]) * m(ix[1], iy[2]) * m(ix[2], iy[0]) +
m(ix[0], iy[2]) * m(ix[1], iy[0]) * m(ix[2], iy[1]) -
m(ix[0], iy[2]) * m(ix[1], iy[1]) * m(ix[2], iy[0]) -
m(ix[0], iy[1]) * m(ix[1], iy[0]) * m(ix[2], iy[2]) -
m(ix[0], iy[0]) * m(ix[1], iy[2]) * m(ix[2], iy[1]);
if ((x + y) % 2 == 1)
cf(x, y) *= -1;
}
}
}
// --------------------------------------------------------------------
Quaternion Normalize(Quaternion q)
{
......@@ -175,7 +434,7 @@ double LargestDepressedQuarticSolution(double a, double b, double c)
Quaternion AlignPoints(const std::vector<Point>& pa, const std::vector<Point>& pb)
{
// First calculate M, a 3x3 Matrix containing the sums of products of the coordinates of A and B
Matrix<double> M(3, 3, 0);
Matrix M(3, 3, 0);
for (uint32_t i = 0; i < pa.size(); ++i)
{
......@@ -188,7 +447,7 @@ Quaternion AlignPoints(const std::vector<Point>& pa, const std::vector<Point>& p
}
// Now calculate N, a symmetric 4x4 Matrix
SymmetricMatrix<double> N(4);
SymmetricMatrix N(4);
N(0, 0) = M(0, 0) + M(1, 1) + M(2, 2);
N(0, 1) = M(1, 2) - M(2, 1);
......@@ -225,6 +484,7 @@ Quaternion AlignPoints(const std::vector<Point>& pa, const std::vector<Point>& p
M(1, 2) * M(2, 0) * M(0, 1) +
M(2, 1) * M(1, 0) * M(0, 2));
// E is the determinant of N:
double E =
(N(0,0) * N(1,1) - N(0,1) * N(0,1)) * (N(2,2) * N(3,3) - N(2,3) * N(2,3)) +
(N(0,1) * N(0,2) - N(0,0) * N(2,1)) * (N(2,1) * N(3,3) - N(2,3) * N(1,3)) +
......@@ -237,44 +497,20 @@ Quaternion AlignPoints(const std::vector<Point>& pa, const std::vector<Point>& p
double lm = LargestDepressedQuarticSolution(C, D, E);
// calculate t = (N - λI)
Matrix<double> li = IdentityMatrix<double>(4) * lm;
Matrix<double> t = N - li;
Matrix t = N - IdentityMatrix(4) * lm;
// calculate a Matrix of cofactors for t
Matrix<double> cf(4, 4);
// calculate a Matrix of cofactors for t, since N is symmetric, t must be symmetric as well and so will be cf
SymmetricMatrix cf(4);
cofactors(t, cf);
const uint32_t ixs[4][3] =
int maxR = 0;
for (int r = 1; r < 4; ++r)
{
{ 1, 2, 3 },
{ 0, 2, 3 },
{ 0, 1, 3 },
{ 0, 1, 2 }
};
uint32_t maxR = 0;
for (uint32_t r = 0; r < 4; ++r)
{
const uint32_t* ir = ixs[r];
for (uint32_t c = 0; c < 4; ++c)
{
const uint32_t* ic = ixs[c];
cf(r, c) =
t(ir[0], ic[0]) * t(ir[1], ic[1]) * t(ir[2], ic[2]) +
t(ir[0], ic[1]) * t(ir[1], ic[2]) * t(ir[2], ic[0]) +
t(ir[0], ic[2]) * t(ir[1], ic[0]) * t(ir[2], ic[1]) -
t(ir[0], ic[2]) * t(ir[1], ic[1]) * t(ir[2], ic[0]) -
t(ir[0], ic[1]) * t(ir[1], ic[0]) * t(ir[2], ic[2]) -
t(ir[0], ic[0]) * t(ir[1], ic[2]) * t(ir[2], ic[1]);
}
if (r > maxR and cf(r, 0) > cf(maxR, 0))
if (cf(r, 0) > cf(maxR, 0))
maxR = r;
}
// NOTE the negation of the y here, why? Maybe I swapped r/c above?
Quaternion q(cf(maxR, 0), cf(maxR, 1), -cf(maxR, 2), cf(maxR, 3));
Quaternion q(cf(maxR, 0), cf(maxR, 1), cf(maxR, 2), cf(maxR, 3));
q = Normalize(q);
return q;
......
test/unit-test.cpp
......@@ -35,6 +35,9 @@
#include "cif++/BondMap.hpp"
#include "cif++/CifValidator.hpp"
namespace tt = boost::test_tools;
std::filesystem::path gTestDir = std::filesystem::current_path(); // filled in first test
// --------------------------------------------------------------------
......@@ -1700,3 +1703,57 @@ BOOST_AUTO_TEST_CASE(reading_file_1)
cif::File file;
BOOST_CHECK_THROW(file.load(is), std::runtime_error);
}
// 3d tests
using namespace mmcif;
BOOST_AUTO_TEST_CASE(t1)
{
// std::random_device rnd;
// std::mt19937 gen(rnd());
// std::uniform_real_distribution<float> dis(0, 1);
// Quaternion q{ dis(gen), dis(gen), dis(gen), dis(gen) };
// q = Normalize(q);
// Quaternion q{ 0.1, 0.2, 0.3, 0.4 };
Quaternion q{ 0.5, 0.5, 0.5, 0.5 };
q = Normalize(q);
const auto &&[angle0, axis0] = QuaternionToAngleAxis(q);
std::vector<Point> p1{
{ 16.979, 13.301, 44.555 },
{ 18.150, 13.525, 43.680 },
{ 18.656, 14.966, 43.784 },
{ 17.890, 15.889, 44.078 },
{ 17.678, 13.270, 42.255 },
{ 16.248, 13.734, 42.347 },
{ 15.762, 13.216, 43.724 }
};
auto p2 = p1;
Point c1 = CenterPoints(p1);
for (auto &p : p2)
p.rotate(q);
Point c2 = CenterPoints(p2);
auto q2 = AlignPoints(p1, p2);
const auto &&[angle, axis] = QuaternionToAngleAxis(q2);
BOOST_TEST(std::fmod(360 + angle, 360) == std::fmod(360 - angle0, 360), tt::tolerance(0.01));
for (auto &p : p1)
p.rotate(q2);
float rmsd = RMSd(p1, p2);
BOOST_TEST(rmsd < 1e-5);
// std::cout << "rmsd: " << RMSd(p1, p2) << std::endl;
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment