RInside Version 0.2.16
rinside_interactive0.cpp
Go to the documentation of this file.
1 // -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; tab-width: 8; -*-
2 //
3 // Example of a planetary motion solver with interactive console
4 //
5 // Copyright (C) 2009 Dirk Eddelbuettel
6 // Copyright (C) 2010 - 2017 Dirk Eddelbuettel and Romain Francois
7 // Copyright (C) 2017 Dirk Eddelbuettel, Romain Francois and Łukasz Łaniewski-Wołłk
8 //
9 // GPL'ed
10 
11 #include <RInside.h> // for the embedded R via RInside
12 
13 class Wrapper;
14 
15 // A planet.
16 struct Planet {
17  double x,y;
18  double m;
19  double vx, vy;
20 };
21 
22 // A gravity simulator
23 class Solver {
24  typedef std::vector<Planet> Planets;
25  Planets tab;
26  double dt;
27  double G;
28  void Iteration() {
29  for (Planets::iterator a=tab.begin(); a != tab.end(); a++) {
30  for (Planets::iterator b=tab.begin(); b != a; b++) {
31  double x = a->x - b->x;
32  double y = a->y - b->y;
33  double r = sqrt(x*x + y*y);
34  double f = a->m * b->m * G / (r*r+1);
35  double fx = f * x/r;
36  double fy = f * y/r;
37  a->vx -= dt * fx / a->m;
38  a->vy -= dt * fy / a->m;
39  b->vx += dt * fx / b->m;
40  b->vy += dt * fy / b->m;
41  }
42  }
43  for (Planets::iterator a=tab.begin(); a != tab.end(); a++) {
44  a->x += dt * a->vx;
45  a->y += dt * a->vy;
46  }
47  }
48 public:
49  Solver(int n): tab(n), dt(1.0e-4), G(1.0) {
50  double v=0;
51  for (Planets::iterator a=tab.begin(); a != tab.end(); a++) {
52  a->x = sin(v);
53  a->y = cos(v);
54  a->m = 1;
55  v += 3.0/n;
56  }
57  }
58  void Iterate(int n) {
59  for (int i=0;i<n;i++) Iteration();
60  }
61  friend class Wrapper;
62 };
63 
64 // A nice wrapper for the solver
65 class Wrapper {
66  Solver * s;
67 public:
68  Wrapper(Solver * s_) : s(s_) {}
69  Rcpp::DataFrame getData() {
70  Rcpp::NumericVector x,y,m,vx,vy;
71  for (Solver::Planets::iterator a=s->tab.begin(); a != s->tab.end(); a++) {
72  x.push_back( a->x);
73  y.push_back( a->y);
74  m.push_back( a->m);
75  vx.push_back(a->vx);
76  vy.push_back(a->vy);
77  }
78  return Rcpp::DataFrame::create(Rcpp::Named("x") = x,
79  Rcpp::Named("y") = y,
80  Rcpp::Named("mass") = m,
81  Rcpp::Named("Vx") = vx,
82  Rcpp::Named("Vy") = vy);
83  }
84  void setData(Rcpp::DataFrame tab) {
85  if ((size_t)tab.nrows() != s->tab.size()) {
86  return;
87  }
88  Rcpp::NumericVector x = tab["x"];
89  Rcpp::NumericVector y = tab["y"];
90  Rcpp::NumericVector m = tab["mass"];
91  Rcpp::NumericVector vx = tab["Vy"];
92  Rcpp::NumericVector vy = tab["Vy"];
93  for (int i=0;i<tab.nrows();i++) {
94  s->tab[i].x = x[i];
95  s->tab[i].y = y[i];
96  s->tab[i].m = m[i];
97  s->tab[i].vx = vx[i];
98  s->tab[i].vy = vy[i];
99  }
100  }
101  double& G() {
102  return s->G;
103  }
104  double& dt() {
105  return s->dt;
106  }
107 };
108 
109 // The function which is called when running Solver$...
110 SEXP Dollar(Rcpp::XPtr<Wrapper> obj, std::string name) {
111  if (name == "data") {
112  return obj->getData();
113  } else if (name == "G") {
114  return Rcpp::NumericVector(1,obj->G());
115  } else if (name == "dt") {
116  return Rcpp::NumericVector(1,obj->dt());
117  } else {
118  return NULL;
119  }
120 }
121 
122 // The function which is called when assigning to Solver$...
123 Rcpp::XPtr<Wrapper> DollarAssign(Rcpp::XPtr<Wrapper> obj, std::string name, SEXP v) {
124  if (name == "data") {
125  obj->setData(v);
126  } else if (name == "G") {
127  obj->G() = Rcpp::NumericVector(v)[0];
128  } else if (name == "dt") {
129  obj->dt() = Rcpp::NumericVector(v)[0];
130  }
131  return obj;
132 }
133 
134 // The function listing the elements of Solver
135 Rcpp::CharacterVector Names(Rcpp::XPtr<Wrapper> obj) {
136  Rcpp::CharacterVector ret;
137  ret.push_back("data");
138  ret.push_back("G");
139  ret.push_back("dt");
140  return ret;
141 }
142 
143 int main(int argc, char *argv[]) {
144  Solver S(70); // Creating the gravity simulator
145 
146  RInside R(argc, argv, false, false, true); // Create an embedded R instance
147  Rcpp::XPtr<Wrapper> wr(new Wrapper(&S)); // Wrapping the solver
148  wr.attr("class") = "Solver";
149  R["Solver"] = wr; // Adding the wrapped solver
150  R["$.Solver"] = Rcpp::InternalFunction(& Dollar); // Adding the functions
151  R["$<-.Solver"] = Rcpp::InternalFunction(& DollarAssign);
152  R["names.Solver"] = Rcpp::InternalFunction(& Names);
153  char type;
154  do {
155  std::cout << "Want to go interactive? [y/n]";
156  std::cin >> type;
157  } while ( !std::cin.fail() && type!='y' && type!='n' );
158  if (type == 'y') { // Running an interactive R session
159  std::cout << "Running an interactive R session. You can explore (and modify) the data in the 'Solver' object." << std::endl;
160  std::cout << "[ You can finish the with Ctrl+D ]" << std::endl;
161  R.parseEval("options(prompt = 'R console > ')");
162  R.parseEval("X11()");
163  R.repl() ;
164  std::cout << "Finishing the interactive mode" << std::endl;
165  R.parseEval("dev.off()");
166  }
167  R.parseEval("X11()");
168  for (int i=0;i<2000;i++) { // Running the some 200'000 iterations in non-interactive mode
169  S.Iterate(100);
170  R.parseEval("P = Solver$data;");
171  R.parseEval("plot(P$x,P$y,cex=P$mass,asp=1,xlim=c(-5,5),ylim=c(-5,5),xlab='X',ylab='Y')");
172  }
173  exit(0);
174 }
175 
SEXP Dollar(Rcpp::XPtr< Wrapper > obj, std::string name)
Wrapper(Solver *s_)
void repl()
Definition: RInside.cpp:424
int parseEval(const std::string &line, SEXP &ans)
Definition: RInside.cpp:326
std::vector< Planet > Planets
Rcpp::CharacterVector Names(Rcpp::XPtr< Wrapper > obj)
Rcpp::DataFrame getData()
int main(int argc, char *argv[])
void Iterate(int n)
void setData(Rcpp::DataFrame tab)
Rcpp::XPtr< Wrapper > DollarAssign(Rcpp::XPtr< Wrapper > obj, std::string name, SEXP v)