API documentation
Attribute.cpp
Go to the documentation of this file.
00001 00008 /* 00009 * JEMRIS Copyright (C) 2007-2010 Tony Stöcker, Kaveh Vahedipour 00010 * Forschungszentrum Jülich, Germany 00011 * 00012 * This program is free software; you can redistribute it and/or modify 00013 * it under the terms of the GNU General Public License as published by 00014 * the Free Software Foundation; either version 2 of the License, or 00015 * (at your option) any later version. 00016 * 00017 * This program is distributed in the hope that it will be useful, 00018 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00019 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00020 * GNU General Public License for more details. 00021 * 00022 * You should have received a copy of the GNU General Public License 00023 * along with this program; if not, write to the Free Software 00024 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 00025 */ 00026 00027 00028 #include "Attribute.h" 00029 #include "Prototype.h" 00030 #include "AtomicSequence.h" 00031 #include "ginac_functions.h" 00032 00033 /***********************************************************/ 00034 Attribute::~Attribute (){ 00035 if (m_datatype == typeid( double*).name() ) delete (( double*) m_backup); 00036 if (m_datatype == typeid( int*).name() ) delete (( int*) m_backup); 00037 if (m_datatype == typeid( long*).name() ) delete (( long*) m_backup); 00038 if (m_datatype == typeid(unsigned*).name() ) delete ((unsigned*) m_backup); 00039 if (m_datatype == typeid( bool*).name() ) delete (( bool*) m_backup); 00040 if (m_datatype == typeid( string*).name() ) delete (( string*) m_backup); 00041 }; 00042 00043 /***********************************************************/ 00044 void Attribute::AttachObserver (Attribute* attrib){ 00045 00046 if ( !IsObservable() ) return; 00047 m_symbol_name = m_prototype->GetName()+"x"+m_name; 00048 for (unsigned int i=0; i<m_observers.size(); i++) if ( attrib == m_observers.at(i) ) return; 00049 m_observers.push_back(attrib); 00050 attrib->AttachSubject(this); 00051 00052 }; 00053 00054 /***********************************************************/ 00055 void Attribute::AttachSubject (Attribute* attrib){ 00056 00057 if ( !attrib->IsObservable() ) return; 00058 for (unsigned int i=0; i<m_subjects.size(); i++) if ( attrib == m_subjects.at(i) ) return; 00059 00060 m_subjects.push_back(attrib); 00061 attrib->AttachObserver(this); 00062 00063 }; 00064 00065 00066 /***********************************************************/ 00067 void Attribute::UpdatePrototype (Prototype* prot){ 00068 prot->Prepare(PREP_UPDATE); 00069 if (prot->GetType() == MOD_PULSE) ((AtomicSequence*) prot->GetParent())->CollectTPOIs(); 00070 }; 00071 00072 /***********************************************************/ 00073 bool Attribute::SetMember (string expr, const vector<Attribute*>& obs_attribs, bool verbose){ 00074 00075 //set my own symbol 00076 m_symbol_name = m_prototype->GetName()+"x"+m_name; 00077 00078 //attribute represents a string 00079 //if (GetTypeID()==typeid(string*).name()) { WriteMember(expr); return true; } 00080 00081 //attribute represents a PulseAxis 00082 if (GetTypeID()==typeid(PulseAxis*).name()) { 00083 if (expr=="RF") { WriteMember(AXIS_RF); return true; } 00084 if (expr=="GX") { WriteMember(AXIS_GX); return true; } 00085 if (expr=="GY") { WriteMember(AXIS_GY); return true; } 00086 if (expr=="GZ") { WriteMember(AXIS_GZ); return true; } 00087 else { WriteMember(AXIS_VOID); return true; } 00088 } 00089 00090 //GiNaC expressions 00091 Prototype::ReplaceString(expr,"step","csgn"); 00092 if (expr.find("I", 0)!=string::npos) m_complex = true; 00093 00094 m_subjects.clear(); 00095 m_symlist.remove_all(); 00096 //GiNaC::symbol d(m_sym_diff); 00097 //loop over all possibly observed subjects 00098 for (unsigned int i=0; i<obs_attribs.size() ; i++) { 00099 //convert string "a1","a2", ... to the matching symbol name 00100 Attribute* subject_attrib = obs_attribs.at(i); 00101 string SymbolName = subject_attrib->GetPrototype()->GetName() + "x" + subject_attrib->GetName(); 00102 stringstream key; key << "a" << i+1; 00103 if (!Prototype::ReplaceString(expr,key.str(),SymbolName)) continue; 00104 //still here? the attribute was in the expression, so it is an observed subject 00105 AttachSubject( subject_attrib ); 00106 m_symlist.append( get_symbol(SymbolName) ); 00107 } 00108 00109 //cout << "!!! " << GetPrototype()->GetName() << " : " << expr << " , " << m_symlist << endl; 00110 m_formula = expr; 00111 00112 //stop for strings 00113 if (GetTypeID()==typeid(string*).name()) { 00114 EvalExpression (); 00115 return true; 00116 } 00117 00118 //build GiNaC expression (maybe not successful at first call, if subjects still missing) 00119 try { 00120 m_expression = GiNaC::ex(m_formula,m_symlist); 00121 //differentiation of expression? 00122 if (m_diff>0) { 00123 m_expression = m_expression.diff(get_symbol(m_sym_diff),m_diff); 00124 stringstream se; se << m_expression; 00125 m_formula = se.str(); 00126 } 00127 //test the expression evaluation once 00128 EvalExpression (); 00129 00130 } catch (exception &p) { 00131 00132 if ( verbose ) { 00133 00134 cout << "Warning in " << m_prototype->GetName() << ": attribute " << GetName() 00135 << " can not evaluate its GiNaC expression E = " << expr 00136 << ". Reason: " << p.what() << endl; 00137 00138 } 00139 return false; 00140 } 00141 return true; 00142 }; 00143 00144 /***********************************************************/ 00145 void Attribute::EvalExpression () { 00146 00147 if (m_formula.empty()) return; 00148 00149 //strings: simply replace attribute symbol by its value 00150 if (GetTypeID()==typeid(string*).name()) { 00151 string expr = m_formula; 00152 for (unsigned int i=0; i<m_subjects.size() ; i++) { 00153 Attribute* a = m_subjects.at(i); 00154 stringstream key; 00155 if (a->GetTypeID()==typeid( double*).name()) key << a->GetMember <double>() ; 00156 if (a->GetTypeID()==typeid( int*).name()) key << a->GetMember <int>() ; 00157 if (a->GetTypeID()==typeid( long*).name()) key << a->GetMember <long>() ; 00158 if (a->GetTypeID()==typeid(unsigned*).name()) key << a->GetMember<unsigned>() ; 00159 if (a->GetTypeID()==typeid( bool*).name()) key << a->GetMember <bool>() ; 00160 string SymbolName = a->GetPrototype()->GetName() + "x" + a->GetName(); 00161 Prototype::ReplaceString(expr,SymbolName,key.str()); 00162 } 00163 WriteMember(expr); 00164 //if (GetName()=="Filename") cout << " Eval: " << expr << endl; 00165 return; 00166 } 00167 00168 //collect symbols and corresponding member-values from observed attributes 00169 GiNaC::lst numlist; 00170 for (unsigned int i=0; i<m_subjects.size() ; i++) { 00171 Attribute* a = m_subjects.at(i); 00172 if (a->GetTypeID()==typeid( double*).name()) { numlist.append(a->GetMember <double>() ); continue; } 00173 if (a->GetTypeID()==typeid( int*).name()) { numlist.append(a->GetMember <int>() ); continue; } 00174 if (a->GetTypeID()==typeid( long*).name()) { numlist.append(a->GetMember <long>() ); continue; } 00175 if (a->GetTypeID()==typeid(unsigned*).name()) { numlist.append(a->GetMember<unsigned>() ); continue; } 00176 if (a->GetTypeID()==typeid( bool*).name()) { numlist.append(a->GetMember <bool>() ); continue; } 00177 } 00178 00179 00180 //numeric evaluation of GiNaC expression 00181 GiNaC::ex e = m_expression.subs(m_symlist,numlist); 00182 m_static_vector = m_prototype->GetVector(); // static pointer to evaluate the Vector function 00183 e = GiNaC::evalf(e); 00184 double d = 0.0; 00185 if (GiNaC::is_a<GiNaC::numeric>(e)) { 00186 if ( m_complex ) m_imaginary = -0.5 * GiNaC::ex_to<GiNaC::numeric>( (e-e.conjugate())*GiNaC::I ).to_double(); 00187 d = GiNaC::ex_to<GiNaC::numeric>( e ).to_double();//default is real-part 00188 } 00189 00190 00191 //overwrite private member 00192 if (m_datatype==typeid( double*).name() ) WriteMember((double) d ); 00193 if (m_datatype==typeid( int*).name() ) WriteMember((int) d ); 00194 if (m_datatype==typeid( long*).name() ) WriteMember((long) d ); 00195 if (m_datatype==typeid(unsigned*).name() ) WriteMember((unsigned) d ); 00196 if (m_datatype==typeid( bool*).name() ) WriteMember((bool) d ); 00197 00198 // for dynamic change of runtime compiled attributes after notification! 00199 if (m_ginac_excomp && m_cur_fp==m_num_fp && m_num_fp>0) { 00200 //cout << "!!! " << GetName() << " : NFP =" << GetNumberFunctionPointers() << endl; 00201 //cout << "!!! " << GetName() << " : CFP =" << GetCurrentFunctionPointer() << endl; 00202 //cout << "!!! " << GetName() << " : size=" << m_compiled.size() << endl; 00203 00204 if (m_compiled.size() == m_cur_fp ) m_compiled.push_back(false); 00205 } 00206 }; 00207 00208 /***********************************************************/ 00209 double Attribute::EvalCompiledExpression (double const val, string const attrib ) { 00210 00211 //cout << GetPrototype()->GetName() << " ?? at pointer num " << m_cur_fp << " -> compiled = " << m_compiled.at(m_cur_fp) << endl; 00212 if (!m_compiled.at(m_cur_fp)) { 00213 //substitute all attributes with numbers in GiNaC expression, except the attribute 00214 //which serves as the free parameter for runtime compilation 00215 GiNaC::lst symlist; 00216 GiNaC::lst numlist; 00217 for (unsigned int i=0; i<m_subjects.size() ; i++) { 00218 Attribute* a = m_subjects.at(i); 00219 if (a->GetName() == attrib) continue; 00220 symlist.append( get_symbol(a->GetSymbol()) ); 00221 if (a->GetTypeID()==typeid( double*).name()) { numlist.append(a->GetMember <double>() ); continue; } 00222 if (a->GetTypeID()==typeid( int*).name()) { numlist.append(a->GetMember <int>() ); continue; } 00223 if (a->GetTypeID()==typeid( long*).name()) { numlist.append(a->GetMember <long>() ); continue; } 00224 if (a->GetTypeID()==typeid(unsigned*).name()) { numlist.append(a->GetMember<unsigned>() ); continue; } 00225 if (a->GetTypeID()==typeid( bool*).name()) { numlist.append(a->GetMember <bool>() ); continue; } 00226 } 00227 00228 GiNaC::ex e = GiNaC::evalf((symlist.nops()==0)?m_expression:m_expression.subs(symlist,numlist)); 00229 00230 //add function pointers 00231 m_fp.push_back(NULL); 00232 m_fpi.push_back(NULL); 00233 //compile the GiNaC expression 00234 try { 00235 //fairly easy for real valued expressions 00236 if (!m_complex) { 00237 compile_ex(e, get_symbol(GetPrototype()->GetAttribute(attrib)->GetSymbol()), m_fp.at(m_num_fp)); 00238 } 00239 //more work to do, since GiNaC::realsymbol does not behave as expected (and it is therefore not used at all) 00240 else { 00241 stringstream se; se << e; string formula = se.str(); 00242 string sym = GetPrototype()->GetAttribute(attrib)->GetSymbol(); 00243 string asym = "abs(VarForEvalCompiledExpression)"; 00244 Prototype::ReplaceString(formula,sym,asym); 00245 GiNaC::lst symlist; 00246 symlist.append( get_symbol("VarForEvalCompiledExpression") ); 00247 GiNaC::ex ea = GiNaC::ex(formula,symlist); 00248 symlist.remove_all(); 00249 symlist.append( get_symbol(sym) ); 00250 00251 GiNaC::ex ear = ea.real_part(); 00252 stringstream ser; ser << ear; formula = ser.str(); 00253 if ( Prototype::ReplaceString(formula,asym,sym) ) { 00254 ear = GiNaC::ex(formula,symlist); 00255 compile_ex(ear, get_symbol(GetPrototype()->GetAttribute(attrib)->GetSymbol()), m_fp.at(m_num_fp)); 00256 } 00257 00258 GiNaC::ex eai = ea.imag_part(); 00259 stringstream sei; sei << eai; formula = sei.str(); 00260 if ( Prototype::ReplaceString(formula,asym,sym) ) { 00261 eai = GiNaC::ex(formula,symlist); 00262 compile_ex(eai, get_symbol(GetPrototype()->GetAttribute(attrib)->GetSymbol()), m_fpi.at(m_num_fp)); 00263 } 00264 } 00265 //cout << " compiling expression " << e << " of attribute " << GetName() << " in module " << GetPrototype()->GetName() << endl; 00266 m_num_fp++; 00267 } 00268 catch (exception &p) { 00269 cout << " Warning: attribute " << GetName() << " of module " << GetPrototype()->GetName() << endl << endl 00270 << " function Attribute::EvalCompiledExpression" << endl 00271 << " No external runtime compiler available: " << p.what() << endl 00272 << " Falling back to (slow) analytic evaluation!" << endl << endl 00273 << " Hint: if you have a shell and gcc on your system, create the one-liner " << endl << endl 00274 << " #!/bin/sh" << endl 00275 << " gcc -x c -fPIC -shared -o $1.so $1" << endl << endl 00276 << " name it \"ginac-excompiler\", and put it somewhere in your search path." << endl << endl; 00277 m_ginac_excomp = false; 00278 } 00279 m_compiled.at(m_cur_fp) = true; //even if compilation failed, as we don't have to try a second time! 00280 } 00281 00282 //if compilation failed, invoke slow analytic evaluation 00283 if (!m_ginac_excomp ) { 00284 *((double*) GetPrototype()->GetAttribute(attrib)-> GetAddress()) = val; 00285 EvalExpression(); 00286 return *((double*) GetAddress()); 00287 } 00288 00289 //invoke fast runtime compiled routines 00290 if (m_fpi.at(m_cur_fp) != NULL ) m_imaginary = m_fpi.at(m_cur_fp)(val); 00291 if ( m_fp.at(m_cur_fp) != NULL ) return m_fp.at(m_cur_fp)(val); 00292 00293 return 0.0; 00294 00295 }; 00296 00297 /***********************************************************/ 00298 double Attribute::EvalCompiledNLGExpression (double const x, double const y ,double const z, double const g ) { 00299 00300 //cout << GetPrototype()->GetName() << " ?? at pointer num " << m_cur_fp << " -> compiled = " << m_compiled.at(m_cur_fp) << endl; 00301 if (m_nlgfp == NULL && m_ginac_excomp) { 00302 //substitute all attributes with numbers in GiNaC expression, except the attribute 00303 //which serves as the free parameter for runtime compilation 00304 GiNaC::lst symlist; 00305 GiNaC::lst numlist; 00306 for (unsigned int i=0; i<m_subjects.size() ; i++) { 00307 Attribute* a = m_subjects.at(i); 00308 if (a->GetName() == "NLG_posX") continue; 00309 if (a->GetName() == "NLG_posY") continue; 00310 if (a->GetName() == "NLG_posZ") continue; 00311 if (a->GetName() == "NLG_value") continue; 00312 symlist.append( get_symbol(a->GetSymbol()) ); 00313 if (a->GetTypeID()==typeid( double*).name()) { numlist.append(a->GetMember <double>() ); continue; } 00314 if (a->GetTypeID()==typeid( int*).name()) { numlist.append(a->GetMember <int>() ); continue; } 00315 if (a->GetTypeID()==typeid( long*).name()) { numlist.append(a->GetMember <long>() ); continue; } 00316 if (a->GetTypeID()==typeid(unsigned*).name()) { numlist.append(a->GetMember<unsigned>() ); continue; } 00317 if (a->GetTypeID()==typeid( bool*).name()) { numlist.append(a->GetMember <bool>() ); continue; } 00318 } 00319 00320 GiNaC::ex e = GiNaC::evalf((symlist.nops()==0)?m_expression:m_expression.subs(symlist,numlist)); 00321 00322 try { 00323 compile_ex (e, 00324 get_symbol(GetPrototype()->GetAttribute("NLG_posX")->GetSymbol()), 00325 get_symbol(GetPrototype()->GetAttribute("NLG_posY")->GetSymbol()), 00326 get_symbol(GetPrototype()->GetAttribute("NLG_posZ")->GetSymbol()), 00327 get_symbol(GetPrototype()->GetAttribute("NLG_value")->GetSymbol()), 00328 m_nlgfp); 00329 //cout << " compiling attribute " << GetName() << " of module " << GetPrototype()->GetName() << endl; 00330 } 00331 catch (exception &p) { 00332 cout << " Warning: attribute " << GetName() << " of module " << GetPrototype()->GetName() << endl << endl 00333 << " function Attribute::EvalCompiledNLGExpression" << endl 00334 << " No external runtime compiler available: " << p.what() << endl 00335 << " Falling back to (slow) analytic evaluation!" << endl << endl 00336 << " Hint: if you have a shell and gcc on your system, create the one-liner " << endl << endl 00337 << " #!/bin/sh" << endl 00338 << " gcc -x c -fPIC -shared -o $1.so $1" << endl << endl 00339 << " name it \"ginac-excompiler\", and put it somewhere in your search path." << endl << endl; 00340 m_ginac_excomp = false; 00341 } 00342 } 00343 00344 //if compilation failed, invoke slow analytic evaluation 00345 if (!m_ginac_excomp ) { 00346 *((double*) GetPrototype()->GetAttribute("NLG_posX")-> GetAddress()) = x; 00347 *((double*) GetPrototype()->GetAttribute("NLG_posY")-> GetAddress()) = y; 00348 *((double*) GetPrototype()->GetAttribute("NLG_posZ")-> GetAddress()) = z; 00349 *((double*) GetPrototype()->GetAttribute("NLG_value")-> GetAddress()) = g; 00350 EvalExpression(); 00351 return *((double*) GetAddress()); 00352 } 00353 00354 //invoke fast runtime compiled routines 00355 if ( m_nlgfp != NULL ) return m_nlgfp(x,y,z,g); 00356 00357 return 0.0; 00358 00359 }; 00360 00361
