Clipper
resol_targetfn.h
1
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
52namespace clipper {
53
54
56
64 template<class T> class TargetFn_meanFnth : public TargetFn_base
65 {
66 public:
68 TargetFn_meanFnth( const HKL_data<T>& hkl_data_, const ftype& n );
70 Rderiv rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const;
72 FNtype type() const { return QUADRATIC; }
73 private:
74 ftype power;
75 const HKL_data<T>* hkl_data;
76 };
77
78
80
88 template<class T> class TargetFn_meanEnth : public TargetFn_base
89 {
90 public:
92 TargetFn_meanEnth( const HKL_data<T>& hkl_data_, const ftype& n );
94 Rderiv rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const;
96 FNtype type() const { return QUADRATIC; }
97 private:
98 ftype power;
99 const HKL_data<T>* hkl_data;
100 };
101
102
104
108 template<class T1, class T2> class TargetFn_scaleF1F2 : public TargetFn_base
109 {
110 public:
112 TargetFn_scaleF1F2( const HKL_data<T1>& hkl_data1_, const HKL_data<T2>& hkl_data2_ );
114 Rderiv rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const;
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
124
132 template<class T1, class T2> class TargetFn_scaleLogF1F2 : public TargetFn_base
133 {
134 public:
136 TargetFn_scaleLogF1F2( const HKL_data<T1>& hkl_data1_, const HKL_data<T2>& hkl_data2_ );
138 Rderiv rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const;
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
151 template<class T1, class T2> class TargetFn_scaleI1I2 : public TargetFn_base
152 {
153 public:
155 TargetFn_scaleI1I2( const HKL_data<T1>& hkl_data1_, const HKL_data<T2>& hkl_data2_ );
157 Rderiv rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const;
159 FNtype type() const { return QUADRATIC; }
160 private:
161 const HKL_data<T1>* hkl_data1;
162 const HKL_data<T2>* hkl_data2;
163 };
164
165
167
175 template<class T1, class T2> class TargetFn_scaleLogI1I2 : public TargetFn_base
176 {
177 public:
179 TargetFn_scaleLogI1I2( const HKL_data<T1>& hkl_data1_, const HKL_data<T2>& hkl_data2_ );
181 Rderiv rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const;
183 FNtype type() const { return QUADRATIC; }
184 private:
185 const HKL_data<T1>* hkl_data1;
186 const HKL_data<T2>* hkl_data2;
187 };
188
189
191
197 template<class T> class TargetFn_scaleEsq : public TargetFn_base
198 {
199 public:
201 TargetFn_scaleEsq( const HKL_data<T>& hkl_data_ );
203 Rderiv rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const;
205 FNtype type() const { return QUADRATIC; }
206 private:
207 const HKL_data<T>* hkl_data;
208 };
209
210
212
229 template<class T> class TargetFn_sigmaa_omegaa : public TargetFn_base
230 {
231 public:
233 TargetFn_sigmaa_omegaa( const HKL_data<T>& eo, const HKL_data<T>& ec );
235 Rderiv rderiv( const HKL_info::HKL_reference_index& ih, const ftype& omegaa ) const;
237 static ftype sigmaa( const ftype& omegaa )
238 {
239 ftype omeg = (omegaa>0.05) ? omegaa : (0.05*exp(omegaa/0.05-1.0));
240 return 0.5 * ( sqrt( 4.0*omeg*omeg + 1.0 ) - 1.0 ) / omeg;
241 }
242 private:
243 const HKL_data<T>* eo_data;
244 const HKL_data<T>* ec_data;
245 };
246
247
249
261 template<class T> class TargetFn_sigmaa : public TargetFn_base
262 {
263 public:
265 TargetFn_sigmaa( const HKL_data<T>& eo, const HKL_data<T>& ec );
267 Rderiv rderiv( const HKL_info::HKL_reference_index& ih, const ftype& sigmaa0 ) const;
269 static ftype sigmaa( const ftype& sigm ) { return sigm; }
270 private:
271 const HKL_data<T>* eo_data;
272 const HKL_data<T>* ec_data;
273 };
274
275
276
277 // implementations for template functions
278
279 // mean F^nth
280
281 template<class T> TargetFn_meanFnth<T>::TargetFn_meanFnth( const HKL_data<T>& hkl_data_, const ftype& n )
282 {
283 power = n;
284 hkl_data = &hkl_data_;
285 }
286
288 {
289 Rderiv result;
290 const HKL_data<T>& data = *hkl_data;
291 if ( !data[ih].missing() ) {
292 ftype d = fh - pow( ftype(data[ih].f()) / sqrt(ih.hkl_class().epsilon()),
293 power );
294 result.r = d * d;
295 result.dr = 2.0 * d;
296 result.dr2 = 2.0;
297 } else {
298 result.r = result.dr = result.dr2 = 0.0;
299 }
300 return result;
301 }
302
303 // mean E^nth
304
305 template<class T> TargetFn_meanEnth<T>::TargetFn_meanEnth( const HKL_data<T>& hkl_data_, const ftype& n )
306 {
307 power = n;
308 hkl_data = &hkl_data_;
309 }
310
312 {
313 Rderiv result;
314 const HKL_data<T>& data = *hkl_data;
315 if ( !data[ih].missing() ) {
316 ftype d = fh - pow( ftype(data[ih].E()), power );
317 result.r = d * d;
318 result.dr = 2.0 * d;
319 result.dr2 = 2.0;
320 } else {
321 result.r = result.dr = result.dr2 = 0.0;
322 }
323 return result;
324 }
325
326 // F1-F2 scaling
327
328 template<class T1, class T2> TargetFn_scaleF1F2<T1,T2>::TargetFn_scaleF1F2( const HKL_data<T1>& hkl_data1_, const HKL_data<T2>& hkl_data2_ )
329 {
330 hkl_data1 = &hkl_data1_;
331 hkl_data2 = &hkl_data2_;
332 }
333
334 template<class T1, class T2> TargetFn_base::Rderiv TargetFn_scaleF1F2<T1,T2>::rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const
335 {
336 Rderiv result;
337 const T1& ft1 = (*hkl_data1)[ih];
338 const T2& ft2 = (*hkl_data2)[ih];
339 if ( !ft1.missing() && !ft2.missing() ) {
340 const ftype eps = ih.hkl_class().epsilon();
341 const ftype f1 = pow( ft1.f(), 2 ) / eps;
342 const ftype f2 = pow( ft2.f(), 2 ) / eps;
343 const ftype d = fh*f1 - f2;
344 result.r = d * d / f1;
345 result.dr = 2.0 * d;
346 result.dr2 = 2.0 * f1;
347 } else {
348 result.r = result.dr = result.dr2 = 0.0;
349 }
350 return result;
351 }
352
353 // Log F1-F2 scaling
354
355 template<class T1, class T2> TargetFn_scaleLogF1F2<T1,T2>::TargetFn_scaleLogF1F2( const HKL_data<T1>& hkl_data1_, const HKL_data<T2>& hkl_data2_ )
356 {
357 hkl_data1 = &hkl_data1_;
358 hkl_data2 = &hkl_data2_;
359 }
360
362 {
363 Rderiv result;
364 result.r = result.dr = result.dr2 = 0.0;
365 const T1& ft1 = (*hkl_data1)[ih];
366 const T2& ft2 = (*hkl_data2)[ih];
367 if ( !ft1.missing() && !ft2.missing() )
368 if ( ft1.f() > 1.0e-6 && ft2.f() > 1.0e-6 ) {
369 const ftype eps = ih.hkl_class().epsilon();
370 const ftype f1 = pow( ft1.f(), 2 ) / eps;
371 const ftype f2 = pow( ft2.f(), 2 ) / eps;
372 const ftype w = 1.0; // f1 * f2;
373 const ftype d = fh + log(f1) - log(f2);
374 result.r = w * d * d;
375 result.dr = 2.0 * w * d;
376 result.dr2 = 2.0 * w;
377 }
378 return result;
379 }
380
381 // I1-I2 scaling
382
383 template<class T1, class T2> TargetFn_scaleI1I2<T1,T2>::TargetFn_scaleI1I2( const HKL_data<T1>& hkl_data1_, const HKL_data<T2>& hkl_data2_ )
384 {
385 hkl_data1 = &hkl_data1_;
386 hkl_data2 = &hkl_data2_;
387 }
388
389 template<class T1, class T2> TargetFn_base::Rderiv TargetFn_scaleI1I2<T1,T2>::rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const
390 {
391 Rderiv result;
392 const T1& ft1 = (*hkl_data1)[ih];
393 const T2& ft2 = (*hkl_data2)[ih];
394 if ( !ft1.missing() && !ft2.missing() ) {
395 const ftype eps = ih.hkl_class().epsilon();
396 const ftype f1 = ft1.I() / eps;
397 const ftype f2 = ft2.I() / eps;
398 const ftype d = fh*f1 - f2;
399 result.r = d * d / f1;
400 result.dr = 2.0 * d;
401 result.dr2 = 2.0 * f1;
402 } else {
403 result.r = result.dr = result.dr2 = 0.0;
404 }
405 return result;
406 }
407
408 // Log I1-I2 scaling
409
410 template<class T1, class T2> TargetFn_scaleLogI1I2<T1,T2>::TargetFn_scaleLogI1I2( const HKL_data<T1>& hkl_data1_, const HKL_data<T2>& hkl_data2_ )
411 {
412 hkl_data1 = &hkl_data1_;
413 hkl_data2 = &hkl_data2_;
414 }
415
417 {
418 Rderiv result;
419 result.r = result.dr = result.dr2 = 0.0;
420 const T1& ft1 = (*hkl_data1)[ih];
421 const T2& ft2 = (*hkl_data2)[ih];
422 if ( !ft1.missing() && !ft2.missing() )
423 if ( ft1.I() > 1.0e-6 && ft2.I() > 1.0e-6 ) {
424 const ftype eps = ih.hkl_class().epsilon();
425 const ftype f1 = ft1.I() / eps;
426 const ftype f2 = ft2.I() / eps;
427 const ftype w = 1.0; // f1 * f2;
428 const ftype d = fh + log(f1) - log(f2);
429 result.r = w * d * d;
430 result.dr = 2.0 * w * d;
431 result.dr2 = 2.0 * w;
432 }
433 return result;
434 }
435
436 // E^2 scaling
437
438 template<class T> TargetFn_scaleEsq<T>::TargetFn_scaleEsq( const HKL_data<T>& hkl_data_ )
439 {
440 hkl_data = &hkl_data_;
441 }
442
444 {
445 Rderiv result;
446 const HKL_data<T>& data = *hkl_data;
447 const ftype two(2.0);
448 if ( !data[ih].missing() ) {
449 ftype fsq = pow( ftype(data[ih].E()), two );
450 ftype d = fsq * fh - 1.0;
451 result.r = d * d / fsq;
452 result.dr = two * d;
453 result.dr2 = two * fsq;
454 } else {
455 result.r = result.dr = result.dr2 = 0.0;
456 }
457 return result;
458 }
459
460
461 // sigmaa (omegaa)
462
464 {
465 eo_data = &eo;
466 ec_data = &ec;
467 }
468
470 {
471 Rderiv result;
472
473 const HKL_data<T>& eodata = *eo_data;
474 const HKL_data<T>& ecdata = *ec_data;
475 if ( eodata[ih].missing() || ecdata[ih].missing() ) {
476 result.r = result.dr = result.dr2 = 0.0;
477 } else {
478 ftype eo = eodata[ih].E();
479 ftype ec = ecdata[ih].E();
480 ftype omeg = (omegaa>0.05) ? omegaa : (0.05*exp(omegaa/0.05-1.0));
481 ftype sigmaa = 0.5*(sqrt(4*omeg*omeg+1)-1)/omeg;
482 ftype dx = 2.0 * eo * ec;
483 ftype x = dx * omeg;
484 ftype t0 = 1.0/(1-sigmaa*sigmaa) + 0.5*log((1-sigmaa*sigmaa));
485 ftype t1 = sigmaa;
486 ftype t2 = pow(1-sigmaa*sigmaa,2)/(1+sigmaa*sigmaa);
487 if ( ih.hkl_class().centric() ) {
488 result.r = 1.0*t0 - log( cosh( x/2 ) );
489 result.dr = 1.0*t1 - dx*0.5*tanh( x/2 );
490 result.dr2 = 1.0*t2 - dx*dx*0.25*(1.0 - pow(tanh(x/2),2) );
491 } else {
492 result.r = 2.0*t0 - Util::sim_integ( x );
493 result.dr = 2.0*t1 - dx*Util::sim( x );
494 result.dr2 = 2.0*t2 - dx*dx*Util::sim_deriv( x );
495 }
496 if ( omegaa < 0.05 ) {
497 ftype dy = exp( omegaa/0.05 ) / exp( 1.0 );
498 ftype dy2 = exp( omegaa/0.05 ) / ( 0.05*exp( 1.0 ) );
499 result.dr2 = result.dr*dy2 + result.dr2*dy*dy;
500 result.dr = result.dr*dy;
501 }
502 }
503 return result;
504 }
505
506
507 // sigmaa (norm)
508
509 template<class T> TargetFn_sigmaa<T>::TargetFn_sigmaa( const HKL_data<T>& eo, const HKL_data<T>& ec )
510 {
511 eo_data = &eo;
512 ec_data = &ec;
513 }
514
515 template<class T> TargetFn_base::Rderiv TargetFn_sigmaa<T>::rderiv( const HKL_info::HKL_reference_index& ih, const ftype& sigmaa0 ) const
516 {
517 Rderiv result;
518
519 const HKL_data<T>& eodata = *eo_data;
520 const HKL_data<T>& ecdata = *ec_data;
521 if ( eodata[ih].missing() || ecdata[ih].missing() ) {
522 result.r = result.dr = result.dr2 = 0.0;
523 } else {
524 ftype eo = eodata[ih].E();
525 ftype ec = ecdata[ih].E();
526 ftype sigmaa = (sigmaa0>0.99)?(0.99):((sigmaa0<0.01)?0.01:sigmaa0);
527 ftype dx = 2.0 * eo * ec;
528 ftype x = dx * sigmaa/(1-sigmaa*sigmaa);
529 ftype t0 = 1.0/(1-sigmaa*sigmaa) + 0.5*log((1-sigmaa*sigmaa));
530 ftype t1 = sigmaa;
531 ftype t2 = pow(1-sigmaa*sigmaa,2)/(1+sigmaa*sigmaa);
532 if ( ih.hkl_class().centric() ) {
533 result.r = 1.0*t0 - log( cosh( x/2 ) );
534 result.dr = 1.0*t1 - dx*0.5*tanh( x/2 );
535 result.dr2 = 1.0*t2 - dx*dx*0.25*(1.0 - pow(tanh(x/2),2) );
536 } else {
537 result.r = 2.0*t0 - Util::sim_integ( x );
538 result.dr = 2.0*t1 - dx*Util::sim( x );
539 result.dr2 = 2.0*t2 - dx*dx*Util::sim_deriv( x );
540 }
541 ftype ds = (1+sigmaa*sigmaa)/pow(1-sigmaa*sigmaa,2);
542 ftype ds2 = 2*sigmaa*(3+sigmaa*sigmaa)/pow(1-sigmaa*sigmaa,3);
543 result.dr2 = result.dr*ds2 + result.dr2*ds*ds;
544 result.dr = result.dr*ds;
545 }
546 return result;
547 }
548
549
550} // namespace clipper
551
552#endif
bool centric() const
is centric?
Definition: coords.h:101
ftype epsilon() const
get epsilon
Definition: coords.h:94
HKL_data<>
Definition: hkl_data.h:235
HKL reference with index-like behaviour.
Definition: hkl_info.h:152
const HKL_class & hkl_class() const
return the reflection class for the reflection
Definition: hkl_info.h:162
object holding the residual function and first two derivatives
Definition: resol_fn.h:119
ftype r
the value of the function
Definition: resol_fn.h:121
ftype dr
first derivative w.r.t basis fn
Definition: resol_fn.h:122
ftype dr2
second derivative w.r.t basis fn
Definition: resol_fn.h:123
abstract base class for least-squares resolution function target functions
Definition: resol_fn.h:115
FNtype
enumeration of function types: optionally used to improve convergence
Definition: resol_fn.h:127
simple mean |E|n target
Definition: resol_targetfn.h:89
TargetFn_meanEnth(const HKL_data< T > &hkl_data_, const ftype &n)
constructor: takes the datalist against which to calc target, and power
Definition: resol_targetfn.h:305
FNtype type() const
the type of the function: optionally used to improve convergence
Definition: resol_targetfn.h:96
Rderiv rderiv(const HKL_info::HKL_reference_index &ih, const ftype &fh) const
return the value and derivatives of the target function
Definition: resol_targetfn.h:311
simple mean |F|n target
Definition: resol_targetfn.h:65
FNtype type() const
the type of the function: optionally used to improve convergence
Definition: resol_targetfn.h:72
Rderiv rderiv(const HKL_info::HKL_reference_index &ih, const ftype &fh) const
return the value and derivatives of the target function
Definition: resol_targetfn.h:287
TargetFn_meanFnth(const HKL_data< T > &hkl_data_, const ftype &n)
constructor: takes the datalist against which to calc target, and power
Definition: resol_targetfn.h:281
|E|2 scaling target
Definition: resol_targetfn.h:198
Rderiv rderiv(const HKL_info::HKL_reference_index &ih, const ftype &fh) const
return the value and derivatives of the target function
Definition: resol_targetfn.h:443
TargetFn_scaleEsq(const HKL_data< T > &hkl_data_)
constructor: takes the datalist against which to calc target
Definition: resol_targetfn.h:438
FNtype type() const
the type of the function: optionally used to improve convergence
Definition: resol_targetfn.h:205
|F|2 scaling target
Definition: resol_targetfn.h:109
Rderiv rderiv(const HKL_info::HKL_reference_index &ih, const ftype &fh) const
return the value and derivatives of the target function
Definition: resol_targetfn.h:334
TargetFn_scaleF1F2(const HKL_data< T1 > &hkl_data1_, const HKL_data< T2 > &hkl_data2_)
constructor: takes the datalist against which to calc target
Definition: resol_targetfn.h:328
FNtype type() const
the type of the function: optionally used to improve convergence
Definition: resol_targetfn.h:116
Definition: resol_targetfn.h:152
FNtype type() const
the type of the function: optionally used to improve convergence
Definition: resol_targetfn.h:159
Rderiv rderiv(const HKL_info::HKL_reference_index &ih, const ftype &fh) const
return the value and derivatives of the target function
Definition: resol_targetfn.h:389
TargetFn_scaleI1I2(const HKL_data< T1 > &hkl_data1_, const HKL_data< T2 > &hkl_data2_)
constructor: takes the datalist against which to calc target
Definition: resol_targetfn.h:383
log |F|2 scaling target
Definition: resol_targetfn.h:133
Rderiv rderiv(const HKL_info::HKL_reference_index &ih, const ftype &fh) const
return the value and derivatives of the target function
Definition: resol_targetfn.h:361
FNtype type() const
the type of the function: optionally used to improve convergence
Definition: resol_targetfn.h:140
TargetFn_scaleLogF1F2(const HKL_data< T1 > &hkl_data1_, const HKL_data< T2 > &hkl_data2_)
constructor: takes the datalist against which to calc target
Definition: resol_targetfn.h:355
log |I| scaling target
Definition: resol_targetfn.h:176
TargetFn_scaleLogI1I2(const HKL_data< T1 > &hkl_data1_, const HKL_data< T2 > &hkl_data2_)
constructor: takes the datalist against which to calc target
Definition: resol_targetfn.h:410
Rderiv rderiv(const HKL_info::HKL_reference_index &ih, const ftype &fh) const
return the value and derivatives of the target function
Definition: resol_targetfn.h:416
FNtype type() const
the type of the function: optionally used to improve convergence
Definition: resol_targetfn.h:183
Definition: resol_targetfn.h:230
Rderiv rderiv(const HKL_info::HKL_reference_index &ih, const ftype &omegaa) const
return the value and derivatives of the target function
Definition: resol_targetfn.h:469
TargetFn_sigmaa_omegaa(const HKL_data< T > &eo, const HKL_data< T > &ec)
constructor: takes the datalist against which to calc target
Definition: resol_targetfn.h:463
static ftype sigmaa(const ftype &omegaa)
convert omegaa to sigmaa
Definition: resol_targetfn.h:237
Definition: resol_targetfn.h:262
TargetFn_sigmaa(const HKL_data< T > &eo, const HKL_data< T > &ec)
constructor: takes the datalist against which to calc target
Definition: resol_targetfn.h:509
Rderiv rderiv(const HKL_info::HKL_reference_index &ih, const ftype &sigmaa0) const
return the value and derivatives of the target function
Definition: resol_targetfn.h:515
static ftype sigmaa(const ftype &sigm)
convert function to sigmaa
Definition: resol_targetfn.h:269
static ftype sim_integ(const ftype &x)
Integral of Sim function: log(I0(X))
Definition: clipper_util.cpp:112
static ftype sim_deriv(const ftype &x)
Derivative of Sim function: d/dx( I1(X)/I0(x) )
Definition: clipper_util.cpp:119
static ftype sim(const ftype &x)
Sim function: I1(X)/I0(X)
Definition: clipper_util.cpp:84
ftype64 ftype
ftype definition for floating point representation
Definition: clipper_precision.h:58