开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 数值算法/人工智能 > 查看源码
VFUNC.CPP
资源名称：BCMATH.rar　[点击查看]
上传用户：zengbais
上传日期：2022-08-08
资源大小：49k
文件大小：14k
源码类别：

数值算法/人工智能

开发平台：

C++ Builder

VFUNC.CPP：源码内容
  1. #include <values.h>
  2. #include <math.h>
  3. #include <stdio.h>
  4. #include "vfunc.h"
  7. matrix valgo::cal(matrix& x) { // 基类的基本算法
  8. matrix xx(x);
  9. if(xshift != 0.0) xx -= xshift;
  10. if(xfactor != 1.0) xx *=xfactor;
  11. calculate(xx);
  12. if(yfactor != 1.0) xx *= yfactor;
  13. if(addconst != 0.0) xx += addconst;
  14. return xx;
  15. }
  17. valgo * valgo::clone() // 克隆自己，必须被继承子类改写
  18. {
  19. return new valgo(*this);
  20. }
  22. valgo * valgo::mul(DOUBLE a) // 乘a
  23. {
  24. valgo * alg;
  25. alg = this;
  26. if(refnum>1) { refnum--;
  27. alg = clone(); // 如引用数大于1，则产生新的对象
  28. }
  29. alg->yfactor *= a;
  30. alg->addconst *= a;
  31. return alg;
  32. }
  34. valgo * valgo::add(DOUBLE a) // 加a
  35. {
  36. valgo * alg;
  37. alg = this;
  38. if(refnum>1) { refnum--;
  39. alg = clone(); // 如引用数大于1，则产生新的对象
  40. }
  41. alg->addconst += a;
  42. return alg;
  43. }
  45. valgo * valgo::neg() // 取负
  46. {
  47. valgo * alg;
  48. alg = this;
  49. if(refnum>1) { refnum--;
  50. alg = clone(); // 如引用数大于1，则产生新的对象
  51. }
  52. alg->addconst = -alg->addconst;
  53. alg->yfactor = -alg->yfactor;
  54. return alg;
  55. }
  57. valgo * valgo::setxfactor(DOUBLE x) // 设置x轴因子
  58. {
  59. valgo * alg;
  60. alg = this;
  61. if(refnum>1) { refnum--;
  62. alg = clone(); // 如引用数大于1，则产生新的对象
  63. }
  64. alg->xfactor = x;
  65. return alg;
  66. }
  68. valgo * valgo::xroom(DOUBLE x) // 将xfactor扩大x倍
  69. {
  70. valgo * alg;
  71. alg = setxfactor(x*xfactor);
  72. if(x!=0.0)
  73. alg->setxshift(alg->xshift/x);
  74. return alg;
  75. }
  77. valgo * valgo::setxshift(DOUBLE x) // 设置xshift的值
  78. {
  79. valgo * alg;
  80. alg = this;
  81. if(refnum>1) { refnum--;
  82. alg = clone(); // 如引用数大于1，则产生新的对象
  83. }
  84. alg->xshift = x;
  85. return alg;
  86. }
  88. valgo * valgo::xshiftas(DOUBLE x) // 从当前开始平移x
  89. {
  90. return setxshift(xshift+x);
  91. }
  93. valgo * valgojoin::clone() // 克隆自己
  94. {
  95. return new valgojoin(*this);
  96. }
  98. matrix& valgojoin::calculate(matrix& x) // 实施结合算法
  99. {
  100. if(leftalgo==0 || rightalgo==0)
  101. throw TMESSAGE("empty valgo pointor!");
  102. if(met == ccom) { // 复合函数
  103. matrix x1;
  104. x1 = rightalgo->cal(x);
  105. x = leftalgo->cal(x1);
  106. return x;
  107. }
  108. matrix xr;
  109. xr = rightalgo->cal(x);
  110. x = leftalgo->cal(x);
  111. if(met == cadd) // 返回各结合运算值
  112. x += xr;
  113. else if(met == csub)
  114. x -= xr;
  115. else if(met == cmul)
  116. x *= xr;
  117. else if(met == cdiv)
  118. x /= xr;
  119. return x;
  120. }
  122. valgo * valgofun::clone() // 克隆自己
  123. {
  124. return new valgofun(*this);
  125. }
  127. matrix& valgofun::calculate(matrix& x) // 实施函数算法
  128. {
  129. if(f)
  130. return x=f(x);
  131. return (x*=0.0);
  132. }
  134. valgo * valgofun1::clone() // 克隆自己
  135. {
  136. return new valgofun1(*this);
  137. }
  139. matrix& valgofun1::calculate(matrix& x) // 实施函数算法
  140. {
  141. DOUBLE xx;
  142. xx = x(0);
  143. if(f)
  144. return (x=f(xx));
  145. return (x*=0.0);
  146. }
  148. matrix& valgodiff::calculate(matrix& x) // 实施函数算法
  149. {
  150. DOUBLE t;
  151. t = x(0);
  152. return calcul(t);
  153. }
  155. valgo * valgodiff::clone() // 克隆自己
  156. {
  157. return new valgodiff(*this);
  158. }
  160. matrix& valgodiff::calcul(DOUBLE t, DOUBLE eps) // 标量算法
  161. {
  162. DOUBLE h,t00;
  163. matrix y00(y0);
  164. t00 = t0;
  165. h = t-t0;
  166. if(t0==t)
  167. return y0;
  168. if(tnow==t)
  169. return ynow;
  170. if(t>tnow) {
  171. t00 = tnow;
  172. h = t - tnow;
  173. y00 = ynow;
  174. }
  175. matrix y(y00);
  176. size_t n = y.rownum;
  177. size_t m,i,j,k;
  178. DOUBLE hh,p,dt,x,tt,q,a[4];
  179. matrix g(n),b(n),c(n),d(n),e(n);
  180. hh=h; m=1; p=1.0+eps; x=t00;
  181. c = y;
  182. while (p>=eps)
  183. { a[0]=hh/2.0; a[1]=a[0]; a[2]=hh; a[3]=hh;
  184.   g = y; y = c;
  185.   dt=h/m; t00=x;
  186.   for (j=0; j<m; j++)
  187.  { d = f(t00,y);  // grkt2f(t,y,n,d);
  188. b = y; e = y;
  189. for (k=0; k<=2; k++)
  190.   {
  191. y = e+(d*a[k]);
  192. b += d*(a[k+1]/3.0);
  193. tt=t00+a[k];
  194. d = f(tt,y);
  195.   }
  196. y = b + d*(hh/6.0);
  197. t00+=dt;
  198.  }
  199.   p=0.0;
  200.   for (i=0; i<n; i++)
  201.  { q=fabs(y(i)-g(i));
  202. if (q>p) p=q;
  203.  }
  204.   hh=hh/2.0; m=m+m;
  205. } // end of while(p>eps)
  206. tnow = t;
  207. ynow = y;
  208. return ynow;
  209. }
  211. vfunc::vfunc()// 缺省构造函数，产生函数f(x)=x;
  212. {
  213. alg = new valgo();
  214. }
  216. vfunc::vfunc(vfunc & fn) // 拷贝构造函数
  217. {
  218. alg = fn.alg;
  219. alg->refnum++;
  220. }
  222. vfunc::vfunc(matrix& (*fun)(matrix&)) // 函数指针的构造函数
  223. {
  224. alg = new valgofun(fun);
  225. }
  227. vfunc::vfunc(matrix (*fun)(DOUBLE)) // 标量到向量的函数指针构造函数
  228. {
  229. alg = new valgofun1(fun);
  230. }
  232. vfunc::vfunc(DOUBLE a) // 常函数构造函数
  233. {
  234. alg = new valgo(a);
  235. }
  237. vfunc& vfunc::operator=(vfunc& fn) // 赋值运算符
  238. {
  239. if(this == &fn) return (*this); // 如果等于自己，干脆啥也别做
  240. if(alg) {
  241. alg->refnum--; // 原算法引用数减1
  242. if(!alg->refnum) // 如结果为0则删除原算法
  243. delete alg;
  244. }
  245. alg = fn.alg; // 将fn的算法移过来
  246. if(alg) alg->refnum++; // 引用数增加
  247. return (*this);
  248. }
  250. vfunc& vfunc::operator=(matrix& (*fn)(matrix&)) // 用函数指针的赋值运算符
  251. {
  252. if(alg) {
  253. alg->refnum--; // 原算法引用数减1
  254. if(!alg->refnum) // 如结果为0则删除原算法
  255. delete alg;
  256. }
  257. alg = new valgofun(fn);
  258. return (*this);
  259. }
  261. vfunc& vfunc::operator=(DOUBLE a) // 常函数的赋值运算符
  262. {
  263. if(alg) {
  264. alg->refnum--; // 原算法引用数减1
  265. if(!alg->refnum) // 如结果为0则删除原算法
  266. delete alg;
  267. }
  268. alg = new valgo(a);
  269. return (*this);
  270. }
  272. vfunc& vfunc::operator+=(vfunc& fn) // 自身加一个函数
  273. {
  274. valgo * a = new valgojoin(alg, fn.alg, cadd);
  275. alg->refnum--; // 因为联合算法对两个算法都加了引用数，因此再减回来
  276. alg = a; // 将新指针赋给alg
  277. return (*this);
  278. }
  280. vfunc& vfunc::operator+=(matrix& (*f)(matrix&)) // 自身加一个函数指针
  281. {
  282. vfunc fn(f); // 将函数指针包装为函数类
  283. operator+=(fn); // 实施加操作
  284. return (*this);
  285. }
  287. vfunc vfunc::operator+(vfunc& fn) // 相加产生新函数
  288. {
  289. valgo * a = new valgojoin(alg, fn.alg, cadd);
  290. vfunc f(a);
  291. return f;
  292. }
  294. vfunc vfunc::operator+(DOUBLE a) // 与常数相加产生新函数
  295. {
  296. vfunc f(*this);
  297. f += a;
  298. return f;
  299. }
  301. vfunc vfunc::operator+(matrix& (*f)(matrix&)) // 加一个函数指针产生新函数
  302. {
  303. vfunc ff(*this);
  304. ff += f;
  305. return ff;
  306. }
  308. vfunc& vfunc::neg() // 自身取负
  309. {
  310. alg=alg->neg();
  311. return (*this);
  312. }
  314. vfunc vfunc::operator-() // 产生负函数
  315. {
  316. vfunc f(*this);
  317. f.neg();
  318. return f;
  319. }
  321. vfunc& vfunc::operator-=(vfunc& fn) // 自身减一个函数
  322. {
  323. valgo * a = new valgojoin(alg, fn.alg, csub);
  324. alg->refnum--; // 因为联合算法对两个算法都加了引用数，因此再减回来
  325. alg = a; // 将新指针赋给alg
  326. return (*this);
  327. }
  329. vfunc& vfunc::operator-=(matrix& (*f)(matrix&)) // 自身减一个函数指针
  330. {
  331. vfunc fn(f); // 将函数指针包装为函数类
  332. operator-=(fn); // 实施减操作
  333. return (*this);
  334. }
  336. vfunc vfunc::operator-(vfunc& fn) // 相减产生新函数
  337. {
  338. valgo * a = new valgojoin(alg, fn.alg, csub);
  339. vfunc f(a);
  340. return f;
  341. }
  343. vfunc vfunc::operator-(DOUBLE a) // 与常数相减产生新函数
  344. {
  345. vfunc f(*this);
  346. f -= a;
  347. return f;
  348. }
  350. vfunc operator-(DOUBLE a, vfunc& f) // 常数减函数
  351. {
  352. return (-f)+a;
  353. }
  355. vfunc vfunc::operator-(matrix& (*f)(matrix&)) // 减一个函数指针产生新函数
  356. {
  357. vfunc ff(*this);
  358. ff -= f;
  359. return ff;
  360. }
  362. vfunc operator-(matrix& (*f)(matrix&),vfunc& fn) // 函数指针减函数
  363. {
  364. vfunc ff(f);
  365. ff -= fn;
  366. return ff;
  367. }
  369. vfunc& vfunc::operator*=(vfunc& fn) // 自身乘一个函数
  370. {
  371. valgo * a = new valgojoin(alg, fn.alg, cmul);
  372. alg->refnum--; // 因为联合算法对两个算法都加了引用数，因此再减回来
  373. alg = a; // 将新指针赋给alg
  374. return (*this);
  375. }
  377. vfunc& vfunc::operator*=(matrix& (*f)(matrix&)) // 自身乘一个函数指针
  378. {
  379. vfunc fn(f); // 将函数指针包装为函数类
  380. operator*=(fn); // 实施乘操作
  381. return (*this);
  382. }
  384. vfunc vfunc::operator*(vfunc& fn) // 相乘产生新函数
  385. {
  386. valgo * a = new valgojoin(alg, fn.alg, cmul);
  387. vfunc f(a);
  388. return f;
  389. }
  391. vfunc vfunc::operator*(DOUBLE a) // 与常数相乘产生新函数
  392. {
  393. vfunc f(*this);
  394. f *= a;
  395. return f;
  396. }
  398. vfunc vfunc::operator*(matrix& (*f)(matrix&)) // 乘一个函数指针产生新函数
  399. {
  400. vfunc ff(*this);
  401. ff *= f;
  402. return ff;
  403. }
  405. vfunc& vfunc::operator/=(vfunc& fn) // 自身除以一个函数
  406. {
  407. valgo * a = new valgojoin(alg, fn.alg, cdiv);
  408. alg->refnum--; // 因为联合算法对两个算法都加了引用数，因此再减回来
  409. alg = a; // 将新指针赋给alg
  410. return (*this);
  411. }
  413. vfunc& vfunc::operator/=(matrix& (*f)(matrix&)) // 自身除以一个函数指针
  414. {
  415. vfunc fn(f); // 将函数指针包装为函数类
  416. operator/=(fn); // 实施除法操作
  417. return (*this);
  418. }
  420. vfunc vfunc::operator/(vfunc& fn) // 相除产生新函数
  421. {
  422. valgo * a = new valgojoin(alg, fn.alg, cdiv);
  423. vfunc f(a);
  424. return f;
  425. }
  427. vfunc vfunc::operator/(DOUBLE a) // 与常数相除产生新函数
  428. {
  429. vfunc f(*this);
  430. f /= a;
  431. return f;
  432. }
  434. vfunc operator/(DOUBLE a, vfunc& f) // 常数除以函数
  435. {
  436. vfunc ff(a);
  437. return ff/f;
  438. }
  440. vfunc vfunc::operator/(matrix& (*f)(matrix&)) // 除以一个函数指针产生新函数
  441. {
  442. vfunc ff(*this);
  443. ff /= f;
  444. return ff;
  445. }
  447. vfunc operator/(matrix& (*f)(matrix&),vfunc& fn) // 函数指针除以函数
  448. {
  449. vfunc ff(f);
  450. ff /= fn;
  451. return ff;
  452. }
  454. void vfunc::setxfactor(DOUBLE a) // 设置x因子为a
  455. {
  456. alg = alg->setxfactor(a);
  457. }
  459. void vfunc::xroom(DOUBLE a)   // x方向扩大a倍
  460. {
  461. alg = alg->xroom(a);
  462. }
  464. vfunc vfunc::operator()(vfunc & fn) // 复合函数，产生新的函数
  465. {
  466. valgo * a = new valgojoin(alg, fn.alg, ccom);
  467. vfunc f(a);
  468. return f;
  469. }
  471. void vfunc::setxshift(DOUBLE a) // 设置函数沿x轴平移a
  472. {
  473. alg = alg->setxshift(a);
  474. }
  476. void vfunc::shiftxas(DOUBLE a) // 函数沿x轴右移a
  477. {
  478. alg = alg->xshiftas(a);
  479. }
  481. linemodel::linemodel(matrix& va, matrix & vh, matrix & q,
  482.  matrix & r, matrix & vx) :a(va),h(vh),w(q),v(r),x(vx),
  483. y(vh.rownum,1)
  484. // 构造函数，va初始状态转移矩阵，vh初始观测矩阵，q当前模型噪声协
  485. // 方差阵，r当前观测噪声协方差阵，vx初始状态变量
  486. {
  487. if(va.rownum != vx.rownum ||
  488. va.rownum != va.colnum ||
  489. vh.colnum != va.rownum ||
  490. q.rownum != va.rownum ||
  491. r.rownum != vh.rownum) throw TMESSAGE("Data not meet need!");
  492. y = h*x+v();
  493. };
  495. void linemodel::setdata(matrix& va, matrix & vh, matrix & q, matrix & r)
  496. {
  497. if(va.rownum != va.colnum ||
  498. va.rownum != x.rownum ||
  499. vh.colnum != va.rownum ||
  500. q.rownum != va.rownum ||
  501. r.rownum != vh.rownum) throw TMESSAGE("Data not meet need!");
  502. a = va;
  503. h = vh;
  504. w.setdata(q);
  505. v.setdata(r);
  506. }
  508. void linemodel::seta(matrix& va)
  509. {
  510. if(va.rownum != va.colnum ||
  511. va.rownum != x.rownum) throw TMESSAGE("a not meet need!");
  512. a = va;
  513. }
  515. void linemodel::seth(matrix& vh)
  516. {
  517. if(vh.rownum != y.rownum ||
  518. vh.colnum != x.rownum ) throw TMESSAGE("h not meet nead!");
  519. h = vh;
  520. }
  522. void linemodel::setq(matrix& q)
  523. {
  524. if(q.rownum != x.rownum ||
  525. !q.issym) throw TMESSAGE("q not meet nead!");
  526. w.setdata(q);
  527. }
  529. void linemodel::setr(matrix& r)
  530. {
  531. if(r.rownum != y.rownum ||
  532. !r.issym) throw TMESSAGE("r not meet nead!");
  533. v.setdata(r);
  534. }
  536. matrix & linemodel::next() // 计算下一级x值并返回新的x对应的y值
  537. {
  538. x = a*x+w(); // 产生新的状态变量值
  539. return (y=h*x+v()); // 产生并返回新的观测值
  540. }
  542. kalman::kalman(matrix &va,matrix& vq,matrix& vr,matrix& vh,matrix& vx,
  543. matrix& vp):
  544. // 构造函数，va为状态转移矩阵，vh观测矩阵，vq模型噪声协方差阵，
  545. // vr当前观测噪声协方差阵，vx初始状态变量估值，vp初始估值协方差阵
  546. a(va),q(vq),r(vr),h(vh),x(vx),p(vp),y(vh.rownum,1)
  547. {
  548. if( va.rownum != va.colnum ||
  549.  !vq.issym || !vr.issym || !vp.issym ||
  550.  vq.rownum != va.rownum || vr.rownum != vh.rownum ||
  551.  vx.rownum != va.rownum || vp.rownum != vx.rownum ||
  552.  vh.colnum != vx.rownum ) throw TMESSAGE("data not meet need!");
  553. if( !vq.ispositive() || !vr.ispositive() )
  554.  throw TMESSAGE("var matrix not positive!");
  555. y = 0.0;
  556. }
  558. void kalman::setdata(matrix &va,matrix& vq,matrix& vr,matrix& vh)
  559. // 为时变系统随时设置系统参数，va为状态转移矩阵，vh观测矩阵，vq模型噪
  560. // 声协方差阵，vr当前观测噪声协方差阵
  561. {
  562. if(va.rownum != va.colnum || va.rownum != x.rownum ||
  563. va.rownum != vq.rownum || !vq.issym ||
  564. vr.rownum != y.rownum || vh.rownum != y.rownum ||
  565. vh.colnum != x.rownum ) throw TMESSAGE("data not meet need!");
  566. if( !vq.ispositive() || !vr.ispositive() )
  567.  throw TMESSAGE("var matrix not positive!");
  568. a = va;
  569. q = vq;
  570. h = vh;
  571. r = vr;
  572. }
  574. void kalman::seta(matrix& va)
  575. {
  576. if(va.rownum != va.colnum || va.rownum != x.rownum)
  577. throw TMESSAGE("data not meet need!");
  578. a = va;
  579. }
  581. void kalman::seth(matrix& vh)
  582. {
  583. if(vh.rownum != y.rownum || vh.colnum != x.rownum)
  584. throw TMESSAGE("data not meet need!");
  585. h = vh;
  586. }
  588. void kalman::setq(matrix& vq)
  589. {
  590. if(vq.rownum != x.rownum || !vq.issym)
  591. throw TMESSAGE("data not meet need!");
  592. if(!vq.ispositive())
  593. throw TMESSAGE("q is not positive!");
  594. q = vq;
  595. }
  597. void kalman::setr(matrix& vr)
  598. {
  599. if(vr.rownum != y.rownum || !vr.issym)
  600. throw TMESSAGE("data not meet need!");
  601. if(!vr.ispositive())
  602. throw TMESSAGE("r is not positive!");
  603. r = vr;
  604. }
  606. matrix& kalman::next(matrix& y) // 根据下一个观测值获得新的状态变量的估值
  607. {
  608. x = a*x; // 预测
  609. p = a*p*a.t()+q; // 预测误差协方差
  610. matrix k=p*h.t()*((h*p*h.t()+r).inv()); // 新息的增益
  611. x += k*(y-h*x); // 由新息校正x得现时估计
  612. p = (unit(p.rownum)-k*h)*p; // 获得最后的误差协方差，注意unit(n)是
  613. // n阶单位矩阵
  614. return x; // 返回x的估计值
  615. }
  618. regress::regress(matrix& x, matrix& y) // x为mXn维矩阵，y为n个观测值
  619. {
  620. if(x.colnum != y.rownum || x.colnum < x.rownum)
  621. throw TMESSAGE("dimension error");
  622. matrix c(x.rownum+1,x.colnum);
  623. size_t i,j;
  624. for(i=0; i<x.rownum; i++)
  625. for(j=0; j<x.colnum; j++)
  626. c.set(i,j,x(i,j));
  627. for(j=0; j<c.colnum; j++)
  628. c.set(x.rownum,j,1.0);
  629. matrix y1(c*y);
  630. c *= c.t();
  631. a = c.cholesky(y1);
  632. v = matrix(x.rownum);
  633. // 计算偏差平方和
  634. q = 0;
  635. DOUBLE yy;
  636. matrix dd(y.rownum);
  637. for(i=0; i<y.rownum; i++) {
  638. yy = y(i);
  639. for(j=0; j<x.rownum; j++)
  640. yy -= a(j)*x(j,i);
  641. yy -= a(a.rownum-1);
  642. dd.set(i,yy);
  643. q += yy*yy;
  644. }
  645. s = sqrt(q/(y.rownum));
  646. yy = 0.0;
  647. for(i=0; i<y.rownum; i++)
  648. yy += y(i);
  649. DOUBLE ye;
  650. ye = yy/(y.rownum);
  651. r = 0.0;
  652. for(i=0; i<y.rownum; i++)
  653. r += (y(i)-ye)*(y(i)-ye);
  654. r = sqrt(1.0-q/r);
  655. for(j=0;j<v.rownum;j++) {
  656. yy = 0;
  657. for(i=0;i<y.rownum;i++)
  658. yy += (dd(i)+a(j)*x(j,i))*(dd(i)+a(j)*x(j,i));
  659. v.set(j,sqrt(1.0-q/yy));
  660. }
  661. u = 0.0;
  662. for(i=0; i<y.rownum; i++) {
  663. yy = ye;
  664. for(j=0; j<x.rownum; j++)
  665. yy -= a(j)*x(j,i);
  666. yy -= a(x.rownum);
  667. u += yy*yy;
  668. }
  669. }
  671. DOUBLE regress::operator()(matrix& x) // 回归后的线性函数，x是自变量
  672. {
  673. if(x.rownum != a.rownum-1) throw TMESSAGE("wrong dimension!");
  674. DOUBLE y = a(a.rownum-1);
  675. for(size_t i=0; i<x.rownum; i++)
  676. y += x(i)*a(i);
  677. return y;
  678. }