1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
|
#include <clipper/clipper.h>
#include <clipper/clipper-ccp4.h>
//#include <clipper/clipper-ccp4.h>
#include <iostream>
using namespace clipper;
using namespace clipper::data32;
int main()
{
CCP4MTZfile file;
// import an mtz
CSpacegroup myspgr( "" );
CCell mycell( myspgr, "" );
CResolution myreso( mycell, "" );
CHKL_info mydata( myreso, "GerE native and MAD.." );
CHKL_data<F_sigF> myfsig( mydata );
CHKL_data<F_sigF_ano> myanom( mydata );
file.open_read("testfile.mtz");
myspgr.init( file.spacegroup() );
mycell.init( file.cell() );
myreso.init( Resolution( 10.0 ) );
file.import_hkl_info( mydata, false );
//myfsig.HKL_data<F_sigF>::debug();
//myanom.HKL_data<F_sigF_ano>::debug();
new CMTZcrystal(mydata, "xtal", MTZcrystal("xtal", "proj", mycell ));
new CMTZdataset(mydata, "xtal/dset", MTZdataset("dset", 1.76) );
file.import_chkl_data( myfsig, "*/native/[FP SIGFP]" );
file.import_chkl_data( myanom, "*/*/[F(+)SEinfl SIGF(+)SEinfl F(-)SEinfl SIGF(-)SEinfl]", "xtal/dset/myanom");
std::vector<String> v1 = file.column_paths(); for ( int i = 0; i < v1.size(); i++ ) std::cout << "Column: " << i << " " << v1[i] << "\n";
std::vector<String> v2 = file.assigned_paths(); for ( int i = 0; i < v2.size(); i++ ) std::cout << "Import: " << i << " " << v2[i] << "\n";
file.close_read();
new Container(mydata, "newproj");
new CMTZcrystal(mydata, "newproj/newcryst", MTZcrystal("newcryst", "newproj", Cell(Cell_descr(10.0,20.0,30.0))));
new CMTZdataset(mydata, "newproj/newcryst/newdset", MTZdataset("newdset", 1.76) );
mydata.Container::debug();
std::cout << mycell.format() << "\n";
for (int i=0; i<mydata.num_reflections(); i+=100) {
std::cout << i << " " << myfsig[i].f() << " " << myfsig[i].sigf() << " " << myanom[i].f_pl() << " " << myanom[i].f_mi() << " " << myanom[i].cov() << "\n";
}
for (HKL_info::HKL_reference_index i = myfsig.first_data(); !i.last(); myfsig.next_data( i )) {
std::cout << i.index() << " " << myfsig[i].f() << " " << myfsig[i].sigf() << "\n";
}
file.open_append("testfile.mtz","junk.mtz");
file.export_chkl_data( myanom, "Gere:native/native/anom");
file.close_append();
file.open_write("junk2.mtz");
file.export_hkl_info( mydata );
file.export_chkl_data( myanom, "Gere:native/native/anom");
file.close_write();
// now test NaN
std::cout << "\n----------------------------------------\n";
std::cout << "NaN: " << sizeof(clipper::ftype32) << "=" << sizeof(clipper::itype32) << "\n";
std::cout << "NaN: " << sizeof(clipper::ftype64) << "=" << sizeof(clipper::itype64) << "\n";
std::cout << "Nan: ";
{float x; Util::set_null(x);std::cout << x << Util::isnan(x) << Util::is_nan(x) << Util::is_null(x);}
{double x; Util::set_null(x);std::cout << x << Util::isnan(x) << Util::is_nan(x) << Util::is_null(x);}
{float x = Util::nanf();std::cout << x << Util::isnan(x) << Util::is_nan(x) << Util::is_null(x);}
{double x = Util::nand();std::cout << x << Util::isnan(x) << Util::is_nan(x) << Util::is_null(x);}
{float x = Util::nand();std::cout << x << Util::isnan(x) << Util::is_nan(x) << Util::is_null(x);}
{double x = Util::nanf();std::cout << x << Util::isnan(x) << Util::is_nan(x) << Util::is_null(x);}
std::cout << "\n";
// now test the resolution function evaluator
std::cout << "\n----------------------------------------\n";
CHKL_info rfl( mycell, "GerE native and MAD..(2)" );
CHKL_data<F_sigF> fsigdata( rfl, "" );
// read data to higher resolution
file.open_read("testfile.mtz");
rfl.init( file.spacegroup(), file.cell(), Resolution(3.0) );
file.import_hkl_info( rfl, false );
file.import_chkl_data( fsigdata, "*/native/[FP SIGFP]" );
file.close_read();
Range<ftype> slim = fsigdata.invresolsq_range();
std::cout << slim.min() << " " << slim.max() << "\n";
TargetFn_meanFnth<F_sigF> targetfn( fsigdata, 2.0 );
BasisFn_binner basisfn( fsigdata, 10 );
std::vector<ftype> p1( 10, 1.0 );
ResolutionFn rfn1( rfl, basisfn, targetfn, p1 );
BasisFn_gaussian basisfn2;
std::vector<ftype> q( 2, 1.0 );
ResolutionFn_nonlinear rfn2( rfl, basisfn2, targetfn, q );
BasisFn_expcubic basisfn3;
std::vector<ftype> q1( 4, 1.0 );
ResolutionFn_nonlinear rfn3( rfl, basisfn3, targetfn, q1 );
BasisFn_spline basisfn4( fsigdata, 10 );
std::vector<ftype> p2( 10, 1.0 );
ResolutionFn rfn4( rfl, basisfn4, targetfn, p2 );
CHKL_data<E_sigE> esigdata( fsigdata );
esigdata.compute( fsigdata, Compute_EsigE_from_FsigF() );
TargetFn_meanEnth<E_sigE> targetfn_e( esigdata, 2.0 );
BasisFn_spline basisfn5( fsigdata, 10 );
ResolutionFn rfn5( rfl, basisfn5, targetfn_e, p2 );
for ( HKL_info::HKL_reference_index ih = rfl.first(); !ih.last(); ih.next() ) {
ftype eps = ih.hkl_class().epsilon();
if ( fabs( rfn4.f(ih) - rfn5.f(ih) ) > 0.001 ) std::cout << "err\n";
if (ih.index()%100 == 0) std::cout << ih.invresolsq() << " " << ih.hkl().format() << " " << rfn1.f(ih)*eps << " " << rfn2.f(ih)*eps << " " << rfn3.f(ih)*eps << " " << rfn4.f(ih)*eps << "\n";
}
}
|
