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
00148
00154 template<class T> class TargetFn_scaleEsq : public TargetFn_base
00155 {
00156 public:
00158 TargetFn_scaleEsq( const HKL_data<T>& hkl_data_ );
00160 Rderiv rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const;
00162 FNtype type() const { return QUADRATIC; }
00163 private:
00164 const HKL_data<T>* hkl_data;
00165 };
00166
00167
00169
00186 template<class T> class TargetFn_sigmaa_omegaa : public TargetFn_base
00187 {
00188 public:
00190 TargetFn_sigmaa_omegaa( const HKL_data<T>& eo, const HKL_data<T>& ec );
00192 Rderiv rderiv( const HKL_info::HKL_reference_index& ih, const ftype& omegaa ) const;
00194 static ftype sigmaa( const ftype& omegaa )
00195 {
00196 ftype omeg = (omegaa>0.05) ? omegaa : (0.05*exp(omegaa/0.05-1.0));
00197 return 0.5 * ( sqrt( 4.0*omeg*omeg + 1.0 ) - 1.0 ) / omeg;
00198 }
00199 private:
00200 const HKL_data<T>* eo_data;
00201 const HKL_data<T>* ec_data;
00202 };
00203
00204
00206
00218 template<class T> class TargetFn_sigmaa : public TargetFn_base
00219 {
00220 public:
00222 TargetFn_sigmaa( const HKL_data<T>& eo, const HKL_data<T>& ec );
00224 Rderiv rderiv( const HKL_info::HKL_reference_index& ih, const ftype& sigmaa0 ) const;
00226 static ftype sigmaa( const ftype& sigm ) { return sigm; }
00227 private:
00228 const HKL_data<T>* eo_data;
00229 const HKL_data<T>* ec_data;
00230 };
00231
00232
00233
00234
00235
00236
00237
00238 template<class T> TargetFn_meanFnth<T>::TargetFn_meanFnth( const HKL_data<T>& hkl_data_, const ftype& n )
00239 {
00240 power = n;
00241 hkl_data = &hkl_data_;
00242 }
00243
00244 template<class T> TargetFn_base::Rderiv TargetFn_meanFnth<T>::rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const
00245 {
00246 Rderiv result;
00247 const HKL_data<T>& data = *hkl_data;
00248 if ( !data[ih].missing() ) {
00249 ftype d = fh - pow( ftype(data[ih].f()) / sqrt(ih.hkl_class().epsilon()),
00250 power );
00251 result.r = d * d;
00252 result.dr = 2.0 * d;
00253 result.dr2 = 2.0;
00254 } else {
00255 result.r = result.dr = result.dr2 = 0.0;
00256 }
00257 return result;
00258 }
00259
00260
00261
00262 template<class T> TargetFn_meanEnth<T>::TargetFn_meanEnth( const HKL_data<T>& hkl_data_, const ftype& n )
00263 {
00264 power = n;
00265 hkl_data = &hkl_data_;
00266 }
00267
00268 template<class T> TargetFn_base::Rderiv TargetFn_meanEnth<T>::rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const
00269 {
00270 Rderiv result;
00271 const HKL_data<T>& data = *hkl_data;
00272 if ( !data[ih].missing() ) {
00273 ftype d = fh - pow( ftype(data[ih].E()), power );
00274 result.r = d * d;
00275 result.dr = 2.0 * d;
00276 result.dr2 = 2.0;
00277 } else {
00278 result.r = result.dr = result.dr2 = 0.0;
00279 }
00280 return result;
00281 }
00282
00283
00284
00285 template<class T1, class T2> TargetFn_scaleF1F2<T1,T2>::TargetFn_scaleF1F2( const HKL_data<T1>& hkl_data1_, const HKL_data<T2>& hkl_data2_ )
00286 {
00287 hkl_data1 = &hkl_data1_;
00288 hkl_data2 = &hkl_data2_;
00289 }
00290
00291 template<class T1, class T2> TargetFn_base::Rderiv TargetFn_scaleF1F2<T1,T2>::rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const
00292 {
00293 Rderiv result;
00294 const T1& ft1 = (*hkl_data1)[ih];
00295 const T2& ft2 = (*hkl_data2)[ih];
00296 if ( !ft1.missing() && !ft2.missing() ) {
00297 const ftype eps = ih.hkl_class().epsilon();
00298 const ftype f1 = pow( ft1.f(), 2 ) / eps;
00299 const ftype f2 = pow( ft2.f(), 2 ) / eps;
00300 const ftype d = fh*f1 - f2;
00301 result.r = d * d / f1;
00302 result.dr = 2.0 * d;
00303 result.dr2 = 2.0 * f1;
00304 } else {
00305 result.r = result.dr = result.dr2 = 0.0;
00306 }
00307 return result;
00308 }
00309
00310
00311
00312 template<class T1, class T2> TargetFn_scaleLogF1F2<T1,T2>::TargetFn_scaleLogF1F2( const HKL_data<T1>& hkl_data1_, const HKL_data<T2>& hkl_data2_ )
00313 {
00314 hkl_data1 = &hkl_data1_;
00315 hkl_data2 = &hkl_data2_;
00316 }
00317
00318 template<class T1, class T2> TargetFn_base::Rderiv TargetFn_scaleLogF1F2<T1,T2>::rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const
00319 {
00320 Rderiv result;
00321 result.r = result.dr = result.dr2 = 0.0;
00322 const T1& ft1 = (*hkl_data1)[ih];
00323 const T2& ft2 = (*hkl_data2)[ih];
00324 if ( !ft1.missing() && !ft2.missing() )
00325 if ( ft1.f() > 1.0e-6 && ft2.f() > 1.0e-6 ) {
00326 const ftype eps = ih.hkl_class().epsilon();
00327 const ftype f1 = pow( ft1.f(), 2 ) / eps;
00328 const ftype f2 = pow( ft2.f(), 2 ) / eps;
00329 const ftype w = 1.0;
00330 const ftype d = fh + log(f1) - log(f2);
00331 result.r = w * d * d;
00332 result.dr = 2.0 * w * d;
00333 result.dr2 = 2.0 * w;
00334 }
00335 return result;
00336 }
00337
00338
00339
00340 template<class T> TargetFn_scaleEsq<T>::TargetFn_scaleEsq( const HKL_data<T>& hkl_data_ )
00341 {
00342 hkl_data = &hkl_data_;
00343 }
00344
00345 template<class T> TargetFn_base::Rderiv TargetFn_scaleEsq<T>::rderiv( const HKL_info::HKL_reference_index& ih, const ftype& fh ) const
00346 {
00347 Rderiv result;
00348 const HKL_data<T>& data = *hkl_data;
00349 const ftype two(2.0);
00350 if ( !data[ih].missing() ) {
00351 ftype fsq = pow( ftype(data[ih].E()), two );
00352 ftype d = fsq * fh - 1.0;
00353 result.r = d * d / fsq;
00354 result.dr = two * d;
00355 result.dr2 = two * fsq;
00356 } else {
00357 result.r = result.dr = result.dr2 = 0.0;
00358 }
00359 return result;
00360 }
00361
00362
00363
00364
00365 template<class T> TargetFn_sigmaa_omegaa<T>::TargetFn_sigmaa_omegaa( const HKL_data<T>& eo, const HKL_data<T>& ec )
00366 {
00367 eo_data = &eo;
00368 ec_data = &ec;
00369 }
00370
00371 template<class T> TargetFn_base::Rderiv TargetFn_sigmaa_omegaa<T>::rderiv( const HKL_info::HKL_reference_index& ih, const ftype& omegaa ) const
00372 {
00373 Rderiv result;
00374
00375 const HKL_data<T>& eodata = *eo_data;
00376 const HKL_data<T>& ecdata = *ec_data;
00377 if ( eodata[ih].missing() || ecdata[ih].missing() ) {
00378 result.r = result.dr = result.dr2 = 0.0;
00379 } else {
00380 ftype eo = eodata[ih].E();
00381 ftype ec = ecdata[ih].E();
00382 ftype omeg = (omegaa>0.05) ? omegaa : (0.05*exp(omegaa/0.05-1.0));
00383 ftype sigmaa = 0.5*(sqrt(4*omeg*omeg+1)-1)/omeg;
00384 ftype dx = 2.0 * eo * ec;
00385 ftype x = dx * omeg;
00386 ftype t0 = 1.0/(1-sigmaa*sigmaa) + 0.5*log((1-sigmaa*sigmaa));
00387 ftype t1 = sigmaa;
00388 ftype t2 = pow(1-sigmaa*sigmaa,2)/(1+sigmaa*sigmaa);
00389 if ( ih.hkl_class().centric() ) {
00390 result.r = 1.0*t0 - log( cosh( x/2 ) );
00391 result.dr = 1.0*t1 - dx*0.5*tanh( x/2 );
00392 result.dr2 = 1.0*t2 - dx*dx*0.25*(1.0 - pow(tanh(x/2),2) );
00393 } else {
00394 result.r = 2.0*t0 - Util::sim_integ( x );
00395 result.dr = 2.0*t1 - dx*Util::sim( x );
00396 result.dr2 = 2.0*t2 - dx*dx*Util::sim_deriv( x );
00397 }
00398 if ( omegaa < 0.05 ) {
00399 ftype dy = exp( omegaa/0.05 ) / exp( 1.0 );
00400 ftype dy2 = exp( omegaa/0.05 ) / ( 0.05*exp( 1.0 ) );
00401 result.dr2 = result.dr*dy2 + result.dr2*dy*dy;
00402 result.dr = result.dr*dy;
00403 }
00404 }
00405 return result;
00406 }
00407
00408
00409
00410
00411 template<class T> TargetFn_sigmaa<T>::TargetFn_sigmaa( const HKL_data<T>& eo, const HKL_data<T>& ec )
00412 {
00413 eo_data = &eo;
00414 ec_data = &ec;
00415 }
00416
00417 template<class T> TargetFn_base::Rderiv TargetFn_sigmaa<T>::rderiv( const HKL_info::HKL_reference_index& ih, const ftype& sigmaa0 ) const
00418 {
00419 Rderiv result;
00420
00421 const HKL_data<T>& eodata = *eo_data;
00422 const HKL_data<T>& ecdata = *ec_data;
00423 if ( eodata[ih].missing() || ecdata[ih].missing() ) {
00424 result.r = result.dr = result.dr2 = 0.0;
00425 } else {
00426 ftype eo = eodata[ih].E();
00427 ftype ec = ecdata[ih].E();
00428 ftype sigmaa = Util::min( Util::max( sigmaa0, 0.01 ), 0.99 );
00429 ftype dx = 2.0 * eo * ec;
00430 ftype x = dx * sigmaa/(1-sigmaa*sigmaa);
00431 ftype t0 = 1.0/(1-sigmaa*sigmaa) + 0.5*log((1-sigmaa*sigmaa));
00432 ftype t1 = sigmaa;
00433 ftype t2 = pow(1-sigmaa*sigmaa,2)/(1+sigmaa*sigmaa);
00434 if ( ih.hkl_class().centric() ) {
00435 result.r = 1.0*t0 - log( cosh( x/2 ) );
00436 result.dr = 1.0*t1 - dx*0.5*tanh( x/2 );
00437 result.dr2 = 1.0*t2 - dx*dx*0.25*(1.0 - pow(tanh(x/2),2) );
00438 } else {
00439 result.r = 2.0*t0 - Util::sim_integ( x );
00440 result.dr = 2.0*t1 - dx*Util::sim( x );
00441 result.dr2 = 2.0*t2 - dx*dx*Util::sim_deriv( x );
00442 }
00443 ftype ds = (1+sigmaa*sigmaa)/pow(1-sigmaa*sigmaa,2);
00444 ftype ds2 = 2*sigmaa*(3+sigmaa*sigmaa)/pow(1-sigmaa*sigmaa,3);
00445 result.dr2 = result.dr*ds2 + result.dr2*ds*ds;
00446 result.dr = result.dr*ds;
00447 }
00448 return result;
00449 }
00450
00451
00452 }
00453
00454 #endif