Swarm-NG  1.1
montecarlo_mcmc_outputs.cpp
Go to the documentation of this file.
1 /*************************************************************************
2  * Copyright (C) 2009-2010 by Eric Ford & 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 <fstream>
34 #include <math.h>
35 #include <signal.h>
36 
37 #include "swarm/swarm.h"
38 #include "swarm/snapshot.hpp"
39 #include "swarm/log/log.hpp"
40 #include "random.hpp"
41 #include "kepler.hpp"
42 
43 #include "swarm/log/host_array_writer.hpp"
44 
45 #define SYNC cudaThreadSynchronize()
46 
47 
48 using namespace swarm;
49 using namespace std;
50 
51 config cfg;
52 
57 defaultEnsemble generate_ensemble_with_initial_conditions_keplerian_from_file(const config& cfg)
58 {
59  defaultEnsemble ens = defaultEnsemble::create( cfg.require("nbod",0), cfg.require("nsys",0) );
60  std::cerr << "# nsystems = " << ens.nsys() << " nbodies/system = " << ens.nbod() << "\n";
61 
62  // std::cerr << "# Set initial time for all systems = ";
63  double time_init = cfg.optional("time_init", 0.0);
64  // std::cerr << time_init << ".\n";
65 
66  const bool use_jacobi = cfg.optional("use_jacobi", 0);
67 
68  // std::cerr << "# Writing stable system ids\n";
69  ifstream jacobi_input( cfg.optional("input_mcmc_keplerian", string("mcmc.out")).c_str() );
70 
71  for(unsigned int sysid=0;sysid<ens.nsys();++sysid)
72  {
73  assert(jacobi_input.good());
74  string id;
75  double mass_star;
76  jacobi_input >> id >> mass_star;
77  ens[sysid].id() = sysid;
78  ens[sysid].time() = time_init;
79  ens[sysid].set_active();
80  double x=0, y=0, z=0, vx=0, vy=0, vz=0;
81  ens.set_body(sysid, 0, mass_star, x, y, z, vx, vy, vz);
82  ens[sysid][0].attribute(0) = 0.005; // solar radius
83  double mass_enclosed = mass_star;
84  for(unsigned int bod=1;bod<ens.nbod();++bod)
85  {
86 
87  double mass_planet, a, e, i, O, w, M;
88  jacobi_input >> mass_planet >> a >> e >> i >> w >> O >> M;
89  // a = pow((period/365.25)*(period/365.25)*mass_enclosed,1.0/3.0);
90  i *= M_PI/180.;
91  O *= M_PI/180.;
92  w *= M_PI/180.;
93  M *= M_PI/180.;
94 
95  mass_enclosed += mass_planet;
96  double mass = use_jacobi ? mass_enclosed : mass_star+mass_planet;
97  if(cfg.count("verbose")&&(sysid<10))
98  std::cout << "# Drawing sysid= " << sysid << " bod= " << bod << ' ' << mass_planet << " " << a << ' ' << e << ' ' << i*180./M_PI << ' ' << O*180./M_PI << ' ' << w*180./M_PI << ' ' << M*180./M_PI << '\n';
99 
100  calc_cartesian_for_ellipse(x,y,z,vx,vy,vz, a, e, i, O, w, M, mass);
101  // printf("%d %d: %lg (%lg %lg %lg) (%lg %lg %lg) \n", sysid, bod, mass_planet, x,y,z,vx,vy,vz);
102 
103  if(use_jacobi)
104  {
105  double bx, by, bz, bvx, bvy, bvz;
106  ens.get_barycenter(sysid,bx,by,bz,bvx,bvy,bvz,bod-1);
107  x += bx; y += by; z += bz;
108  vx += bvx; vy += bvy; vz += bvz;
109  }
110 
111  // assign body a mass, position and velocity
112  ens.set_body(sysid, bod, mass_planet, x, y, z, vx, vy, vz);
113  ens[sysid][bod].attribute(0) = 0.;
114 
115  if(cfg.count("verbose")&&(sysid<10))
116  {
117  double x_t = ens.x(sysid,bod);
118  double y_t = ens.y(sysid,bod);
119  double z_t = ens.z(sysid,bod);
120  double vx_t = ens.vx(sysid,bod);
121  double vy_t = ens.vy(sysid,bod);
122  double vz_t = ens.vz(sysid,bod);
123 
124  std::cout << " x= " << x << "=" << x_t << " ";
125  std::cout << " y= " << y << "=" << y_t << " ";
126  std::cout << " z= " << z << "=" << z_t << " ";
127  std::cout << "vx= " << vx << "=" << vx_t << " ";
128  std::cout << "vy= " << vy << "=" << vy_t << " ";
129  std::cout << "vz= " << vz << "=" << vz_t << "\n";
130  }
131 
132  } // end loop over bodies
133 
134  // Shift into barycentric frame
135  ens.get_barycenter(sysid,x,y,z,vx,vy,vz);
136  for(unsigned int bod=0;bod<ens.nbod();++bod)
137  {
138  ens.set_body(sysid, bod, ens.mass(sysid,bod),
139  ens.x(sysid,bod)-x, ens.y(sysid,bod)-y, ens.z(sysid,bod)-z,
140  ens.vx(sysid,bod)-vx, ens.vy(sysid,bod)-vy, ens.vz(sysid,bod)-vz);
141  ens[sysid][bod].attribute(0) = 0.;
142  } // end loop over bodies
143 
144  } // end loop over systems
145  return ens;
146 }
147 
148 
153 defaultEnsemble generate_ensemble_with_initial_conditions_cartesian_from_file(const config& cfg)
154 {
155  defaultEnsemble ens = defaultEnsemble::create( cfg.require("nbod",0), cfg.require("nsys",0) );
156  std::cerr << "# nsystems = " << ens.nsys() << " nbodies/system = " << ens.nbod() << "\n";
157 
158  // std::cerr << "# Set initial time for all systems = ";
159  double time_init = cfg.optional("time_init", 0.0);
160  // std::cerr << time_init << ".\n";
161 
162  const bool use_jacobi = cfg.optional("use_jacobi", 0);
163 
164  // std::cerr << "# Writing stable system ids\n";
165  ifstream mcmc_input( cfg.optional("input_mcmc_cartesian", string("mcmc.out")).c_str() );
166 
167  assert(mcmc_input.good());
168 
169  const int lines_to_skip = cfg.optional("lines_to_skip", 0);
170  for(int i=0;i<lines_to_skip;++i)
171  {
172  assert(mcmc_input.good());
173  char junk[1024];
174  mcmc_input.getline(junk,1024);
175  }
176 
177  for(unsigned int sysid=0;sysid<ens.nsys();++sysid)
178  {
179  assert(mcmc_input.good());
180 #if 1
181  std::vector<double> masses(ens.nbod(),0.0);
182  for(unsigned int bod=0;bod<ens.nbod();++bod)
183  {
184  mcmc_input >> masses[bod];
185  } // end loop over bodies
186 
187  ens[sysid].id() = sysid;
188  ens[sysid].time() = time_init;
189  ens[sysid].set_active();
190  double x=0, y=0, z=0, vx=0, vy=0, vz=0;
191  ens.set_body(sysid, 0, masses[0], x, y, z, vx, vy, vz);
192  ens[sysid][0].attribute(0) = 0.005; // solar radius
193  double mass_enclosed = masses[0];
194  for(unsigned int bod=1;bod<ens.nbod();++bod)
195  {
196  mcmc_input >> x >> y >> z >> vx >> vy >> vz;
197 
198  mass_enclosed += masses[bod];
199  double mass = use_jacobi ? mass_enclosed : masses[0]+masses[bod];
200  if(cfg.count("verbose")&&(sysid<10))
201  std::cout << "# Drawing sysid= " << sysid << " bod= " << bod << ' ' << masses[bod] << " " << x << ' ' << y << ' ' << z << ' ' << vx << ' ' << vy << ' ' << vz << '\n';
202 
203  if(use_jacobi)
204  {
205  double bx, by, bz, bvx, bvy, bvz;
206  ens.get_barycenter(sysid,bx,by,bz,bvx,bvy,bvz,bod-1);
207  x += bx; y += by; z += bz;
208  vx += bvx; vy += bvy; vz += bvz;
209  }
210 
211  // assign body a mass, position and velocity
212  ens.set_body(sysid, bod, masses[bod], x, y, z, vx, vy, vz);
213  ens[sysid][bod].attribute(0) = 0.0;
214 
215  if(cfg.count("verbose")&&(sysid<10))
216  {
217  double x_t = ens.x(sysid,bod);
218  double y_t = ens.y(sysid,bod);
219  double z_t = ens.z(sysid,bod);
220  double vx_t = ens.vx(sysid,bod);
221  double vy_t = ens.vy(sysid,bod);
222  double vz_t = ens.vz(sysid,bod);
223 
224  std::cout << " x= " << x << "=" << x_t << " ";
225  std::cout << " y= " << y << "=" << y_t << " ";
226  std::cout << " z= " << z << "=" << z_t << " ";
227  std::cout << "vx= " << vx << "=" << vx_t << " ";
228  std::cout << "vy= " << vy << "=" << vy_t << " ";
229  std::cout << "vz= " << vz << "=" << vz_t << "\n";
230  }
231  }
232 #else
233  std::vector<double> masses(ens.nbod(),0.0);
234  double x=0, y=0, z=0, vx=0, vy=0, vz=0;
235  double mass_enclosed = 0.0;
236  ens[sysid].id() = sysid;
237  ens[sysid].time() = time_init;
238  ens[sysid].set_active();
239  for(unsigned int bod=0;bod<ens.nbod();++bod)
240  {
241  mcmc_input >> masses[bod];
242  mcmc_input >> x >> y >> z >> vx >> vy >> vz;
243 
244  mass_enclosed += masses[bod];
245  double mass = use_jacobi ? mass_enclosed : masses[0]+masses[bod];
246  if(cfg.count("verbose")&&(sysid<10))
247  std::cout << "# Drawing sysid= " << sysid << " bod= " << bod << ' ' << masses[bod] << " " << x << ' ' << y << ' ' << z << ' ' << vx << ' ' << vy << ' ' << vz << '\n';
248 
249  if(use_jacobi)
250  {
251  double bx, by, bz, bvx, bvy, bvz;
252  ens.get_barycenter(sysid,bx,by,bz,bvx,bvy,bvz,bod-1);
253  x += bx; y += by; z += bz;
254  vx += bvx; vy += bvy; vz += bvz;
255  }
256 
257  // assign body a mass, position and velocity
258  ens.set_body(sysid, bod, masses[bod], x, y, z, vx, vy, vz);
259 
260  if(cfg.count("verbose")&&(sysid<10))
261  {
262  double x_t = ens.x(sysid,bod);
263  double y_t = ens.y(sysid,bod);
264  double z_t = ens.z(sysid,bod);
265  double vx_t = ens.vx(sysid,bod);
266  double vy_t = ens.vy(sysid,bod);
267  double vz_t = ens.vz(sysid,bod);
268 
269  std::cout << " x= " << x << "=" << x_t << " ";
270  std::cout << " y= " << y << "=" << y_t << " ";
271  std::cout << " z= " << z << "=" << z_t << " ";
272  std::cout << "vx= " << vx << "=" << vx_t << " ";
273  std::cout << "vy= " << vy << "=" << vy_t << " ";
274  std::cout << "vz= " << vz << "=" << vz_t << "\n";
275  }
276  } // end loop over bodies
277 #endif
278 
279  // Shift into barycentric frame
280  ens.get_barycenter(sysid,x,y,z,vx,vy,vz);
281  for(unsigned int bod=0;bod<ens.nbod();++bod)
282  {
283  ens.set_body(sysid, bod, ens.mass(sysid,bod),
284  ens.x(sysid,bod)-x, ens.y(sysid,bod)-y, ens.z(sysid,bod)-z,
285  ens.vx(sysid,bod)-vx, ens.vy(sysid,bod)-vy, ens.vz(sysid,bod)-vz);
286  } // end loop over bodies
287 
288  } // end loop over systems
289  return ens;
290 }
291 
293 void print_system(const swarm::ensemble& ens, const int systemid, std::ostream &os = std::cout)
294 {
295  enum {
296  JACOBI, BARYCENTRIC, ASTROCENTRIC
297  } COORDINATE_SYSTEM = BARYCENTRIC;
298  const bool use_jacobi = cfg.optional("use_jacobi_output", 0);
299  if(use_jacobi) COORDINATE_SYSTEM = JACOBI;
300 
301  std::streamsize cout_precision_old = os.precision();
302  os.precision(10);
303  os << "sys_idx= " << systemid << " sys_id= " << ens[systemid].id() << " time= " << ens.time(systemid) << "\n";
304  double star_mass = ens.mass(systemid,0);
305  double mass_effective = star_mass;
306  double bx, by, bz, bvx, bvy, bvz;
307  switch(COORDINATE_SYSTEM)
308  {
309  case JACOBI:
310  ens.get_barycenter(systemid,bx,by,bz,bvx,bvy,bvz,0);
311  break;
312 
313  case BARYCENTRIC:
314  ens.get_barycenter(systemid,bx,by,bz,bvx,bvy,bvz);
315  break;
316 
317  case ASTROCENTRIC:
318  ens.get_body(systemid,0,star_mass,bx,by,bz,bvx,bvy,bvz);
319  break;
320  }
321 
322  for(unsigned int bod=1;bod<ens.nbod();++bod) // Skip star since printing orbits
323  {
324  // std::cout << "pos= (" << ens.x(systemid, bod) << ", " << ens.y(systemid, bod) << ", " << ens.z(systemid, bod) << ") vel= (" << ens.vx(systemid, bod) << ", " << ens.vy(systemid, bod) << ", " << ens.vz(systemid, bod) << ").\n";
325  double mass = ens.mass(systemid,bod);
326 
327  switch(COORDINATE_SYSTEM)
328  {
329  case JACOBI:
330  ens.get_barycenter(systemid,bx,by,bz,bvx,bvy,bvz,bod-1);
331  mass_effective += mass;
332  break;
333 
334  case BARYCENTRIC:
335  mass_effective = star_mass + mass;
336  break;
337 
338  case ASTROCENTRIC:
339  mass_effective = star_mass + mass;
340  break;
341  }
342  double x = ens.x(systemid,bod)-bx;
343  double y = ens.y(systemid,bod)-by;
344  double z = ens.z(systemid,bod)-bz;
345  double vx = ens.vx(systemid,bod)-bvx;
346  double vy = ens.vy(systemid,bod)-bvy;
347  double vz = ens.vz(systemid,bod)-bvz;
348 
349  double a, e, i, O, w, M;
350  calc_keplerian_for_cartesian(a,e,i,O,w,M, x,y,z,vx,vy,vz, mass_effective);
351  i *= 180/M_PI;
352  O *= 180/M_PI;
353  w *= 180/M_PI;
354  M *= 180/M_PI;
355  // os << " b= " << bod << " m= " << mass << " a= " << a << " e= " << e << " i= " << i << " Omega= " << O << " omega= " << w << " M= " << M << "\n";
356  os << ens[systemid].id() << " " << bod << " " << mass << " " << a << " " << e << " " << i << " " << O << " " << w << " " << M << "\n";
357  }
358 
359  os.precision(cout_precision_old);
360  os << std::flush;
361 }
362 
364 void print_selected_systems(swarm::ensemble& ens, std::vector<unsigned int> systemindices, std::ostream &os = std::cout)
365 {
366  for(unsigned int i=0; i<systemindices.size(); ++i)
367  print_system(ens,systemindices[i], os);
368 }
370 void print_selected_systems_for_demo(swarm::ensemble& ens, unsigned int nprint, std::ostream &os = std::cout)
371 {
372  if(nprint>ens.nsys()) nprint = ens.nsys();
373  for(unsigned int systemid = 0; systemid< nprint; ++systemid)
374  print_system(ens,systemid,os);
375 }
376 
378 void write_stable_systems(defaultEnsemble &ens, defaultEnsemble &ens_init)
379 {
380  // find the stable ones and output the initial conditions for the stable
381  // ones in keplerian coordinates
382  // std::cerr << "# Finding stable system ids\n";
383  std::vector<unsigned int> stable_system_indices, unstable_system_indices;
384  for(int i = 0; i < ens.nsys() ; i++ )
385  {
386  if(ens[i].is_disabled())
387  unstable_system_indices.push_back(i);
388  else
389  stable_system_indices.push_back(i);
390  }
391 
392  // std::cerr << "# Writing stable system ids\n";
393  if(cfg.count("stable_init_output"))
394  {
395  ofstream stable_init_output( cfg.optional("stable_init_output", string("stable_init_output.txt")).c_str() );
396  print_selected_systems(ens_init,stable_system_indices, stable_init_output);
397  }
398 
399  if(cfg.count("stable_final_output"))
400  {
401  ofstream stable_final_output( cfg.optional("stable_final_output", string("stable_final_output.txt")).c_str() );
402  print_selected_systems(ens,stable_system_indices, stable_final_output);
403  }
404 
405  if(cfg.count("unstable_init_output"))
406  {
407  ofstream unstable_init_output( cfg.optional("unstable_init_output", string("unstable_init_output.txt")).c_str() );
408  print_selected_systems(ens_init,unstable_system_indices, unstable_init_output);
409  }
410 
411  if(cfg.count("unstable_final_output"))
412  {
413  ofstream unstable_final_output( cfg.optional("unstable_final_output", string("unstable_final_output.txt")).c_str() );
414  print_selected_systems(ens,unstable_system_indices, unstable_final_output);
415  }
416 
417 }
419 std::vector<std::vector<double> > calc_semimajor_axes(defaultEnsemble& ens)
420 {
421  std::vector<std::vector<double> > semimajor_axes(ens.nsys(),std::vector<double>(ens.nbod(),0.));
422 
423  for(int sys_idx = 0; sys_idx < ens.nsys() ; sys_idx++)
424  {
425  int sys_id = ens[sys_idx].id();
426  assert(sys_id>=0);
427  assert(sys_id<ens.nsys());
428  double star_mass = ens.mass(sys_idx,0);
429  double mass_effective = star_mass;
430  double bx, by, bz, bvx, bvy, bvz;
431  ens.get_barycenter(sys_idx,bx,by,bz,bvx,bvy,bvz,0);
432  for(unsigned int bod=1;bod<ens.nbod();++bod) // Skip star since calculating orbits
433  {
434  double mass = ens.mass(sys_idx,bod);
435  mass_effective += mass;
436  ens.get_barycenter(sys_idx,bx,by,bz,bvx,bvy,bvz,bod-1);
437  double x = ens.x(sys_idx,bod)-bx;
438  double y = ens.y(sys_idx,bod)-by;
439  double z = ens.z(sys_idx,bod)-bz;
440  double vx = ens.vx(sys_idx,bod)-bvx;
441  double vy = ens.vy(sys_idx,bod)-bvy;
442  double vz = ens.vz(sys_idx,bod)-bvz;
443  double a, e, i, O, w, M;
444  calc_keplerian_for_cartesian(a,e,i,O,w,M, x,y,z,vx,vy,vz, mass_effective);
445  semimajor_axes[sys_id][bod-1] = a;
446  }
447  }
448  return semimajor_axes;
449 }
450 
451 void disable_unstable_systems(defaultEnsemble& ens, const std::vector<std::vector<double> >& semimajor_axes_init, const double deltaa_threshold )
452 {
453  for(int sys_idx = 0; sys_idx < ens.nsys() ; sys_idx++)
454  {
455  if(ens[sys_idx].is_disabled() ) continue;
456  bool disable = false;
457  int sys_id = ens[sys_idx].id();
458  double star_mass = ens.mass(sys_idx,0);
459  double mass_effective = star_mass;
460  double bx, by, bz, bvx, bvy, bvz;
461  ens.get_barycenter(sys_idx,bx,by,bz,bvx,bvy,bvz,0);
462  for(unsigned int bod=1;bod<ens.nbod();++bod) // Skip star since calculating orbits
463  {
464  // std::cout << "body= " << bod << ": ";
465  // std::cout << "pos= (" << ens.x(sys_idx, bod) << ", " << ens.y(sys_idx, bod) << ", " << ens.z(sys_idx, bod) << ") vel= (" << ens.vx(sys_idx, bod) << ", " << ens.vy(sys_idx, bod) << ", " << ens.vz(sys_idx, bod) << ").\n";
466  double mass = ens.mass(sys_idx,bod);
467  mass_effective += mass;
468  ens.get_barycenter(sys_idx,bx,by,bz,bvx,bvy,bvz,bod-1);
469  double x = ens.x(sys_idx,bod)-bx;
470  double y = ens.y(sys_idx,bod)-by;
471  double z = ens.z(sys_idx,bod)-bz;
472  double vx = ens.vx(sys_idx,bod)-bvx;
473  double vy = ens.vy(sys_idx,bod)-bvy;
474  double vz = ens.vz(sys_idx,bod)-bvz;
475  double a, e, i, O, w, M;
476  calc_keplerian_for_cartesian(a,e,i,O,w,M, x,y,z,vx,vy,vz, mass_effective);
477  i *= 180/M_PI;
478  O *= 180/M_PI;
479  w *= 180/M_PI;
480  M *= 180/M_PI;
481 
482  if(!((e>=0.)&&(e<1.0))) { disable = true; }
483  if(!((a>0.)&&(a<10.0))) { disable = true; }
484  assert(sys_id>=0);
485  assert(sys_id<semimajor_axes_init.size());
486  assert(bod-1>=0);
487  assert(bod-1<semimajor_axes_init[sys_id].size());
488  double da = a - semimajor_axes_init[sys_id][bod-1];
489  if(fabs(da) >= deltaa_threshold * semimajor_axes_init[sys_id][bod-1])
490  { disable = true; }
491 
492  if(disable)
493  {
494  if(cfg.count("verbose"))
495  std::cout << "# Disabling idx=" << sys_idx << " id=" << sys_id << " b=" << bod << " ainit= " << semimajor_axes_init[sys_id][bod-1] << " a= " << a << " e= " << e << " i= " << i << " Omega= " << O << " omega= " << w << " M= " << M << "\n";
496  break;
497  }
498  }
499  if(disable) ens[sys_idx].set_disabled();
500  }
501 }
503 bool needs_shrinking( const defaultEnsemble& ens )
504 {
505  // This is the ratio we use when we are shrinking
506  // if the ratio of active ones to all is less
507  // than this number then we trim the fat
508  const double critical_ratio = 0.75;
509 
510  int count_disabled = number_of_disabled_systems( ens ) ;
511  double ratio = double(count_disabled) / double( ens.nsys() );
512 
513  return ratio > critical_ratio;
514 }
515 
516 
520 void save_snapshot( defaultEnsemble& ens )
521 {
522  if(cfg.count("snapshot")) snapshot::save( ens, cfg["snapshot"] );
523 }
524 
525 
533 defaultEnsemble trim_disabled_systems( const defaultEnsemble& ens )
534 {
535  int nsys = ens.nsys();
536  int active_nsys = nsys - number_of_disabled_systems( ens ) ;
537  // WARNING: Needed to add this to prevent segfaults. TODO: Figure out why.
538  if(active_nsys==0) return ens;
539  int nbod = ens.nbod();
540 
541  defaultEnsemble active_ens = defaultEnsemble::create( nbod, active_nsys );
542 
543  // Copy the active ones to the new ensemble
544  for(int i = 0, j = 0; (i < nsys) && (j < active_nsys); i++)
545  {
546  if( !ens[i].is_disabled() )
547  {
548  ens[i].copyTo( active_ens[j] );
549  j++;
550  }
551  }
552 
553  return active_ens;
554 }
555 
556 
558 void reactivate_systems(defaultEnsemble&ens)
559 {
560  for(int i = 0; i < ens.nsys() ; i++)
561  {
562  if(ens[i].is_inactive())
563  ens.set_active(i);
564  }
565 }
567 volatile bool integration_loop_not_aborted_yet = true;
574 void ctrl_c_trap(int)
575 {
576  fprintf(stderr, "# Break requested... Finishing the current GPU kernel call and will then save the results before exitting.\n");
577  integration_loop_not_aborted_yet = false;
578 }
579 void catch_ctrl_c()
580 {
581  signal(SIGINT, &ctrl_c_trap );
582 }
583 
585 int main(int argc, char* argv[] )
586 {
587  // We keep it simple, later on one can use boost::program_options to
588  // have more options
589  // but now we only use two configuration files. It is because the
590  // initial conditions configuration file can get really big and
591  // it has very little to do with other configuration options.
592  if(argc < 3) cout << "Usage: montecarlo_mcmc_outputs <integration configuration> <initial conditions configuration>" << endl;
593 
594  // First one is the configuration for integration
595  string integ_configfile = argv[1];
596  // the second one is the configuration for generating
597  // initial conditions it is used by \ref generate_ensemble_with_randomized_initial_conditions
598  string initc_configfile = argv[2];
599 
600  cfg = config::load(integ_configfile);
601 
602  // 1.read keplerian coordinates from a file
603  // 2.generate guesses based on the keplerian coordinates
604  // 3.convert keplerian coordinates to an ensemble
605  // The following line that is taken from swarm_tutorial_montecarlo.cpp
606  // does the first three steps. Its pretty amazing.
607  defaultEnsemble ens ;
608  if( cfg.count("input") )
609  { ens = snapshot::load(cfg["input"]); }
610  else if(cfg.count("input_mcmc_keplerian") )
611  { ens = generate_ensemble_with_initial_conditions_keplerian_from_file( config::load(initc_configfile) ); }
612  else if(cfg.count("input_mcmc_cartesian") )
613  { ens = generate_ensemble_with_initial_conditions_cartesian_from_file( config::load(initc_configfile) ); }
614  else
615  {
616  std::cerr << "# Must specify one of input [for binary snapshot], input_mcmc_keplerian or input_mcmc_cartesian.\n";
617  return 255;
618  }
619 
620  // save the ensemble as a snapshot
621  if(cfg.count("initial_snapshot"))
622  { snapshot::save( ens, cfg["initial_snapshot"] ); }
623 
624  defaultEnsemble ens_init = ens.clone() ;
625  std::vector<std::vector<double> > semimajor_axes_init = calc_semimajor_axes(ens);
626 
627  // We want to find the ones that are really stable, so we integrate for
628  // a really long time and over time we get rid of unstable ones.
629  double destination_time = cfg.optional("destination_time", 1.0E6);
630 
631 
632  swarm::init(cfg);
633  Pintegrator integ = integrator::create(cfg);
634  integ->set_ensemble(ens);
635  integ->set_destination_time ( destination_time );
636  // We can set the following two items if we really need
637  // longer integrations before we stop for checking the
638  // ensemble and saving snapshots.
639  // one kernel call to allow for prompt CPU pruning of unstable systems
640  int max_kernel_calls_per_integrate = cfg.optional("max_kernel_calls_per_integrate",1);
641  // 10^2 inner orbits at 200 time steps per inner orbit
642  int max_itterations_per_kernel_call = cfg.optional("max_itterations_per_kernel_call",20000);
643  integ->set_max_attempts( max_kernel_calls_per_integrate );
644  integ->set_max_iterations ( max_itterations_per_kernel_call );
645  SYNC;
646 
647  // integrate ensemble
648  // - drop the unstable ones as you go
649  // - redistribute the systems for better efficiency
650  // - save the results periodically
651  // This can be an infitie loop. but we can always get out of this
652  // the best way to do it is to use Ctrl-C. Further wecan
653  // define Ctrl-C signal handler to get out
654  // of this loop in a safe way. But we really want this loop
655  // to run for a long time
656  reactivate_systems(ens);
657  // EBF Experiment trying to expose host log.
658  swarm::log::ensemble(*(swarm::log::manager::default_log()->get_hostlog()),ens);
659  integ->flush_log();
660 
661  catch_ctrl_c();
662  int num_integrate_calls = 0;
663  while( number_of_active_systems(ens) > 0 && integration_loop_not_aborted_yet ) {
664 
665  // 1. Integrate, we could use core_integrate but the general integrate
666  // saves all the headache. We should only use core_integrate if we are
667  // going to do everything on GPU, but since we are saving a snapshot in
668  // the middle, there's no point. It also has a nice for loop and can
669  // to several kernel calls.
670  integ->integrate();
671  ++num_integrate_calls;
672 
673  // 2. CPU-based tests to identify systems that can be terminated
674  int active_ones = number_of_active_systems(ens);
675  const double deltaa_frac_threshold = cfg.optional("deltaa_frac_threshold", 0.5);
676  disable_unstable_systems( ens, semimajor_axes_init, deltaa_frac_threshold );
677  double max_deltaE = find_max_energy_conservation_error(ens, ens_init );
678  std::cerr << "# Time: " << ens.time_ranges().min << " " << ens.time_ranges().median << " " << ens.time_ranges().max << " Systems: " << ens.nsys() << ", " << active_ones << ", ";
679  active_ones = number_of_active_systems(ens);
680  std::cerr << active_ones << " max|dE/E|= " << max_deltaE << "\n";
681 
682  // EBF Experiment trying to expose host log.
683  if(num_integrate_calls%10==0)
684  {
685  swarm::log::ensemble_enabled(*(swarm::log::manager::default_log()->get_hostlog()),ens);
686  integ->flush_log();
687  }
688 
689  // 3. Now we need to get rid of the inactive ones. There
690  // should be some criteria, whatever it is we are
691  // going to put it in a function and call it \ref needs_shrinking
692  if ( needs_shrinking( ens ) )
693  {
694  // Now this function will take care of trimming for us
695  // We need to create a new ensemble of a smaller size
696  // thats why we need to call set_ensemble again. because
697  // the GPU ensemble should get recreated for us.
698  // we need better memory management to safely allow
699  // the following statement. Right now we don't have a
700  // very good automatic memory management
701  // std::cerr << "# Removing disabled systems\n";
702  ens = trim_disabled_systems( ens );
703  // std::cerr << "# Testing if ens has any systems left.\n";
704  if(ens.nsys()>0)
705  {
706  // std::cerr << "# Setting ensemble for integrator\n";
707  // WARNING: This appears to be the cause of SEGFAULT
708  // if the ensemble is empty.
709  integ->set_ensemble( ens );
710  }
711  // std::cerr << "# Time: " << ens.time_ranges() << " Total Systems: " << ens.nsys() << " Active Systems: " << active_ones << ", ";
712  // active_ones = number_of_active_systems(ens);
713  // std::cerr << active_ones << "\n";
714  }
715 
716  // std::cerr << std::endl;
717 
718  // 4. We need to save a snapshot in case system crashes or someone
719  // gets tired of it and hits Ctrl-C
720  // std::cerr << "# Saving snapshot\n";
721  save_snapshot( ens );
722  write_stable_systems(ens,ens_init);
723  }
724 
725  // Now we are at the end of the system, before we examine
726  // the output we need to
727  // std::cerr << "# Removing disabled systems\n";
728  if(number_of_active_systems(ens)>0)
729  {
730  ens = trim_disabled_systems( ens );
731  // std::cerr << "# Saving snapshot\n";
732  save_snapshot( ens );
733  }
734 
735  write_stable_systems(ens,ens_init);
736 
737 #if ACCESS_HOST_ARRAY_OF_TRANSIT_TIMES
738 
739  static const int max_num_doubles_per_event = 2;
740  typedef swarm::event_record<max_num_doubles_per_event> event_record_type;
741  typedef std::vector<event_record_type> event_log_one_system_type;
742  typedef std::vector<event_log_one_system_type> event_log_one_code_type;
743  std::vector<event_log_one_code_type> event_log;
744 
745  event_log_one_code_type data = (static_cast<host_array_writer* >(swarm::log::manager::default_log()->get_writer().get()))->get_event_log_all_systems(0);
746  std::vector<std::vector<std::vector<double> > > transit_times_model(ens_init.nsys(), std::vector<std::vector<double> >(ens_init.nbod()));
747  std::vector<std::vector<std::vector<double> > > transit_durations_model(ens_init.nsys(), std::vector<std::vector<double> >(ens_init.nbod()));
748 
749  for(int sysid=0;sysid<data.size();++sysid)
750  {
751  for(int e=0;e<data[sysid].size();++e)
752  {
753  event_record_type er = data[sysid][e];
754  int bodid = er.bodid1;
755  double time = er.time;
756  double b = er.data[0];
757  double vproj = er.data[1];
758  transit_times_model[sysid][bodid].push_back(time);
759  transit_durations_model[sysid][bodid].push_back(sqrt(1.-b*b)/vproj);
760  }
761  }
762 
763  std::cout.precision(10);
764 
765  for(int sysid=0;sysid<transit_times_model.size();++sysid)
766  {
767  for(int bodid=0;bodid<transit_times_model[sysid].size();++bodid)
768  {
769  for(int trid=0;trid<transit_times_model[sysid][bodid].size();++trid)
770  {
771  std::cout << "s= " << sysid << " b= " << bodid << " n= " << trid << " t= " << transit_times_model[sysid][bodid][trid] << " D= " << transit_durations_model[sysid][bodid][trid] << " \n";
772  }
773  }
774  }
775 
776 #endif
777 
778 }
779