source: src/ellipsoid.cpp@ 042f82

Action_Thermostats Add_AtomRandomPerturbation Add_FitFragmentPartialChargesAction Add_RotateAroundBondAction Add_SelectAtomByNameAction Added_ParseSaveFragmentResults AddingActions_SaveParseParticleParameters Adding_Graph_to_ChangeBondActions Adding_MD_integration_tests Adding_ParticleName_to_Atom Adding_StructOpt_integration_tests AtomFragments Automaking_mpqc_open AutomationFragmentation_failures Candidate_v1.5.4 Candidate_v1.6.0 Candidate_v1.6.1 ChangeBugEmailaddress ChangingTestPorts ChemicalSpaceEvaluator CombiningParticlePotentialParsing Combining_Subpackages Debian_Package_split Debian_package_split_molecuildergui_only Disabling_MemDebug Docu_Python_wait EmpiricalPotential_contain_HomologyGraph EmpiricalPotential_contain_HomologyGraph_documentation Enable_parallel_make_install Enhance_userguide Enhanced_StructuralOptimization Enhanced_StructuralOptimization_continued Example_ManyWaysToTranslateAtom Exclude_Hydrogens_annealWithBondGraph FitPartialCharges_GlobalError Fix_BoundInBox_CenterInBox_MoleculeActions Fix_ChargeSampling_PBC Fix_ChronosMutex Fix_FitPartialCharges Fix_FitPotential_needs_atomicnumbers Fix_ForceAnnealing Fix_IndependentFragmentGrids Fix_ParseParticles Fix_ParseParticles_split_forward_backward_Actions Fix_PopActions Fix_QtFragmentList_sorted_selection Fix_Restrictedkeyset_FragmentMolecule Fix_StatusMsg Fix_StepWorldTime_single_argument Fix_Verbose_Codepatterns Fix_fitting_potentials Fixes ForceAnnealing_goodresults ForceAnnealing_oldresults ForceAnnealing_tocheck ForceAnnealing_with_BondGraph ForceAnnealing_with_BondGraph_continued ForceAnnealing_with_BondGraph_continued_betteresults ForceAnnealing_with_BondGraph_contraction-expansion FragmentAction_writes_AtomFragments FragmentMolecule_checks_bonddegrees GeometryObjects Gui_Fixes Gui_displays_atomic_force_velocity ImplicitCharges IndependentFragmentGrids IndependentFragmentGrids_IndividualZeroInstances IndependentFragmentGrids_IntegrationTest IndependentFragmentGrids_Sole_NN_Calculation JobMarket_RobustOnKillsSegFaults JobMarket_StableWorkerPool JobMarket_unresolvable_hostname_fix MoreRobust_FragmentAutomation ODR_violation_mpqc_open PartialCharges_OrthogonalSummation PdbParser_setsAtomName PythonUI_with_named_parameters QtGui_reactivate_TimeChanged_changes Recreated_GuiChecks Rewrite_FitPartialCharges RotateToPrincipalAxisSystem_UndoRedo SaturateAtoms_findBestMatching SaturateAtoms_singleDegree StoppableMakroAction Subpackage_CodePatterns Subpackage_JobMarket Subpackage_LinearAlgebra Subpackage_levmar Subpackage_mpqc_open Subpackage_vmg Switchable_LogView ThirdParty_MPQC_rebuilt_buildsystem TrajectoryDependenant_MaxOrder TremoloParser_IncreasedPrecision TremoloParser_MultipleTimesteps TremoloParser_setsAtomName Ubuntu_1604_changes stable
Last change on this file since 042f82 was 042f82, checked in by Frederik Heber <heber@…>, 15 years ago

fixed indentation from tabs to two spaces.

  • Property mode set to 100644
