00001
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045 #ifndef CLIPPER_RESOL_TARGETFN
00046 #define CLIPPER_RESOL_TARGETFN
00047
00048 #include "resol_basisfn.h"
00049 #include "hkl_datatypes.h"
00050
00051
00052 namespace clipper {
00053
00054
00056
00064 template<class T> class TargetFn_meanFnth : public TargetFn_base
00065 {
00066 public:
00068 TargetFn_meanFnth( const HKL_data<T>& hkl_data_, const ftype& n );
00070 Rderiv rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const;
00072 FNtype type() const { return QUADRATIC; }
00073 private:
00074 ftype power;
00075 const HKL_data<T>* hkl_data;
00076 };
00077
00078
00080
00088 template<class T> class TargetFn_meanEnth : public TargetFn_base
00089 {
00090 public:
00092 TargetFn_meanEnth( const HKL_data<T>& hkl_data_, const ftype& n );
00094 Rderiv rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const;
00096 FNtype type() const { return QUADRATIC; }
00097 private:
00098 ftype power;
00099 const HKL_data<T>* hkl_data;
00100 };
00101
00102
00104
00108 template<class T1, class T2> class TargetFn_scaleF1F2 : public TargetFn_base
00109 {
00110 public:
00112 TargetFn_scaleF1F2( const HKL_data<T1>& hkl_data1_, const HKL_data<T2>& hkl_data2_ );
00114 Rderiv rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const;
00116 FNtype type() const { return QUADRATIC; }
00117 private:
00118 const HKL_data<T1>* hkl_data1;
00119 const HKL_data<T2>* hkl_data2;
00120 };
00121
00122
00124
00132 template<class T1, class T2> class TargetFn_scaleLogF1F2 : public TargetFn_base
00133 {
00134 public:
00136 TargetFn_scaleLogF1F2( const HKL_data<T1>& hkl_data1_, const HKL_data<T2>& hkl_data2_ );
00138 Rderiv rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const;
00140 FNtype type() const { return QUADRATIC; }
00141 private:
00142 const HKL_data<T1>* hkl_data1;
00143 const HKL_data<T2>* hkl_data2;
00144 };
00145
00146
00151 template<class T1, class T2> class TargetFn_scaleI1I2 : public TargetFn_base
00152 {
00153 public:
00155 TargetFn_scaleI1I2( const HKL_data<T1>& hkl_data1_, const HKL_data<T2>& hkl_data2_ );
00157 Rderiv rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const;
00159 FNtype type() const { return QUADRATIC; }
00160 private:
00161 const HKL_data<T1>* hkl_data1;
00162 const HKL_data<T2>* hkl_data2;
00163 };
00164
00165
00167
00175 template<class T1, class T2> class TargetFn_scaleLogI1I2 : public TargetFn_base
00176 {
00177 public:
00179 TargetFn_scaleLogI1I2( const HKL_data<T1>& hkl_data1_, const HKL_data<T2>& hkl_data2_ );
00181 Rderiv rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const;
00183 FNtype type() const { return QUADRATIC; }
00184 private:
00185 const HKL_data<T1>* hkl_data1;
00186 const HKL_data<T2>* hkl_data2;
00187 };
00188
00189
00191
00197 template<class T> class TargetFn_scaleEsq : public TargetFn_base
00198 {
00199 public:
00201 TargetFn_scaleEsq( const HKL_data<T>& hkl_data_ );
00203 Rderiv rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const;
00205 FNtype type() const { return QUADRATIC; }
00206 private:
00207 const HKL_data<T>* hkl_data;
00208 };
00209
00210
00212
00229 template<class T> class TargetFn_sigmaa_omegaa : public TargetFn_base
00230 {
00231 public:
00233 TargetFn_sigmaa_omegaa( const HKL_data<T>& eo, const HKL_data<T>& ec );
00235 Rderiv rderiv( const HKL_info::HKL_reference_index& ih, const ftype& omegaa ) const;
00237 static ftype sigmaa( const ftype& omegaa )
00238 {
00239 ftype omeg = (omegaa>0.05) ? omegaa : (0.05*exp(omegaa/0.05-1.0));
00240 return 0.5 * ( sqrt( 4.0*omeg*omeg + 1.0 ) - 1.0 ) / omeg;
00241 }
00242 private:
00243 const HKL_data<T>* eo_data;
00244 const HKL_data<T>* ec_data;
00245 };
00246
00247
00249
00261 template<class T> class TargetFn_sigmaa : public TargetFn_base
00262 {
00263 public:
00265 TargetFn_sigmaa( const HKL_data<T>& eo, const HKL_data<T>& ec );
00267 Rderiv rderiv( const HKL_info::HKL_reference_index& ih, const ftype& sigmaa0 ) const;
00269 static ftype sigmaa( const ftype& sigm ) { return sigm; }
00270 private:
00271 const HKL_data<T>* eo_data;
00272 const HKL_data<T>* ec_data;
00273 };
00274
00275
00276
00277
00278
00279
00280
00281 template<class T> TargetFn_meanFnth<T>::TargetFn_meanFnth( const HKL_data<T>& hkl_data_, const ftype& n )
00282 {
00283 power = n;
00284 hkl_data = &hkl_data_;
00285 }
00286
00287 template<class T> TargetFn_base::Rderiv TargetFn_meanFnth<T>::rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const
00288 {
00289 Rderiv result;
00290 const HKL_data<T>& data = *hkl_data;
00291 if ( !data[ih].missing() ) {
00292 ftype d = fh - pow( ftype(data[ih].f()) / sqrt(ih.hkl_class().epsilon()),
00293 power );
00294 result.r = d * d;
00295 result.dr = 2.0 * d;
00296 result.dr2 = 2.0;
00297 } else {
00298 result.r = result.dr = result.dr2 = 0.0;
00299 }
00300 return result;
00301 }
00302
00303
00304
00305 template<class T> TargetFn_meanEnth<T>::TargetFn_meanEnth( const HKL_data<T>& hkl_data_, const ftype& n )
00306 {
00307 power = n;
00308 hkl_data = &hkl_data_;
00309 }
00310
00311 template<class T> TargetFn_base::Rderiv TargetFn_meanEnth<T>::rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const
00312 {
00313 Rderiv result;
00314 const HKL_data<T>& data = *hkl_data;
00315 if ( !data[ih].missing() ) {
00316 ftype d = fh - pow( ftype(data[ih].E()), power );
00317 result.r = d * d;
00318 result.dr = 2.0 * d;
00319 result.dr2 = 2.0;
00320 } else {
00321 result.r = result.dr = result.dr2 = 0.0;
00322 }
00323 return result;
00324 }
00325
00326
00327
00328 template<class T1, class T2> TargetFn_scaleF1F2<T1,T2>::TargetFn_scaleF1F2( const HKL_data<T1>& hkl_data1_, const HKL_data<T2>& hkl_data2_ )
00329 {
00330 hkl_data1 = &hkl_data1_;
00331 hkl_data2 = &hkl_data2_;
00332 }
00333
00334 template<class T1, class T2> TargetFn_base::Rderiv TargetFn_scaleF1F2<T1,T2>::rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const
00335 {
00336 Rderiv result;
00337 const T1& ft1 = (*hkl_data1)[ih];
00338 const T2& ft2 = (*hkl_data2)[ih];
00339 if ( !ft1.missing() && !ft2.missing() ) {
00340 const ftype eps = ih.hkl_class().epsilon();
00341 const ftype f1 = pow( ft1.f(), 2 ) / eps;
00342 const ftype f2 = pow( ft2.f(), 2 ) / eps;
00343 const ftype d = fh*f1 - f2;
00344 result.r = d * d / f1;
00345 result.dr = 2.0 * d;
00346 result.dr2 = 2.0 * f1;
00347 } else {
00348 result.r = result.dr = result.dr2 = 0.0;
00349 }
00350 return result;
00351 }
00352
00353
00354
00355 template<class T1, class T2> TargetFn_scaleLogF1F2<T1,T2>::TargetFn_scaleLogF1F2( const HKL_data<T1>& hkl_data1_, const HKL_data<T2>& hkl_data2_ )
00356 {
00357 hkl_data1 = &hkl_data1_;
00358 hkl_data2 = &hkl_data2_;
00359 }
00360
00361 template<class T1, class T2> TargetFn_base::Rderiv TargetFn_scaleLogF1F2<T1,T2>::rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const
00362 {
00363 Rderiv result;
00364 result.r = result.dr = result.dr2 = 0.0;
00365 const T1& ft1 = (*hkl_data1)[ih];
00366 const T2& ft2 = (*hkl_data2)[ih];
00367 if ( !ft1.missing() && !ft2.missing() )
00368 if ( ft1.f() > 1.0e-6 && ft2.f() > 1.0e-6 ) {
00369 const ftype eps = ih.hkl_class().epsilon();
00370 const ftype f1 = pow( ft1.f(), 2 ) / eps;
00371 const ftype f2 = pow( ft2.f(), 2 ) / eps;
00372 const ftype w = 1.0;
00373 const ftype d = fh + log(f1) - log(f2);
00374 result.r = w * d * d;
00375 result.dr = 2.0 * w * d;
00376 result.dr2 = 2.0 * w;
00377 }
00378 return result;
00379 }
00380
00381
00382
00383 template<class T1, class T2> TargetFn_scaleI1I2<T1,T2>::TargetFn_scaleI1I2( const HKL_data<T1>& hkl_data1_, const HKL_data<T2>& hkl_data2_ )
00384 {
00385 hkl_data1 = &hkl_data1_;
00386 hkl_data2 = &hkl_data2_;
00387 }
00388
00389 template<class T1, class T2> TargetFn_base::Rderiv TargetFn_scaleI1I2<T1,T2>::rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const
00390 {
00391 Rderiv result;
00392 const T1& ft1 = (*hkl_data1)[ih];
00393 const T2& ft2 = (*hkl_data2)[ih];
00394 if ( !ft1.missing() && !ft2.missing() ) {
00395 const ftype eps = ih.hkl_class().epsilon();
00396 const ftype f1 = ft1.I() / eps;
00397 const ftype f2 = ft2.I() / eps;
00398 const ftype d = fh*f1 - f2;
00399 result.r = d * d / f1;
00400 result.dr = 2.0 * d;
00401 result.dr2 = 2.0 * f1;
00402 } else {
00403 result.r = result.dr = result.dr2 = 0.0;
00404 }
00405 return result;
00406 }
00407
00408
00409
00410 template<class T1, class T2> TargetFn_scaleLogI1I2<T1,T2>::TargetFn_scaleLogI1I2( const HKL_data<T1>& hkl_data1_, const HKL_data<T2>& hkl_data2_ )
00411 {
00412 hkl_data1 = &hkl_data1_;
00413 hkl_data2 = &hkl_data2_;
00414 }
00415
00416 template<class T1, class T2> TargetFn_base::Rderiv TargetFn_scaleLogI1I2<T1,T2>::rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const
00417 {
00418 Rderiv result;
00419 result.r = result.dr = result.dr2 = 0.0;
00420 const T1& ft1 = (*hkl_data1)[ih];
00421 const T2& ft2 = (*hkl_data2)[ih];
00422 if ( !ft1.missing() && !ft2.missing() )
00423 if ( ft1.I() > 1.0e-6 && ft2.I() > 1.0e-6 ) {
00424 const ftype eps = ih.hkl_class().epsilon();
00425 const ftype f1 = ft1.I() / eps;
00426 const ftype f2 = ft2.I() / eps;
00427 const ftype w = 1.0;
00428 const ftype d = fh + log(f1) - log(f2);
00429 result.r = w * d * d;
00430 result.dr = 2.0 * w * d;
00431 result.dr2 = 2.0 * w;
00432 }
00433 return result;
00434 }
00435
00436
00437
00438 template<class T> TargetFn_scaleEsq<T>::TargetFn_scaleEsq( const HKL_data<T>& hkl_data_ )
00439 {
00440 hkl_data = &hkl_data_;
00441 }
00442
00443 template<class T> TargetFn_base::Rderiv TargetFn_scaleEsq<T>::rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const
00444 {
00445 Rderiv result;
00446 const HKL_data<T>& data = *hkl_data;
00447 const ftype two(2.0);
00448 if ( !data[ih].missing() ) {
00449 ftype fsq = pow( ftype(data[ih].E()), two );
00450 ftype d = fsq * fh - 1.0;
00451 result.r = d * d / fsq;
00452 result.dr = two * d;
00453 result.dr2 = two * fsq;
00454 } else {
00455 result.r = result.dr = result.dr2 = 0.0;
00456 }
00457 return result;
00458 }
00459
00460
00461
00462
00463 template<class T> TargetFn_sigmaa_omegaa<T>::TargetFn_sigmaa_omegaa( const HKL_data<T>& eo, const HKL_data<T>& ec )
00464 {
00465 eo_data = &eo;
00466 ec_data = &ec;
00467 }
00468
00469 template<class T> TargetFn_base::Rderiv TargetFn_sigmaa_omegaa<T>::rderiv( const HKL_info::HKL_reference_index& ih, const ftype& omegaa ) const
00470 {
00471 Rderiv result;
00472
00473 const HKL_data<T>& eodata = *eo_data;
00474 const HKL_data<T>& ecdata = *ec_data;
00475 if ( eodata[ih].missing() || ecdata[ih].missing() ) {
00476 result.r = result.dr = result.dr2 = 0.0;
00477 } else {
00478 ftype eo = eodata[ih].E();
00479 ftype ec = ecdata[ih].E();
00480 ftype omeg = (omegaa>0.05) ? omegaa : (0.05*exp(omegaa/0.05-1.0));
00481 ftype sigmaa = 0.5*(sqrt(4*omeg*omeg+1)-1)/omeg;
00482 ftype dx = 2.0 * eo * ec;
00483 ftype x = dx * omeg;
00484 ftype t0 = 1.0/(1-sigmaa*sigmaa) + 0.5*log((1-sigmaa*sigmaa));
00485 ftype t1 = sigmaa;
00486 ftype t2 = pow(1-sigmaa*sigmaa,2)/(1+sigmaa*sigmaa);
00487 if ( ih.hkl_class().centric() ) {
00488 result.r = 1.0*t0 - log( cosh( x/2 ) );
00489 result.dr = 1.0*t1 - dx*0.5*tanh( x/2 );
00490 result.dr2 = 1.0*t2 - dx*dx*0.25*(1.0 - pow(tanh(x/2),2) );
00491 } else {
00492 result.r = 2.0*t0 - Util::sim_integ( x );
00493 result.dr = 2.0*t1 - dx*Util::sim( x );
00494 result.dr2 = 2.0*t2 - dx*dx*Util::sim_deriv( x );
00495 }
00496 if ( omegaa < 0.05 ) {
00497 ftype dy = exp( omegaa/0.05 ) / exp( 1.0 );
00498 ftype dy2 = exp( omegaa/0.05 ) / ( 0.05*exp( 1.0 ) );
00499 result.dr2 = result.dr*dy2 + result.dr2*dy*dy;
00500 result.dr = result.dr*dy;
00501 }
00502 }
00503 return result;
00504 }
00505
00506
00507
00508
00509 template<class T> TargetFn_sigmaa<T>::TargetFn_sigmaa( const HKL_data<T>& eo, const HKL_data<T>& ec )
00510 {
00511 eo_data = &eo;
00512 ec_data = &ec;
00513 }
00514
00515 template<class T> TargetFn_base::Rderiv TargetFn_sigmaa<T>::rderiv( const HKL_info::HKL_reference_index& ih, const ftype& sigmaa0 ) const
00516 {
00517 Rderiv result;
00518
00519 const HKL_data<T>& eodata = *eo_data;
00520 const HKL_data<T>& ecdata = *ec_data;
00521 if ( eodata[ih].missing() || ecdata[ih].missing() ) {
00522 result.r = result.dr = result.dr2 = 0.0;
00523 } else {
00524 ftype eo = eodata[ih].E();
00525 ftype ec = ecdata[ih].E();
00526 ftype sigmaa = Util::min( Util::max( sigmaa0, 0.01 ), 0.99 );
00527 ftype dx = 2.0 * eo * ec;
00528 ftype x = dx * sigmaa/(1-sigmaa*sigmaa);
00529 ftype t0 = 1.0/(1-sigmaa*sigmaa) + 0.5*log((1-sigmaa*sigmaa));
00530 ftype t1 = sigmaa;
00531 ftype t2 = pow(1-sigmaa*sigmaa,2)/(1+sigmaa*sigmaa);
00532 if ( ih.hkl_class().centric() ) {
00533 result.r = 1.0*t0 - log( cosh( x/2 ) );
00534 result.dr = 1.0*t1 - dx*0.5*tanh( x/2 );
00535 result.dr2 = 1.0*t2 - dx*dx*0.25*(1.0 - pow(tanh(x/2),2) );
00536 } else {
00537 result.r = 2.0*t0 - Util::sim_integ( x );
00538 result.dr = 2.0*t1 - dx*Util::sim( x );
00539 result.dr2 = 2.0*t2 - dx*dx*Util::sim_deriv( x );
00540 }
00541 ftype ds = (1+sigmaa*sigmaa)/pow(1-sigmaa*sigmaa,2);
00542 ftype ds2 = 2*sigmaa*(3+sigmaa*sigmaa)/pow(1-sigmaa*sigmaa,3);
00543 result.dr2 = result.dr*ds2 + result.dr2*ds*ds;
00544 result.dr = result.dr*ds;
00545 }
00546 return result;
00547 }
00548
00549
00550 }
00551
00552 #endif