00001 /*************************************************************************** 00002 ssa.h 00003 ----------- 00004 00005 begin : Tue Apr 24 17:56:51 CEST 2001 00006 author : (C) 2001 by Michael Peeters 00007 email : Michael.Peeters@vub.ac.be 00008 ***************************************************************************/ 00009 00010 /*************************************************************************** 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 ***************************************************************************/ 00018 00019 #ifndef SSA_H 00020 #define SSA_H 00021 00022 #include "utility.h" 00023 #include <stdexcept> 00024 #include "numerictraits.h" 00025 #include "numerictypes.h" 00026 #include "jacobian.h" 00027 #include "rootscan.h" 00028 00029 namespace MODEL { 00030 00035 template <integer dims, typename nelem=number, class NT = NumericTraits<nelem,dims> > 00036 class SSA 00037 { 00038 public: 00039 typedef typename NT::number numT; 00040 typedef typename NT::vect vect; 00041 typedef typename NT::matrix matrix; 00042 typedef typename Jacobian<2*dims>::matrix system; 00043 typedef typename Jacobian<2*dims>::vect input; 00044 typedef typename NumericTraits<complex,dims>::vect response; 00045 typedef response cinput; 00046 typedef typename NT::vf vf; 00047 00048 public: 00050 SSA(vf& my_sys, const string& parm); 00051 00056 counter set_stat_param(const number& value); 00057 00059 void set_stat_point(const counter& po=1); 00060 00063 void set_stat_point(const vect& right_there); 00064 00066 const ScanList<dims,NT>& get_stat_points(void); 00067 00069 response calc_point(const number& omega); 00070 00072 ScanList<dims,complex> calc_response(const number& from, const 00073 number& to, const counter& 00074 n); 00075 00078 ScanList<dims> calc_response_norm(const number& from, const 00079 number& to, const counter& 00080 n); 00081 00085 void set_dep(const vect& input) 00086 { 00087 linpar = 0.; 00088 for (integer i=0;i<dims;i++) 00089 linpar[i]=input[i]; 00090 } 00091 00095 void set_dep(const cinput& input) 00096 { 00097 linpar = 0.; 00098 for (integer i=0;i<dims;i++) 00099 { 00100 linpar[i]=input[i].real(); 00101 linpar[dims+i]=input[i].imag(); 00102 } 00103 } 00104 00105 private: 00107 void calc_dep(void); 00108 00110 void put_omega(const number& omega); 00111 00112 private: 00115 bool all_is_well; 00116 00118 Jacobian<dims>* J; 00119 00121 matrix j; // small j: nice and confusing 00122 00124 system A; 00125 00127 RootScan<dims,nelem>* rs; 00128 00130 ScanList<dims,nelem> statp; 00131 00133 vect here; 00134 00136 ParameterP p; 00137 00140 input linpar; 00141 00142 private: 00143 // to avoid silly errors 00144 NO_COPY(SSA); 00145 00146 }; 00147 00148 // ------------------------------------------------------------ 00149 00150 template <integer dims, typename nelem, class NT > 00151 SSA<dims,nelem,NT>::SSA(vf& my_sys, const string& parm) 00152 { 00153 all_is_well=false; 00154 00155 // get the jacobian, copy the function as is 00156 J=new Jacobian<dims>(my_sys,1E-6,true); 00157 00158 // get the parameter 00159 p=&(J->get_function().get_parameter(parm)); 00160 00161 // init the scanner 00162 rs=new RootScan<dims,nelem>(J->get_function()); 00163 } 00164 00165 template <integer dims, typename nelem, class NT > 00166 counter SSA<dims,nelem,NT>::set_stat_param(const number& value) 00167 { 00168 *p=value; 00169 00170 statp=rs->scan(0.,1.,1); // these are dummy values 00171 statp=statp.select(RootScan<dims,nelem>::stable); 00172 00173 // statp.print_raw(cout); 00174 00175 // cerr << "SSA: started number of solutions " << statp.get_data()[0.].size() << endl; 00176 00177 return statp.get_data()[0.].size(); 00178 // watch out here: you might find only one (wrong point), because 00179 // we are NOT scanning. This is a problem in RootScan. 00180 } 00181 00182 template <integer dims, typename nelem, class NT > 00183 void SSA<dims,nelem,NT>::set_stat_point(const counter& po) 00184 { 00185 // go to correct point 00186 if(statp.get_data()[0.].size()>0) 00187 { 00188 // take the first parameter value in the list 00189 ScanList<dims,nelem>::parpointlist::iterator 00190 getit=statp.get_data()[0].begin(); 00191 00192 // go to the right value 00193 for(counter i=po;i>1;i--) getit++; 00194 00195 // the first value is the point... 00196 here=(*getit)[0]; 00197 00198 // now that we know where we, calculate the approximations 00199 calc_dep(); 00200 00201 // now we can do something 00202 all_is_well=true; 00203 } 00204 00205 } 00206 template <integer dims, typename nelem, class NT > 00207 void SSA<dims,nelem,NT>::set_stat_point(const vect& right_there) 00208 { 00209 here=right_there; 00210 00211 // now that we know where we, calculate the approximations 00212 calc_dep(); 00213 00214 // now we can do something 00215 all_is_well=true; 00216 } 00217 00218 template <integer dims, typename nelem, class NT > 00219 SSA<dims,nelem,NT>::response SSA<dims,nelem, NT>::calc_point(const number& omega) 00220 { 00221 if(all_is_well) { 00222 00223 put_omega(omega); 00224 LUSolve<2*dims> solA(A,true); // make a copy 00225 00226 // copy into correct form 00227 input t1 (solA(linpar)); 00228 response t2; 00229 for (integer i=0;i<dims;i++) 00230 { 00231 t2[i]=t1[i]+I*t1[i+dims]; 00232 } 00233 00234 return t2; 00235 } 00236 00237 } 00238 00239 template <integer dims, typename nelem, class NT > 00240 ScanList<dims,complex> 00241 SSA<dims,nelem, NT>::calc_response(const number& from, const 00242 number& to, const counter& 00243 n) 00244 { 00245 if(all_is_well) { 00246 00247 ScanList<dims,complex> res; 00248 00249 number scaler=pow(to/from,1./number(n)); 00250 for (number omega=from;omega<to;omega*=scaler) 00251 { 00252 res.add_point(calc_point(omega)); 00253 res.add_param(omega); 00254 } 00255 return res; 00256 } 00257 00258 } 00259 00260 template <integer dims, typename nelem, class NT > 00261 ScanList<dims> 00262 SSA<dims,nelem, NT>::calc_response_norm(const number& from, const 00263 number& to, const counter& 00264 n) 00265 { 00266 if(all_is_well) { 00267 00268 ScanList<dims> res; 00269 00270 number scaler=pow(to/from,1./number(n)); 00271 for (number omega=from;omega<to;omega*=scaler) 00272 { 00273 response t1=calc_point(omega); 00274 vect t2; 00275 for(integer i=0;i<dims;i++) 00276 t2[i]=norm(t1[i]); 00277 00278 res.add_point(t2); 00279 res.add_param(omega); 00280 } 00281 return res; 00282 } 00283 00284 } 00285 00286 template <integer dims, typename nelem, class NT > 00287 void SSA<dims,nelem, NT>::calc_dep(void) 00288 { 00289 // generate the Jacobian accurately (slow) 00290 j=J->calculate_accurate(here); 00291 00292 // And copy it to the larger system 00293 for (int i=0;i<dims;i++) 00294 { 00295 for (int k=0;k<dims;k++) 00296 { 00297 A[i][k] = -j[i][k]; 00298 A[dims+i][dims+k] =-j[i][k]; 00299 } 00300 } 00301 00302 // find the (linear) dependecies of the equations to the 00303 // (tiny/small/microscopic input). This is done in a way very 00304 // similar to the Jacobian. 00305 // it is store in linpar 00306 linpar=0.; // just to be safe 00307 00308 // value now 00309 vf& f=J->get_function(); 00310 vect fu ( f(here) ); 00311 00312 number pdp=*p; 00313 00315 const number epsilon = 1; 00316 number dp=epsilon*abs(pdp); 00317 00318 if(dp==0.0) dp = epsilon; // avoid numerical error 00319 pdp += dp; // " 00320 dp = pdp - *p; // " 00321 00322 // get value at new point 00323 number oldp=*p; 00324 *p=pdp; 00325 00326 Parameter herep = f.get_parameter("current"); 00327 00328 vect fudp(f(here)); 00329 00330 *p=oldp; // restore parameter 00331 00332 for(integer i=0;i<dims;i++) 00333 { 00334 linpar[i]=(fudp[i]-fu[i])/dp; 00335 } 00336 } 00337 00338 template <integer dims, typename nelem, class NT > 00339 void SSA<dims,nelem, NT>::put_omega(const number& omega) 00340 { 00341 for (int i=0;i<dims;i++) 00342 { 00343 A[i][dims+i] = -omega; 00344 A[dims+i][i] =omega; 00345 } 00346 } 00347 00348 } 00349 // end MODEL 00350 #endif 00351 00352 // CVS Test: hello Guy. 00353 00354 /********************************************************************* 00355 $Id: ssa.h,v 1.1.2.2 2001/08/30 14:48:27 mpeeters Exp $ 00356 ********************************************************************** 00357 00358 $Log: ssa.h,v $ 00359 Revision 1.1.2.2 2001/08/30 14:48:27 mpeeters 00360 Testing CVS Branching. 00361 00362 Revision 1.1.2.1 2001/08/30 14:01:22 mpeeters 00363 Fixed stupid warning about extra token (due to #endif SSHA_H). 00364 00365 Revision 1.1 2001/05/22 10:54:55 mpeeters 00366 Moved sources and headers for libModel to model/ subdirectory, in an attempt to rationalize the source tree. This should make things "netter". 00367 00368 Revision 1.1 2001/05/21 11:56:43 mpeeters 00369 Commited lots of small changes (too many to tell) at the previous commit (accidentally). In this commit, we add the small signal analysis stuff. Whoohoo. 00370 00371 *********************************************************************/
More Info? Michael Peeters. Also, check our research website: www.alna.vub.ac.be
Last update: June 2002.
Looking for Open Source? Check out
!