← Back to branch summary

~mok0/clipper/trunk

« back to all changes in this revision

Viewing changes to clipper/core/resol_targetfn.h

  • Committer: Morten Kjeldgaard
  • Date: 2011-02-09 21:25:56 UTC
  • Revision ID: mok@bioxray.dk-20110209212556-3zbon5ukh4zf3x00
Tags: 2.0.3
version 2.0.3

Show diffs side-by-side

added added

removed removed

Lines of Context:
clipper/core/resol_targetfn.h
 
1
/*! \file lib/resol_targetfn.h
 
2
    Header file for resolution function generator
 
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_RESOL_TARGETFN
 
46
#define CLIPPER_RESOL_TARGETFN
 
47
 
 
48
#include "resol_basisfn.h"
 
49
#include "hkl_datatypes.h"
 
50
 
 
51
 
 
52
namespace clipper {
 
53
 
 
54
 
 
55
  //! simple mean |F|<sup>n</sup> target
 
56
  /*! This class implements the target function for calculating mean
 
57
    |F|<sup>n</sup> as a function of position in reciprocal space. It
 
58
    includes the appropriate multiplicity correction, and so can be
 
59
    applied to any type with an 'f' member with the same dimensions as
 
60
    an |F| or |U| (or an uncorrected |E|).
 
61
 
 
62
    \Note This function should not be used to scale F's to E's.
 
63
    See TargetFn_scaleEsq. */
 
64
  template<class T> class TargetFn_meanFnth : public TargetFn_base
 
65
  {
 
66
  public:
 
67
    //! constructor: takes the datalist against which to calc target, and power
 
68
    TargetFn_meanFnth( const HKL_data<T>& hkl_data_, const ftype& n );
 
69
    //! return the value and derivatives of the target function
 
70
    Rderiv rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const;
 
71
    //! the type of the function: optionally used to improve convergence
 
72
    FNtype type() const { return QUADRATIC; }
 
73
  private:
 
74
    ftype power;
 
75
    const HKL_data<T>* hkl_data;
 
76
  };
 
77
 
 
78
 
 
79
  //! simple mean |E|<sup>n</sup> target
 
80
  /*! This class implements the target function for calculating mean
 
81
    |E|<sup>n</sup> as a function of position in reciprocal space. It
 
82
    includes the appropriate multiplicity correction, and so can be
 
83
    applied to any type with an 'E' member with the same dimensions as
 
84
    an |E| (or corrected |F| or |U|).
 
85
 
 
86
    \Note This function should not be used to scale F's to E's.
 
87
    See TargetFn_scaleEsq. */
 
88
  template<class T> class TargetFn_meanEnth : public TargetFn_base
 
89
  {
 
90
  public:
 
91
    //! constructor: takes the datalist against which to calc target, and power
 
92
    TargetFn_meanEnth( const HKL_data<T>& hkl_data_, const ftype& n );
 
93
    //! return the value and derivatives of the target function
 
94
    Rderiv rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const;
 
95
    //! the type of the function: optionally used to improve convergence
 
96
    FNtype type() const { return QUADRATIC; }
 
97
  private:
 
98
    ftype power;
 
99
    const HKL_data<T>* hkl_data;
 
100
  };
 
101
 
 
102
 
 
103
  //! |F|<sup>2</sup> scaling target
 
104
  /*! This class implements the target function for calculating the
 
105
    scale factor to scale one set of F's to another. The resulting
 
106
    scale is the square of the factor that scales the first set of
 
107
    data to match the second. */
 
108
  template<class T1, class T2> class TargetFn_scaleF1F2 : public TargetFn_base
 
109
  {
 
110
  public:
 
111
    //! constructor: takes the datalist against which to calc target
 
112
    TargetFn_scaleF1F2( const HKL_data<T1>& hkl_data1_, const HKL_data<T2>& hkl_data2_ );
 
113
    //! return the value and derivatives of the target function
 
114
    Rderiv rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const;
 
115
    //! the type of the function: optionally used to improve convergence
 
116
    FNtype type() const { return QUADRATIC; }
 
117
  private:
 
118
    const HKL_data<T1>* hkl_data1;
 
119
    const HKL_data<T2>* hkl_data2;
 
120
  };
 
121
 
 
122
 
 
123
  //! log |F|<sup>2</sup> scaling target
 
