Swarm-NG  1.1
montecarlo.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 #define SYNC cudaThreadSynchronize()
44 
45 
46 using namespace swarm;
47 using namespace std;
48 
49 config cfg;
51 void inspect(defaultEnsemble &ens, const int sys, const int bod )
52 {
53  fprintf(stderr,"%d %d: %lg (%lg %lg %lg) (%lg %lg %lg) \n", sys, bod,
54  ens[sys][bod].mass(),
55  ens[sys][bod][0].pos(),
56  ens[sys][bod][1].pos(),
57  ens[sys][bod][2].pos(),
58  ens[sys][bod][0].vel(),
59  ens[sys][bod][1].vel(),
60  ens[sys][bod][2].vel()
61  );
62 }
64 void generate_initial_conditions_for_system(const config& cfg, defaultEnsemble &ens, const int sysidx, const int sysid)
65 {
66  double time_init = cfg.optional("time_init", 0.0);
67  const bool use_jacobi = cfg.optional("use_jacobi", 0);
68  bool keplerian_ephemeris = cfg.optional("use_keplerian", 1);
69  bool transit_ephemeris = cfg.optional("use_transit", 0);
70  if(transit_ephemeris) keplerian_ephemeris = 0;
71  assert(!(keplerian_ephemeris && transit_ephemeris)); // find a better way
72 
73  ens[sysidx].id() = sysid;
74  ens[sysidx].time() = time_init;
75  ens[sysidx].set_active();
76  // set sun to unit mass and at origin
77  double mass_star = cfg.optional("mass_star", 1.);
78  double x=0, y=0, z=0, vx=0, vy=0, vz=0;
79  ens.set_body(sysidx, 0, mass_star, x, y, z, vx, vy, vz);
80  if(ens[sysidx][0].num_attributes()>=1)
81  ens[sysidx][0].attribute(0) = cfg.optional("radius_star", 1.) * 0.00464912633;
82 
83  double mass_enclosed = mass_star;
84  for(unsigned int bod=1;bod<ens.nbod();++bod)
85  {
86  double mass_planet = draw_value_from_config(cfg,"mass",bod,0.,mass_star);
87  mass_enclosed += mass_planet;
88 
89  if(ens[sysidx][bod].num_attributes()>=1)
90  ens[sysidx][bod].attribute(0) = draw_value_from_config(cfg,"radius",bod,0.,200.) * 4.26349283e-5;
91 
92 
93  double a, e, i, O, w, M;
94  if(transit_ephemeris)
95  {
96  double period = draw_value_from_config(cfg,"period",bod,0.2,365250.);
97  double epoch = draw_value_from_config(cfg,"epoch",bod,0.2,365250.);
98  a = pow((period/365.25)*(period/365.25)*mass_enclosed,1.0/3.0);
99  e = draw_value_from_config(cfg,"ecc",bod,0.,1.);
100  i = draw_value_from_config(cfg,"inc",bod,-180.,180.);
101  O = draw_value_from_config(cfg,"node",bod,-720.,720.);
102  w = draw_value_from_config(cfg,"omega",bod,-720.,720.);
103 
104  i *= M_PI/180.;
105  O *= M_PI/180.;
106  w *= M_PI/180.;
107 
108  double T = (0.5*M_PI-w)+2.0*M_PI*((epoch/period)-floor(epoch/period));
109  double E = atan2(sqrt(1.-e)*(1.+e)*sin(T),e+cos(T));
110  M = E-e*sin(E);
111  }
112  else if (keplerian_ephemeris)
113  {
114  a = draw_value_from_config(cfg,"a",bod,0.001,10000.);
115  e = draw_value_from_config(cfg,"ecc",bod,0.,1.);
116  i = draw_value_from_config(cfg,"inc",bod,-180.,180.);
117  O = draw_value_from_config(cfg,"node",bod,-720.,720.);
118  w = draw_value_from_config(cfg,"omega",bod,-720.,720.);
119  M = draw_value_from_config(cfg,"meananom",bod,-720.,720.);
120 
121  i *= M_PI/180.;
122  O *= M_PI/180.;
123  w *= M_PI/180.;
124  M *= M_PI/180.;
125  }
126 
127  double mass = use_jacobi ? mass_enclosed : mass_star+mass_planet;
128  if(cfg.count("verbose")&&(sysid<10))
129  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';
130 
131  calc_cartesian_for_ellipse(x,y,z,vx,vy,vz, a, e, i, O, w, M, mass);
132 
133  // printf("%d %d: %lg (%lg %lg %lg) (%lg %lg %lg) \n", sysidx, bod, mass_planet, x,y,z,vx,vy,vz);
134 
135  if(use_jacobi)
136  {
137  double bx, by, bz, bvx, bvy, bvz;
138  ens.get_barycenter(sysidx,bx,by,bz,bvx,bvy,bvz,bod-1);
139  x += bx; y += by; z += bz;
140  vx += bvx; vy += bvy; vz += bvz;
141  }
142 
143  // assign body a mass, position and velocity
144  ens.set_body(sysidx, bod, mass_planet, x, y, z, vx, vy, vz);
145 
146  if(cfg.count("verbose")&&(sysid<10))
147  {
148  double x_t = ens.x(sysidx,bod);
149  double y_t = ens.y(sysidx,bod);
150  double z_t = ens.z(sysidx,bod);
151  double vx_t = ens.vx(sysidx,bod);
152  double vy_t = ens.vy(sysidx,bod);
153  double vz_t = ens.vz(sysidx,bod);
154 
155  std::cout << " x= " << x << "=" << x_t << " ";
156  std::cout << " y= " << y << "=" << y_t << " ";
157  std::cout << " z= " << z << "=" << z_t << " ";
158  std::cout << "vx= " << vx << "=" << vx_t << " ";
159  std::cout << "vy= " << vy << "=" << vy_t << " ";
160  std::cout << "vz= " << vz << "=" << vz_t << "\n";
161  }
162 
163  } // end loop over bodies
164 
165  // Shift into barycentric frame
166  ens.get_barycenter(sysidx,x,y,z,vx,vy,vz);
167  for(unsigned int bod=0;bod<ens.nbod();++bod)
168  {
169  ens.set_body(sysidx, bod, ens.mass(sysidx,bod),
170  ens.x(sysidx,bod)-x, ens.y(sysidx,bod)-y, ens.z(sysidx,bod)-z,
171  ens.vx(sysidx,bod)-vx, ens.vy(sysidx,bod)-vy, ens.vz(sysidx,bod)-vz);
172  } // end loop over bodies
173 }
174 
189 defaultEnsemble generate_ensemble_with_randomized_initial_conditions(const config& cfg)
190 {
191  defaultEnsemble ens = defaultEnsemble::create( cfg.require("nbod",0), cfg.require("nsys",0) );
192  std::cerr << "# nsystems = " << ens.nsys() << " nbodies/system = " << ens.nbod() << "\n";
193 
194  // std::cerr << "# Set initial time for all systems = ";
195  double time_init = cfg.optional("time_init", 0.0);
196  // std::cerr << time_init << ".\n";
197 
198  const bool use_jacobi = cfg.optional("use_jacobi", 0);
199  bool keplerian_ephemeris = cfg.optional("use_keplerian", 1);
200  bool transit_ephemeris = cfg.optional("use_transit", 0);
201  if(transit_ephemeris) keplerian_ephemeris = 0;
202  assert(!(keplerian_ephemeris && transit_ephemeris)); // find a better way
203 
204  for(unsigned int sys=0;sys<ens.nsys();++sys)
205  {
206  generate_initial_conditions_for_system(cfg,ens,sys,sys);
207  } // end loop over systems
208  return ens;
209 }
210 
211 
212 void print_system(const swarm::ensemble& ens, const int systemid, std::ostream &os = std::cout)
213 {
214  enum {
215  JACOBI, BARYCENTRIC, ASTROCENTRIC
216  } COORDINATE_SYSTEM = BARYCENTRIC;
217  const bool use_jacobi = cfg.optional("use_jacobi", 0);
218  if(use_jacobi) COORDINATE_SYSTEM = JACOBI;
219 
220  std::streamsize cout_precision_old = os.precision();
221  os.precision(10);
222  os << "sys_idx= " << systemid << " sys_id= " << ens[systemid].id() << " time= " << ens.time(systemid) << "\n";
223  double star_mass = ens.mass(systemid,0);
224  double mass_effective = star_mass;
225  double bx, by, bz, bvx, bvy, bvz;
226  switch(COORDINATE_SYSTEM)
227  {
228  case JACOBI:
229  ens.get_barycenter(systemid,bx,by,bz,bvx,bvy,bvz,0);
230  break;
231 
232  case BARYCENTRIC:
233  ens.get_barycenter(systemid,bx,by,bz,bvx,bvy,bvz);
234  break;
235 
236  case ASTROCENTRIC:
237  ens.get_body(systemid,0,star_mass,bx,by,bz,bvx,bvy,bvz);
238  break;
239  }
240 
241  for(unsigned int bod=1;bod<ens.nbod();++bod) // Skip star since printing orbits
242  {
243  // 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";
244  double mass = ens.mass(systemid,bod);
245 
246  switch(COORDINATE_SYSTEM)
247  {
248  case JACOBI:
249  ens.get_barycenter(systemid,bx,by,bz,bvx,bvy,bvz,bod-1);
250  mass_effective += mass;
251  break;
252 
253  case BARYCENTRIC:
254  mass_effective = star_mass + mass;
255  break;
256 
257  case ASTROCENTRIC:
258  mass_effective = star_mass + mass;
259  break;
260  }
261  double x = ens.x(systemid,bod)-bx;
262  double y = ens.y(systemid,bod)-by;
263  double z = ens.z(systemid,bod)-bz;
264  double vx = ens.vx(systemid,bod)-bvx;
265  double vy = ens.vy(systemid,bod)-bvy;
266  double vz = ens.vz(systemid,bod)-bvz;
267 
268  double a, e, i, O, w, M;
269  calc_keplerian_for_cartesian(a,e,i,O,w,M, x,y,z,vx,vy,vz, mass_effective);
270  i *= 180/M_PI;
271  O *= 180/M_PI;
272  w *= 180/M_PI;
273  M *= 180/M_PI;
274  // os << " b= " << bod << " m= " << mass << " a= " << a << " e= " << e << " i= " << i << " Omega= " << O << " omega= " << w << " M= " << M << "\n";
275  os << ens[systemid].id() << " " << bod << " " << mass << " " << a << " " << e << " " << i << " " << O << " " << w << " " << M << "\n";
276  }
277 
278  os.precision(cout_precision_old);
279  os << std::flush;
280 }
281 
282 
283 void print_selected_systems(swarm::ensemble& ens, std::vector<unsigned int> systemindices, std::ostream &os = std::cout)
284 {
285  for(unsigned int i=0; i<systemindices.size(); ++i)
286  print_system(ens,systemindices[i], os);
287 }
288 
289 void print_selected_systems_for_demo(swarm::ensemble& ens, unsigned int nprint, std::ostream &os = std::cout)
290 {
291  if(nprint>ens.nsys()) nprint = ens.nsys();
292  for(unsigned int systemid = 0; systemid< nprint; ++systemid)
293  print_system(ens,systemid,os);
294 }
295 
296 
297 void write_stable_systems(defaultEnsemble &ens, defaultEnsemble &ens_init)
298 {
299  // find the stable ones and output the initial conditions for the stable
300  // ones in keplerian coordinates
301  // std::cerr << "# Finding stable system ids\n";
302  std::vector<unsigned int> stable_system_indices, unstable_system_indices;
303  for(int i = 0; i < ens.nsys() ; i++ )
304  {
305  if(ens[i].is_disabled())
306  unstable_system_indices.push_back(i);
307  else
308  stable_system_indices.push_back(i);
309  }
310 
311  // std::cerr << "# Writing stable system ids\n";
312  if(cfg.count("stable_init_output"))
313  {
314  ofstream stable_init_output( cfg.optional("stable_init_output", string("stable_init_output.txt")).c_str() );
315  print_selected_systems(ens_init,stable_system_indices, stable_init_output);
316  }
317 
318  if(cfg.count("stable_final_output"))
319  {
320  ofstream stable_final_output( cfg.optional("stable_final_output", string("stable_final_output.txt")).c_str() );
321  print_selected_systems(ens,stable_system_indices, stable_final_output);
322  }
323 
324  if(cfg.count("unstable_init_output"))
325  {
326  ofstream unstable_init_output( cfg.optional("unstable_init_output", string("unstable_init_output.txt")).c_str() );
327  print_selected_systems(ens_init,unstable_system_indices, unstable_init_output);
328  }
329 
330  if(cfg.count("unstable_final_output"))
331  {
332  ofstream unstable_final_output( cfg.optional("unstable_final_output", string("unstable_final_output.txt")).c_str() );
333  print_selected_systems(ens,unstable_system_indices, unstable_final_output);
334  }
335 
336 }
337 
339 std::vector<std::vector<double> > calc_semimajor_axes(defaultEnsemble& ens)
340 {
341  std::vector<std::vector<double> > semimajor_axes(ens.nsys(),std::vector<double>(ens.nbod(),0.));
342 
343  for(int sys_idx = 0; sys_idx < ens.nsys() ; sys_idx++)
344  {
345  int sys_id = ens[sys_idx].id();
346  assert(sys_id>=0);
347  assert(sys_id<ens.nsys());
348  double star_mass = ens.mass(sys_idx,0);
349  double mass_effective = star_mass;
350  double bx, by, bz, bvx, bvy, bvz;
351  ens.get_barycenter(sys_idx,bx,by,bz,bvx,bvy,bvz,0);
352  for(unsigned int bod=1;bod<ens.nbod();++bod) // Skip star since calculating orbits
353  {
354  double mass = ens.mass(sys_idx,bod);
355  mass_effective += mass;
356  ens.get_barycenter(sys_idx,bx,by,bz,bvx,bvy,bvz,bod-1);
357  double x = ens.x(sys_idx,bod)-bx;
358  double y = ens.y(sys_idx,bod)-by;
359  double z = ens.z(sys_idx,bod)-bz;
360  double vx = ens.vx(sys_idx,bod)-bvx;
361  double vy = ens.vy(sys_idx,bod)-bvy;
362  double vz = ens.vz(sys_idx,bod)-bvz;
363  double a, e, i, O, w, M;
364  calc_keplerian_for_cartesian(a,e,i,O,w,M, x,y,z,vx,vy,vz, mass_effective);
365  semimajor_axes[sys_id][bod-1] = a;
366  }
367  }
368  return semimajor_axes;
369 }
371 void disable_unstable_systems(defaultEnsemble& ens, const std::vector<std::vector<double> >& semimajor_axes_init, const double deltaa_threshold )
372 {
373  for(int sys_idx = 0; sys_idx < ens.nsys() ; sys_idx++)
374  {
375  if(ens[sys_idx].is_disabled() ) continue;
376  bool disable = false;
377  int sys_id = ens[sys_idx].id();
378  double star_mass = ens.mass(sys_idx,0);
379  double mass_effective = star_mass;
380  double bx, by, bz, bvx, bvy, bvz;
381  ens.get_barycenter(sys_idx,bx,by,bz,bvx,bvy,bvz,0);
382  for(unsigned int bod=1;bod<ens.nbod();++bod) // Skip star since calculating orbits
383  {
384  // std::cout << "body= " << bod << ": ";
385  // 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";
386  double mass = ens.mass(sys_idx,bod);
387  mass_effective += mass;
388  ens.get_barycenter(sys_idx,bx,by,bz,bvx,bvy,bvz,bod-1);
389  double x = ens.x(sys_idx,bod)-bx;
390  double y = ens.y(sys_idx,bod)-by;
391  double z = ens.z(sys_idx,bod)-bz;
392  double vx = ens.vx(sys_idx,bod)-bvx;
393  double vy = ens.vy(sys_idx,bod)-bvy;
394  double vz = ens.vz(sys_idx,bod)-bvz;
395  double a, e, i, O, w, M;
396  calc_keplerian_for_cartesian(a,e,i,O,w,M, x,y,z,vx,vy,vz, mass_effective);
397  i *= 180/M_PI;
398  O *= 180/M_PI;
399  w *= 180/M_PI;
400  M *= 180/M_PI;
401 
402  if(!((e>=0.)&&(e<1.0))) { disable = true; }
403  if(!((a>0.)&&(a<10.0))) { disable = true; }
404  assert(sys_id>=0);
405  assert(sys_id<semimajor_axes_init.size());
406  assert(bod-1>=0);
407  assert(bod-1<semimajor_axes_init[sys_id].size());
408  double da = a - semimajor_axes_init[sys_id][bod-1];
409  if(fabs(da) >= deltaa_threshold * semimajor_axes_init[sys_id][bod-1])
410  { disable = true; }
411 
412  if(disable)
413  {
414  if(cfg.count("verbose"))
415  std::cout << "# Disabling idx=" << sys_idx << " id=" << sys_id << " b=" << bod << " a= " << a << " e= " << e << " i= " << i << " Omega= " << O << " omega= " << w << " M= " << M << "\n";
416  break;
417  }
418  }
419  if(disable) ens[sys_idx].set_disabled();
420  }
421 }
422 
423 bool needs_shrinking( const defaultEnsemble& ens )
424 {
425  // This is the ratio we use when we are shrinking
426  // if the ratio of active ones to all is less
427  // than this number then we trim the fat
428  const double critical_ratio = 0.75;
429 
430  int count_disabled = number_of_disabled_systems( ens ) ;
431  double ratio = double(count_disabled) / double( ens.nsys() );
432 
433  return ratio > critical_ratio;
434 }
435 
436 
440 void save_snapshot( defaultEnsemble& ens )
441 {
442  if(cfg.count("snapshot")) snapshot::save( ens, cfg["snapshot"] );
443 }
444 
445 
453 defaultEnsemble trim_disabled_systems( const defaultEnsemble& ens )
454 {
455  int nsys = ens.nsys();
456  int active_nsys = nsys - number_of_disabled_systems( ens ) ;
457  // WARNING: Needed to add this to prevent segfaults. TODO: Figure out why.
458  if(active_nsys==0) return ens;
459  int nbod = ens.nbod();
460 
461  defaultEnsemble active_ens = defaultEnsemble::create( nbod, active_nsys );
462 
463  // Copy the active ones to the new ensemble
464  for(int i = 0, j = 0; (i < nsys) && (j < active_nsys); i++)
465  {
466  if( !ens[i].is_disabled() )
467  {
468  ens[i].copyTo( active_ens[j] );
469  j++;
470  }
471  }
472 
473  return active_ens;
474 }
475 
476 
477 void reactivate_systems(defaultEnsemble&ens)
478 {
479  for(int i = 0; i < ens.nsys() ; i++)
480  {
481  if(ens[i].is_inactive())
482  ens.set_active(i);
483  }
484 }
485 
486 volatile bool integration_loop_not_aborted_yet = true;
493 void ctrl_c_trap(int)
494 {
495  fprintf(stderr, "# Break requested... Finishing the current GPU kernel call and will then save the results before exitting.\n");
496  integration_loop_not_aborted_yet = false;
497 }
498 void catch_ctrl_c()
499 {
500  signal(SIGINT, &ctrl_c_trap );
501 }
502 
503 int main(int argc, char* argv[] )
504 {
505  // We keep it simple, later on one can use boost::program_options to
506  // have more options
507  // but now we only use two configuration files. It is because the
508  // initial conditions configuration file can get really big and
509  // it has very little to do with other configuration options.
510  if(argc < 3) cout << "Usage: montecarlo <integration configuration> <initial conditions configuration>" << endl;
511 
512  // First one is the configuration for integration
513  string integ_configfile = argv[1];
514  // the second one is the configuration for generating
515  // initial conditions it is used by \ref generate_ensemble_with_randomized_initial_conditions
516  string initc_configfile = argv[2];
517 
518  cfg = config::load(integ_configfile);
519 
520  // 1.read keplerian coordinates from a file
521  // 2.generate guesses based on the keplerian coordinates
522  // 3.convert keplerian coordinates to an ensemble
523  // The following line that is taken from swarm_tutorial_montecarlo.cpp
524  // does the first three steps. Its pretty amazing.
525  defaultEnsemble ens ;
526  if( cfg.count("input") )
527  { ens = snapshot::load(cfg["input"]); }
528  else
529  { ens = generate_ensemble_with_randomized_initial_conditions( config::load(initc_configfile) ); }
530 
531  // save the ensemble as a snapshot
532  if(cfg.count("initial_snapshot"))
533  { snapshot::save( ens, cfg["initial_snapshot"] ); }
534 
535  defaultEnsemble ens_init = ens.clone() ;
536  std::vector<std::vector<double> > semimajor_axes_init = calc_semimajor_axes(ens);
537 
538  // We want to find the ones that are really stable, so we integrate for
539  // a really long time and over time we get rid of unstable ones.
540  double destination_time = cfg.optional("destination_time", 1.0E6);
541 
542 
543  swarm::init(cfg);
544  Pintegrator integ = integrator::create(cfg);
545  integ->set_ensemble(ens);
546  integ->set_destination_time ( destination_time );
547  // We can set the following two items if we really need
548  // longer integrations before we stop for checking the
549  // ensemble and saving snapshots.
550  // one kernel call to allow for prompt CPU pruning of unstable systems
551  int max_kernel_calls_per_integrate = cfg.optional("max_kernel_calls_per_integrate",1);
552  // 10^2 inner orbits at 200 time steps per inner orbit
553  int max_itterations_per_kernel_call = cfg.optional("max_itterations_per_kernel_call",20000);
554  integ->set_max_attempts( max_kernel_calls_per_integrate );
555  integ->set_max_iterations ( max_itterations_per_kernel_call );
556  SYNC;
557 
558  // integrate ensemble
559  // - drop the unstable ones as you go
560  // - redistribute the systems for better efficiency
561  // - save the results periodically
562  // This can be an infitie loop. but we can always get out of this
563  // the best way to do it is to use Ctrl-C. Further wecan
564  // define Ctrl-C signal handler to get out
565  // of this loop in a safe way. But we really want this loop
566  // to run for a long time
567  reactivate_systems(ens);
568  // EBF Experiment trying to expose host log.
569  swarm::log::ensemble(*(swarm::log::manager::default_log()->get_hostlog()),ens);
570  integ->flush_log();
571 
572  catch_ctrl_c();
573  while( number_of_active_systems(ens) > 0 && integration_loop_not_aborted_yet ) {
574 
575  // 1. Integrate, we could use core_integrate but the general integrate
576  // saves all the headache. We should only use core_integrate if we are
577  // going to do everything on GPU, but since we are saving a snapshot in
578  // the middle, there's no point. It also has a nice for loop and can
579  // to several kernel calls.
580  integ->integrate();
581 
582  // 2. CPU-based tests to identify systems that can be terminated
583  int active_ones = number_of_active_systems(ens);
584  const double deltaa_frac_threshold = cfg.optional("deltaa_frac_threshold", 0.5);
585  disable_unstable_systems( ens, semimajor_axes_init, deltaa_frac_threshold );
586  std::cerr << "# Time: " << ens.time_ranges().min << " " << ens.time_ranges().median << " " << ens.time_ranges().max << " Systems: " << ens.nsys() << ", " << active_ones << ", ";
587  active_ones = number_of_active_systems(ens);
588  std::cerr << active_ones << "\n";
589 
590  // EBF Experiment trying to expose host log.
591  swarm::log::ensemble_enabled(*(swarm::log::manager::default_log()->get_hostlog()),ens);
592  integ->flush_log();
593 
594  // 3. Now we need to get rid of the inactive ones. There
595  // should be some criteria, whatever it is we are
596  // going to put it in a function and call it \ref needs_shrinking
597  if ( needs_shrinking( ens ) )
598  {
599  // Now this function will take care of trimming for us
600  // We need to create a new ensemble of a smaller size
601  // thats why we need to call set_ensemble again. because
602  // the GPU ensemble should get recreated for us.
603  // we need better memory management to safely allow
604  // the following statement. Right now we don't have a
605  // very good automatic memory management
606  // std::cerr << "# Removing disabled systems\n";
607  ens = trim_disabled_systems( ens );
608  // std::cerr << "# Testing if ens has any systems left.\n";
609  if(ens.nsys()>0)
610  {
611  // std::cerr << "# Setting ensemble for integrator\n";
612  // WARNING: This appears to be the cause of SEGFAULT
613  // if the ensemble is empty.
614  integ->set_ensemble( ens );
615  }
616  // std::cerr << "# Time: " << ens.time_ranges() << " Total Systems: " << ens.nsys() << " Active Systems: " << active_ones << ", ";
617  // active_ones = number_of_active_systems(ens);
618  // std::cerr << active_ones << "\n";
619  }
620 
621  // std::cerr << std::endl;
622 
623  // 4. We need to save a snapshot in case system crashes or someone
624  // gets tired of it and hits Ctrl-C
625  // std::cerr << "# Saving snapshot\n";
626  save_snapshot( ens );
627  write_stable_systems(ens,ens_init);
628  }
629 
630  // Now we are at the end of the system, before we examine
631  // the output we need to
632  // std::cerr << "# Removing disabled systems\n";
633  if(number_of_active_systems(ens)>0)
634  {
635  ens = trim_disabled_systems( ens );
636  // std::cerr << "# Saving snapshot\n";
637  save_snapshot( ens );
638  }
639 
640  write_stable_systems(ens,ens_init);
641 
642 }
643