project(gromacs-fitng CXX)
+set(CMAKE_CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+set(CMAKE_CXX_EXTENSIONS OFF)
+
if (NOT CMAKE_BUILD_TYPE)
set(CMAKE_BUILD_TYPE "Release" CACHE STRING "Choose the type of build, options are: Debug Release RelWithDebInfo MinSizeRel." FORCE)
endif()
set(GROMACS_SUFFIX ${GMX_SUFFIX})
endif()
-find_package(GROMACS 2016 REQUIRED)
+if (GMX_OPENMP)
+ find_package(OpenMP REQUIRED)
+else()
+ find_package(OpenMP)
+endif()
+
+find_package(GROMACS 2019 REQUIRED)
gromacs_check_double(GMX_DOUBLE)
gromacs_check_compiler(CXX)
-include(gmxManageOpenMP)
add_definitions(${GROMACS_DEFINITIONS})
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2018,2019, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \file
+ * \brief
+ * Declares gmx::TopologyInformation.
+ *
+ * \author Teemu Murtola <teemu.murtola@gmail.com>
+ * \author Mark Abraham <mark.j.abraham@gmail.com>
+ * \inlibraryapi
+ * \ingroup module_trajectoryanalysis
+ */
+#ifndef GMX_TRAJECTORYANALYSIS_TOPOLOGYINFORMATION_H
+#define GMX_TRAJECTORYANALYSIS_TOPOLOGYINFORMATION_H
+
+#include <memory>
+#include <string>
+#include <vector>
+
+#include "gromacs/math/vectypes.h"
+#include "gromacs/topology/atoms.h"
+#include "gromacs/topology/topology.h"
+#include "gromacs/utility/classhelpers.h"
+
+//! Forward declaration
+typedef struct gmx_rmpbc* gmx_rmpbc_t;
+
+namespace gmx
+{
+
+template<typename T>
+class ArrayRef;
+
+class TopologyInformation;
+class TrajectoryAnalysisRunnerCommon;
+
+/*! \libinternal
+ * \brief Topology information available to a trajectory analysis module.
+ *
+ * This class is used to provide topology information to trajectory
+ * analysis modules and to manage memory for them. Having a single
+ * wrapper object instead of passing each item separately makes
+ * TrajectoryAnalysisModule interface simpler, and also reduces the
+ * need to change existing code if additional information is added.
+ *
+ * It is intended that eventually most clients of this class will be
+ * analysis tools ported to the new analysis framework, but we will
+ * use this infrastructure also from the legacy analysis tools during
+ * the transition period. That will make it easier to put those tools
+ * under tests, and eventually port them.
+ *
+ * Methods in this class do not throw if not explicitly stated.
+ *
+ * The main data content is constant once loaded, but some content is
+ * constructed only when required (e.g. atoms_ and
+ * expandedTopology_). Their data members are mutable, so that the
+ * lazy construction idiom works properly. Some clients wish to modify
+ * the t_atoms, so there is an efficient mechanism for them to get a
+ * copy they can modify without disturbing this class. (The
+ * implementation releases the cached lazily constructed atoms_, but
+ * from the point of view of the caller, that is a copy.) The getters
+ * and copy functions are const for callers, which correctly expresses
+ * that the topology information is not being changed, merely copied
+ * and presented in a different form.
+ *
+ * \ingroup module_trajectoryanalysis
+ */
+class TopologyInformation
+{
+public:
+ //! Returns true if a topology file was loaded.
+ bool hasTopology() const { return hasLoadedMtop_; }
+ //! Returns true if a full topology file was loaded.
+ bool hasFullTopology() const { return bTop_; }
+ /*! \brief Builder function to fill the contents of
+ * TopologyInformation in \c topInfo from \c filename.
+ *
+ * Different tools require, might need, would benefit from, or
+ * do not need topology information. This functions implements
+ * the two-phase construction that is currently needed to
+ * support that.
+ *
+ * Any coordinate or run input file format will work, but the
+ * kind of data available from the getter methods afterwards
+ * will vary. For example, the mtop() available after reading
+ * a plain structure file will have a single molecule block and
+ * molecule type, regardless of contents.
+ *
+ * After reading, this object can return many kinds of primary
+ * and derived data structures to its caller.
+ *
+ * \todo This should throw upon error but currently does
+ * not. */
+ void fillFromInputFile(const std::string& filename);
+ /*! \brief Returns the loaded topology, or nullptr if not loaded. */
+ gmx_mtop_t* mtop() const { return mtop_.get(); }
+ //! Returns the loaded topology fully expanded, or nullptr if no topology is available.
+ const gmx_localtop_t* expandedTopology() const;
+ /*! \brief Returns a read-only handle to the fully expanded
+ * atom data arrays, which might be valid but empty if no
+ * topology is available. */
+ const t_atoms* atoms() const;
+ /*! \brief Copies the fully expanded atom data arrays, which
+ * might be valid but empty if no topology is available. */
+ AtomsDataPtr copyAtoms() const;
+ //! Returns the ePBC field from the topology.
+ int ePBC() const { return ePBC_; }
+ /*! \brief
+ * Gets the configuration positions from the topology file.
+ *
+ * If TrajectoryAnalysisSettings::efUseTopX has not been specified,
+ * this method should not be called.
+ *
+ * \throws APIError if topology position coordinates are not available
+ */
+ ArrayRef<const RVec> x() const;
+ /*! \brief
+ * Gets the configuration velocities from the topology file.
+ *
+ * If TrajectoryAnalysisSettings::efUseTopV has not been specified,
+ * this method should not be called.
+ *
+ * \throws APIError if topology velocity coordinates are not available
+ */
+ ArrayRef<const RVec> v() const;
+ /*! \brief
+ * Gets the configuration box from the topology file.
+ *
+ * \param[out] box Box size from the topology file, must not be nullptr.
+ */
+ void getBox(matrix box) const;
+ /*! \brief Returns a name for the topology.
+ *
+ * If a full topology was read from a a file, returns the name
+ * it contained, otherwise the empty string. */
+ const char* name() const;
+
+ TopologyInformation();
+ ~TopologyInformation();
+
+private:
+ //! The topology structure, or nullptr if no topology loaded.
+ std::unique_ptr<gmx_mtop_t> mtop_;
+ //! Whether a topology has been loaded.
+ bool hasLoadedMtop_;
+ //! The fully expanded topology structure, nullptr if not yet constructed.
+ mutable ExpandedTopologyPtr expandedTopology_;
+ //! The fully expanded atoms data structure, nullptr if not yet constructed.
+ mutable AtomsDataPtr atoms_;
+ //! true if full tpx file was loaded, false otherwise.
+ bool bTop_;
+ //! Position coordinates from the topology (can be nullptr).
+ std::vector<RVec> xtop_;
+ //! Velocity coordinates from the topology (can be nullptr).
+ std::vector<RVec> vtop_;
+ //! The box loaded from the topology file.
+ matrix boxtop_{};
+ //! The ePBC field loaded from the topology file.
+ int ePBC_;
+
+ // TODO This type is probably movable if we need that.
+ GMX_DISALLOW_COPY_AND_ASSIGN(TopologyInformation);
+
+ /*! \brief
+ * Needed to initialize the data.
+ */
+ friend class TrajectoryAnalysisRunnerCommon;
+};
+
+//! Convenience overload useful for implementing legacy tools.
+gmx_rmpbc_t gmx_rmpbc_init(const gmx::TopologyInformation& topInfo);
+
+} // namespace gmx
+
+#endif
set(GROMACS_CXX_FLAGS "${GROMACS_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
set( CMAKE_CXX_FLAGS "-fopenmp")
-add_executable(fitng fitng.cpp)
-#add_library(fittingn STATIC fittingn.h fittingn.cpp)
+
+add_executable(fitng fitng.cpp newfit.cpp)
include_directories(
${GROMACS_INCLUDE_DIRS}
- ${CMAKE_SOURCE_DIR}/include)
+ ${CMAKE_SOURCE_DIR}/include
+ ${CMAKE_SOURCE_DIR}/src)
set_target_properties(fitng PROPERTIES
COMPILE_FLAGS "${GROMACS_CXX_FLAGS}")
-target_link_libraries(fitng ${GROMACS_LIBRARIES})
+target_link_libraries(fitng ${GROMACS_LIBRARIES} ${OpenMP_CXX_LIBRARIES})
#include <omp.h>
-#include <gromacs/trajectoryanalysis.h>
-#include <gromacs/pbcutil/pbc.h>
-#include <gromacs/utility/smalloc.h>
-#include <gromacs/math/do_fit.h>
#include <gromacs/fileio/trxio.h>
#include <gromacs/trajectoryanalysis/topologyinformation.h>
-using namespace gmx;
-
-using gmx::RVec;
-
-double F (double aix, double aiy, double aiz, double bix_plus_x, double biy_plus_y, double biz_plus_z, double A, double B, double C) {
- return sqrt( pow(aix + biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * bix_plus_x, 2) +
- pow(aiy - biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * bix_plus_x, 2) +
- pow(aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y, 2) );
-}
-
-void searchF0xyzabc (double &F, double &Fx, double &Fy, double &Fz, double &Fa, double &Fb, double &Fc, double aix, double aiy, double aiz, double bix_plus_x, double biy_plus_y, double biz_plus_z, double A, double B, double C) {
- double t01, t02, t03, t04, t05, t06, t07;
- t01 = pow(aix + biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * bix_plus_x, 2);
- t02 = pow(aiy - biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * bix_plus_x, 2);
- t03 = pow(aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y, 2);
- t04 = sqrt(t01 + t02 + t03);
- t05 = (aix + biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * bix_plus_x);
- t06 = (aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y);
- t07 = (aiy - biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * bix_plus_x);
-
- //#pragma omp parallel sections
- //{
- //#pragma omp section
- F += t04;
- //#pragma omp section
- Fx += -(2 * cos(B) * cos(C) * t05 - 2 * sin(B) * t06 + 2 * cos(B) * sin(C) * t07) / (2 * t04);
- //#pragma omp section
- Fy += -(2 * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) * t07 - 2 * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) * t05 + 2 * cos(B) * sin(A) * t06) / (2 * t04);
- //#pragma omp section
- Fz += -(2 * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) * t05 - 2 * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) * t07 + 2 * cos(A) * cos(B) * t06) / (2 * t04);
- //#pragma omp section
- Fa += -(2 * (cos(A) * cos(B) * biy_plus_y - cos(B) * sin(A) * biz_plus_z) * t06 -
- 2 * (biy_plus_y * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) + biz_plus_z * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C))) * t07 +
- 2 * (biy_plus_y * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) + biz_plus_z * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B))) * t05) / (2 * t04);
- //#pragma omp section
- Fb += -(2 * (cos(A) * cos(B) * sin(C) * biz_plus_z - sin(B) * sin(C) * bix_plus_x + cos(B) * sin(A) * sin(C) * biy_plus_y) * t07 +
- 2 * (cos(A) * cos(B) * cos(C) * biz_plus_z - cos(C) * sin(B) * bix_plus_x + cos(B) * cos(C) * sin(A) * biy_plus_y) * t05 -
- 2 * (cos(B) * bix_plus_x + sin(A) * sin(B) * biy_plus_y + cos(A) * sin(B) * biz_plus_z) * t06) / (2 * t04);
- //#pragma omp section
- Fc += (2 * (biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) - biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) + cos(B) * sin(C) * bix_plus_x) * t05 -
- 2 * (biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) + cos(B) * cos(C) * bix_plus_x) * t07) / (2 * t04);
- //}
-}
-
-/*void searchL (double &fLd, double &fLdd, double aix, double aiy, double aiz, double bix, double biy, double biz, double x, double y, double z, double A, double B, double C, double L, double fx, double fy, double fz, double fa, double fb, double fc) {
- double AmLfa, BmLfb, CmLfc, BIXpXmLfx, BIYpYmLfy, BIZpZmLfz, t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11, t12, t13, t14, t15, t16, t17, t18, t19, t20;
- AmLfa = A - L * fa;
- BmLfb = B - L * fb;
- CmLfc = C - L * fc;
- BIXpXmLfx = bix + x - L * fx;
- BIYpYmLfy = biy + y - L * fy;
- BIZpZmLfz = biz + z - L * fz;
- t1 = (cos(AmLfa) * cos(CmLfc) + sin(AmLfa) * sin(BmLfb) * sin(CmLfc));
- t2 = (cos(CmLfc) * sin(AmLfa) - cos(AmLfa) * sin(BmLfb) * sin(CmLfc));
- t3 = (cos(AmLfa) * sin(CmLfc) - cos(CmLfc) * sin(AmLfa) * sin(BmLfb));
- t4 = (sin(AmLfa) * sin(CmLfc) + cos(AmLfa) * cos(CmLfc) * sin(BmLfb));
- t5 = pow(aiz + sin(BmLfb) * BIXpXmLfx - cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - cos(BmLfb) * sin(AmLfa) * BIYpYmLfy, 2);
- t6 = pow(aix + t3 * BIYpYmLfy - t4 * BIZpZmLfz - cos(BmLfb) * cos(CmLfc) * BIXpXmLfx, 2);
- t7 = pow(aiy - t1 * BIYpYmLfy + t2 * BIZpZmLfz - cos(BmLfb) * sin(CmLfc) * BIXpXmLfx, 2);
-
- t8 = fc * sin(AmLfa) * sin(CmLfc);
- t9 = fa * cos(AmLfa) * cos(CmLfc);
- t10 = fb * cos(AmLfa) * cos(BmLfb) * sin(CmLfc);
- t11 = fc * cos(AmLfa) * cos(CmLfc) * sin(BmLfb);
- t12 = fa * sin(AmLfa) * sin(BmLfb) * sin(CmLfc);
-
- t13 = fa * cos(AmLfa) * sin(BmLfb) * sin(CmLfc);
- t14 = fc * cos(AmLfa) * sin(CmLfc);
- t15 = fa * cos(CmLfc) * sin(AmLfa);
- t16 = fb * cos(BmLfb) * sin(AmLfa) * sin(CmLfc);
- t17 = fc * cos(CmLfc) * sin(AmLfa) * sin(BmLfb);
-
- t18 = fx * cos(BmLfb) * sin(CmLfc);
- t19 = fc * cos(BmLfb) * cos(CmLfc);
- t20 = fb * sin(BmLfb) * sin(CmLfc);
-
- fLd += (2 * (aiy - t1 * BIYpYmLfy + t2 * BIZpZmLfz - cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
- ( BIZpZmLfz * (t8 - t9 + t10 + t11 - t12) + fy * t1 - fz * t2 + BIYpYmLfy * (t13 - t14 - t15 + t16 + t17) + t18 + t19 * BIXpXmLfx - t20 * BIXpXmLfx) -
- 2 * (aiz + sin(BmLfb) * BIXpXmLfx - cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - cos(BmLfb) * sin(AmLfa) * BIYpYmLfy) *
- (fx * sin(BmLfb) - fy * cos(BmLfb) * sin(AmLfa) + fb * cos(BmLfb) * BIXpXmLfx - fz * cos(AmLfa) * cos(BmLfb) -
- fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy + fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy) +
- 2 * (aix + t3 * BIYpYmLfy - t4 * BIZpZmLfz - cos(BmLfb) * cos(CmLfc) * BIXpXmLfx) *
- (BIZpZmLfz * (fa * cos(AmLfa) * sin(CmLfc) + fc * cos(CmLfc) * sin(AmLfa) + fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) - fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) +
- BIYpYmLfy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) + t9 * sin(BmLfb) +
- fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) - fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) - fy * t3 +
- fz * t4 + fx * cos(BmLfb) * cos(CmLfc) - fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx)) /
- (2 * sqrt( t5 + t6 + pow(aiy - t1 * BIYpYmLfy + t2 * BIZpZmLfz - cos(BmLfb) * sin(CmLfc) * BIXpXmLfx, 2)));
-
-
- fLdd += ( (2 * aiz + 2 * sin(BmLfb) * BIXpXmLfx - 2 * cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - 2 * cos(BmLfb) * sin(AmLfa) * BIYpYmLfy) *
- (2 * fb * fx * cos(BmLfb) - pow(fb, 2) * sin(BmLfb) * BIXpXmLfx - 2 * fa * fy * cos(AmLfa) * cos(BmLfb) + 2 * fa * fz * cos(BmLfb) * sin(AmLfa) + 2 * fb * fz * cos(AmLfa) * sin(BmLfb) +
- 2 * fb * fy * sin(AmLfa) * sin(BmLfb) + pow(fa, 2) * cos(AmLfa) * cos(BmLfb) * BIZpZmLfz + pow(fb, 2) * cos(AmLfa) * cos(BmLfb) * BIZpZmLfz +
- pow(fa, 2) * cos(BmLfb) * sin(AmLfa) * BIYpYmLfy + pow(fb, 2) * cos(BmLfb) * sin(AmLfa) * BIYpYmLfy + 2 * fa * fb * cos(AmLfa) * sin(BmLfb) * BIYpYmLfy -
- 2 * fa * fb * sin(AmLfa) * sin(BmLfb) * BIZpZmLfz) +
- (BIZpZmLfz * (fa * cos(AmLfa) * sin(CmLfc) + fc * cos(CmLfc) * sin(AmLfa) +
- fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) - fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) +
- BIYpYmLfy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) + t9 * sin(BmLfb) +
- fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) - fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) - fy * (cos(AmLfa) * sin(CmLfc) -
- cos(CmLfc) * sin(AmLfa) * sin(BmLfb)) + fz * t4 + fx * cos(BmLfb) * cos(CmLfc) -
- fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
- (2 * BIZpZmLfz * (fa * cos(AmLfa) * sin(CmLfc) +
- fc * cos(CmLfc) * sin(AmLfa) + fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) -
- fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) + 2 * BIYpYmLfy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) +
- t9 * sin(BmLfb) + fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) - fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) -
- 2 * fy * t3 + 2 * fz * t4 +
- 2 * fx * cos(BmLfb) * cos(CmLfc) - 2 * fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - 2 * fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) +
- (2 * aix + 2 * t3 * BIYpYmLfy - 2 * (sin(AmLfa) * sin(CmLfc) +
- cos(AmLfa) * cos(CmLfc) * sin(BmLfb)) * BIZpZmLfz - 2 * cos(BmLfb) * cos(CmLfc) * BIXpXmLfx) *
- (BIZpZmLfz * (pow(fa, 2) * sin(AmLfa) * sin(CmLfc) +
- pow(fc, 2) * sin(AmLfa) * sin(CmLfc) - 2 * fa * fc * cos(AmLfa) * cos(CmLfc) + pow(fa, 2) * cos(AmLfa) * cos(CmLfc) * sin(BmLfb) +
- pow(fb, 2) * cos(AmLfa) * cos(CmLfc) * sin(BmLfb) + pow(fc, 2) * cos(AmLfa) * cos(CmLfc) * sin(BmLfb) + 2 * fa * fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) +
- 2 * fb * fc * cos(AmLfa) * cos(BmLfb) * sin(CmLfc) - 2 * fa * fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) + BIYpYmLfy * (pow(fa, 2) * cos(CmLfc) * sin(AmLfa) * sin(BmLfb) -
- pow(fc, 2) * cos(AmLfa) * sin(CmLfc) - 2 * fa * fc * cos(CmLfc) * sin(AmLfa) - pow(fa, 2) * cos(AmLfa) * sin(CmLfc) + pow(fb, 2) * cos(CmLfc) * sin(AmLfa) * sin(BmLfb) +
- pow(fc, 2) * cos(CmLfc) * sin(AmLfa) * sin(BmLfb) - 2 * fa * fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) + 2 * fa * fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc) +
- 2 * fb * fc * cos(BmLfb) * sin(AmLfa) * sin(CmLfc)) - 2 * fz * (fa * cos(AmLfa) * sin(CmLfc) + fc * cos(CmLfc) * sin(AmLfa) +
- fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) - fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) -
- 2 * fy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) + t9 * sin(BmLfb) + fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) -
- fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) + 2 * fb * fx * cos(CmLfc) * sin(BmLfb) + 2 * fc * t18 + pow(fb, 2) * cos(BmLfb) * cos(CmLfc) * BIXpXmLfx +
- pow(fc, 2) * cos(BmLfb) * cos(CmLfc) * BIXpXmLfx - 2 * fb * fc * sin(BmLfb) * sin(CmLfc) * BIXpXmLfx) +
- (BIZpZmLfz * (t8 - t9 + t10 + t11 - t12) + fy * t1 - fz * t2 + BIYpYmLfy * (t13 - t14 - t15 + t16 + t17) + t18 + t19 * BIXpXmLfx - t20 * BIXpXmLfx) *
- (2 * BIZpZmLfz * (t8 - t9 + t10 + t11 - t12) + 2 * fy * t1 - 2 * fz * t2 + 2 * BIYpYmLfy * (t13 - t14 - t15 + t16 + t17) + 2 * t18 + 2 * t19 * BIXpXmLfx - 2 * t20 * BIXpXmLfx) +
- (fx * sin(BmLfb) - fy * cos(BmLfb) * sin(AmLfa) + fb * cos(BmLfb) * BIXpXmLfx - fz * cos(AmLfa) * cos(BmLfb) - fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy +
- fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy) *
- (2 * fx * sin(BmLfb) - 2 * fy * cos(BmLfb) * sin(AmLfa) + 2 * fb * cos(BmLfb) * BIXpXmLfx - 2 * fz * cos(AmLfa) * cos(BmLfb) - 2 * fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy +
- 2 * fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + 2 * fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + 2 * fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy) +
- (2 * aiy - 2 * t1 * BIYpYmLfy + 2 * (cos(CmLfc) * sin(AmLfa) - cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) * BIZpZmLfz - 2 * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
- (BIZpZmLfz * (pow(fa, 2) * cos(AmLfa) * sin(BmLfb) * sin(CmLfc) -
- pow(fc, 2) * cos(CmLfc) * sin(AmLfa) - 2 * fa * t14 - pow(fa, 2) * cos(CmLfc) * sin(AmLfa) + pow(fb, 2) * cos(AmLfa) * sin(BmLfb) * sin(CmLfc) +
- pow(fc, 2) * cos(AmLfa) * sin(BmLfb) * sin(CmLfc) - 2 * fb * fc * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) + 2 * fa * t16 +
- 2 * fa * t17) + BIYpYmLfy * (pow(fa, 2) * cos(AmLfa) * cos(CmLfc) + pow(fc, 2) * cos(AmLfa) * cos(CmLfc) - 2 * fa * t8 +
- pow(fa, 2) * sin(AmLfa) * sin(BmLfb) * sin(CmLfc) + pow(fb, 2) * sin(AmLfa) * sin(BmLfb) * sin(CmLfc) + pow(fc, 2) * sin(AmLfa) * sin(BmLfb) * sin(CmLfc) -
- 2 * fa * t10 - 2 * fa * t11 - 2 * fb * t19 * sin(AmLfa)) - 2 * fz * (t8 - t9 + t10 + t11 - t12) - 2 * fy * (t13 - t14 - t15 + t16 + t17) -
- 2 * fc * fx * cos(BmLfb) * cos(CmLfc) + 2 * fb * fx * sin(BmLfb) * sin(CmLfc) +
- pow(fb, 2) * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx + pow(fc, 2) * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx + 2 * fb * fc * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx)) /
- (2 * sqrt( t5 + t6 + t7)) -
- ( (2 * (aiy - t1 * BIYpYmLfy + t2 * BIZpZmLfz - cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
- (BIZpZmLfz * (t8 - t9 + t10 + t11 - t12) + fy * t1 - fz * t2 + BIYpYmLfy * (t13 - t14 - t15 + t16 + t17) + t18 + t19 * BIXpXmLfx - t20 * BIXpXmLfx) -
- 2 * (aiz + sin(BmLfb) * BIXpXmLfx - cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - cos(BmLfb) * sin(AmLfa) * BIYpYmLfy) *
- (fx * sin(BmLfb) - fy * cos(BmLfb) * sin(AmLfa) + fb * cos(BmLfb) * BIXpXmLfx - fz * cos(AmLfa) * cos(BmLfb) - fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy +
- fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy) +
- 2 * (aix + t3 * BIYpYmLfy - t4 * BIZpZmLfz - cos(BmLfb) * cos(CmLfc) * BIXpXmLfx) *
- (BIZpZmLfz * (fa * cos(AmLfa) * sin(CmLfc) + fc * cos(CmLfc) * sin(AmLfa) + fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) - fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) +
- BIYpYmLfy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) + t9 * sin(BmLfb) + fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) - fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) -
- fy * t3 + fz * t4 + fx * cos(BmLfb) * cos(CmLfc) - fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx)) *
- ( (2 * aix + 2 * t3 * BIYpYmLfy - 2 * t4 * BIZpZmLfz - 2 * cos(BmLfb) * cos(CmLfc) * BIXpXmLfx) *
- (BIZpZmLfz * (fa * cos(AmLfa) * sin(CmLfc) + fc * cos(CmLfc) * sin(AmLfa) + fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) - fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) +
- BIYpYmLfy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) + t9 * sin(BmLfb) + fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) - fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) -
- fy * t3 + fz * t4 + fx * cos(BmLfb) * cos(CmLfc) - fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) +
- (2 * aiy - 2 * t1 * BIYpYmLfy + 2 * t2 * BIZpZmLfz - 2 * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
- (BIZpZmLfz * (t8 - t9 + t10 + t11 - t12) + fy * (cos(AmLfa) * cos(CmLfc) + sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) - fz * t2 + BIYpYmLfy * (t13 - t14 - t15 + t16 + t17) +
- t18 + t19 * BIXpXmLfx - t20 * BIXpXmLfx) - (2 * aiz + 2 * sin(BmLfb) * BIXpXmLfx - 2 * cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - 2 * cos(BmLfb) * sin(AmLfa) * BIYpYmLfy) *
- (fx * sin(BmLfb) - fy * cos(BmLfb) * sin(AmLfa) + fb * cos(BmLfb) * BIXpXmLfx - fz * cos(AmLfa) * cos(BmLfb) - fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy +
- fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy))) /
- (4 * pow ( t5 + t6 + t7, 3 / 2));
-
-}*/
-
-void ApplyFit (std::vector< RVec > &b, double x, double y, double z, double A, double B, double C) {
- #pragma omp parallel
- {
- #pragma omp for schedule(dynamic)
- for (int i = 0; i < b.size(); i++) {
- RVec temp = b[i];
- b[i][0] = (temp[2] + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - (temp[1] + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) + cos(B) * cos(C) * (temp[0] + x);
- b[i][1] = (temp[1] + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) - (temp[2] + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) + cos(B) * sin(C) * (temp[0] + x);
- b[i][2] = cos(A) * cos(B) * (temp[2] + z) - sin(B) * (temp[0] + x) + cos(B) * sin(A) * (temp[1] + y);
- }
- }
-}
-
-void CalcMid (std::vector< RVec > a, std::vector< RVec > b, RVec &midA, RVec &midB, std::vector< std::pair< int, int > > pairs) {
- midA[0] = 0;
- midA[1] = 0;
- midA[2] = 0;
-
- midB[0] = 0;
- midB[1] = 0;
- midB[2] = 0;
-
- for (int i = 0; i < pairs.size(); i++) {
- rvec_inc(midA, a[pairs[i].first]);
- rvec_inc(midB, b[pairs[i].second]);
- }
- midA[0] /= pairs.size();
- midA[1] /= pairs.size();
- midA[2] /= pairs.size();
-
- midB[0] /= pairs.size();
- midB[1] /= pairs.size();
- midB[2] /= pairs.size();
-}
-
-/*float*/double FrameMidLength(std::vector< RVec > a) {
- RVec b;
- b[0] = 0;
- b[1] = 0;
- b[2] = 0;
- for (int i = 0; i < a.size(); i++) {
- b += a[i];
- }
- return b.norm();
-}
-
-void reset_coord(std::vector< RVec > &frame, RVec move) {
- for (int i = 0; i < frame.size(); i++) {
- frame[i] -= move;
- }
-}
-
-void MyFitNew (std::vector< RVec > a, std::vector< RVec > &b, std::vector< std::pair< int, int > > pairs, double FtCnst) {
- double f1 = 0, f2 = 0, fx = 0, fy = 0, fz = 0, fa = 0, fb = 0, fc = 0, l = 1;
- RVec ma, mb;
- CalcMid(a, b, ma, mb, pairs);
-
- reset_coord(a, ma);
- reset_coord(b, mb);
-
- rvec_dec(ma, mb);
-
- double x = ma[0], y = ma[1], z = ma[2], A = 0, B = 0, C = 0;
- for (int i = 0; i < pairs.size(); i++) {
- f1 += F(a[pairs[i].first][0], a[pairs[i].first][1], a[pairs[i].first][2], b[pairs[i].second][0] + x, b[pairs[i].second][1] + y, b[pairs[i].second][2] + z, A, B, C);
- }
- if (FtCnst == 0) {
- FtCnst = 0.0001;
- }
- if (f1 == 0) {
- return;
- } else {
- //int count = 0;
- while (f1 - f2 < 0 || f1 - f2 > FtCnst) {
- f1 = 0; fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0;
- for (int i = 0; i < pairs.size(); i++) {
- searchF0xyzabc(f1, fx, fy, fz, fa, fb, fc, a[pairs[i].first][0], a[pairs[i].first][1], a[pairs[i].first][2], b[pairs[i].second][0] + x, b[pairs[i].second][1] + y, b[pairs[i].second][2] + z, A, B, C);
- }
- while (true) {
- f2 = 0;
- for (int i = 0; i < pairs.size(); i++) {
- f2 += F(a[pairs[i].first][0], a[pairs[i].first][1], a[pairs[i].first][2], b[pairs[i].second][0] + (x - l * fx), b[pairs[i].second][1] + (y - l * fy), b[pairs[i].second][2] + (z - l * fz), A - l * fa, B - l * fb, C - l * fc);
- }
- //count++;
- if (f2 < f1) {
- x -= l * fx; y -= l * fy; z -= l * fz; A -= l * fa; B -= l * fb; C -= l * fc;
- fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0;
- break;
- } else {
- l *= 0.99; // may be put somewhere l > 0
- }
- }
- }
- //std::cout << " " << count << "\n";
- ApplyFit(b, x, y, z, A, B, C);
- }
-}
-
-void printPDBtraj(const char* Fname, std::vector< std::vector < RVec > > trj, std::vector< std::vector< std::pair< int, int > > > prs, std::vector< int > ndx) {
- FILE* a = fopen (Fname, "a");
- for (int i = 1; i < trj.size(); i++) {
- fprintf(a, "TITLE test t= %9.5f step= %5d\nMODEL %d\n", (/*float*/double)(i + 1), i + 1, i + 1);
- for (int j = 0; j < prs.size(); j++) {
- for (int k = 0; k < prs[j].size(); k++) {
- fprintf(a, "ATOM %5d CA CA %5d %8.3f%8.3f%8.3f%6.2f%6.2f %4s%2s \n", ndx[prs[j][k].first], ndx[prs[j][k].first], (trj[i][prs[j][k].first][0] * 10), (trj[i][prs[j][k].first][1] * 10), (trj[i][prs[j][k].first][2] * 10), 1.0, 20.0, " ", " ");
- }
- }
- fprintf(a, "TER\nENDMDL\n");
- }
- //fprintf(a, "END\n");
- fclose(a);
-}
-
-bool myfunction (const std::vector< std::pair< int, int > > i, const std::vector< std::pair< int, int > > j) {
- return i.size() < j.size();
-}
-
-void domain_reading(SelectionList sList, std::vector< int > ind, std::vector< std::vector < std::pair< int, int > > > &d_ind_num) {
- d_ind_num.resize(sList.size() - 1);
- for (int i = 1; i < sList.size(); i++) {
- d_ind_num[i - 1].resize(0);
- for (ArrayRef< const int >::iterator ai = sList[i].atomIndices().begin(); (ai < sList[i].atomIndices().end()); ai++) {
- d_ind_num[i - 1].push_back(std::make_pair(*ai, 0));
- }
- }
-
- for (int i = 0; i < d_ind_num.size(); i++) {
- int k;
- for (int j = 0; j < d_ind_num[i].size(); j++) {
- k = 0;
- while (ind[k] != d_ind_num[i][j].first) {
- k++;
- }
- d_ind_num[i][j].second = k;
- }
- }
-
- // its a point of interest as we can either take biggest domains of the given set and work from those or we take them in given order
- //std::sort(d_ind_num.begin(), d_ind_num.end(), myfunction);
-
- for (int i = 0; i < d_ind_num.size(); i++) {
- if (d_ind_num[i].size() >= 4) {
- for (int k = 0; k < d_ind_num[i].size(); k++) {
- for (int j = i + 1; j < d_ind_num.size(); j++) {
- for (int x = 0; x < d_ind_num[j].size(); x++) {
- if (d_ind_num[j][x].first == d_ind_num[i][k].first) {
- d_ind_num[j].erase(d_ind_num[j].begin() + x);
- x--;
- }
- }
- }
- }
- }
- }
-}
-
-void domain_expansion(std::vector< int > indx, std::vector< std::vector < std::pair< int, int > > > &d_ind_num, std::vector< RVec > frame) {
- std::vector< bool > flag;
- flag.resize(indx.size(), 1);
- for (int i = 0; i < d_ind_num.size(); i++) {
- for (int j = 0; j < d_ind_num[i].size(); j++) {
- flag[d_ind_num[i][j].second] = 0;
- }
- }
- /*float*/double dist;
- int x3 = -1;
- for (int i = 0; i < indx.size(); i++) {
- if (flag[i]) {
- dist = 0;
- for (int x1 = 0; x1 < d_ind_num.size(); x1++) {
- for (int x2 = 0; x2 < d_ind_num[x1].size(); x2++) {
- if (dist == 0) {
- dist = (frame[i] - frame[d_ind_num[0][0].first]).norm();
- x3 = 0;
- }
- if ((frame[i] - frame[d_ind_num[x1][x2].first]).norm() < dist) {
- x3 = x1;
- dist = (frame[i] - frame[d_ind_num[x1][x2].first]).norm();
- }
- }
- }
- d_ind_num[x3].push_back(std::make_pair(indx[i], i));
- }
- }
-}
-
-void MakeDomainFitPairs(std::vector< std::vector< std::pair< int, int > > > &p, std::vector< std::vector< std::pair< int, int > > > d) {
- p.resize(d.size());
- for (int i = 0; i < d.size(); i++) {
- p[i].resize(0);
- for (int j = 0; j < d[i].size(); j++) {
- p[i].push_back(std::make_pair(d[i][j].second, d[i][j].second));
- }
- }
-}
-
-void TrajFitting(std::vector< std::vector < RVec > > &trj, double FC, std::vector< std::vector< std::pair< int, int > > > prs) {
-
- // its needed to be thought through on whats shared and whats private
- //#pragma omp parallel for shared(trj, prs, FC)/* private(clone)*/
- for (int i = 1; i < trj.size()/*2*/; i++) {
- std::vector< std::vector < RVec > > clone;
- clone.resize(0);
- clone.resize(prs.size(), trj[i]);
- //std::chrono::time_point<std::chrono::system_clock> now1, now2;
- for (int j = 0; j < prs.size(); j++) {
- //now1 = std::chrono::system_clock::now();
- MyFitNew(trj[0], clone[j], prs[j], FC);
- //now2 = std::chrono::system_clock::now();
- //std::cout << " " << std::chrono::duration_cast<std::chrono::milliseconds>(now2 - now1).count() << "\n";
- for (int k = 0; k < prs[j].size(); k++) {
- trj[i][prs[j][k].first] = clone[j][prs[j][k].first];
- }
- }
- clone.resize(0);
- //std::cout << "frame " << i << " fitted\n";
- }
-}
+#include <newfit.h>
class Fitng : public TrajectoryAnalysisModule
{
virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
const TopologyInformation &top);
- virtual void initAfterFirstFrame(const TrajectoryAnalysisSettings &settings,
- const t_trxframe &fr);
-
//! Call for each frame of the trajectory
// virtual void analyzeFrame(const t_trxframe &fr, t_pbc *pbc);
virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
virtual void writeOutput();
private:
- SelectionList sel_;
-
- std::vector< int > index;
- const int *ind;
- std::vector< std::vector< std::pair< int, int > > > domains_index_number;
- std::vector< std::vector < RVec > > trajectory;
- std::vector< RVec > reference;
-
- std::vector< t_trxframe > saved_traj;
- //Selection selec;
- std::vector< std::vector< std::pair< int, int > > > pairs;
- int frames = 0;
- int method = 1;
- double FitConst = 0;
- std::string OutPutTrjName;
- t_trxstatus *op;
+ Selection sel_;
+ std::vector< int > index;
+ std::vector < RVec > trajectoryFrame;
+ std::vector< RVec > reference;
+ std::vector< std::pair< size_t, size_t > > fitPairs;
+ double fitConst {0.000'001};
+ std::string outputTrjName;
+ t_trxframe tempFrame;
+ t_trxstatus *op;
};
Fitng::Fitng(): TrajectoryAnalysisModule()
}
void
-Fitng::initOptions(IOptionsContainer *options,
- TrajectoryAnalysisSettings *settings)
+Fitng::initOptions( IOptionsContainer *options,
+ TrajectoryAnalysisSettings *settings)
{
static const char *const desc[] = {
"[THISMODULE] to be done"
};
// Add the descriptive text (program help text) to the options
settings->setHelpText(desc);
- // Add option for selecting a subset of atoms
- //options->addOption(SelectionOption("select")
- // .store(&selec).required()
- // .description("Atoms that are considered as part of the excluded volume"));
// Add option for selection list
- options->addOption(SelectionOption("select_residue_and_domains").storeVector(&sel_)
- .required().dynamicMask().multiValue()
- .description("Selected group and domains based on it"));
+ options->addOption(SelectionOption("select")
+ .store(&sel_).required()
+ .description("Atoms that are considered as part of the excluded volume"));
// Add option for Fit constant
- options->addOption(DoubleOption("FitConst")
- .store(&FitConst)
- .description("Fitting untill diff <= FitConst, by default == 0.00001"));
+ options->addOption(DoubleOption("fitConst")
+ .store(&fitConst)
+ .description("Fitting untill diff <= FitConst, by default == 0.000'001"));
// Add option for output pdb traj (meh edition)
- options->addOption(StringOption("trj_out")
- .store(&OutPutTrjName)
- .description("transformed trajectory"));
- // Add option for trajectory Fit transformation
- options->addOption(IntegerOption("method")
- .store(&method)
- .description("1 - fitting all -> all | 2 - keep only domains, all -> all | 3 - keep only domains, domain -> domain | 4 - extend domains to full coverage, domain -> domain"));
+ options->addOption(StringOption("out")
+ .store(&outputTrjName)
+ .description("transformed trajectory"));
// Control input settings
settings->setFlags(TrajectoryAnalysisSettings::efNoUserPBC);
settings->setFlag(TrajectoryAnalysisSettings::efUseTopX);
settings->setPBC(true);
}
-// /home/toluk/Develop/samples/reca_test_ground/
-// -s '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.xtc' -sf '/home/toluk/Develop/samples/reca_test_ground/SL0' -n '/home/toluk/Develop/samples/reca_test_ground/t0.ndx' -pdb_out '/home/toluk/Develop/samples/reca_test_ground/FullCAFit.pdb' -e 1000
-// -s '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.xtc' -sf '/home/toluk/Develop/samples/reca_test_ground/SL2' -n '/home/toluk/Develop/samples/reca_test_ground/t2.ndx' -pdb_out '/home/toluk/Develop/samples/reca_test_ground/CAdomainsFit.pdb' -e 1000 -method 2
-// -s '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.xtc' -sf '/home/toluk/Develop/samples/reca_test_ground/SL2' -n '/home/toluk/Develop/samples/reca_test_ground/t2.ndx' -pdb_out '/home/toluk/Develop/samples/reca_test_ground/CAdomainsOnlyFit.pdb' -e 1000 -method 3
-// -s '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.xtc' -sf '/home/toluk/Develop/samples/reca_test_ground/SL2' -n '/home/toluk/Develop/samples/reca_test_ground/t2.ndx' -pdb_out '/home/toluk/Develop/samples/reca_test_ground/CAdomainsWExtraOnlyFit.pdb' -e 1000 -method 4
-
-// -s '/home/titov_ai/Develop/samples/reca_rd/reca_rd.mono.tpr' -f '/home/titov_ai/Develop/samples/reca_rd/reca_rd.mono.xtc' -sf '/home/titov_ai/Develop/samples/reca_rd/SLCa' -n '/home/titov_ai/Develop/samples/reca_rd/CA.ndx' -pdb_out '/home/titov_ai/Develop/samples/reca_rd/test.xtc' -method 1
-
-// -s '/home/titov_ai/BioStore/titov_ai/trp_cage/kv1/trp-cage.pdb' -f '/home/titov_ai/BioStore/titov_ai/trp_cage/kv1/test.xtc' -sf '/home/titov_ai/BioStore/titov_ai/trp_cage/kv1/SLfull' -n '/home/titov_ai/BioStore/titov_ai/trp_cage/kv1/full.ndx' -pdb_out '/home/titov_ai/BioStore/titov_ai/trp_cage/kv1/testfit.xtc' -method 1
-
-// -f '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/local_test_trjconvfit_time.xtc' -s '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/trp-cage.pdb' -sf '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/SLfull' -n '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/full.ndx' -pdb_out '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/local_test_prefit_time.xtc' -method 1 -e 500 -FitConst 0.000001
-// -s '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/trp-cage.pdb' -f '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/test_Done.xtc' -sf '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/SLfull' -n '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/full.ndx' -pdb_out '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/testfit.xtc' -method 1
-
-// -s '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/trp-cage.md_npt.tpr' -f '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/test_Done.xtc' -sf '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/SLfull' -n '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/full.ndx' -trj_out '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/testfit.xtc' -method 1 -h
-
-
-// -s '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/trp-cage.md_npt.tpr' -f '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/test1k.xtc' -sf '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/SLfull' -n '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/full.ndx' -trj_out '/home/titov_ai/BioStore/titov_ai/trp_cage/new_try/kv1/local_test/testfit.xtc' -method 1 -h
-
-
void
-Fitng::initAnalysis(const TrajectoryAnalysisSettings &settings, const TopologyInformation &top)
+Fitng::initAnalysis(const TrajectoryAnalysisSettings &settings,
+ const TopologyInformation &top)
{
index.resize(0);
- for (ArrayRef< const int >::iterator ai = sel_[0].atomIndices().begin(); (ai < sel_[0].atomIndices().end()); ai++) {
+ for (ArrayRef< const int >::iterator ai {sel_.atomIndices().begin()}; ai < sel_.atomIndices().end(); ++ai) {
index.push_back(*ai);
}
reference.resize(0);
if (top.hasFullTopology()) {
- for (int i = 0; i < index.size(); i++) {
+ for (size_t i {0}; i < index.size(); ++i) {
reference.push_back(top.x().at(index[i]));
}
}
- //std::cout << reference.size() << "\n important \n";
- saved_traj.resize(0);
+ op = open_trx(outputTrjName.c_str(), "w+");
+ trajectoryFrame.reserve(index.size());
- if (sel_.size() > 1 && method >= 2) {
- domain_reading(sel_, index, domains_index_number);
+ fitPairs.resize(index.size());
+ for (size_t i {0}; i < index.size(); ++i) {
+ fitPairs[i] = std::make_pair(i, i);
}
}
void
-Fitng::initAfterFirstFrame(const TrajectoryAnalysisSettings &settings, const t_trxframe &fr)
+Fitng::analyzeFrame(int frnr,
+ const t_trxframe &fr,
+ t_pbc *pbc,
+ TrajectoryAnalysisModuleData *pdata)
{
- trajectory.resize(0);
- std::vector < RVec > a;
- a.resize(0);
- trajectory.push_back(reference);
- /*for (int i = 0; i < index.size(); i++) {
- trajectory[0].push_back(fr.x[index[i]]);
- }*/
- trajectory.push_back(a);
- switch (method) {
- case 1:
- pairs.resize(1);
- for (int i = 0; i < index.size(); i++) {
- pairs[0].push_back(std::make_pair(i, i));
- }
- break;
- case 2:
- pairs.resize(1);
- for (int i = 0; i < domains_index_number.size(); i++) {
- for (int j = 0; j < domains_index_number[i].size(); j++) {
- pairs[0].push_back(std::make_pair(domains_index_number[i][j].second, domains_index_number[i][j].second));
- }
- }
- break;
- case 3:
- MakeDomainFitPairs(pairs, domains_index_number);
- break;
- case 4:
- domain_expansion(index, domains_index_number, trajectory[0]);
- MakeDomainFitPairs(pairs, domains_index_number);
- break;
+ trajectoryFrame.resize(0);
+ for (size_t i {0}; i < index.size(); ++i) {
+ trajectoryFrame.push_back(fr.x[index[i]]);
}
- op = open_trx(OutPutTrjName.c_str(), "w+");
-}
+ MyFitNew(reference, trajectoryFrame, fitPairs, fitConst);
-void
-Fitng::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
- TrajectoryAnalysisModuleData *pdata)
-{
- trajectory[1].resize(0);
- for (int i = 0; i < index.size(); i++) {
- trajectory[1].push_back(fr.x[index[i]]);
+ tempFrame = fr;
+ for (size_t i {0}; i < index.size(); ++i) {
+ copy_rvec(trajectoryFrame[i].as_vec(), tempFrame.x[index[i]]);
}
- TrajFitting(trajectory, FitConst, pairs);
- t_trxframe c = fr;
- for (int i = 0; i < index.size(); i++) {
- copy_rvec(trajectory[1][i].as_vec(), c.x[index[i]]);
- }
- //trajectory.push_back(trajectory[1]);
- if (write_trxframe_indexed(op, static_cast<const t_trxframe *>(&c), index.size(), static_cast<const int *>(&index.front())/*saved_traj[i].index*/, NULL) != 0) {
- std::cout << "\nFileOutPut error\nframe number = " << frnr << "\n";
+ if (write_trxframe_indexed(op, static_cast<const t_trxframe *>(&tempFrame), index.size(), static_cast<const int *>(&index.front()), NULL) != 0) {
+ std::cout << "\nFileOutPut error | frame number = " << frnr << "\n";
}
-
- //printPDBtraj((OutPutTrjName + ".pdb").c_str(), trajectory, pairs, index);
}
-
void
Fitng::finishAnalysis(int nframes)
{
- std::cout << "\nGame Over\n";
+ std::cout << "\nGAME OVER\n";
}
void
Fitng::writeOutput()
{
- //printPDBtraj((OutPutTrjName + ".pdb").c_str(), trajectory, pairs, index);
+
}
+++ /dev/null
-#ifndef FITNG_H
-#define FITNG_H
-
-#endif // FITNG_H
+++ /dev/null
-#include <gromacs/trajectoryanalysis.h>
-
-#include <vector>
-#include <math.h>
-
-using namespace gmx;
-using gmx::RVec;
-
-double F (double aix, double aiy, double aiz, double bix_plus_x, double biy_plus_y, double biz_plus_z, double A, double B, double C) {
- return sqrt( pow(aix + biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * bix_plus_x, 2) +
- pow(aiy - biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * bix_plus_x, 2) +
- pow(aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y, 2) );
-}
-
-void searchF0xyzabc (double &F, double &Fx, double &Fy, double &Fz, double &Fa, double &Fb, double &Fc, double aix, double aiy, double aiz, double bix_plus_x, double biy_plus_y, double biz_plus_z, double A, double B, double C) {
- double t01, t02, t03, t04, t05, t06, t07;
- t01 = pow(aix + biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * bix_plus_x, 2);
- t02 = pow(aiy - biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * bix_plus_x, 2);
- t03 = pow(aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y, 2);
- t04 = sqrt(t01 + t02 + t03);
- t05 = (aix + biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * bix_plus_x);
- t06 = (aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y);
- t07 = (aiy - biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * bix_plus_x);
-
- //#pragma omp parallel sections
- //{
- //#pragma omp section
- F += t04;
- //#pragma omp section
- Fx += -(2 * cos(B) * cos(C) * t05 - 2 * sin(B) * t06 + 2 * cos(B) * sin(C) * t07) / (2 * t04);
- //#pragma omp section
- Fy += -(2 * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) * t07 - 2 * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) * t05 + 2 * cos(B) * sin(A) * t06) / (2 * t04);
- //#pragma omp section
- Fz += -(2 * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) * t05 - 2 * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) * t07 + 2 * cos(A) * cos(B) * t06) / (2 * t04);
- //#pragma omp section
- Fa += -(2 * (cos(A) * cos(B) * biy_plus_y - cos(B) * sin(A) * biz_plus_z) * t06 -
- 2 * (biy_plus_y * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) + biz_plus_z * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C))) * t07 +
- 2 * (biy_plus_y * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) + biz_plus_z * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B))) * t05) / (2 * t04);
- //#pragma omp section
- Fb += -(2 * (cos(A) * cos(B) * sin(C) * biz_plus_z - sin(B) * sin(C) * bix_plus_x + cos(B) * sin(A) * sin(C) * biy_plus_y) * t07 +
- 2 * (cos(A) * cos(B) * cos(C) * biz_plus_z - cos(C) * sin(B) * bix_plus_x + cos(B) * cos(C) * sin(A) * biy_plus_y) * t05 -
- 2 * (cos(B) * bix_plus_x + sin(A) * sin(B) * biy_plus_y + cos(A) * sin(B) * biz_plus_z) * t06) / (2 * t04);
- //#pragma omp section
- Fc += (2 * (biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) - biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) + cos(B) * sin(C) * bix_plus_x) * t05 -
- 2 * (biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) + cos(B) * cos(C) * bix_plus_x) * t07) / (2 * t04);
- //}
-}
-
-void searchL (double &fLd, double &fLdd, double aix, double aiy, double aiz, double bix, double biy, double biz, double x, double y, double z, double A, double B, double C, double L, double fx, double fy, double fz, double fa, double fb, double fc) {
- double AmLfa, BmLfb, CmLfc, BIXpXmLfx, BIYpYmLfy, BIZpZmLfz, t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11, t12, t13, t14, t15, t16, t17, t18, t19, t20;
- AmLfa = A - L * fa;
- BmLfb = B - L * fb;
- CmLfc = C - L * fc;
- BIXpXmLfx = bix + x - L * fx;
- BIYpYmLfy = biy + y - L * fy;
- BIZpZmLfz = biz + z - L * fz;
- t1 = (cos(AmLfa) * cos(CmLfc) + sin(AmLfa) * sin(BmLfb) * sin(CmLfc));
- t2 = (cos(CmLfc) * sin(AmLfa) - cos(AmLfa) * sin(BmLfb) * sin(CmLfc));
- t3 = (cos(AmLfa) * sin(CmLfc) - cos(CmLfc) * sin(AmLfa) * sin(BmLfb));
- t4 = (sin(AmLfa) * sin(CmLfc) + cos(AmLfa) * cos(CmLfc) * sin(BmLfb));
- t5 = pow(aiz + sin(BmLfb) * BIXpXmLfx - cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - cos(BmLfb) * sin(AmLfa) * BIYpYmLfy, 2);
- t6 = pow(aix + t3 * BIYpYmLfy - t4 * BIZpZmLfz - cos(BmLfb) * cos(CmLfc) * BIXpXmLfx, 2);
- t7 = pow(aiy - t1 * BIYpYmLfy + t2 * BIZpZmLfz - cos(BmLfb) * sin(CmLfc) * BIXpXmLfx, 2);
-
- t8 = fc * sin(AmLfa) * sin(CmLfc);
- t9 = fa * cos(AmLfa) * cos(CmLfc);
- t10 = fb * cos(AmLfa) * cos(BmLfb) * sin(CmLfc);
- t11 = fc * cos(AmLfa) * cos(CmLfc) * sin(BmLfb);
- t12 = fa * sin(AmLfa) * sin(BmLfb) * sin(CmLfc);
-
- t13 = fa * cos(AmLfa) * sin(BmLfb) * sin(CmLfc);
- t14 = fc * cos(AmLfa) * sin(CmLfc);
- t15 = fa * cos(CmLfc) * sin(AmLfa);
- t16 = fb * cos(BmLfb) * sin(AmLfa) * sin(CmLfc);
- t17 = fc * cos(CmLfc) * sin(AmLfa) * sin(BmLfb);
-
- t18 = fx * cos(BmLfb) * sin(CmLfc);
- t19 = fc * cos(BmLfb) * cos(CmLfc);
- t20 = fb * sin(BmLfb) * sin(CmLfc);
-
- fLd += (2 * (aiy - t1 * BIYpYmLfy + t2 * BIZpZmLfz - cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
- ( BIZpZmLfz * (t8 - t9 + t10 + t11 - t12) + fy * t1 - fz * t2 + BIYpYmLfy * (t13 - t14 - t15 + t16 + t17) + t18 + t19 * BIXpXmLfx - t20 * BIXpXmLfx) -
- 2 * (aiz + sin(BmLfb) * BIXpXmLfx - cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - cos(BmLfb) * sin(AmLfa) * BIYpYmLfy) *
- (fx * sin(BmLfb) - fy * cos(BmLfb) * sin(AmLfa) + fb * cos(BmLfb) * BIXpXmLfx - fz * cos(AmLfa) * cos(BmLfb) -
- fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy + fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy) +
- 2 * (aix + t3 * BIYpYmLfy - t4 * BIZpZmLfz - cos(BmLfb) * cos(CmLfc) * BIXpXmLfx) *
- (BIZpZmLfz * (fa * cos(AmLfa) * sin(CmLfc) + fc * cos(CmLfc) * sin(AmLfa) + fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) - fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) +
- BIYpYmLfy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) + t9 * sin(BmLfb) +
- fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) - fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) - fy * t3 +
- fz * t4 + fx * cos(BmLfb) * cos(CmLfc) - fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx)) /
- (2 * sqrt( t5 + t6 + pow(aiy - t1 * BIYpYmLfy + t2 * BIZpZmLfz - cos(BmLfb) * sin(CmLfc) * BIXpXmLfx, 2)));
-
-
- fLdd += ( (2 * aiz + 2 * sin(BmLfb) * BIXpXmLfx - 2 * cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - 2 * cos(BmLfb) * sin(AmLfa) * BIYpYmLfy) *
- (2 * fb * fx * cos(BmLfb) - pow(fb, 2) * sin(BmLfb) * BIXpXmLfx - 2 * fa * fy * cos(AmLfa) * cos(BmLfb) + 2 * fa * fz * cos(BmLfb) * sin(AmLfa) + 2 * fb * fz * cos(AmLfa) * sin(BmLfb) +
- 2 * fb * fy * sin(AmLfa) * sin(BmLfb) + pow(fa, 2) * cos(AmLfa) * cos(BmLfb) * BIZpZmLfz + pow(fb, 2) * cos(AmLfa) * cos(BmLfb) * BIZpZmLfz +
- pow(fa, 2) * cos(BmLfb) * sin(AmLfa) * BIYpYmLfy + pow(fb, 2) * cos(BmLfb) * sin(AmLfa) * BIYpYmLfy + 2 * fa * fb * cos(AmLfa) * sin(BmLfb) * BIYpYmLfy -
- 2 * fa * fb * sin(AmLfa) * sin(BmLfb) * BIZpZmLfz) +
- (BIZpZmLfz * (fa * cos(AmLfa) * sin(CmLfc) + fc * cos(CmLfc) * sin(AmLfa) +
- fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) - fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) +
- BIYpYmLfy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) + t9 * sin(BmLfb) +
- fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) - fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) - fy * (cos(AmLfa) * sin(CmLfc) -
- cos(CmLfc) * sin(AmLfa) * sin(BmLfb)) + fz * t4 + fx * cos(BmLfb) * cos(CmLfc) -
- fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
- (2 * BIZpZmLfz * (fa * cos(AmLfa) * sin(CmLfc) +
- fc * cos(CmLfc) * sin(AmLfa) + fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) -
- fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) + 2 * BIYpYmLfy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) +
- t9 * sin(BmLfb) + fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) - fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) -
- 2 * fy * t3 + 2 * fz * t4 +
- 2 * fx * cos(BmLfb) * cos(CmLfc) - 2 * fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - 2 * fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) +
- (2 * aix + 2 * t3 * BIYpYmLfy - 2 * (sin(AmLfa) * sin(CmLfc) +
- cos(AmLfa) * cos(CmLfc) * sin(BmLfb)) * BIZpZmLfz - 2 * cos(BmLfb) * cos(CmLfc) * BIXpXmLfx) *
- (BIZpZmLfz * (pow(fa, 2) * sin(AmLfa) * sin(CmLfc) +
- pow(fc, 2) * sin(AmLfa) * sin(CmLfc) - 2 * fa * fc * cos(AmLfa) * cos(CmLfc) + pow(fa, 2) * cos(AmLfa) * cos(CmLfc) * sin(BmLfb) +
- pow(fb, 2) * cos(AmLfa) * cos(CmLfc) * sin(BmLfb) + pow(fc, 2) * cos(AmLfa) * cos(CmLfc) * sin(BmLfb) + 2 * fa * fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) +
- 2 * fb * fc * cos(AmLfa) * cos(BmLfb) * sin(CmLfc) - 2 * fa * fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) + BIYpYmLfy * (pow(fa, 2) * cos(CmLfc) * sin(AmLfa) * sin(BmLfb) -
- pow(fc, 2) * cos(AmLfa) * sin(CmLfc) - 2 * fa * fc * cos(CmLfc) * sin(AmLfa) - pow(fa, 2) * cos(AmLfa) * sin(CmLfc) + pow(fb, 2) * cos(CmLfc) * sin(AmLfa) * sin(BmLfb) +
- pow(fc, 2) * cos(CmLfc) * sin(AmLfa) * sin(BmLfb) - 2 * fa * fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) + 2 * fa * fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc) +
- 2 * fb * fc * cos(BmLfb) * sin(AmLfa) * sin(CmLfc)) - 2 * fz * (fa * cos(AmLfa) * sin(CmLfc) + fc * cos(CmLfc) * sin(AmLfa) +
- fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) - fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) -
- 2 * fy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) + t9 * sin(BmLfb) + fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) -
- fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) + 2 * fb * fx * cos(CmLfc) * sin(BmLfb) + 2 * fc * t18 + pow(fb, 2) * cos(BmLfb) * cos(CmLfc) * BIXpXmLfx +
- pow(fc, 2) * cos(BmLfb) * cos(CmLfc) * BIXpXmLfx - 2 * fb * fc * sin(BmLfb) * sin(CmLfc) * BIXpXmLfx) +
- (BIZpZmLfz * (t8 - t9 + t10 + t11 - t12) + fy * t1 - fz * t2 + BIYpYmLfy * (t13 - t14 - t15 + t16 + t17) + t18 + t19 * BIXpXmLfx - t20 * BIXpXmLfx) *
- (2 * BIZpZmLfz * (t8 - t9 + t10 + t11 - t12) + 2 * fy * t1 - 2 * fz * t2 + 2 * BIYpYmLfy * (t13 - t14 - t15 + t16 + t17) + 2 * t18 + 2 * t19 * BIXpXmLfx - 2 * t20 * BIXpXmLfx) +
- (fx * sin(BmLfb) - fy * cos(BmLfb) * sin(AmLfa) + fb * cos(BmLfb) * BIXpXmLfx - fz * cos(AmLfa) * cos(BmLfb) - fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy +
- fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy) *
- (2 * fx * sin(BmLfb) - 2 * fy * cos(BmLfb) * sin(AmLfa) + 2 * fb * cos(BmLfb) * BIXpXmLfx - 2 * fz * cos(AmLfa) * cos(BmLfb) - 2 * fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy +
- 2 * fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + 2 * fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + 2 * fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy) +
- (2 * aiy - 2 * t1 * BIYpYmLfy + 2 * (cos(CmLfc) * sin(AmLfa) - cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) * BIZpZmLfz - 2 * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
- (BIZpZmLfz * (pow(fa, 2) * cos(AmLfa) * sin(BmLfb) * sin(CmLfc) -
- pow(fc, 2) * cos(CmLfc) * sin(AmLfa) - 2 * fa * t14 - pow(fa, 2) * cos(CmLfc) * sin(AmLfa) + pow(fb, 2) * cos(AmLfa) * sin(BmLfb) * sin(CmLfc) +
- pow(fc, 2) * cos(AmLfa) * sin(BmLfb) * sin(CmLfc) - 2 * fb * fc * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) + 2 * fa * t16 +
- 2 * fa * t17) + BIYpYmLfy * (pow(fa, 2) * cos(AmLfa) * cos(CmLfc) + pow(fc, 2) * cos(AmLfa) * cos(CmLfc) - 2 * fa * t8 +
- pow(fa, 2) * sin(AmLfa) * sin(BmLfb) * sin(CmLfc) + pow(fb, 2) * sin(AmLfa) * sin(BmLfb) * sin(CmLfc) + pow(fc, 2) * sin(AmLfa) * sin(BmLfb) * sin(CmLfc) -
- 2 * fa * t10 - 2 * fa * t11 - 2 * fb * t19 * sin(AmLfa)) - 2 * fz * (t8 - t9 + t10 + t11 - t12) - 2 * fy * (t13 - t14 - t15 + t16 + t17) -
- 2 * fc * fx * cos(BmLfb) * cos(CmLfc) + 2 * fb * fx * sin(BmLfb) * sin(CmLfc) +
- pow(fb, 2) * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx + pow(fc, 2) * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx + 2 * fb * fc * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx)) /
- (2 * sqrt( t5 + t6 + t7)) -
- ( (2 * (aiy - t1 * BIYpYmLfy + t2 * BIZpZmLfz - cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
- (BIZpZmLfz * (t8 - t9 + t10 + t11 - t12) + fy * t1 - fz * t2 + BIYpYmLfy * (t13 - t14 - t15 + t16 + t17) + t18 + t19 * BIXpXmLfx - t20 * BIXpXmLfx) -
- 2 * (aiz + sin(BmLfb) * BIXpXmLfx - cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - cos(BmLfb) * sin(AmLfa) * BIYpYmLfy) *
- (fx * sin(BmLfb) - fy * cos(BmLfb) * sin(AmLfa) + fb * cos(BmLfb) * BIXpXmLfx - fz * cos(AmLfa) * cos(BmLfb) - fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy +
- fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy) +
- 2 * (aix + t3 * BIYpYmLfy - t4 * BIZpZmLfz - cos(BmLfb) * cos(CmLfc) * BIXpXmLfx) *
- (BIZpZmLfz * (fa * cos(AmLfa) * sin(CmLfc) + fc * cos(CmLfc) * sin(AmLfa) + fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) - fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) +
- BIYpYmLfy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) + t9 * sin(BmLfb) + fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) - fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) -
- fy * t3 + fz * t4 + fx * cos(BmLfb) * cos(CmLfc) - fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx)) *
- ( (2 * aix + 2 * t3 * BIYpYmLfy - 2 * t4 * BIZpZmLfz - 2 * cos(BmLfb) * cos(CmLfc) * BIXpXmLfx) *
- (BIZpZmLfz * (fa * cos(AmLfa) * sin(CmLfc) + fc * cos(CmLfc) * sin(AmLfa) + fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) - fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) +
- BIYpYmLfy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) + t9 * sin(BmLfb) + fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) - fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) -
- fy * t3 + fz * t4 + fx * cos(BmLfb) * cos(CmLfc) - fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) +
- (2 * aiy - 2 * t1 * BIYpYmLfy + 2 * t2 * BIZpZmLfz - 2 * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
- (BIZpZmLfz * (t8 - t9 + t10 + t11 - t12) + fy * (cos(AmLfa) * cos(CmLfc) + sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) - fz * t2 + BIYpYmLfy * (t13 - t14 - t15 + t16 + t17) +
- t18 + t19 * BIXpXmLfx - t20 * BIXpXmLfx) - (2 * aiz + 2 * sin(BmLfb) * BIXpXmLfx - 2 * cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - 2 * cos(BmLfb) * sin(AmLfa) * BIYpYmLfy) *
- (fx * sin(BmLfb) - fy * cos(BmLfb) * sin(AmLfa) + fb * cos(BmLfb) * BIXpXmLfx - fz * cos(AmLfa) * cos(BmLfb) - fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy +
- fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy))) /
- (4 * pow ( t5 + t6 + t7, 3 / 2));
-
-}
-
-
-void ApplyFit (std::vector< RVec > &b, double x, double y, double z, double A, double B, double C) {
- RVec temp;
- //#pragma omp parallel
- //{
- //#pragma omp for schedule(dynamic)
- for (int i = 0; i < b.size(); i++) {
- temp = b[i];
- b[i][0] = (temp[2] + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - (temp[1] + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) + cos(B) * cos(C) * (temp[0] + x);
- b[i][1] = (temp[1] + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) - (temp[2] + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) + cos(B) * sin(C) * (temp[0] + x);
- b[i][2] = cos(A) * cos(B) * (temp[2] + z) - sin(B) * (temp[0] + x) + cos(B) * sin(A) * (temp[1] + y);
- }
- //}
-}
-
-void CalcMid (std::vector< RVec > a, std::vector< RVec > b, RVec &midA, RVec &midB, std::vector< std::pair< int, int > > pairs) {
- midA[0] = 0;
- midA[1] = 0;
- midA[2] = 0;
-
- midB[0] = 0;
- midB[1] = 0;
- midB[2] = 0;
-
- for (int i = 0; i < pairs.size(); i++) {
- rvec_inc(midA, a[pairs[i].first]);
- rvec_inc(midB, b[pairs[i].second]);
- }
- midA[0] /= pairs.size();
- midA[1] /= pairs.size();
- midA[2] /= pairs.size();
-
- midB[0] /= pairs.size();
- midB[1] /= pairs.size();
- midB[2] /= pairs.size();
-}
-
-void MyFitNew (std::vector< RVec > a, std::vector< RVec > &b, std::vector< std::pair< int, int > > pairs) {
- double f1 = 0, f2 = 0, fx = 0, fy = 0, fz = 0, fa = 0, fb = 0, fc = 0, l = 1;
- RVec ma, mb;
- CalcMid(a, b, ma, mb, pairs);
- rvec_dec(ma, mb);
- double x = ma[0], y = ma[1], z = ma[2], A = 0, B = 0, C = 0;
- for (int i = 0; i < pairs.size(); i++) {
- f1 += F(a[pairs[i].first][0], a[pairs[i].first][1], a[pairs[i].first][2], b[pairs[i].second][0] + x, b[pairs[i].second][1] + y, b[pairs[i].second][2] + z, A, B, C);
- }
- if (f1 == 0) {
- return;
- } else {
- double fLd, fLdd;
- //l = 0;
- int count = 0;
- while (f1 - f2 < 0 || f1 - f2 > 0.00001) {
- f1 = 0; fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0;
- for (int i = 0; i < pairs.size(); i++) {
- searchF0xyzabc(f1, fx, fy, fz, fa, fb, fc, a[pairs[i].first][0], a[pairs[i].first][1], a[pairs[i].first][2], b[pairs[i].second][0] + x, b[pairs[i].second][1] + y, b[pairs[i].second][2] + z, A, B, C);
- }
- /*do {
- fLd = 0;
- fLdd = 0;
- for (int i = 0; i < pairs.size(); i++) {
- searchL(fLd, fLdd, a[pairs[i].first][0], a[pairs[i].first][1], a[pairs[i].first][2], b[pairs[i].second][0], b[pairs[i].second][1], b[pairs[i].second][2], x, y, z, A, B, C, l, fx, fy, fz, fa, fb, fc);
- }
- l -= fLd / fLdd;
- } while (std::abs(fLd) > 1);*/
- //std::cout << " Ltempo= " << Ltempo << " iter= " << Lco << "\n";
- while (true) {
- f2 = 0;
- for (int i = 0; i < pairs.size(); i++) {
- f2 += F(a[pairs[i].first][0], a[pairs[i].first][1], a[pairs[i].first][2], b[pairs[i].second][0] + (x - l * fx), b[pairs[i].second][1] + (y - l * fy), b[pairs[i].second][2] + (z - l * fz), A - l * fa, B - l * fb, C - l * fc);
- }
- count++;
- if (f2 < f1) {
- x -= l * fx; y -= l * fy; z -= l * fz; A -= l * fa; B -= l * fb; C -= l * fc;
- fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0;
- break;
- } else {
- l *= 0.85;
- }
- }
- count++;
- }
- //std::cout << "iteration count: " << count << "\n";
- ApplyFit(b, x, y, z, A, B, C);
- }
-}
+++ /dev/null
-#ifndef FITTINGN_H
-#define FITTINGN_H
-
-// где-то косяк в Смаке, без этой строчки undefined reference to
-
-#include <type_traits>
-#include <iostream>
-#include <string>
-
-using namespace gmx;
-using gmx::RVec;
-
-
-//template< class t1, class = typename std::enable_if<std::is_floating_point< t1 >::value >::type >
-//double F (t1 aix, t1 aiy, t1 aiz, t1 bix_plus_x, t1 biy_plus_y, t1 biz_plus_z, t1 A, t1 B, t1 C);
-
-namespace ftn {
-
- double F (double aix, double aiy, double aiz, double bix_plus_x, double biy_plus_y, double biz_plus_z, double A, double B, double C);
-
- void searchF0xyzabc (double &F, double &Fx, double &Fy, double &Fz, double &Fa, double &Fb, double &Fc, double aix, double aiy, double aiz, double bix_plus_x, double biy_plus_y, double biz_plus_z, double A, double B, double C);
-
- void searchL (double &fLd, double &fLdd, double aix, double aiy, double aiz, double bix, double biy, double biz, double x, double y, double z, double A, double B, double C, double L, double fx, double fy, double fz, double fa, double fb, double fc);
-
- void ApplyFit (std::vector< RVec > &b, double x, double y, double z, double A, double B, double C);
-
- void CalcMid (std::vector< RVec > a, std::vector< RVec > b, RVec &midA, RVec &midB, std::vector< std::pair< int, int > > pairs);
-
- void MyFitNew (std::vector< RVec > a, std::vector< RVec > &b, std::vector< std::pair< int, int > > pairs);
-
-}
-
-#endif // FITTINGN_H
--- /dev/null
+#include "newfit.h"
+
+// вычисление модуля расстояния между двумя точками, с учётом геометрического преобразования
+double F (const DVec &ai, const DVec &biR, const DVec &angl) {
+ const double cosA {cos(angl[0])}, cosB {cos(angl[1])}, cosC {cos(angl[2])}, sinA {sin(angl[0])}, sinB {sin(angl[1])}, sinC {sin(angl[2])};
+ return sqrt( pow(ai[0] + biR[1] * (cosA * sinC - cosC * sinA * sinB) - biR[2] * (sinA * sinC + cosA * cosC * sinB) - cosB * cosC * biR[0], 2) +
+ pow(ai[1] - biR[1] * (cosA * cosC + sinA * sinB * sinC) + biR[2] * (cosC * sinA - cosA * sinB * sinC) - cosB * sinC * biR[0], 2) +
+ pow(ai[2] + sinB * biR[0] - cosA * cosB * biR[2] - cosB * sinA * biR[1], 2) );
+}
+
+template< typename T1, typename T2 >
+void checkIfNan(T1 &a, const T2 &b) {
+ if (!std::isnan(b)) {
+ a += b;
+ } else {
+ a += 0.;
+ }
+}
+
+// вычисление функции F и её частных производных
+void searchF0xyzabc (long double &F, double &Fx, double &Fy, double &Fz, double &Fa, double &Fb, double &Fc, const RVec &ai, const RVec &biR, const DVec &angl) {
+ const double cosA {cos(angl[0])}, cosB {cos(angl[1])}, cosC {cos(angl[2])}, sinA {sin(angl[0])}, sinB {sin(angl[1])}, sinC {sin(angl[2])};
+ double t01 {pow(ai[0] + biR[1] * (cosA * sinC - cosC * sinA * sinB) - biR[2] * (sinA * sinC + cosA * cosC * sinB) - cosB * cosC * biR[0], 2)};
+ double t02 {pow(ai[1] - biR[1] * (cosA * cosC + sinA * sinB * sinC) + biR[2] * (cosC * sinA - cosA * sinB * sinC) - cosB * sinC * biR[0], 2)};
+ double t03 {pow(ai[2] + sinB * biR[0] - cosA * cosB * biR[2] - cosB * sinA * biR[1], 2)};
+ double t04 {sqrt(t01 + t02 + t03)};
+ double t05 {(ai[0] + biR[1] * (cosA * sinC - cosC * sinA * sinB) - biR[2] * (sinA * sinC + cosA * cosC * sinB) - cosB * cosC * biR[0])};
+ double t06 {(ai[2] + sinB * biR[0] - cosA * cosB * biR[2] - cosB * sinA * biR[1])};
+ double t07 {(ai[1] - biR[1] * (cosA * cosC + sinA * sinB * sinC) + biR[2] * (cosC * sinA - cosA * sinB * sinC) - cosB * sinC * biR[0])};
+
+ checkIfNan(F, t04);
+ checkIfNan(Fx, -(2 * cosB * cosC * t05 - 2 * sinB * t06 + 2 * cosB * sinC * t07) / (2 * t04));
+ checkIfNan(Fy, -(2 * (cosA * cosC + sinA * sinB * sinC) * t07 - 2 * (cosA * sinC - cosC * sinA * sinB) * t05 + 2 * cosB * sinA * t06) / (2 * t04));
+ checkIfNan(Fz, -(2 * (sinA * sinC + cosA * cosC * sinB) * t05 - 2 * (cosC * sinA - cosA * sinB * sinC) * t07 + 2 * cosA * cosB * t06) / (2 * t04));
+ checkIfNan(Fa, -(2 * (cosA * cosB * biR[1] - cosB * sinA * biR[2]) * t06 -
+ 2 * (biR[1] * (cosC * sinA - cosA * sinB * sinC) + biR[2] * (cosA * cosC + sinA * sinB * sinC)) * t07 +
+ 2 * (biR[1] * (sinA * sinC + cosA * cosC * sinB) + biR[2] * (cosA * sinC - cosC * sinA * sinB)) * t05) / (2 * t04));
+ checkIfNan(Fb, -(2 * (cosA * cosB * sinC * biR[2] - sinB * sinC * biR[0] + cosB * sinA * sinC * biR[1]) * t07 +
+ 2 * (cosA * cosB * cosC * biR[2] - cosC * sinB * biR[0] + cosB * cosC * sinA * biR[1]) * t05 -
+ 2 * (cosB * biR[0] + sinA * sinB * biR[1] + cosA * sinB * biR[2]) * t06) / (2 * t04));
+ checkIfNan(Fc, (2 * (biR[1] * (cosA * cosC + sinA * sinB * sinC) - biR[2] * (cosC * sinA - cosA * sinB * sinC) + cosB * sinC * biR[0]) * t05 -
+ 2 * (biR[2] * (sinA * sinC + cosA * cosC * sinB) - biR[1] * (cosA * sinC - cosC * sinA * sinB) + cosB * cosC * biR[0]) * t07) / (2 * t04));
+}
+
+// применения геометрического преобразования: смещение на вектор (x, y, z) и повороты на градусы A, B, C относительно осей (x, y, z)
+//
+// А порядок поворотов какой?
+void ApplyFit (std::vector< RVec > &b, const DVec &R, const DVec &angl) {
+ DVec t(0., 0., 0.);
+ const double cosA {cos(angl[0])}, cosB {cos(angl[1])}, cosC {cos(angl[2])}, sinA {sin(angl[0])}, sinB {sin(angl[1])}, sinC {sin(angl[2])};
+ for (auto &i : b) {
+ t = i.toDVec();
+ i[0] = static_cast< float >((t[2] + R[2]) * (sinA * sinC + cosA * cosC * sinB) - (t[1] + R[1]) * (cosA * sinC - cosC * sinA * sinB) + cosB * cosC * (t[0] + R[0]));
+ i[1] = static_cast< float >((t[1] + R[1]) * (cosA * cosC + sinA * sinB * sinC) - (t[2] + R[2]) * (cosC * sinA - cosA * sinB * sinC) + cosB * sinC * (t[0] + R[0]));
+ i[2] = static_cast< float >(cosA * cosB * (t[2] + R[2]) - sinB * (t[0] + R[0]) + cosB * sinA * (t[1] + R[1]));
+ }
+}
+
+// подсчёт геометрических центров множеств a и b на основе пар pairs
+void CalcMid (const std::vector< RVec > &a, const std::vector< RVec > &b, DVec &midA, DVec &midB, const std::vector< std::pair< size_t, size_t > > &pairs) {
+ for (const auto &i : pairs) {
+ midA += a[i.first].toDVec();
+ midB += b[i.second].toDVec();
+ }
+ midA[0] /= pairs.size();
+ midA[1] /= pairs.size();
+ midA[2] /= pairs.size();
+
+ midB[0] /= pairs.size();
+ midB[1] /= pairs.size();
+ midB[2] /= pairs.size();
+}
+
+// функция фитирования фрейма b на фрейм a на основе пар pairs с "точностью" FtCnst
+void MyFitNew (const std::vector< RVec > &a, std::vector< RVec > &b, const std::vector< std::pair< size_t, size_t > > &pairs, long double FtCnst) {
+ double fx {0.}, fy {0.}, fz {0.}, fa {0.}, fb {0.}, fc {0.};
+ long double f1 {0.}, f2 {0.}, l {1};
+ DVec tempA(0., 0., 0.), tempB(0., 0., 0.), tempC(0., 0., 0.), tempD(0., 0., 0.);
+ CalcMid(a, b, tempA, tempB, pairs);
+ tempA -= tempB;
+ // поиск стартового отклонения множеств
+ for (const auto &i : pairs) {
+ f1 += F(a[i.first].toDVec(), b[i.second].toDVec() + tempA, tempC);
+ }
+ if (FtCnst == 0) {
+ FtCnst = 0.000'001; // experimental
+ }
+ if (f1 == 0) {
+ ApplyFit(b, tempA, tempC);
+ return;
+ } else {
+ // поиск оптимального геометрического преобразования методом градиентного спуска
+ while (f1 < f2 || f1 - f2 > FtCnst) {
+ f1 = 0.; f2 = 0.; fx = 0.; fy = 0.; fz = 0.; fa = 0.; fb = 0.; fc = 0.; l = 1.;
+ // поиск производных и отклонения при заданных параметрах
+ for (const auto &i : pairs) {
+ searchF0xyzabc(f1, fx, fy, fz, fa, fb, fc, a[i.first], b[i.second] + tempA.toRVec(), tempC);
+ }
+ while (true) {
+ f2 = 0;
+ // поиск отклонения при новых параметрах
+ tempB[0] = tempA[0] - l * fx;
+ tempB[1] = tempA[1] - l * fy;
+ tempB[2] = tempA[2] - l * fz;
+ tempD[0] = tempC[0] - l * fa;
+ tempD[1] = tempC[1] - l * fb;
+ tempD[2] = tempC[2] - l * fc;
+ for (const auto &i : pairs) {
+ f2 += F(a[i.first].toDVec(), b[i.second].toDVec() + tempB, tempD);
+ }
+ if (f2 < f1) {
+ tempA = tempB;
+ tempC = tempD;
+ break;
+ } else {
+ // по факту существует минимальное значение переменной типа double
+ // в случае его достижения дважды останавливаем цикл т.к. упёрлись в предел допустимой точности
+ // таким образом гарантируем выход из цикла
+ if (l == DBL_MIN || l == 0.) { //DBL_TRUE_MIN
+ f2 = f1;
+ break;
+ }
+ l *= 0.1;
+ }
+ }
+ }
+ ApplyFit(b, tempA, tempC);
+ }
+}
--- /dev/null
+#ifndef NEWFIT_H
+#define NEWFIT_H
+
+#include <gromacs/trajectoryanalysis.h>
+#include <gromacs/trajectoryanalysis/topologyinformation.h>
+//#include "/home/titov_ai/Develop/gromacs-original/src/gromacs/trajectoryanalysis/topologyinformation.h"
+#include <gromacs/math/vectypes.h>
+#include <gromacs/math/vec.h>
+
+#include <iostream>
+#include <cfloat>
+#include <omp.h>
+#include <vector>
+#include <cstdio>
+
+using gmx::RVec;
+using gmx::DVec;
+using namespace gmx;
+
+// вычисление модуля расстояния между двумя точками, с учётом геометрического преобразования
+double F (const DVec &ai, const DVec &bi_plusR, const DVec &angl);
+
+template< typename T1, typename T2 > void checkIfNan(T1 &a, const T2 &b);
+
+// вычисление функции F и её частных производных
+void searchF0xyzabc (long double &F, double &Fx, double &Fy, double &Fz, double &Fa, double &Fb, double &Fc, const RVec &ai, const RVec &biR, const DVec &angl);
+
+// применения геометрического преобразования: смещение на вектор (x, y, z) и повороты на градусы A, B, C относительно осей (x, y, z)
+void ApplyFit (std::vector< RVec > &b, const DVec &R, const DVec &angl);
+
+// подсчёт геометрических центров множеств a и b на основе пар pairs
+void CalcMid (const std::vector< RVec > &a, const std::vector< RVec > &b, DVec &midA, DVec &midB, const std::vector< std::pair< size_t, size_t > > &pairs);
+
+// функция фитирования фрейма b на фрейм a на основе пар pairs с "точностью" FtCnst
+void MyFitNew (const std::vector< RVec > &a, std::vector< RVec > &b, const std::vector< std::pair< size_t, size_t > > &pairs, long double FtCnst);
+
+#endif // NEWFIT_H
biod.pnpi.spb.ru / alexxy / gromacs-fitng.git / commitdiff