124
  /*! This class implements the target function for calculating the
 
125
    scale factor to scale the weighted log of one set of F's to
 
126
    another. The resulting scale is the square of the factor that
 
127
    scales the first set of data to match the second. The log scaling
 
128
    target is used in conjunction with the log-Gaussian basis
 
129
    functions for a fast and robust approximation to iso/aniso
 
130
    Gaussian scaling.
 
131
  */
 
132
  template<class T1, class T2> class TargetFn_scaleLogF1F2 : public TargetFn_base
 
133
  {
 
134
  public:
 
135
    //! constructor: takes the datalist against which to calc target
 
136
    TargetFn_scaleLogF1F2( const HKL_data<T1>& hkl_data1_, const HKL_data<T2>& hkl_data2_ );
 
137
    //! return the value and derivatives of the target function
 
138
    Rderiv rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const;
 
139
    //! the type of the function: optionally used to improve convergence
 
140
    FNtype type() const { return QUADRATIC; }
 
141
  private:
 
142
    const HKL_data<T1>* hkl_data1;
 
143
    const HKL_data<T2>* hkl_data2;
 
144
  };
 
145
 
 
146
 
 
147
  //! |E|<sup>2</sup> scaling target
 
148
  /*! This class implements the target function for calculating the
 
149
    scale factor to normalise to <|E|<sup>2</sup>> = 1. Note that this
 
150
    is not the same as dividing by <|E|<sup>2</sup>>, except in a few
 
151
    special cases, e.g. a simple resolution bins calculation. The
 
152
    resulting targen function is the square of the value by which |E|
 
153
    should be multiplied to acheive the correct normalisation. */
 
154
  template<class T> class TargetFn_scaleEsq : public TargetFn_base
 
155
  {
 
156
  public:
 
157
    //! constructor: takes the datalist against which to calc target
 
158
    TargetFn_scaleEsq( const HKL_data<T>& hkl_data_ );
 
159
    //! return the value and derivatives of the target function
 
160
    Rderiv rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const;
 
161
    //! the type of the function: optionally used to improve convergence
 
162
    FNtype type() const { return QUADRATIC; }
 
163
  private:
 
164
    const HKL_data<T>* hkl_data;
 
165
  };
 
166
 
 
167
 
 
168
  //! \deprecated simple sigma_a target function
 
169
  /*! This class implements the target function for calculating sigma_a.
 
170
    Required is a datalist containing Eo, Ec.
 
171
 
 
172
    It actually refines omegaa = sigmaa/(1-sigmaa^2). This has better
 
173
    proerties for refinement. To get sigmaa use
 
174
    \code sigmaa = ( sqrt( 4*omegaa^2 + 1 ) - 1 ) / ( 2*omegaa ) \endcode
 
175
    This is available as a static function:
 
176
    \code sigmaa = targetfn.sigmaa( omegaa ) \endcode
 
177
 
 
178
    This version simplifies terms in |Eo|^2 and |Ec|^2 which should
 
179
    average out to 1 if the normalisation scheme is consistent with
 
180
    the sigmaa calc.
 
181
 
 
182
    Convergence is good for calculations using the 'binner' basis
 
183
    function, however the smooth basis function have convergence
 
184
    problems. This is still under investigation.
 
185
  */
 
186
  template<class T> class TargetFn_sigmaa_omegaa : public TargetFn_base
 
187
  {
 
188
  public:
 
189
    //! constructor: takes the datalist against which to calc target
 
190
    TargetFn_sigmaa_omegaa( const HKL_data<T>& eo, const HKL_data<T>& ec );
 
191
    //! return the value and derivatives of the target function
 
192
    Rderiv rderiv( const HKL_info::HKL_reference_index& ih, const ftype& omegaa ) const;
 
193
    //! convert omegaa to sigmaa
 
194
    static ftype sigmaa( const ftype& omegaa )
 
195
    {
 
196
      ftype omeg = (omegaa>0.05) ? omegaa : (0.05*exp(omegaa/0.05-1.0));
 
197
      return 0.5 * ( sqrt( 4.0*omeg*omeg + 1.0 ) - 1.0 ) / omeg;
 
198
    }
 
199
  private:
 
200
    const HKL_data<T>* eo_data;
 
201
    const HKL_data<T>* ec_data;
 
202
  };
 
203
 
 
204
 
 
205
  //! \deprecated simple sigma_a target function
 
