1 /*! \file lib/map_utils.h
2     Header file for map statistics type
3 */
4 //C Copyright (C) 2000-2006 Kevin Cowtan and University of York
5 //L
6 //L  This library is free software and is distributed under the terms
7 //L  and conditions of version 2.1 of the GNU Lesser General Public
8 //L  Licence (LGPL) with the following additional clause:
9 //L
10 //L     `You may also combine or link a "work that uses the Library" to
11 //L     produce a work containing portions of the Library, and distribute
12 //L     that work under terms of your choice, provided that you give
13 //L     prominent notice with each copy of the work that the specified
14 //L     version of the Library is used in it, and that you include or
15 //L     provide public access to the complete corresponding
16 //L     machine-readable source code for the Library including whatever
17 //L     changes were used in the work. (i.e. If you make changes to the
18 //L     Library you must distribute those, but you do not need to
19 //L     distribute source or object code to those portions of the work
20 //L     not covered by this licence.)'
21 //L
22 //L  Note that this clause grants an additional right and does not impose
23 //L  any additional restriction, and so does not affect compatibility
24 //L  with the GNU General Public Licence (GPL). If you wish to negotiate
25 //L  other terms, please contact the maintainer.
26 //L
27 //L  You can redistribute it and/or modify the library under the terms of
28 //L  the GNU Lesser General Public License as published by the Free Software
29 //L  Foundation; either version 2.1 of the License, or (at your option) any
30 //L  later version.
31 //L
32 //L  This library is distributed in the hope that it will be useful, but
33 //L  WITHOUT ANY WARRANTY; without even the implied warranty of
34 //L  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
35 //L  Lesser General Public License for more details.
36 //L
37 //L  You should have received a copy of the CCP4 licence and/or GNU
38 //L  Lesser General Public License along with this library; if not, write
39 //L  to the CCP4 Secretary, Daresbury Laboratory, Warrington WA4 4AD, UK.
40 //L  The GNU Lesser General Public can also be obtained by writing to the
41 //L  Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
42 //L  MA 02111-1307 USA
43 
44 
45 #ifndef CLIPPER_MAP_UTILS
46 #define CLIPPER_MAP_UTILS
47 
48 
49 #include "derivs.h"
50 
51 
52 namespace clipper
53 {
54 
55   //! Generic map statistics class
56   /*! This class is used to calculate and store the mean and standard
57     deviation of a generic map object of scalar types (e.g. Xmap,
58     NXmap). If the map contains NaN values, those points are excluded
59     for the calculation. In the case of an Xmap, the proper
60     multiplicty corrections are applied to give statistics for a whole
61     unit cell */
62   class Map_stats
63   {
64   public:
Map_stats()65     Map_stats() {}                                //!< null constructor
66     template<class M> Map_stats( const M& map );  //!< Constructor: from Xmap
mean()67     const ftype& mean() const { return mean_; }        //!< Mean of map
std_dev()68     const ftype& std_dev() const { return std_dev_; }  //!< Std deviation of map
min()69     const ftype& min() const { return min_; }          //!< Minimum of map
max()70     const ftype& max() const { return max_; }          //!< Maximum of map
range()71     const Range<> range() const { return Range<>( min_, max_ ); } //!< Range
72   private:
73     ftype mean_, std_dev_, min_, max_;
74   };
75 
76 
77   //! Generic map sorting class
78   /*! This class is used to sort a vector of integer indices into a
79     map. This includes sorting the whole map to get highest or lowest
80     density first, or sorting some subset, e.g. a list of peak
81     indices. Integer indices are used because they are the most
82     compact way of referencing a unique map location. e.g.
83     \code
84     clipper::Xmap<float>::Map_reference_index ix;
85     std::vector<int> index;
86     for ( ix = xmap.first(); !ix.last(); ix.next() )
87       index.push_back( ix.index() );
88     Map_index_sort::sort_decreasing( xmap, index );
89     \endcode
90   */
91   class Map_index_sort
92   {
93   public:
94     //! Sort a list into increasing order
95     template<class M> static void sort_increasing( const M& map, std::vector<int>& index );
96     //! Sort a list into decreasing order
97     template<class M> static void sort_decreasing( const M& map, std::vector<int>& index );
98     //! Internal helper class used for sorting
99   private:
100     template<class M> class Compare_density {
101     public:
Compare_density(const M & m)102       Compare_density( const M& m ) { p = &m; }
operator()103       bool operator() ( const int& i1, const int& i2 ) const { return p->get_data(i1) < p->get_data(i2); }
104     private:
105       const M* p;
106     };
107   };
108 
109 
110   // template implementations
111 
112   /*! For float and double maps
113     \params map The map for which moments are to be calculated. */
Map_stats(const M & map)114   template<class M> Map_stats::Map_stats( const M& map )
115   {
116     ftype64 w, x, s, sx, sxx;
117     s = sx = sxx = 0.0;
118     min_ =  1.0e12;
119     max_ = -1.0e12;
120     for ( typename M::Map_reference_index im = map.first();
121 	  !im.last(); im.next() ) {
122       w = 1.0 / ftype64( map.multiplicity( im.coord() ) );
123       x = ftype64( map[im] );
124       if ( !Util::is_nan(x) ) {
125 	s += w;
126 	sx += w*x;
127 	sxx += w*x*x;
128 	if ( x < min_ ) min_ = x;
129 	if ( x > max_ ) max_ = x;
130       }
131     }
132     sx /= s;
133     sxx /= s;
134     mean_ = sx;
135     std_dev_ = sqrt( sxx - sx*sx );
136   }
137 
138 
139 } // namespace clipper
140 
141 #endif
142