JEMRIS  2.8.1
open-source MRI simulations
ginac_functions.h
Go to the documentation of this file.
1 
5 /*
6  * JEMRIS Copyright (C)
7  * 2006-2015 Tony Stoecker
8  * 2007-2015 Kaveh Vahedipour
9  * 2009-2015 Daniel Pflugfelder
10  *
11  *
12  * This program is free software; you can redistribute it and/or modify
13  * it under the terms of the GNU General Public License as published by
14  * the Free Software Foundation; either version 2 of the License, or
15  * (at your option) any later version.
16  *
17  * This program is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU General Public License for more details.
21  *
22  * You should have received a copy of the GNU General Public License
23  * along with this program; if not, write to the Free Software
24  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
25  */
26 
27 #ifndef GINAC_FUNCTIONS_H_
28 #define GINAC_FUNCTIONS_H_
29 
30 #include <ginac/ginac.h>
31 #include <ginac/power.h>
32 #include <vector>
33 #ifndef __MINGW32__
34 #include <dlfcn.h>
35 #endif
36 #include <fstream>
37 #include <ios>
38 #include <sstream>
39 #include <stdexcept>
40 #include <string>
41 #include <stdlib.h>
42 #include "config.h"
43 
44 #if defined(_MSC_VER)
45 #define MKTEMP(path,n) _mktemp_s(path, n)
46 #elif defined(__MINGW32__)
47 #include <io.h>
48 #define MKTEMP(path,n) _mktemp(path)
49 #else
50 #ifdef HAVE_MKSTEMPS
51 extern "C" int mkstemps (char *path, int len);
52 #define MKTEMP(path,n) mkstemps(path,n)
53 #endif
54 #endif
55 
56 using namespace std;
57 
58 using namespace GiNaC ;
59 
60 
61 
68 const symbol & get_symbol(const string & sym_name)
69 {
70  static map<string, symbol> symbol_list;
71  map<string, symbol>::iterator it = symbol_list.find(sym_name);
72  if (it != symbol_list.end())
73  return it->second;
74  else
75  return symbol_list.insert(pair<string, symbol>(sym_name, symbol(sym_name))).first->second;
76 }
77 
78 
83 /*
84 static ex sqrt(const ex &x){ return pow(x,0.5); }
85 static ex sqrt_evalf(const ex & x){ return (ex_to<numeric>(pow(x,0.5))).to_double(); }
86 static ex sqrt_deriv(const ex & x, unsigned diff_param) { return 0.5*pow(x,-0.5); }
87 static ex sqrt_real_part(const ex & x) { return real_part(pow(x,0.5)); }
88 static ex sqrt_imag_part(const ex & x) { return imag_part(pow(x,0.5)); }
89 static ex sqrt_conjugate(const ex & x) { return pow(x.conjugate(),0.5); }
90 DECLARE_FUNCTION_1P(sqrt)
91 REGISTER_FUNCTION(sqrt, eval_func (sqrt_evalf ).
92  evalf_func (sqrt_evalf ).
93  derivative_func(sqrt_deriv ).
94  real_part_func (sqrt_real_part).
95  imag_part_func (sqrt_imag_part).
96  conjugate_func (sqrt_conjugate))
97 */
98 
104 static ex sinc_eval(const ex &x){
105  if ( x.is_zero() )
106  return 1;
107  else
108  return (sin(x)/x);
109 }
110 
111 static ex sinc_evalf(const ex & x){
112  //if (is_exactly_a<numeric>(x)) {
113  double d = (ex_to<numeric>(x)).to_double();
114  return ( d==0.0 ? 1 : sin(d)/d );
115  //}
116  //else
117  // return (sin(x)/x);
118  }
119 
120 static ex sinc_deriv(const ex & x, unsigned diff_param) {
121  if ( x.is_zero() )
122  return 0;
123  else
124  return cos(x)/x - sin(x)/(x*x);
125 }
126 
127 //real(sinc(a+I*b)) = (a*cosh(b)*sin(a) + b*cos(a)*sinh(b))/(a^2 + b^2)
128 static ex sinc_real_part(const ex & x)
129 {
130  if ( x.is_zero() )
131  return 1;
132  else
133  return ( ( ( real_part(x)*cosh(imag_part(x))*sin(real_part(x)) ) +
134  ( imag_part(x)*sinh(imag_part(x))*cos(real_part(x)) ) ) /
135  ( pow(real_part(x),2) + pow(imag_part(x),2) ) );
136 }
137 
138 //imag(sinc(a+I*b)) = (a*cos(a)*sinh(b) - b*cosh(b)*sin(a))/(a^2 + b^2)
139 static ex sinc_imag_part(const ex & x)
140 {
141  if ( x.is_zero() )
142  return 0;
143  else
144  return ( ( ( real_part(x)*sinh(imag_part(x))*cos(real_part(x)) ) -
145  ( imag_part(x)*cosh(imag_part(x))*sin(real_part(x)) ) ) /
146  ( pow(real_part(x),2) + pow(imag_part(x),2) ) );
147 }
148 
149 static ex sinc_conjugate(const ex & x) {
150  if ( x.is_zero() )
151  return 1;
152  else
153  return sin(x.conjugate())/x.conjugate();
154 }
155 
156 DECLARE_FUNCTION_1P(sinc)
157 REGISTER_FUNCTION(sinc, eval_func (sinc_eval ).
158  evalf_func (sinc_evalf ).
159  derivative_func(sinc_deriv ).
160  real_part_func (sinc_real_part).
161  imag_part_func (sinc_imag_part).
162  conjugate_func (sinc_conjugate))
163 
164 
170 static ex floor_evalf(const ex &x){
171  if (!is_a<numeric>(x) ) return x;
172  ex xn = ex_to<numeric>(x);
173  return ((int) ex_to<numeric>(xn).to_double() );
174 }
175 DECLARE_FUNCTION_1P(floor)
176 REGISTER_FUNCTION(floor, evalf_func(floor_evalf))
177 
178 
183 static ex mod_evalf(const ex &x, const ex &y){
184  if (!is_a<numeric>(x) || !is_a<numeric>(y)) return 0;
185  ex xn = ex_to<numeric>(x);
186  ex yn = ex_to<numeric>(y);
187  return (xn - floor_evalf(xn/yn)*yn);
188 }
189 DECLARE_FUNCTION_2P(mod)
190 REGISTER_FUNCTION(mod, evalf_func(mod_evalf))
191 
192 
197 static ex equal_evalf(const ex &x, const ex &y){
198  if (!is_a<numeric>(x) || !is_a<numeric>(y)) return 0;
199  ex xn = ex_to<numeric>(x);
200  ex yn = ex_to<numeric>(y);
201  int b = ((int) ( ex_to<numeric>(xn).to_double()
202  == ex_to<numeric>(yn).to_double() ) );
203  return b;
204 }
205 DECLARE_FUNCTION_2P(equal)
206 REGISTER_FUNCTION(equal, evalf_func(equal_evalf))
207 
208 
213 static ex gt_evalf(const ex &x, const ex &y){
214  if (!is_a<numeric>(x) || !is_a<numeric>(y)) return 0;
215  ex xn = ex_to<numeric>(x);
216  ex yn = ex_to<numeric>(y);
217  int b = ((int) ( ex_to<numeric>(xn).to_double()
218  > ex_to<numeric>(yn).to_double() ) );
219  return b;
220 }
221 DECLARE_FUNCTION_2P(gt)
222 REGISTER_FUNCTION(gt, evalf_func(gt_evalf))
223 
224 
229 static ex lt_evalf(const ex &x, const ex &y){
230  if (!is_a<numeric>(x) || !is_a<numeric>(y)) return 0;
231  ex xn = ex_to<numeric>(x);
232  ex yn = ex_to<numeric>(y);
233  int b = ((int) ( ex_to<numeric>(xn).to_double()
234  < ex_to<numeric>(yn).to_double() ) );
235  return b;
236 }
237 DECLARE_FUNCTION_2P(lt)
238 REGISTER_FUNCTION(lt, evalf_func(lt_evalf))
239 
240 
245 static ex ite_evalf(const ex &a, const ex &b, const ex &x, const ex &y){
246  if (!is_a<numeric>(x) || !is_a<numeric>(y)) return 0;
247  ex xn = ex_to<numeric>(x);
248  ex yn = ex_to<numeric>(y);
249  if ( equal_evalf(a,b)==1 ) return xn;
250  else return yn;
251 }
252 DECLARE_FUNCTION_4P(ite)
253 REGISTER_FUNCTION(ite, evalf_func(ite_evalf))
254 
255 
256 
261 static vector<double>* m_static_vector;
262 static ex Vector_evalf(const ex &i){
263  if (!is_a<numeric>(i) ) return 0;
264  unsigned int in = ((int) ex_to<numeric>(i).to_double() );
265  if ( (*m_static_vector).size() > in )
266  return (*m_static_vector).at(in);
267  else
268  return 0.0;
269 }
270 DECLARE_FUNCTION_1P(Vector)
271 REGISTER_FUNCTION(Vector, evalf_func(Vector_evalf))
272 
273 
274 
279 /* copy-paste of the GiNaC excompiler helper class (Removed all comments; see GiNaC sources for more info) */
280 class excompiler
281 {
282  struct filedesc
283  {
284  void* module;
285  std::string name;
286  bool clean_up;
287  };
288  std::vector<filedesc> filelist;
289 public:
290  ~excompiler()
291  {
292  for (std::vector<filedesc>::const_iterator it = filelist.begin(); it != filelist.end(); ++it) {
293  clean_up(it);
294  }
295  }
296  void add_opened_module(void* module, const std::string& name, bool clean_up)
297  {
298  filedesc fd;
299  fd.module = module;
300  fd.name = name;
301  fd.clean_up = clean_up;
302  filelist.push_back(fd);
303  }
304  void clean_up(const std::vector<filedesc>::const_iterator it)
305  {
306 #ifndef __MINGW32__
307  dlclose(it->module);
308 #endif
309  if (it->clean_up) {
310  remove(it->name.c_str());
311  }
312  }
313  void create_src_file(std::string& filename, std::ofstream& ofs)
314  {
315 #ifndef __MINGW32__
316  if (filename.empty()) {
317  const char* filename_pattern = "./GiNaCXXXXXX";
318  char* new_filename = new char[strlen(filename_pattern)+1];
319  strcpy(new_filename, filename_pattern);
320  #ifndef HAVE_MKSTEMPS
321  if (!mkstemp(new_filename)) {
322  delete[] new_filename;
323  throw std::runtime_error("mktemp failed");
324  }
325  #else
326  if (!MKTEMP(new_filename, 0)) {
327  delete[] new_filename;
328  throw std::runtime_error("mktemps failed");
329  }
330  #endif
331  filename = std::string(new_filename);
332  ofs.open(new_filename, std::ios::out);
333  delete[] new_filename;
334  } else {
335  ofs.open(filename.c_str(), std::ios::out);
336  }
337 
338  if (!ofs) {
339  throw std::runtime_error("could not create source code file for compilation");
340  }
341 
342  ofs << "#include <stddef.h> " << std::endl;
343  ofs << "#include <stdlib.h> " << std::endl;
344  ofs << "#include <math.h> " << std::endl;
345  ofs << std::endl;
346 #endif
347  }
348  void compile_src_file(const std::string filename, bool clean_up)
349  {
350  std::string strcompile = "ginac-excompiler " + filename;
351  if (system(strcompile.c_str())) {
352  throw std::runtime_error("excompiler::compile_src_file: error compiling source file!");
353  }
354  if (clean_up) {
355  remove(filename.c_str());
356  }
357  }
358  void* link_so_file(const std::string filename, bool clean_up)
359  {
360  void* module = NULL;
361 #ifndef __MINGW32__
362  module = dlopen(filename.c_str(), RTLD_NOW);
363 #endif
364  if (module == NULL) {
365  throw std::runtime_error("excompiler::link_so_file: could not open compiled module!");
366  }
367 
368  add_opened_module(module, filename, clean_up);
369 
370 #ifndef __MINGW32__
371  return dlsym(module, "compiled_ex");
372 #endif
373  }
374  void unlink(const std::string filename)
375  {
376  for (std::vector<filedesc>::iterator it = filelist.begin(); it != filelist.end();) {
377  if (it->name == filename) {
378  clean_up(it);
379  it = filelist.erase(it);
380  } else {
381  ++it;
382  }
383  }
384  }
385 };
386 
387 typedef double (*FUNCP_4P) (double, double, double, double);
388 static excompiler global_excompiler;
389 void compile_ex(const ex& expr, const symbol& sym1, const symbol& sym2, const symbol& sym3, const symbol& sym4,
390  FUNCP_4P& fp, const std::string filename = "") {
391 
392  symbol x("x"), y("y"), z("z"), g("g");
393  ex expr_with_xyzg = expr.subs(lst(sym1==x, sym2==y, sym3==z, sym4==g));
394 
395  std::ofstream ofs;
396  std::string unique_filename = filename;
397  global_excompiler.create_src_file(unique_filename, ofs);
398 
399  ofs << "double compiled_ex(double x, double y, double z, double g)" << std::endl;
400  ofs << "{" << std::endl;
401  ofs << "double res = ";
402  expr_with_xyzg.print(GiNaC::print_csrc_double(ofs));
403  ofs << ";" << std::endl;
404  ofs << "return(res); " << std::endl;
405  ofs << "}" << std::endl;
406 
407  ofs.close();
408 
409  global_excompiler.compile_src_file(unique_filename, filename.empty());
410  fp = (FUNCP_4P) global_excompiler.link_so_file(unique_filename+".so", filename.empty());
411 }
412 
413 
414 #endif
static ex mod_evalf(const ex &x, const ex &y)
Mod routine.
Definition: ginac_functions.h:183
static ex equal_evalf(const ex &x, const ex &y)
equal routine.
Definition: ginac_functions.h:197
const symbol & get_symbol(const string &sym_name)
Get a unique GiNaC symbol.
Definition: ginac_functions.h:68
static vector< double > * m_static_vector
Vector element.
Definition: ginac_functions.h:261
static ex lt_evalf(const ex &x, const ex &y)
less_than routine.
Definition: ginac_functions.h:229
static ex ite_evalf(const ex &a, const ex &b, const ex &x, const ex &y)
IfThenElse routine.
Definition: ginac_functions.h:245
static ex sinc_eval(const ex &x)
sqrt-function. (problem on mac: libc++abi.dylib throws exception: "sqrt" function is not recognized i...
Definition: ginac_functions.h:104
REGISTER_FUNCTION(sinc, eval_func(sinc_eval ).evalf_func(sinc_evalf ).derivative_func(sinc_deriv ).real_part_func(sinc_real_part).imag_part_func(sinc_imag_part).conjugate_func(sinc_conjugate)) static ex floor_evalf(const ex &x)
floor routine.
Definition: ginac_functions.h:157
static ex gt_evalf(const ex &x, const ex &y)
greater_than routine.
Definition: ginac_functions.h:213

-- last change 17.06.2016 | Tony Stoecker | Imprint --