Swarm-NG  1.1
swarm.cpp
Go to the documentation of this file.
1 /*************************************************************************
2  * Copyright (C) 2011 by Saleh Dindar and the Swarm-NG Development Team *
3  * *
4  * This program is free software; you can redistribute it and/or modify *
5  * it under the terms of the GNU General Public License as published by *
6  * the Free Software Foundation; either version 3 of the License. *
7  * *
8  * This program is distributed in the hope that it will be useful, *
9  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
11  * GNU General Public License for more details. *
12  * *
13  * You should have received a copy of the GNU General Public License *
14  * along with this program; if not, write to the *
15  * Free Software Foundation, Inc., *
16  * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
17  ************************************************************************/
18 
32 #include <iostream>
33 #include <boost/program_options.hpp>
34 #include <boost/program_options/positional_options.hpp>
35 #include "swarm.h"
36 #include "query.hpp"
37 #include "snapshot.hpp"
38 #include "stopwatch.h"
39 #include "gpu/device_settings.hpp"
40 
41 int DEBUG_LEVEL = 0;
42 
43 using namespace swarm;
44 using namespace std;
45 namespace po = boost::program_options;
46 using boost::bind;
47 using boost::shared_ptr;
48 
49 
50 
51 // Runtime variables
52 string command;
53 config cfg;
54 config base_cfg;
55 defaultEnsemble initial_ens;
56 defaultEnsemble current_ens;
57 defaultEnsemble reference_ens;
58 po::variables_map argvars_map;
59 Pintegrator integ;
60 
61 void stability_test() {
62 
63  defaultEnsemble &ens = current_ens;
64 
65  double begin_time = ens.time_ranges().median;
66  double destination_time = cfg.optional("destination_time", begin_time + 10 * M_PI );
67  double interval = cfg.optional("interval", (destination_time-begin_time) ) ;
68  double logarithmic = cfg.optional("logarithmic", 0 ) ;
69  double allowed_deltaE = cfg.optional("allowed_deltaE", 0.01 );
70 
71  if(destination_time < begin_time ) ERROR("Destination time should be larger than begin time");
72  if(interval < 0) ERROR("Interval cannot be negative");
73  if(interval < 0.001) ERROR("Interval too small");
74  if(logarithmic != 0 && logarithmic <= 1) ERROR("logarithm base should be greater than 1");
75 
76 
77  std::cout << "Time, Energy Conservation Factor (delta E/E), Active Systems" << std::endl;
78  for(double time = begin_time; time < destination_time; ) {
79 
80  if(logarithmic > 1)
81  time = (time > 0) ? time * logarithmic : interval;
82  else
83  time += interval;
84 
85  double effective_time = min(time,destination_time);
86  integ->set_destination_time ( effective_time );
87 
88  DEBUG_OUTPUT(2, "Integrator ensemble" );
89  integ->integrate();
90 
91  SYNC;
92  DEBUG_OUTPUT(2, "Check energy conservation" );
93  ensemble::range_t deltaE_range = energy_conservation_error_range(ens, initial_ens );
94 
95  int active_systems = ens.nsys() - number_of_disabled_systems( ens );
96  std::cout << effective_time << ", " << deltaE_range.max << ", " << deltaE_range.median << ", " << active_systems << std::endl;
97 
98  if(deltaE_range.median > allowed_deltaE){
99  INFO_OUTPUT(0, "At least half of systems are too unstable to conserve energy. dE/E: " << deltaE_range.median << std::endl );
100  exit(1);
101  }
102 
103  }
104 
105 }
106 
107 bool read_input_file(defaultEnsemble& ens, config& cfg){
108 
109  if(cfg.valid("input") ) {
110  INFO_OUTPUT(1, "Loading from binary file " << cfg["input"]);
111  ens = swarm::snapshot::load(cfg["input"]);
112  INFO_OUTPUT(1,", time = " << ens.time_ranges() << endl);
113  return true;
114 
115  }else if(cfg.valid("text_input")) {
116  INFO_OUTPUT(1, "Loading from text file " << cfg["text_input"]);
117  ens = swarm::snapshot::load_text(cfg["text_input"]);
118  INFO_OUTPUT(1,", time = " << ens.time_ranges() << endl);
119  return true;
120 
121  }else
122  return false;
123 
124 }
125 
126 bool read_output_file(defaultEnsemble& ens, config& cfg){
127 
128  if(cfg.valid("output") ) {
129  INFO_OUTPUT(1, "Loading from binary file " << cfg["output"]);
130  ens = swarm::snapshot::load(cfg["output"]);
131  INFO_OUTPUT(1,", time = " << ens.time_ranges() << endl);
132  return true;
133 
134  }else if(cfg.valid("text_output")) {
135  INFO_OUTPUT(1, "Loading from text file " << cfg["text_output"]);
136  ens = swarm::snapshot::load_text(cfg["text_output"]);
137  INFO_OUTPUT(1,", time = " << ens.time_ranges() << endl);
138  return true;
139 
140  }else
141  return false;
142 
143 }
144 
145 void load_generate_ensemble(){
146  // Load/Generate the ensemble
147  if(read_input_file(initial_ens,cfg)) {
148 
149  }else{
150  INFO_OUTPUT(1, "Generating new ensemble: " << cfg["nsys"] << ", " << cfg["nbod"] << endl);
151  srand(time(NULL));
152  initial_ens = generate_ensemble(cfg);
153  }
154 }
155 
156 bool save_ensemble(){
157  // Save the ensemble
158  if(cfg.valid("output")) {
159  INFO_OUTPUT(1, "Saving as binary " << cfg["output"]);
160  INFO_OUTPUT(1, ", time = " << current_ens.time_ranges() << endl);
161  swarm::snapshot::save(current_ens,cfg["output"]);
162  return true;
163 
164  }else if(cfg.valid("text_output")) {
165  INFO_OUTPUT(1, "Saving as text " << cfg["text_output"]);
166  INFO_OUTPUT(1, ", time = " << current_ens.time_ranges() << endl);
167  swarm::snapshot::save_text(current_ens,cfg["text_output"]);
168  return true;
169 
170  }else
171  return false;
172 }
173 
174 void init_cuda(){
175  // Initialize Swarm
176  swarm::init(cfg);
177 }
178 
179 void prepare_integrator () {
180  // Initialize Integrator
181  DEBUG_OUTPUT(2, "Initializing integrator" );
182  double begin_time = initial_ens.time_ranges().average;
183  double destination_time = cfg.optional("destination_time", begin_time + 10 * M_PI );
184  int max_iterations = cfg.optional("max_iterations", integrator::_default_max_iterations );
185  integ = integrator::create(cfg);
186  integ->set_max_iterations(max_iterations);
187  integ->set_ensemble(current_ens);
188  integ->set_destination_time ( destination_time );
189  SYNC;
190 }
191 
192 void generic_integrate () {
193  DEBUG_OUTPUT(2, "Integrator ensemble" );
194  integ->integrate();
195 }
196 
197 void run_integration(){
198  if(!validate_configuration(cfg) ) ERROR( "Invalid configuration" );
199 
200  load_generate_ensemble();
201 
202  DEBUG_OUTPUT(2, "Make a copy of ensemble for energy conservation test" );
203  current_ens = initial_ens.clone();
204 
205  prepare_integrator();
206 
207  double integration_time = watch_time ( cfg.valid("interval") ? stability_test : generic_integrate );
208 
209  save_ensemble();
210 
211  INFO_OUTPUT( 1, "Integration time: " << integration_time << " ms " << std::endl);
212 }
213 
214 
215 void reference_integration() {
216  DEBUG_OUTPUT(2, "Make a copy of ensemble for reference ensemble" );
217  reference_ens = initial_ens.clone();
218  DEBUG_OUTPUT(1, "Reference integration" );
219  prepare_integrator();
220  integ->set_ensemble( reference_ens );
221  integ->integrate();
222 }
223 
224 
225 
226 bool verification_results = true, verify_mode = false;
227 double pos_threshold = 0, vel_threshold = 0, time_threshold = 0;
228 
229 void fail_verify() {
230  if(verify_mode){
231  INFO_OUTPUT(0, "Verify failed" << endl);
232  exit(1);
233  }else{
234  verification_results = false;
235  }
236 }
237 
238 void output_test() {
239  if(!validate_configuration(cfg) ) ERROR( "Invalid configuration" );
240 
241  if(!read_input_file(initial_ens, cfg) ) {
242  ERROR("you should have a tested input file");
243  }
244 
245  DEBUG_OUTPUT(2, "Make a copy of ensemble for energy conservation test" );
246  current_ens = initial_ens.clone();
247 
248  prepare_integrator();
249 
250  double integration_time = watch_time ( cfg.valid("interval") ? stability_test : generic_integrate );
251 
252  if(read_output_file(reference_ens,cfg)){
253 
254  // Compare with reneference ensemble for integrator verification
255  double pos_diff = 0, vel_diff = 0, time_diff = 0;
256  bool comparison = compare_ensembles( current_ens, reference_ens , pos_diff, vel_diff, time_diff );
257 
258  INFO_OUTPUT(1,"\tPosition difference: " << pos_diff << endl
259  << "\tVelocity difference: " << vel_diff << endl
260  << "\tTime difference: " << time_diff << endl );
261 
262  if( !comparison || pos_diff > pos_threshold || vel_diff > vel_threshold || time_diff > time_threshold ){
263  INFO_OUTPUT(0, "Test failed" << endl);
264  exit(1);
265  }else {
266  INFO_OUTPUT(0, "Test success" << endl);
267  }
268 
269 
270  }else{
271  ERROR("You should provide a test output file");
272  }
273 
274  INFO_OUTPUT( 1, "Integration time: " << integration_time << " ms " << std::endl);
275 }
276 
277 void benchmark_item(const string& param, const string& value) {
278  //outputConfigSummary(std::cout,cfg);
279  if(!validate_configuration(cfg) ) ERROR( "Invalid configuration" );
280 
281  if(param == "input" || param == "nsys" || param == "nbod" || cfg.count("reinitialize"))
282  load_generate_ensemble();
283 
284  DEBUG_OUTPUT(2, "Make a copy of ensemble for energy conservation test" );
285  current_ens = initial_ens.clone();
286 
287  double init_time = watch_time( prepare_integrator );
288 
289  double integration_time = watch_time( generic_integrate );
290 
291  double max_deltaE = find_max_energy_conservation_error(current_ens, initial_ens );
292 
293  // Compare with reneference ensemble for integrator verification
294  double pos_diff = 0, vel_diff = 0, time_diff = 0;
295  bool comparison = compare_ensembles( current_ens, reference_ens , pos_diff, vel_diff, time_diff );
296 
298  std::cout << param << ", "
299  << value << ", "
300  << current_ens.time_ranges() << ", "
301  << max_deltaE << ", "
302  << pos_diff << ", "
303  << vel_diff << ", "
304  << time_diff << ", "
305  << integration_time << ", "
306  << init_time
307  << std::endl;
308 
309  if( !comparison || pos_diff > pos_threshold || vel_diff > vel_threshold || time_diff > time_threshold ){
310  fail_verify();
311  }
312 
313 }
314 
319 struct parameter_range {
320  string parameter;
321  vector<string> values;
322  bool interpolate;
323 };
324 
327 #if BOOST_VERSION < 104200
328  void raise_validation_error(const string& s){
329  throw po::validation_error(s);
330  }
331 #else
332  void raise_validation_error(const string& s){
333  throw po::validation_error(po::validation_error::invalid_option_value, s);
334  }
335 #endif
336 
338 void validate(boost::any& v, const std::vector<std::string>& values
339  , parameter_range* , int){
340 
341  // Make sure no previous assignment to 'v' was made.
342  po::validators::check_first_occurrence(v);
343 
344  // Extract the first string from 'values'. If there is more than
345  // one string, it's an error, and exception will be thrown.
346  const std::string& s = po::validators::get_single_string(values);
347 
348  parameter_range p;
349  static boost::regex assign("^(\\w+)=(.+)$");
350  static boost::regex range("^([0-9\\.Ee]+)\\.\\.([0-9\\.Ee]+)$");
351  static boost::regex rangeinc("^([0-9\\.Ee]+)\\.\\.([0-9\\.Ee]+)\\.\\.([0-9\\.Ee]+)$");
352  static boost::regex list("^([a-zA-Z0-9_\\.]+)(?:,([a-zA-Z0-9_\\.]+))+$");
353  static boost::regex item("([a-zA-Z0-9_\\.]+)");
354  boost::smatch match, items;
355  if (boost::regex_match(s, match, assign)) {
356  p.parameter = match[1];
357  string value_range = match[2];
358  boost::smatch range_match, list_match;
359  if (boost::regex_match( value_range , range_match, rangeinc )) {
360  p.values.push_back(range_match[1]);
361  p.values.push_back(range_match[3]);
362  p.values.push_back(range_match[2]);
363  p.interpolate = true;
364  }else if (boost::regex_match( value_range , range_match, range )) {
365  p.values.push_back(range_match[1]);
366  p.values.push_back(range_match[2]);
367  p.interpolate = true;
368  }else if (boost::regex_match( value_range , list_match, list )) {
369  boost::sregex_iterator items( value_range.begin(), value_range.end() , item ), end;
370  for(;items != end; items++)
371  p.values.push_back((*items)[0]);
372  p.interpolate = false;
373  }else {
374  raise_validation_error("Range is invalid");
375  }
376  }else
377  raise_validation_error("Wrong parameter-value pair ");
378  v = boost::any(p);
379 }
380 
381 void benchmark(const parameter_range& pr){
382 
383 
384  // Go through the loop
385  std::cout << "Parameter, Value, Time, Energy Conservation Factor (delta E/E)"
386  ", Position Difference, Velocity Difference, Time difference "
387  ", Integration (ms), Integrator initialize (ms) \n";
388 
389  string param = pr.parameter;
390  vector<string> values = pr.values;
391 
392  if(!(param == "input" || param == "nsys" || param == "nbod" || cfg.count("reinitialize")))
393  load_generate_ensemble();
394  outputConfigSummary(std::cout,cfg);
395  if(verify_mode)
396  reference_integration();
397 
398  if(pr.interpolate) {
399  // Numeric range
400  double from = atof(values[0].c_str()), to = atof(values[1].c_str());
401  double increment = (values.size() == 3) ? atof(values[2].c_str()) : 1;
402 
403  for(double i = from; i <= to ; i+= increment ){
404  std::ostringstream stream;
405  stream << i;
406  string value = stream.str();
407 
408  cfg[param] = value;
409 
410  DEBUG_OUTPUT(1, "=========================================");
411  benchmark_item(param,value);
412  }
413  }else {
414  // Iterate over values
415  for(vector<string>::const_iterator i = values.begin(); i != values.end(); i++){
416  string value = *i;
417 
418  if(param == "config")
419  cfg = config::load(value,base_cfg);
420  else
421  cfg[param] = *i;
422 
423  DEBUG_OUTPUT(1, "=========================================");
424  benchmark_item(param,value);
425  }
426  }
427 }
428 
429 void print_version(){
430  bool amd64 = sizeof(void*) == 8;
431  cout << "Swarm is running as " << (amd64 ? "64bit" : "32bit" ) << endl;
432  exit(0);
433 }
434 
435 void parse_commandline_and_config(int argc, char* argv[]){
436 
437  po::positional_options_description pos;
438  pos.add("command", 1);
439  pos.add("param_value_pair", -1);
440 
441  po::options_description desc("Usage:\n \tswarm [options] COMMAND [PARAMETER VALUE VALUE ...]\n\n"
442  "Possible commands are: \n"
443  "\tintegrate : Integrate a [loaded|generated] ensemble\n"
444  "\tbenchmark : Compare outputs for different methods of integrations\n"
445  "\tverify : Verify an integrator against a reference integrator\n"
446  "\tquery : Query data from a log file\n"
447  "\ttest : Test a configuration against input/output files\n"
448  "\tgenerate : Generate a new ensemble and save it to output file\n"
449  "\tconvert : Read input file and write it to output file (converts to/from text)\n"
450  "\nOptions"
451  );
452 
453 
454  po::options_description integrate("Integation Options");
455  integrate.add_options()
456  ("destination_time,d", po::value<std::string>() , "Destination time to achieve")
457  ("logarithmic,l", po::value<std::string>() , "Produce times in logarithmic scale" )
458  ("interval,n", po::value<std::string>() , "Energy test intervals")
459  ("input,i", po::value<std::string>(), "Input file")
460  ("output,o", po::value<std::string>(), "Output file")
461  ("text_input,I", po::value<std::string>(), "Text Input file")
462  ("text_output,O", po::value<std::string>(), "Text Output file");
463 
464  po::options_description benchmark("Benchmark Options");
465  benchmark.add_options()
466  ("range", po::value<parameter_range>(), "Parameter value range param=v1,v2,v3");
467 
468  po::options_description query("Query Options");
469  query.add_options()
470  ("time,t", po::value<query::time_range_t>(), "range of times to query")
471  ("system,s", po::value<query::sys_range_t>(), "range of systems to query")
472  ("body,b", po::value<query::sys_range_t>(), "range of bodies to query")
473  ("keplerian,k", "output in Keplerian coordinates")
474  ("astrocentric", "output coordinates in astrocentric frame")
475  ("barycentric", "output coordinates in barycentric frame")
476  ("origin", "output coordinates in origin frame [default w/ Cartesian]")
477  ("jacobi", "output coordinates in Jacobi frame [default w/ Keplerian]")
478  ("logfile,f", po::value<std::string>(), "the log file to query");
479 
480  po::options_description positional("Positional Options");
481  positional.add_options()
482  ("command" , po::value<string>() , "Swarm command: integrate, benchmark, verify, query ")
483  ("param_value_pair" , po::value< vector<string> >() , "Parameter value pairs: param=value")
484  ;
485 
486  po::options_description general("General Options");
487  general.add_options()
488  ("cfg,c", po::value<std::string>(), "Integrator configuration file")
489  ("nogpu,g", "Only use CPU for integration, use this if your machine does not have an NVIDIA GPU")
490  ("defaults", "Use default values for required configuration options")
491  ("help,h", "produce help message")
492  ("plugins,p", "list all of the plugins")
493  ("verbose,v", po::value<int>(), "Verbosity level (debug output) ")
494  ("quiet,q", "Suppress all the notifications (for CSV generation)")
495  ("version,V", "Print version message");
496 
497  desc.add(general).add(integrate).add(benchmark).add(query);
498 
499  po::options_description all;
500  all.add(desc).add(positional);
501 
502 
503 
504  po::variables_map &vm = argvars_map;
505  po::store(po::command_line_parser(argc, argv).
506  options(all).positional(pos).run(), vm);
507  po::notify(vm);
508 
510  //
511  if (vm.count("help")) { std::cout << desc << "\n"; exit(1); }
512  if (vm.count("version")) print_version();
513  if (vm.count("verbose") ) DEBUG_LEVEL = vm["verbose"].as<int>();
514  if (vm.count("quiet") ) DEBUG_LEVEL = -1;
515 
516  if (vm.count("plugins")) {
517  cout << swarm::plugin::help;
518  exit(0);
519  }
520 
521  if(vm.count("defaults")){
522  cfg = default_config();
523  }
524  if(vm.count("cfg")){
525  std::string icfgfn = vm["cfg"].as<std::string>();
526  cfg = config::load(icfgfn, cfg);
527  }
528 
529  if(vm.count("command") == 0){
530  std::cout << desc << "\n"; exit(1);
531  }
532 
533  const int cmd_to_config_len = 7;
534  const char* cmd_to_config[cmd_to_config_len] = { "input", "output", "text_input" , "text_output", "destination_time", "logarithmic", "interval" };
535  for(int i = 0; i < cmd_to_config_len; i++)
536  if(vm.count(cmd_to_config[i]))
537  cfg[cmd_to_config[i]] = vm[cmd_to_config[i]].as<std::string>();
538 
539  if(vm.count("nogpu"))
540  cfg["nogpu"] = "1";
541 
542  if(vm.count("param_value_pair") > 0) {
543  vector<string> pairs = vm["param_value_pair"].as< vector<string> >();
544 
545  for(vector<string>::const_iterator i = pairs.begin(); i != pairs.end(); i++){
546  string::size_type pos = i->find_first_of( "=" );
547  if( string::npos != pos ){
548  string parameter = i->substr( 0, pos );
549  string value = i->substr( pos + 1 );
550  cfg[parameter] = value;
551  INFO_OUTPUT(3, parameter << " set to " << value << endl );
552  }else{
553  ERROR( ("Wrong parameter value pair : " + *i ).c_str() );
554  }
555  }
556  }
557 
558 
559 }
560 
561 
562 int main(int argc, char* argv[]){
563 
564  // Command line and Config file
565  parse_commandline_and_config(argc,argv);
566 // outputConfigSummary(std::cout,cfg);
567  base_cfg = cfg;
568  command = argvars_map["command"].as< string >();
569 
570 
571  // Set some variables
572  pos_threshold = cfg.optional("pos_threshold", 1e-10);
573  vel_threshold = cfg.optional("vel_threshold", 1e-10);
574  time_threshold = cfg.optional("time_threshold", 1e-10);
575 
576  // Branch based on COMMAND
577  if(command == "integrate"){
578  init_cuda();
579  run_integration();
580 
581  }else if(command == "test"){
582  init_cuda();
583  output_test();
584  }else if(command == "benchmark" || command == "verify") {
585  init_cuda();
586  if(command == "verify")
587  verify_mode = true;
588 
589  if(argvars_map.count("range") > 0) {
590  parameter_range pr = argvars_map["range"].as<parameter_range>();
591  benchmark(pr);
592  }else
593  cerr << "A parameter-value range is missing" << endl;
594 
595  if(command == "verify") {
596  cout << (verification_results ? "Verify success" : "Verify failed") << endl;
597  return verification_results ? 0 : 1;
598  }
599  }
600 
601 
602  else if(command == "query" ) {
603  using namespace query;
604  if (!argvars_map.count("logfile")) { cerr << "Name of input log file is missing \n"; return 1; }
605 
606  time_range_t T;
607  sys_range_t sys, body_range;
608  if (argvars_map.count("time")) { T = argvars_map["time"].as<time_range_t>(); }
609  if (argvars_map.count("system")) { sys = argvars_map["system"].as<sys_range_t>(); }
610  if (argvars_map.count("body")) { body_range = argvars_map["body"].as<sys_range_t>(); }
611  if (argvars_map.count("keplerian")) { query::set_keplerian_output(); }
612  if (argvars_map.count("origin")) { query::set_coordinate_system(query::origin); }
613  if (argvars_map.count("astrocentric")) { query::set_coordinate_system(query::astrocentric); }
614  if (argvars_map.count("barycentric")) { query::set_coordinate_system(query::barycentric); }
615  if (argvars_map.count("jacobi")) { query::set_coordinate_system(query::jacobi); }
616 
617 
618  std::cout.flush();
619  cerr << "# Systems: " << sys << " Bodies: " << body_range << " Times: " << T << endl;
620  if (argvars_map.count("keplerian"))
621  { std::cerr << "# Output in Keplerian coordinates "; }
622  else { std::cerr << "# Output in Cartesian coordinates "; }
623 #if 0
624  switch(query::planets_coordinate_system)
625  {
626  case astrocentric:
627  std::cerr << "(astrocentric)";
628  break;
629  case barycentric:
630  std::cerr << "(barycentric)";
631  break;
632  case jacobi:
633  std::cerr << "(jacobi)";
634  break;
635  case origin:
636  std::cerr << "(origin)";
637  break;
638  }
639 #endif
640  std::cerr << "\n";
641 
642 
643  std::string datafile(argvars_map["logfile"].as<std::string>());
644 
645  query::execute(datafile, T, sys, body_range);
646  }
647 
648  else if(command == "generate" ) {
649  srand(time(NULL));
650  INFO_OUTPUT(1, "Generating new ensemble: " << cfg["nsys"] << ", " << cfg["nbod"] << endl);
651  current_ens = generate_ensemble(cfg);
652  save_ensemble();
653  }
654 
655  else if(command == "convert" ) {
656  if( read_input_file(current_ens,cfg) && save_ensemble() ) {
657  INFO_OUTPUT(1,"Converted!");
658  }else{
659  INFO_OUTPUT(1,"Convertion failed");
660  exit(1);
661  }
662  }
663 
664  else
665  std::cerr << "Valid commands are: integrate, benchmark, verify, test, query, generate " << std::endl;
666 
667  return 0;
668 }