206
  /*! \par Warning: Convergence of this basis-function can be
 
207
    unreliable under some circumstances. Use
 
208
    clipper::TargetFn_sigmaa_omegaa instead, except for development
 
209
    purposes.
 
210
 
 
211
    This class implements the target function for calculating sigma_a.
 
212
    Required is a datalist containing Eo, Ec.
 
213
 
 
214
    This version simplifies terms in |Eo|^2 and |Ec|^2 which should
 
215
    average out to 1 if the normalisation scheme is consistent with
 
216
    the sigmaa calc.
 
217
  */
 
218
  template<class T> class TargetFn_sigmaa : public TargetFn_base
 
219
  {
 
220
  public:
 
221
    //! constructor: takes the datalist against which to calc target
 
222
    TargetFn_sigmaa( const HKL_data<T>& eo, const HKL_data<T>& ec );
 
223
    //! return the value and derivatives of the target function
 
224
    Rderiv rderiv( const HKL_info::HKL_reference_index& ih, const ftype& sigmaa0 ) const;
 
225
    //! convert function to sigmaa
 
226
    static ftype sigmaa( const ftype& sigm ) { return sigm; }
 
227
  private:
 
228
    const HKL_data<T>* eo_data;
 
229
    const HKL_data<T>* ec_data;
 
230
  };
 
231
 
 
232
 
 
233
 
 
234
  // implementations for template functions
 
235
 
 
236
  // mean F^nth
 
237
 
 
238
  template<class T> TargetFn_meanFnth<T>::TargetFn_meanFnth( const HKL_data<T>& hkl_data_, const ftype& n )
 
239
  {
 
240
    power = n;
 
241
    hkl_data = &hkl_data_;
 
242
  }
 
243
 
 
244
  template<class T> TargetFn_base::Rderiv TargetFn_meanFnth<T>::rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const
 
245
  {
 
246
    Rderiv result;
 
247
    const HKL_data<T>& data = *hkl_data;
 
248
    if ( !data[ih].missing() ) {
 
249
      ftype d = fh - pow( ftype(data[ih].f()) / sqrt(ih.hkl_class().epsilon()),
 
250
                          power );
 
251
      result.r = d * d;
 
252
      result.dr = 2.0 * d;
 
253
      result.dr2 = 2.0;
 
254
    } else {
 
255
      result.r = result.dr = result.dr2 = 0.0;
 
256
    }
 
257
    return result;
 
258
  }
 
259
 
 
260
  // mean E^nth
 
261
 
 
262
  template<class T> TargetFn_meanEnth<T>::TargetFn_meanEnth( const HKL_data<T>& hkl_data_, const ftype& n )
 
263
  {
 
264
    power = n;
 
265
    hkl_data = &hkl_data_;
 
266
  }
 
267
 
 
268
  template<class T> TargetFn_base::Rderiv TargetFn_meanEnth<T>::rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const
 
269
  {
 
270
    Rderiv result;
 
271
    const HKL_data<T>& data = *hkl_data;
 
272
    if ( !data[ih].missing() ) {
 
273
      ftype d = fh - pow( ftype(data[ih].E()), power );
 
274
      result.r = d * d;
 
275
      result.dr = 2.0 * d;
 
276
      result.dr2 = 2.0;
 
277
    } else {
 
278
      result.r = result.dr = result.dr2 = 0.0;
 
279
    }
 
280
    return result;
 
281
  }
 
282
 
 
283
  // F1-F2 scaling
 
284
 
 
285
  template<class T1, class T2> TargetFn_scaleF1F2<T1,T2>::TargetFn_scaleF1F2( const HKL_data<T1>& hkl_data1_, const HKL_data<T2>& hkl_data2_ )
 
286
  {
 
287
    hkl_data1 = &hkl_data1_;
 
288
    hkl_data2 = &hkl_data2_;
 
289
  }
 
290
 
 
291
  template<class T1, class T2> TargetFn_base::Rderiv TargetFn_scaleF1F2<T1,T2>::rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const
 
292
  {
 
293
    Rderiv result;
 
294
    const T1& ft1 = (*hkl_data1)[ih];
 
295
    const T2& ft2 = (*hkl_data2)[ih];
 
296
    if ( !ft1.missing() && !ft2.missing() ) {
 
297
      const ftype eps = ih.hkl_class().epsilon();
 
298
      const ftype f1 = pow( ft1.f(), 2 ) / eps;
 
299
      const ftype f2 = pow( ft2.f(), 2 ) / eps;
 
300
      const ftype d = fh*f1 - f2;
 
301
      result.r = d * d / f1;
 
302
      result.dr = 2.0 * d;
 
303
      result.dr2 = 2.0 * f1;
 
304
    } else {
 
305
      result.r = result.dr = result.dr2 = 0.0;
 
306
    }
 
307
    return result;
 
308
  }
 