File size: 16.3 KB
Line 
1/*
2 * ellipsoid.cpp
3 *
4 * Created on: Jan 20, 2009
5 * Author: heber
6 */
7
8#include "ellipsoid.hpp"
9
10/** Determines squared distance for a given point \a x to surface of ellipsoid.
11 * \param x given point
12 * \param EllipsoidCenter center of ellipsoid
13 * \param EllipsoidLength[3] three lengths of half axis of ellipsoid
14 * \param EllipsoidAngle[3] three rotation angles of ellipsoid
15 * \return squared distance from point to surface
16 */
17double SquaredDistanceToEllipsoid(Vector &x, Vector &EllipsoidCenter, double *EllipsoidLength, double *EllipsoidAngle)
18{
19 Vector helper, RefPoint;
20 double distance = -1.;
21 double Matrix[NDIM*NDIM];
22 double InverseLength[3];
23 double psi,theta,phi; // euler angles in ZX'Z'' convention
24
25 //cout << Verbose(3) << "Begin of SquaredDistanceToEllipsoid" << endl;
26
27 for(int i=0;i<3;i++)
28 InverseLength[i] = 1./EllipsoidLength[i];
29
30 // 1. translate coordinate system so that ellipsoid center is in origin
31 helper.CopyVector(&x);
32 helper.SubtractVector(&EllipsoidCenter);
33 RefPoint.CopyVector(&helper);
34 //cout << Verbose(4) << "Translated given point is at " << RefPoint << "." << endl;
35
36 // 2. transform coordinate system by inverse of rotation matrix and of diagonal matrix
37 psi = EllipsoidAngle[0];
38 theta = EllipsoidAngle[1];
39 phi = EllipsoidAngle[2];
40 Matrix[0] = cos(psi)*cos(phi) - sin(psi)*cos(theta)*sin(phi);
41 Matrix[1] = -cos(psi)*sin(phi) - sin(psi)*cos(theta)*cos(phi);
42 Matrix[2] = sin(psi)*sin(theta);
43 Matrix[3] = sin(psi)*cos(phi) + cos(psi)*cos(theta)*sin(phi);
44 Matrix[4] = cos(psi)*cos(theta)*cos(phi) - sin(psi)*sin(phi);
45 Matrix[5] = -cos(psi)*sin(theta);
46 Matrix[6] = sin(theta)*sin(phi);
47 Matrix[7] = sin(theta)*cos(phi);
48 Matrix[8] = cos(theta);
49 helper.MatrixMultiplication(Matrix);
50 helper.Scale(InverseLength);
51 //cout << Verbose(4) << "Transformed RefPoint is at " << helper << "." << endl;
52
53 // 3. construct intersection point with unit sphere and ray between origin and x
54 helper.Normalize(); // is simply normalizes vector in distance direction
55 //cout << Verbose(4) << "Transformed intersection is at " << helper << "." << endl;
56
57 // 4. transform back the constructed intersection point
58 psi = -EllipsoidAngle[0];
59 theta = -EllipsoidAngle[1];
60 phi = -EllipsoidAngle[2];
61 helper.Scale(EllipsoidLength);
62 Matrix[0] = cos(psi)*cos(phi) - sin(psi)*cos(theta)*sin(phi);
63 Matrix[1] = -cos(psi)*sin(phi) - sin(psi)*cos(theta)*cos(phi);
64 Matrix[2] = sin(psi)*sin(theta);
65 Matrix[3] = sin(psi)*cos(phi) + cos(psi)*cos(theta)*sin(phi);
66 Matrix[4] = cos(psi)*cos(theta)*cos(phi) - sin(psi)*sin(phi);
67 Matrix[5] = -cos(psi)*sin(theta);
68 Matrix[6] = sin(theta)*sin(phi);
69 Matrix[7] = sin(theta)*cos(phi);
70 Matrix[8] = cos(theta);
71 helper.MatrixMultiplication(Matrix);
72 //cout << Verbose(4) << "Intersection is at " << helper << "." << endl;
73
74 // 5. determine distance between backtransformed point and x
75 distance = RefPoint.DistanceSquared(&helper);
76 //cout << Verbose(4) << "Squared distance between intersection and RefPoint is " << distance << "." << endl;
77
78 return distance;
79 //cout << Verbose(3) << "End of SquaredDistanceToEllipsoid" << endl;
80};
81
82/** structure for ellipsoid minimisation containing points to fit to.
83 */
84struct EllipsoidMinimisation {
85 int N; //!< dimension of vector set
86 Vector *x; //!< array of vectors
87};
88
89/** Sum of squared distance to ellipsoid to be minimised.
90 * \param *x parameters for the ellipsoid
91 * \param *params EllipsoidMinimisation with set of data points to minimise distance to and dimension
92 * \return sum of squared distance, \sa SquaredDistanceToEllipsoid()
93 */
94double SumSquaredDistance (const gsl_vector * x, void * params)
95{
96 Vector *set= ((struct EllipsoidMinimisation *)params)->x;
97 int N = ((struct EllipsoidMinimisation *)params)->N;
98 double SumDistance = 0.;
99 double distance;
100 Vector Center;
101 double EllipsoidLength[3], EllipsoidAngle[3];
102
103 // put parameters into suitable ellipsoid form
104 for (int i=0;i<3;i++) {
105 Center.x[i] = gsl_vector_get(x, i+0);
106 EllipsoidLength[i] = gsl_vector_get(x, i+3);
107 EllipsoidAngle[i] = gsl_vector_get(x, i+6);
108 }
109
110 // go through all points and sum distance
111 for (int i=0;i<N;i++) {
112 distance = SquaredDistanceToEllipsoid(set[i], Center, EllipsoidLength, EllipsoidAngle);
113 if (!isnan(distance)) {
114 SumDistance += distance;
115 } else {
116 SumDistance = GSL_NAN;
117 break;
118 }
119 }
120
121 //cout << "Current summed distance is " << SumDistance << "." << endl;
122 return SumDistance;
123};
124
125/** Finds best fitting ellipsoid parameter set in Least square sense for a given point set.
126 * \param *out output stream for debugging
127 * \param *set given point set
128 * \param N number of points in set
129 * \param EllipsoidParamter[3] three parameters in ellipsoid equation
130 * \return true - fit successful, false - fit impossible
131 */
132bool FitPointSetToEllipsoid(ofstream *out, Vector *set, int N, Vector *EllipsoidCenter, double *EllipsoidLength, double *EllipsoidAngle)
133{
134 int status = GSL_SUCCESS;
135 *out << Verbose(2) << "Begin of FitPointSetToEllipsoid " << endl;
136 if (N >= 3) { // check that enough points are given (9 d.o.f.)
137 struct EllipsoidMinimisation par;
138 const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
139 gsl_multimin_fminimizer *s = NULL;
140 gsl_vector *ss, *x;
141 gsl_multimin_function minex_func;
142
143 size_t iter = 0;
144 double size;
145
146 /* Starting point */
147 x = gsl_vector_alloc (9);
148 for (int i=0;i<3;i++) {
149 gsl_vector_set (x, i+0, EllipsoidCenter->x[i]);
150 gsl_vector_set (x, i+3, EllipsoidLength[i]);
151 gsl_vector_set (x, i+6, EllipsoidAngle[i]);
152 }
153 par.x = set;
154 par.N = N;
155
156 /* Set initial step sizes */
157 ss = gsl_vector_alloc (9);
158 for (int i=0;i<3;i++) {
159 gsl_vector_set (ss, i+0, 0.1);
160 gsl_vector_set (ss, i+3, 1.0);
161 gsl_vector_set (ss, i+6, M_PI/20.);
162 }
163
164 /* Initialize method and iterate */
165 minex_func.n = 9;
166 minex_func.f = &SumSquaredDistance;
167 minex_func.params = (void *)&par;
168
169 s = gsl_multimin_fminimizer_alloc (T, 9);
170 gsl_multimin_fminimizer_set (s, &minex_func, x, ss);
171
172 do {
173 iter++;
174 status = gsl_multimin_fminimizer_iterate(s);
175
176 if (status)
177 break;
178
179 size = gsl_multimin_fminimizer_size (s);
180 status = gsl_multimin_test_size (size, 1e-2);
181
182 if (status == GSL_SUCCESS) {
183 for (int i=0;i<3;i++) {
184 EllipsoidCenter->x[i] = gsl_vector_get (s->x,i+0);
185 EllipsoidLength[i] = gsl_vector_get (s->x, i+3);
186 EllipsoidAngle[i] = gsl_vector_get (s->x, i+6);
187 }
188 *out << setprecision(3) << Verbose(4) << "Converged fit at: " << *EllipsoidCenter << ", lengths " << EllipsoidLength[0] << ", " << EllipsoidLength[1] << ", " << EllipsoidLength[2] << ", angles " << EllipsoidAngle[0] << ", " << EllipsoidAngle[1] << ", " << EllipsoidAngle[2] << " with summed distance " << s->fval << "." << endl;
189 }
190
191 } while (status == GSL_CONTINUE && iter < 1000);
192
193 gsl_vector_free(x);
194 gsl_vector_free(ss);
195 gsl_multimin_fminimizer_free (s);
196
197 } else {
198 *out << Verbose(3) << "Not enough points provided for fit to ellipsoid." << endl;
199 return false;
200 }
201 *out << Verbose(2) << "End of FitPointSetToEllipsoid" << endl;
202 if (status == GSL_SUCCESS)
203 return true;
204 else
205 return false;
206};
207
208/** Picks a number of random points from a LC neighbourhood as a fitting set.
209 * \param *out output stream for debugging
210 * \param *T Tesselation containing boundary points
211 * \param *LC linked cell list of all atoms
212 * \param *&x random point set on return (not allocated!)
213 * \param PointsToPick number of points in set to pick
214 */
215void PickRandomNeighbouredPointSet(ofstream *out, class Tesselation *T, class LinkedCell *LC, Vector *&x, size_t PointsToPick)
216{
217 size_t PointsLeft = 0;
218 size_t PointsPicked = 0;
219 int Nlower[NDIM], Nupper[NDIM];
220 set<int> PickedAtomNrs; // ordered list of picked atoms
221 set<int>::iterator current;
222 int index;
223 atom *Candidate = NULL;
224 LinkedAtoms *List = NULL;
225 *out << Verbose(2) << "Begin of PickRandomPointSet" << endl;
226
227 // allocate array
228 if (x == NULL) {
229 x = new Vector[PointsToPick];
230 } else {
231 *out << "WARNING: Given pointer to vector array seems already allocated." << endl;
232 }
233
234 do {
235 for(int i=0;i<NDIM;i++) // pick three random indices
236 LC->n[i] = (rand() % LC->N[i]);
237 *out << Verbose(2) << "INFO: Center cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " ... ";
238 // get random cell
239 List = LC->GetCurrentCell();
240 if (List == NULL) { // set index to it
241 continue;
242 }
243 *out << "with No. " << LC->index << "." << endl;
244
245 *out << Verbose(2) << "LC Intervals:";
246 for (int i=0;i<NDIM;i++) {
247 Nlower[i] = ((LC->n[i]-1) >= 0) ? LC->n[i]-1 : 0;
248 Nupper[i] = ((LC->n[i]+1) < LC->N[i]) ? LC->n[i]+1 : LC->N[i]-1;
249 *out << " [" << Nlower[i] << "," << Nupper[i] << "] ";
250 }
251 *out << endl;
252
253 // count whether there are sufficient atoms in this cell+neighbors
254 PointsLeft=0;
255 for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
256 for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
257 for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
258 List = LC->GetCurrentCell();
259 PointsLeft += List->size();
260 }
261 *out << Verbose(2) << "There are " << PointsLeft << " atoms in this neighbourhood." << endl;
262 if (PointsLeft < PointsToPick) { // ensure that we can pick enough points in its neighbourhood at all.
263 continue;
264 }
265
266 // pre-pick a fixed number of atoms
267 PickedAtomNrs.clear();
268 do {
269 index = (rand() % PointsLeft);
270 current = PickedAtomNrs.find(index); // not present?
271 if (current == PickedAtomNrs.end()) {
272 //*out << Verbose(2) << "Picking atom nr. " << index << "." << endl;
273 PickedAtomNrs.insert(index);
274 }
275 } while (PickedAtomNrs.size() < PointsToPick);
276
277 index = 0; // now go through all and pick those whose from PickedAtomsNr
278 PointsPicked=0;
279 current = PickedAtomNrs.begin();
280 for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
281 for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
282 for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
283 List = LC->GetCurrentCell();
284// *out << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << " containing " << List->size() << " points." << endl;
285 if (List != NULL) {
286// if (List->begin() != List->end())
287// *out << Verbose(2) << "Going through candidates ... " << endl;
288// else
289// *out << Verbose(2) << "Cell is empty ... " << endl;
290 for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
291 if ((current != PickedAtomNrs.end()) && (*current == index)) {
292 Candidate = (*Runner);
293 *out << Verbose(2) << "Current picked node is " << **Runner << " with index " << index << "." << endl;
294 x[PointsPicked++].CopyVector(&(Candidate->x)); // we have one more atom picked
295 current++; // next pre-picked atom
296 }
297 index++; // next atom nr.
298 }
299// } else {
300// *out << Verbose(2) << "List for this index not allocated!" << endl;
301 }
302 }
303 *out << Verbose(2) << "The following points were picked: " << endl;
304 for (size_t i=0;i<PointsPicked;i++)
305 *out << Verbose(2) << x[i] << endl;
306 if (PointsPicked == PointsToPick) // break out of loop if we have all
307 break;
308 } while(1);
309
310 *out << Verbose(2) << "End of PickRandomPointSet" << endl;
311};
312
313/** Picks a number of random points from a set of boundary points as a fitting set.
314 * \param *out output stream for debugging
315 * \param *T Tesselation containing boundary points
316 * \param *&x random point set on return (not allocated!)
317 * \param PointsToPick number of points in set to pick
318 */
319void PickRandomPointSet(ofstream *out, class Tesselation *T, Vector *&x, size_t PointsToPick)
320{
321 size_t PointsLeft = (size_t) T->PointsOnBoundaryCount;
322 size_t PointsPicked = 0;
323 double value, threshold;
324 PointMap *List = &T->PointsOnBoundary;
325 *out << Verbose(2) << "Begin of PickRandomPointSet" << endl;
326
327 // allocate array
328 if (x == NULL) {
329 x = new Vector[PointsToPick];
330 } else {
331 *out << "WARNING: Given pointer to vector array seems already allocated." << endl;
332 }
333
334 if (List != NULL)
335 for (PointMap::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
336 threshold = 1. - (double)(PointsToPick - PointsPicked)/(double)PointsLeft;
337 value = (double)rand()/(double)RAND_MAX;
338 //*out << Verbose(3) << "Current node is " << *Runner->second->node << " with " << value << " ... " << threshold << ": ";
339 if (value > threshold) {
340 x[PointsPicked].CopyVector(&(Runner->second->node->x));
341 PointsPicked++;
342 //*out << "IN." << endl;
343 } else {
344 //*out << "OUT." << endl;
345 }
346 PointsLeft--;
347 }
348 *out << Verbose(2) << "The following points were picked: " << endl;
349 for (size_t i=0;i<PointsPicked;i++)
350 *out << Verbose(3) << x[i] << endl;
351
352 *out << Verbose(2) << "End of PickRandomPointSet" << endl;
353};
354
355/** Finds best fitting ellipsoid parameter set in least square sense for a given point set.
356 * \param *out output stream for debugging
357 * \param *T Tesselation containing boundary points
358 * \param *LCList linked cell list of all atoms
359 * \param N number of unique points in ellipsoid fit, must be greater equal 6
360 * \param number of fits (i.e. parameter sets in output file)
361 * \param *filename name for output file
362 */
363void FindDistributionOfEllipsoids(ofstream *out, class Tesselation *T, class LinkedCell *LCList, int N, int number, const char *filename)
364{
365 ofstream output;
366 Vector *x = NULL;
367 Vector Center;
368 Vector EllipsoidCenter;
369 double EllipsoidLength[3];
370 double EllipsoidAngle[3];
371 double distance, MaxDistance, MinDistance;
372 *out << Verbose(0) << "Begin of FindDistributionOfEllipsoids" << endl;
373
374 // construct center of gravity of boundary point set for initial ellipsoid center
375 Center.Zero();
376 for (PointMap::iterator Runner = T->PointsOnBoundary.begin(); Runner != T->PointsOnBoundary.end(); Runner++)
377 Center.AddVector(&Runner->second->node->x);
378 Center.Scale(1./T->PointsOnBoundaryCount);
379 *out << Verbose(1) << "Center is at " << Center << "." << endl;
380
381 // Output header
382 output.open(filename, ios::trunc);
383 output << "# Nr.\tCenterX\tCenterY\tCenterZ\ta\tb\tc\tpsi\ttheta\tphi" << endl;
384
385 // loop over desired number of parameter sets
386 for (;number >0;number--) {
387 *out << Verbose(1) << "Determining data set " << number << " ... " << endl;
388 // pick the point set
389 x = NULL;
390 //PickRandomPointSet(out, T, LCList, x, N);
391 PickRandomNeighbouredPointSet(out, T, LCList, x, N);
392
393 // calculate some sensible starting values for parameter fit
394 MaxDistance = 0.;
395 MinDistance = x[0].ScalarProduct(&x[0]);
396 for (int i=0;i<N;i++) {
397 distance = x[i].ScalarProduct(&x[i]);
398 if (distance > MaxDistance)
399 MaxDistance = distance;
400 if (distance < MinDistance)
401 MinDistance = distance;
402 }
403 //*out << Verbose(2) << "MinDistance " << MinDistance << ", MaxDistance " << MaxDistance << "." << endl;
404 EllipsoidCenter.CopyVector(&Center); // use Center of Gravity as initial center of ellipsoid
405 for (int i=0;i<3;i++)
406 EllipsoidAngle[i] = 0.;
407 EllipsoidLength[0] = sqrt(MaxDistance);
408 EllipsoidLength[1] = sqrt((MaxDistance+MinDistance)/2.);
409 EllipsoidLength[2] = sqrt(MinDistance);
410
411 // fit the parameters
412 if (FitPointSetToEllipsoid(out, x, N, &EllipsoidCenter, &EllipsoidLength[0], &EllipsoidAngle[0])) {
413 *out << Verbose(1) << "Picking succeeded!" << endl;
414 // output obtained parameter set
415 output << number << "\t";
416 for (int i=0;i<3;i++)
417 output << setprecision(9) << EllipsoidCenter.x[i] << "\t";
418 for (int i=0;i<3;i++)
419 output << setprecision(9) << EllipsoidLength[i] << "\t";
420 for (int i=0;i<3;i++)
421 output << setprecision(9) << EllipsoidAngle[i] << "\t";
422 output << endl;
423 } else { // increase N to pick one more
424 *out << Verbose(1) << "Picking failed!" << endl;
425 number++;
426 }
427 delete[](x); // free allocated memory for point set
428 }
429 // close output and finish
430 output.close();
431
432 *out << Verbose(0) << "End of FindDistributionOfEllipsoids" << endl;
433};
Note: See TracBrowser for help on using the repository browser.