From 40a96e3a1c71c0c1c4da24f238155e6995be37f9 Mon Sep 17 00:00:00 2001
From: Gareth Tribello <gareth.tribello@gmail.com>
Date: Sun, 13 Nov 2016 19:11:10 +0000
Subject: [PATCH] Started adding stuff to calculate projections using PCA
---
regtest/dimred/rt-pca-2/plumed.dat | 3 +-
regtest/dimred/rt-pca/plumed.dat | 3 +-
src/dimred/DimensionalityReductionBase.cpp | 2 +-
src/dimred/DimensionalityReductionBase.h | 4 +-
src/dimred/OutputPCAProjections.cpp | 114 +++++++++++++++++++++
src/dimred/PCA.cpp | 94 ++++++++---------
src/dimred/PCA.h | 58 +++++++++++
src/dimred/ProjectNonLandmarkPoints.cpp | 28 +++--
8 files changed, 236 insertions(+), 70 deletions(-)
create mode 100644 src/dimred/OutputPCAProjections.cpp
create mode 100644 src/dimred/PCA.h
diff --git a/regtest/dimred/rt-pca-2/plumed.dat b/regtest/dimred/rt-pca-2/plumed.dat
index 7b8302858..0ec9d4852 100644
--- a/regtest/dimred/rt-pca-2/plumed.dat
+++ b/regtest/dimred/rt-pca-2/plumed.dat
@@ -2,4 +2,5 @@ c1: READ FILE=colvar.in VALUES=cv1
c2: READ FILE=colvar.in VALUES=cv2
ff: COLLECT_FRAMES ARG=c1,c2
-PCA USE_OUTPUT_DATA_FROM=ff METRIC=EUCLIDEAN NLOW_DIM=2 OFILE=pca-comp.pdb FMT=%8.4f
+pca: PCA USE_OUTPUT_DATA_FROM=ff METRIC=EUCLIDEAN NLOW_DIM=2
+OUTPUT_PCA_PROJECTION USE_OUTPUT_DATA_FROM=pca FILE=pca-comp.pdb FMT=%8.4f
diff --git a/regtest/dimred/rt-pca/plumed.dat b/regtest/dimred/rt-pca/plumed.dat
index ad9156983..40c098e4e 100644
--- a/regtest/dimred/rt-pca/plumed.dat
+++ b/regtest/dimred/rt-pca/plumed.dat
@@ -1,2 +1,3 @@
ff: COLLECT_FRAMES ATOMS=1-22 STRIDE=1
-PCA USE_OUTPUT_DATA_FROM=ff NLOW_DIM=2 OFILE=pca-comp.pdb METRIC=OPTIMAL
+pca: PCA USE_OUTPUT_DATA_FROM=ff NLOW_DIM=2 METRIC=OPTIMAL
+OUTPUT_PCA_PROJECTION USE_OUTPUT_DATA_FROM=pca FILE=pca-comp.pdb
diff --git a/src/dimred/DimensionalityReductionBase.cpp b/src/dimred/DimensionalityReductionBase.cpp
index 0aeea0384..5a6860910 100644
--- a/src/dimred/DimensionalityReductionBase.cpp
+++ b/src/dimred/DimensionalityReductionBase.cpp
@@ -40,7 +40,7 @@ dimredbase(NULL)
{
// Check that some dissimilarity information is available
if( my_input_data ){
- if( !dissimilaritiesWereSet() ) error("dissimilarities have not been calcualted in input actions");
+ if( getName()!="PCA" && !dissimilaritiesWereSet() ) error("dissimilarities have not been calcualted in input actions");
// Now we check if the input was a dimensionality reduction object
dimredbase = dynamic_cast<DimensionalityReductionBase*>( my_input_data );
}
diff --git a/src/dimred/DimensionalityReductionBase.h b/src/dimred/DimensionalityReductionBase.h
index 902bfd5da..61a44e7cc 100644
--- a/src/dimred/DimensionalityReductionBase.h
+++ b/src/dimred/DimensionalityReductionBase.h
@@ -45,11 +45,11 @@ public:
static void registerKeywords( Keywords& keys );
DimensionalityReductionBase( const ActionOptions& );
/// Get the ith data point (this returns the projection)
- void getProjection( const unsigned& idata, std::vector<double>& point, double& weight );
+ virtual void getProjection( const unsigned& idata, std::vector<double>& point, double& weight );
/// Get a reference configuration (this returns the projection)
analysis::DataCollectionObject& getStoredData( const unsigned& idata, const bool& calcdist );
/// Actually perform the analysis
- void performAnalysis();
+ virtual void performAnalysis();
/// Overwrite getArguments so we get arguments from underlying class
std::vector<Value*> getArgumentList();
/// Calculate the projections of points
diff --git a/src/dimred/OutputPCAProjections.cpp b/src/dimred/OutputPCAProjections.cpp
new file mode 100644
index 000000000..78cff11ee
--- /dev/null
+++ b/src/dimred/OutputPCAProjections.cpp
@@ -0,0 +1,114 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2012 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed-code.org for more information.
+
+ This file is part of plumed, version 2.0.
+
+ plumed 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 3 of the License, or
+ (at your option) any later version.
+
+ plumed 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 plumed. If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "analysis/AnalysisBase.h"
+#include "reference/ReferenceAtoms.h"
+#include "reference/ReferenceArguments.h"
+#include "core/ActionRegister.h"
+#include "core/PlumedMain.h"
+#include "core/ActionSet.h"
+#include "core/Atoms.h"
+#include "core/SetupMolInfo.h"
+#include "tools/PDB.h"
+#include "PCA.h"
+
+namespace PLMD {
+namespace dimred {
+
+//+PLUMEDOC ANALYSIS OUTPUT_PCA_PROJECTION
+/*
+This is used to output the projection calculated by principle component analysis
+
+\par Examples
+
+*/
+//+ENDPLUMEDOC
+
+class OutputPCAProjection : public analysis::AnalysisBase {
+private:
+ PDB mypdb;
+ PCA* mypca;
+ std::string fmt;
+ std::string filename;
+public:
+ static void registerKeywords( Keywords& keys );
+ OutputPCAProjection( const ActionOptions& );
+ void performTask( const unsigned& , const unsigned& , MultiValue& ) const { plumed_error(); }
+ void performAnalysis();
+};
+
+PLUMED_REGISTER_ACTION(OutputPCAProjection,"OUTPUT_PCA_PROJECTION")
+
+void OutputPCAProjection::registerKeywords( Keywords& keys ){
+ analysis::AnalysisBase::registerKeywords( keys );
+ keys.add("compulsory","FILE","the name of the file to output to");
+ keys.add("optional","FMT","the format to use in the output file");
+ keys.add("compulsory","STRIDE","0","the frequency with which to perform the required analysis and to output the data. The default value of 0 tells plumed to use all the data");
+}
+
+OutputPCAProjection::OutputPCAProjection( const ActionOptions& ao ):
+Action(ao),
+analysis::AnalysisBase(ao),
+fmt("%f")
+{
+ // Setup the PCA object
+ mypca = dynamic_cast<PCA*>( my_input_data );
+ if( !mypca ) error("input must be a PCA object");
+
+ // Get setup the pdb
+ mypdb.setAtomNumbers( my_input_data->getAtomIndexes() );
+ mypdb.setArgumentNames( (mypca->my_input_data)->getArgumentNames() );
+
+ // Find a moldata object
+ std::vector<SetupMolInfo*> moldat=plumed.getActionSet().select<SetupMolInfo*>();
+ if( moldat.empty() ) warning("PDB output files do not have atom types unless you use MOLDATA");
+
+ parse("FILE",filename); parse("FMT",fmt);
+ if( !getRestart() ){ OFile ofile; ofile.link(*this); ofile.setBackupString("analysis"); ofile.backupAllFiles(filename); }
+ log.printf(" printing data to file named %s \n",filename.c_str() );
+}
+
+void OutputPCAProjection::performAnalysis(){
+ // Find a moldata object
+ std::vector<SetupMolInfo*> moldat=plumed.getActionSet().select<SetupMolInfo*>();
+ if( moldat.size()>1 ) error("you should only have one MOLINFO action in your input file");
+ SetupMolInfo* mymoldat=NULL; if( moldat.size()==1 ) mymoldat=moldat[0];
+ // Output the embedding in plumed pdb format
+ OFile afile; afile.link(*this); afile.setBackupString("analysis");
+ mypdb.setAtomPositions( (mypca->myref)->getReferencePositions() );
+ for(unsigned j=0;j<mypca->getArguments().size();++j) mypdb.setArgumentValue( (mypca->getArguments()[j])->getName(), (mypca->myref)->getReferenceArgument(j) );
+ // And output the first frame
+ afile.open( filename.c_str() ); afile.printf("REMARK TYPE=%s \n", mypca->mtype.c_str() );
+ if( plumed.getAtoms().usingNaturalUnits() ) mypdb.print( 1.0, mymoldat, afile, fmt );
+ else mypdb.print( atoms.getUnits().getLength()/0.1, mymoldat, afile, fmt );
+ // And now ouput the eigenvectors
+ for(unsigned dim=0;dim<mypca->nlow;++dim){
+ afile.printf("REMARK TYPE=DIRECTION \n");
+ mypdb.setAtomPositions( mypca->directions[dim].getReferencePositions() );
+ for(unsigned j=0;j<mypca->getArguments().size();++j) mypdb.setArgumentValue( (mypca->getArguments()[j])->getName(), mypca->directions[dim].getReferenceArgument(j) );
+ if( plumed.getAtoms().usingNaturalUnits() ) mypdb.print( 1.0, mymoldat, afile, fmt );
+ else mypdb.print( atoms.getUnits().getLength()/0.1, mymoldat, afile, fmt );
+ }
+ afile.close();
+}
+
+}
+}
diff --git a/src/dimred/PCA.cpp b/src/dimred/PCA.cpp
index 0f30f9ce8..c67c464cb 100644
--- a/src/dimred/PCA.cpp
+++ b/src/dimred/PCA.cpp
@@ -19,11 +19,9 @@
You should have received a copy of the GNU Lesser General Public License
along with plumed. If not, see <http://www.gnu.org/licenses/>.
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-#include "analysis/AnalysisBase.h"
+#include "PCA.h"
#include "tools/Matrix.h"
-#include "reference/Direction.h"
#include "reference/MetricRegister.h"
-#include "reference/ReferenceConfiguration.h"
#include "reference/ReferenceValuePack.h"
#include "analysis/ReadAnalysisFrames.h"
#include "core/ActionRegister.h"
@@ -93,39 +91,18 @@ PCA ARG=d1,d2,d3,d4,d5,d6 METRIC=EUCLIDEAN STRIDE=5 RUN=1000 NLOW_DIM=2 REWEIGHT
namespace PLMD {
namespace dimred {
-class PCA : public analysis::AnalysisBase {
-private:
- unsigned ndim;
-/// The way we are measuring distances
- std::string mtype;
-/// The position of the reference configuration (the one we align to)
- PDB mypdb;
-/// The eigenvectors for the atomic displacements
- Matrix<Vector> atom_eigv;
-/// The eigenvectors for the displacements in argument space
- Matrix<double> arg_eigv;
- std::string ofilename, fmt;
-public:
- static void registerKeywords( Keywords& keys );
- explicit PCA(const ActionOptions&ao);
- void performAnalysis();
- void performTask( const unsigned& , const unsigned& , MultiValue& ) const { plumed_error(); }
-};
-
PLUMED_REGISTER_ACTION(PCA,"PCA")
void PCA::registerKeywords( Keywords& keys ){
- AnalysisBase::registerKeywords( keys ); keys.use("ARG"); keys.reset_style("ARG","optional");
+ DimensionalityReductionBase::registerKeywords( keys ); keys.use("ARG"); keys.reset_style("ARG","optional");
keys.add("compulsory","METRIC","EUCLIDEAN","the method that you are going to use to measure the distances between points");
keys.add("atoms","ATOMS","the list of atoms that you are going to use in the measure of distance that you are using");
- keys.add("compulsory","NLOW_DIM","number of PCA coordinates required");
- keys.add("compulsory","OFILE","the file on which to output the eigenvectors");
- keys.add("compulsory","FMT","%f","the format to use for output files");
}
PCA::PCA(const ActionOptions&ao):
Action(ao),
-AnalysisBase(ao)
+DimensionalityReductionBase(ao),
+myref(NULL)
{
// Get the input PDB file from the underlying ReadAnalysisFrames object
analysis::ReadAnalysisFrames* myframes = dynamic_cast<analysis::ReadAnalysisFrames*>( my_input_data );
@@ -172,25 +149,19 @@ AnalysisBase(ao)
mypdb.setArgumentNames( argnames ); requestArguments( myargs );
}
}
- parse("NLOW_DIM",ndim);
- if( mypdb.getPositions().size()>0 ) atom_eigv.resize( ndim, mypdb.getPositions().size() );
- if( mypdb.getArgumentNames().size()>0 ) arg_eigv.resize( ndim, mypdb.getArgumentNames().size() );
-
- // Read stuff for output file
- parse("OFILE",ofilename); parse("FMT",fmt);
- if( !getRestart() ){ OFile ofile; ofile.link(*this); ofile.setBackupString("analysis"); ofile.backupAllFiles(ofilename); }
+ if( nlow==0 ) error("dimensionality of output not set");
checkRead();
}
void PCA::performAnalysis(){
// Align everything to the first frame
- getStoredData( 0, false ).transferDataToPDB( mypdb ); ReferenceConfiguration* myconf0=metricRegister().create<ReferenceConfiguration>(mtype, mypdb);
+ my_input_data->getStoredData( 0, false ).transferDataToPDB( mypdb );
+ ReferenceConfiguration* myconf0=metricRegister().create<ReferenceConfiguration>(mtype, mypdb);
MultiValue myval( 1, myconf0->getNumberOfReferenceArguments() + 3*myconf0->getNumberOfReferencePositions() + 9 );
ReferenceValuePack mypack( myconf0->getNumberOfReferenceArguments(), myconf0->getNumberOfReferencePositions(), myval );
for(unsigned i=0;i<myconf0->getNumberOfReferencePositions();++i) mypack.setAtomIndex( i, i );
// Setup some PCA storage
myconf0->setupPCAStorage ( mypack );
-
// Create some arrays to store the average position
std::vector<double> sarg( myconf0->getNumberOfReferenceArguments(), 0 );
std::vector<Vector> spos( myconf0->getNumberOfReferencePositions() );
@@ -199,7 +170,7 @@ void PCA::performAnalysis(){
// Calculate the average displacement from the first frame
double norm=getWeight(0); std::vector<double> args( getNumberOfArguments() );
for(unsigned i=1;i<getNumberOfDataPoints();++i){
- getStoredData( i, false ).transferDataToPDB( mypdb );
+ my_input_data->getStoredData( i, false ).transferDataToPDB( mypdb );
for(unsigned j=0;j<getArguments().size();++j) mypdb.getArgumentValue( getArguments()[j]->getName(), args[j] );
double d = myconf0->calc( mypdb.getPositions(), getPbc(), getArguments(), args, mypack, true );
// Accumulate average displacement of arguments (Here PBC could do fucked up things - really needs Berry Phase ) GAT
@@ -216,9 +187,7 @@ void PCA::performAnalysis(){
unsigned narg=myconf0->getNumberOfReferenceArguments(), natoms=myconf0->getNumberOfReferencePositions();
Matrix<double> covar( narg+3*natoms, narg+3*natoms ); covar=0;
for(unsigned i=0;i<getNumberOfDataPoints();++i){
- // double d = data[i]->calc( spos, getPbc(), getArguments(), sarg, mypack, true );
- // ReferenceConfiguration* myconf = getReferenceConfiguration( i, false );
- getStoredData( i, false ).transferDataToPDB( mypdb );
+ my_input_data->getStoredData( i, false ).transferDataToPDB( mypdb );
for(unsigned j=0;j<getArguments().size();++j) mypdb.getArgumentValue( getArguments()[j]->getName(), args[j] );
double d = myconf0->calc( mypdb.getPositions(), getPbc(), getArguments(), args, mypack, true );
for(unsigned jarg=0;jarg<narg;++jarg){
@@ -252,32 +221,49 @@ void PCA::performAnalysis(){
Matrix<double> eigvec( narg+3*natoms, narg+3*natoms );
diagMat( covar, eigval, eigvec );
- // Open an output file
- OFile ofile; ofile.link(*this); ofile.setBackupString("analysis");
- ofile.open( ofilename );
// Output the reference configuration
mypdb.setAtomPositions( spos );
for(unsigned j=0;j<sarg.size();++j) mypdb.setArgumentValue( getArguments()[j]->getName(), sarg[j] );
- ofile.printf("REMARK TYPE=%s \n",mtype.c_str() );
- mypdb.print( atoms.getUnits().getLength()/0.1, NULL, ofile, fmt );
+ // Reset the reference configuration
+ if( myref ) delete myref;
+ myref = metricRegister().create<ReferenceConfiguration>( mtype, mypdb );
// Store and print the eigenvectors
std::vector<Vector> tmp_atoms( natoms );
- for(unsigned dim=0;dim<ndim;++dim){
+ for(unsigned dim=0;dim<nlow;++dim){
unsigned idim = covar.ncols() - 1 - dim;
- for(unsigned i=0;i<narg;++i){
- arg_eigv(dim,i)=eigvec(idim,i);
- mypdb.setArgumentValue( getArguments()[i]->getName(), eigvec(idim,i) );
- }
+ for(unsigned i=0;i<narg;++i) mypdb.setArgumentValue( getArguments()[i]->getName(), eigvec(idim,i) );
for(unsigned i=0;i<natoms;++i){
- for(unsigned k=0;k<3;++k) tmp_atoms[i][k]=atom_eigv(dim,i)[k]=eigvec(idim,narg+3*i+k);
+ for(unsigned k=0;k<3;++k) tmp_atoms[i][k]=eigvec(idim,narg+3*i+k);
}
mypdb.setAtomPositions( tmp_atoms );
- ofile.printf("REMARK TYPE=DIRECTION \n");
- mypdb.print( atoms.getUnits().getLength()/0.1, NULL, ofile, fmt );
+ // Create a direction object so that we can calculate other PCA components
+ directions.push_back( Direction(ReferenceConfigurationOptions("DIRECTION")));
+ directions[dim].read( mypdb );
}
// Close the output file
- delete myconf0; ofile.close();
+ delete myconf0;
+}
+
+void PCA::getProjection( const unsigned& idata, std::vector<double>& point, double& weight ){
+ if( point.size()!=nlow ) point.resize( nlow );
+ // Retrieve the weight
+ weight = getWeight(idata);
+ // Retrieve the data point
+ getProjection( my_input_data->getStoredData( idata, false ), point );
+}
+
+void PCA::getProjection( analysis::DataCollectionObject& myidata, std::vector<double>& point ){
+ myidata.transferDataToPDB( mypdb ); std::vector<double> args( getArguments().size() );
+ for(unsigned j=0;j<getArguments().size();++j) mypdb.getArgumentValue( getArguments()[j]->getName(), args[j] );
+ // Create some storage space
+ MultiValue myval( 1, 3*mypdb.getPositions().size() + args.size() + 9);
+ ReferenceValuePack mypack( args.size(), mypdb.getPositions().size(), myval );
+ for(unsigned i=0;i<mypdb.getPositions().size();++i) mypack.setAtomIndex( i, i );
+ myref->setupPCAStorage( mypack );
+ // And calculate
+ myref->calculate( mypdb.getPositions(), getPbc(), getArguments(), mypack, true );
+ for(unsigned i=0;i<nlow;++i) point[i]=myref->projectDisplacementOnVector( directions[i], mypdb.getPositions(), getArguments(), args, mypack );
}
}
diff --git a/src/dimred/PCA.h b/src/dimred/PCA.h
new file mode 100644
index 000000000..35dfd3f49
--- /dev/null
+++ b/src/dimred/PCA.h
@@ -0,0 +1,58 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2012-2014 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed-code.org for more information.
+
+ This file is part of plumed, version 2.
+
+ plumed 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 3 of the License, or
+ (at your option) any later version.
+
+ plumed 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 plumed. If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#ifndef __PLUMED_dimred_PCA_h
+#define __PLUMED_dimred_PCA_h
+
+#include "DimensionalityReductionBase.h"
+#include "reference/ReferenceConfiguration.h"
+#include "reference/Direction.h"
+#include "tools/PDB.h"
+
+namespace PLMD {
+namespace dimred {
+
+class PCA : public DimensionalityReductionBase {
+friend class OutputPCAProjection;
+private:
+/// The way we are measuring distances
+ std::string mtype;
+/// The position of the reference configuration (the one we align to)
+ PDB mypdb;
+/// The eigenvectors for the displacements in argument space
+ std::string ofilename, fmt;
+/// The eigenvectors that we are using
+ ReferenceConfiguration* myref;
+ std::vector<Direction> directions;
+public:
+ static void registerKeywords( Keywords& keys );
+ explicit PCA(const ActionOptions&ao);
+ void performAnalysis();
+ void getProjection( const unsigned& idata, std::vector<double>& point, double& weight );
+ void getProjection( analysis::DataCollectionObject& myidata, std::vector<double>& point );
+ void calculateProjections( const Matrix<double>& , Matrix<double>& ){ plumed_error(); }
+ void setTargetDistance( const unsigned& , const double& ){ plumed_error(); }
+ double calculateStress( const std::vector<double>& pp, std::vector<double>& der ){ plumed_error(); }
+};
+
+}
+}
+#endif
diff --git a/src/dimred/ProjectNonLandmarkPoints.cpp b/src/dimred/ProjectNonLandmarkPoints.cpp
index 47e9d1eb8..7bfa29f76 100644
--- a/src/dimred/ProjectNonLandmarkPoints.cpp
+++ b/src/dimred/ProjectNonLandmarkPoints.cpp
@@ -27,6 +27,7 @@
#include "analysis/AnalysisBase.h"
#include "reference/ReferenceConfiguration.h"
#include "DimensionalityReductionBase.h"
+#include "PCA.h"
//+PLUMEDOC DIMRED PROJECT_ALL_ANALYSIS_DATA
/*
@@ -101,18 +102,23 @@ std::vector<Value*> ProjectNonLandmarkPoints::getArgumentList(){
}
void ProjectNonLandmarkPoints::generateProjection( const unsigned& idat, std::vector<double>& point ){
- ConjugateGradient<ProjectNonLandmarkPoints> myminimiser( this );
- unsigned closest=0; double mindist = sqrt( getDissimilarity( mybase->getDataPointIndexInBase(0), idat ) );
- mybase->setTargetDistance( 0, mindist );
- for(unsigned i=1;i<mybase->getNumberOfDataPoints();++i){
- double dist = sqrt( getDissimilarity( mybase->getDataPointIndexInBase(i), idat ) );
- mybase->setTargetDistance( i, dist );
- if( dist<mindist ){ mindist=dist; closest=i; }
+ PCA* ispca = dynamic_cast<PCA*>( mybase );
+ if( ispca ){
+ ispca->getProjection( getStoredData(idat,false), point );
+ } else {
+ ConjugateGradient<ProjectNonLandmarkPoints> myminimiser( this );
+ unsigned closest=0; double mindist = sqrt( getDissimilarity( mybase->getDataPointIndexInBase(0), idat ) );
+ mybase->setTargetDistance( 0, mindist );
+ for(unsigned i=1;i<mybase->getNumberOfDataPoints();++i){
+ double dist = sqrt( getDissimilarity( mybase->getDataPointIndexInBase(i), idat ) );
+ mybase->setTargetDistance( i, dist );
+ if( dist<mindist ){ mindist=dist; closest=i; }
+ }
+ // Put the initial guess near to the closest landmark -- may wish to use grid here again Sandip??
+ Random random; random.setSeed(-1234);
+ for(unsigned j=0;j<nlow;++j) point[j]=mybase->projections(closest,j) + (random.RandU01() - 0.5)*0.01;
+ myminimiser.minimise( cgtol, point, &ProjectNonLandmarkPoints::calculateStress );
}
- // Put the initial guess near to the closest landmark -- may wish to use grid here again Sandip??
- Random random; random.setSeed(-1234);
- for(unsigned j=0;j<nlow;++j) point[j]=mybase->projections(closest,j) + (random.RandU01() - 0.5)*0.01;
- myminimiser.minimise( cgtol, point, &ProjectNonLandmarkPoints::calculateStress );
}
analysis::DataCollectionObject& ProjectNonLandmarkPoints::getStoredData( const unsigned& idat, const bool& calcdist ){
--
GitLab