309
 
 
310
  // Log F1-F2 scaling
 
311
 
 
312
  template<class T1, class T2> TargetFn_scaleLogF1F2<T1,T2>::TargetFn_scaleLogF1F2( const HKL_data<T1>& hkl_data1_, const HKL_data<T2>& hkl_data2_ )
 
313
  {
 
314
    hkl_data1 = &hkl_data1_;
 
315
    hkl_data2 = &hkl_data2_;
 
316
  }
 
317
 
 
318
  template<class T1, class T2> TargetFn_base::Rderiv TargetFn_scaleLogF1F2<T1,T2>::rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const
 
319
  {
 
320
    Rderiv result;
 
321
    result.r = result.dr = result.dr2 = 0.0;
 
322
    const T1& ft1 = (*hkl_data1)[ih];
 
323
    const T2& ft2 = (*hkl_data2)[ih];
 
324
    if ( !ft1.missing() && !ft2.missing() )
 
325
      if ( ft1.f() > 1.0e-6 && ft2.f() > 1.0e-6 ) {
 
326
        const ftype eps = ih.hkl_class().epsilon();
 
327
        const ftype f1 = pow( ft1.f(), 2 ) / eps;
 
328
        const ftype f2 = pow( ft2.f(), 2 ) / eps;
 
329
        const ftype w = sqrt( f1 * f2 );  // f1 * f2;
 
330
        const ftype d = fh + log(f1) - log(f2);
 
331
        result.r   =       w * d * d;
 
332
        result.dr  = 2.0 * w * d;
 
333
        result.dr2 = 2.0 * w;
 
334
    }
 
335
    return result;
 
336
  }
 
337
 
 
338
  // E^2 scaling
 
339
 
 
340
  template<class T> TargetFn_scaleEsq<T>::TargetFn_scaleEsq( const HKL_data<T>& hkl_data_ )
 
341
  {
 
342
    hkl_data = &hkl_data_;
 
343
  }
 
344
 
 
345
  template<class T> TargetFn_base::Rderiv TargetFn_scaleEsq<T>::rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const
 
346
  {
 
347
    Rderiv result;
 
348
    const HKL_data<T>& data = *hkl_data;
 
349
    const ftype two(2.0);
 
350
    if ( !data[ih].missing() ) {
 
351
      ftype fsq = pow( ftype(data[ih].E()), two );
 
352
      ftype d = fsq * fh - 1.0;
 
353
      result.r = d * d / fsq;
 
354
      result.dr = two * d;
 
355
      result.dr2 = two * fsq;
 
356
    } else {
 
357
      result.r = result.dr = result.dr2 = 0.0;
 
358
    }
 
359
    return result;
 
360
  }
 
361
 
 
362
 
 
363
  // sigmaa (omegaa)
 
364
 
 
365
  template<class T> TargetFn_sigmaa_omegaa<T>::TargetFn_sigmaa_omegaa( const HKL_data<T>& eo, const HKL_data<T>& ec )
 
366
  {
 
367
    eo_data = &eo;
 
368
    ec_data = &ec;
 
369
  }
 
370
 
 
371
  template<class T> TargetFn_base::Rderiv TargetFn_sigmaa_omegaa<T>::rderiv( const HKL_info::HKL_reference_index& ih, const ftype& omegaa ) const
 
