Ignore:
Timestamp:
Jul 23, 2009, 2:21:57 PM (16 years ago)
Author:
Frederik Heber <heber@…>
Children:
c3a303
Parents:
c95b69
Message:

fixed indentation from tabs to two spaces.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • molecuilder/src/ellipsoid.cpp

    rc95b69 ra048fa  
    22 * ellipsoid.cpp
    33 *
    4  *      Created on: Jan 20, 2009
    5  *                      Author: heber
     4 *  Created on: Jan 20, 2009
     5 *      Author: heber
    66 */
    77
    … …  
    1717double SquaredDistanceToEllipsoid(Vector &x, Vector &EllipsoidCenter, double *EllipsoidLength, double *EllipsoidAngle)
    1818{
    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;
     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;
    8080};
    8181
    … …  
    8383 */
    8484struct EllipsoidMinimisation {
    85         int N;                  //!< dimension of vector set
    86         Vector *x;      //!< array of vectors
     85  int N;      //!< dimension of vector set
     86  Vector *x;  //!< array of vectors
    8787};
    8888
    … …  
    9494double SumSquaredDistance (const gsl_vector * x, void * params)
    9595{
    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;
     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;
    123123};
    124124
    … …  
    132132bool FitPointSetToEllipsoid(ofstream *out, Vector *set, int N, Vector *EllipsoidCenter, double *EllipsoidLength, double *EllipsoidAngle)
    133133{
    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;
     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;
    206206};
    207207
    … …  
    217217  size_t PointsLeft = 0;
    218218  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;
     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;
    311311};
    312312
    … …  
    321321  size_t PointsLeft = (size_t) T->PointsOnBoundaryCount;
    322322  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;
     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;
    353353};
    354354
    … …  
    363363void FindDistributionOfEllipsoids(ofstream *out, class Tesselation *T, class LinkedCell *LCList, int N, int number, const char *filename)
    364364{
    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 };
     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 TracChangeset for help on using the changeset viewer.