IMP  2.0.1
The Integrative Modeling Platform
compute_chi.cpp
1 /**
2  This is the program for validation of SAXS profiles.
3  */
4 #include <IMP/saxs/Profile.h>
6 
7 #include <fstream>
8 #include <vector>
9 #include <string>
10 
11 #include <boost/lexical_cast.hpp>
12 #include <boost/program_options.hpp>
13 namespace po = boost::program_options;
14 
15 int main(int argc, char **argv)
16 {
17  // output arguments
18  for (int i = 0; i < argc; i++) std::cerr << argv[i] << " ";
19  std::cerr << std::endl;
20 
21  bool use_offset = false;
22 
23  po::options_description
24  desc("Usage: <experimental_profile> <profile_file1> <profile_file2> ...");
25  desc.add_options()
26  ("help", "Any number of input profiles is supported. \
27 The chi value is computed relative to the first profile using its error column")
28  ("input-files", po::value< std::vector<std::string> >(),
29  "input PDB and profile files")
30  ("offset,o", "use offset in fitting (default = false)")
31  ;
32  po::positional_options_description p;
33  p.add("input-files", -1);
34  po::variables_map vm;
35  po::store(
36  po::command_line_parser(argc,argv).options(desc).positional(p).run(), vm);
37  po::notify(vm);
38 
39  std::vector<std::string> files, dat_files;
40  if(vm.count("input-files")) {
41  files = vm["input-files"].as< std::vector<std::string> >();
42  }
43  if(vm.count("help") || files.size() == 0) {
44  std::cout << desc << "\n";
45  return 0;
46  }
47  if(vm.count("offset")) use_offset=true;
48 
49  std::vector<IMP::saxs::Profile *> exp_profiles;
50  for(unsigned int i=0; i<files.size(); i++) {
51  // check if file exists
52  std::ifstream in_file(files[i].c_str());
53  if(!in_file) {
54  std::cerr << "Can't open file " << files[i] << std::endl;
55  exit(1);
56  }
57 
58  IMP::saxs::Profile *profile = new IMP::saxs::Profile(files[i]);
59  if(profile->size() == 0) {
60  std::cerr << "can't parse input file " << files[i] << std::endl;
61  return 1;
62  } else {
63  dat_files.push_back(files[i]);
64  exp_profiles.push_back(profile);
65  std::cout << "Profile read from file " << files[i] << " size = "
66  << profile->size() << std::endl;
67  }
68  }
69 
70  IMP::saxs::Profile* exp_saxs_profile = exp_profiles[0];
71  IMP::Pointer<IMP::saxs::ProfileFitter<IMP::saxs::ChiScore> >saxs_score =
72  new IMP::saxs::ProfileFitter<IMP::saxs::ChiScore>(*exp_saxs_profile);
73  for(unsigned int i=1; i<exp_profiles.size(); i++) {
74  std::string fit_file = "fit" +
75  std::string(boost::lexical_cast<std::string>(i)) + ".dat";
76  float chi = saxs_score->compute_score(*exp_profiles[i], use_offset, "fit");
77  std::cout << "File " << files[i] << " chi=" << chi << std::endl;
78  }
79  return 0;
80 }