372
  {
 
373
    Rderiv result;
 
374
    
 
375
    const HKL_data<T>& eodata = *eo_data;
 
376
    const HKL_data<T>& ecdata = *ec_data;
 
377
    if ( eodata[ih].missing() || ecdata[ih].missing()  ) {
 
378
      result.r = result.dr = result.dr2 = 0.0;
 
379
    } else {
 
380
      ftype eo = eodata[ih].E();
 
381
      ftype ec = ecdata[ih].E();
 
382
      ftype omeg = (omegaa>0.05) ? omegaa : (0.05*exp(omegaa/0.05-1.0));
 
383
      ftype sigmaa = 0.5*(sqrt(4*omeg*omeg+1)-1)/omeg;
 
384
      ftype dx = 2.0 * eo * ec;
 
385
      ftype x = dx * omeg;
 
386
      ftype t0 = 1.0/(1-sigmaa*sigmaa) + 0.5*log((1-sigmaa*sigmaa));
 
387
      ftype t1 = sigmaa;
 
388
      ftype t2 = pow(1-sigmaa*sigmaa,2)/(1+sigmaa*sigmaa);
 
389
      if ( ih.hkl_class().centric() ) {
 
390
        result.r = 1.0*t0 - log( cosh( x/2 ) );
 
391
        result.dr = 1.0*t1 - dx*0.5*tanh( x/2 );
 
392
        result.dr2 = 1.0*t2 - dx*dx*0.25*(1.0 - pow(tanh(x/2),2) );
 
393
      } else {
 
394
        result.r = 2.0*t0 - Util::sim_integ( x );
 
395
        result.dr = 2.0*t1 - dx*Util::sim( x );
 
396
        result.dr2 = 2.0*t2 - dx*dx*Util::sim_deriv( x );
 
397
      }
 
398
      if ( omegaa < 0.05 ) {
 
399
        ftype dy = exp( omegaa/0.05 ) / exp( 1.0 );
 
400
        ftype dy2 = exp( omegaa/0.05 ) / ( 0.05*exp( 1.0 ) );
 
401
        result.dr2 = result.dr*dy2 + result.dr2*dy*dy;
 
402
        result.dr = result.dr*dy;
 
403
      }
 
404
    }
 
405
    return result;
 
406
  }
 
407
 
 
408
 
 
409
  // sigmaa (norm)
 
410
 
 
411
  template<class T> TargetFn_sigmaa<T>::TargetFn_sigmaa( const HKL_data<T>& eo, const HKL_data<T>& ec )
 
412
  {
 
413
    eo_data = &eo;
 
414
    ec_data = &ec;
 
415
  }
 
416
 
 
417
  template<class T> TargetFn_base::Rderiv TargetFn_sigmaa<T>::rderiv( const HKL_info::HKL_reference_index& ih, const ftype& sigmaa0 ) const
 
418
  {
 
419
    Rderiv result;
 
420
    
 
421
    const HKL_data<T>& eodata = *eo_data;
 
422
    const HKL_data<T>& ecdata = *ec_data;
 
423
    if ( eodata[ih].missing() || ecdata[ih].missing()  ) {
 
424
      result.r = result.dr = result.dr2 = 0.0;
 
425
    } else {
 
426
      ftype eo = eodata[ih].E();
 
427
      ftype ec = ecdata[ih].E();
 
428
      ftype sigmaa = Util::min( Util::max( sigmaa0, 0.01 ), 0.99 );
 
429
      ftype dx = 2.0 * eo * ec;
 
430
      ftype x = dx * sigmaa/(1-sigmaa*sigmaa);
 
431
      ftype t0 = 1.0/(1-sigmaa*sigmaa) + 0.5*log((1-sigmaa*sigmaa));
 
432
      ftype t1 = sigmaa;
 
433
      ftype t2 = pow(1-sigmaa*sigmaa,2)/(1+sigmaa*sigmaa);
 
434
      if ( ih.hkl_class().centric() ) {
 
435
        result.r = 1.0*t0 - log( cosh( x/2 ) );
 
436
        result.dr = 1.0*t1 - dx*0.5*tanh( x/2 );
 
437
        result.dr2 = 1.0*t2 - dx*dx*0.25*(1.0 - pow(tanh(x/2),2) );
 
438
      } else {
 
439
        result.r = 2.0*t0 - Util::sim_integ( x );
 
440
        result.dr = 2.0*t1 - dx*Util::sim( x );
 
441
        result.dr2 = 2.0*t2 - dx*dx*Util::sim_deriv( x );
 
442
      }
 
443
      ftype ds = (1+sigmaa*sigmaa)/pow(1-sigmaa*sigmaa,2);
 
444
      ftype ds2 = 2*sigmaa*(3+sigmaa*sigmaa)/pow(1-sigmaa*sigmaa,3);
 
445
      result.dr2 = result.dr*ds2 + result.dr2*ds*ds;
 
446
      result.dr = result.dr*ds;
 
447
    }
 
448
    return result;
 
449
  }
 
450
 
 
451
 
 
452
} // namespace clipper
 
453
 
 
454
#endif