1 /*----------------------------------------------------------------------------
2 linearfit.cc (linear fitting routine)
3 This file is a part of topaz systems
4 Copyright: Hisao Kawaura, All rights reserved
5 1997 - 98
6
7
8 This program is free software; you can redistribute it and/or modify
9 it under the terms of the GNU General Public License as published by
10 the Free Software Foundation; either version 2 of the License, or
11 (at your option) any later version.
12
13 This program is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License for more details.
17
18 You should have received a copy of the GNU General Public License
19 along with this program; if not, write to the Free Software
20 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21
22 ----------------------------------------------------------------------------*/
23
24 #include <stdio.h>
25 #include <stdlib.h>
26 #include <math.h>
27 #include <string>
28
29 #include "linearfit.h"
30 #include "naninf.h"
31 #include "calc.h"
32 #include "graph.h"
33 #include "frame.h"
34 #include "data.h"
35
36 extern int CalcFormulastring(char *s, double *n, int InitFormula, data* dat, int cellno, int seriesno, int Isminmax);
37
getlinearfitfunctionname(int id,std::string * out)38 bool getlinearfitfunctionname(int id, std::string *out)
39 {
40 switch(id)
41 {
42 case FIT_POLY00:
43 *out = std::string("0-poly.");
44 break;
45 case FIT_POLY01:
46 *out = std::string("1-poly.");
47 break;
48 case FIT_POLY02:
49 *out = std::string("2-poly.");
50 break;
51 case FIT_POLY03:
52 *out = std::string("3-poly.");
53 break;
54 case FIT_POLY04:
55 *out = std::string("4-poly.");
56 break;
57 case FIT_POLY05:
58 *out = std::string("5-poly.");
59 break;
60 case FIT_POLY06:
61 *out = std::string("6-poly.");
62 break;
63 case FIT_POLY07:
64 *out = std::string("7-poly.");
65 break;
66 case FIT_POLY08:
67 *out = std::string("8-poly.");
68 break;
69 case FIT_POLY09:
70 *out = std::string("9-poly.");
71 break;
72 case FIT_POLY10:
73 *out = std::string("10-poly.");
74 break;
75 case FIT_EXP:
76 *out = std::string("A*exp(B*x)");
77 break;
78 case FIT_LOG:
79 *out = std::string("A*ln(x) + B");
80 break;
81 case FIT_POW1:
82 *out = std::string("A*x^B");
83 break;
84 case FIT_POW2:
85 *out = std::string("A*B^x");
86 break;
87 default:
88 *out = std::string("");
89 return false;
90 }
91 return true;
92 }
93
CalcLinearLeastSq(doublearray * inx,doublearray * iny,int n,int m,int mode,double * factor,double * cor,double x0,double y0,bool Fixed,int StartLine,int EndLine)94 bool CalcLinearLeastSq(doublearray *inx,
95 doublearray *iny,
96 int n, // number of input data point
97 int m, // polynomial order
98 int mode, // fitting function type
99 double *factor, // fitting factor
100 double *cor, // correlation coefficient
101 double x0, // fixed point
102 double y0,
103 bool Fixed, // go through given fixed point or not
104 int StartLine, // start line for evaluation
105 int EndLine // end line for evaluation
106 )
107 {
108 doublearray *r[12], *a[12], *inxrow = 0, *inyrow = 0;
109 int i, j, k;
110 int normalcount;
111 double dbtemp = 0.0, sum, pivot, p, t;
112 bool ErrorFlag = false;
113 double rss=0.0,ave=0.0,var=0.0;
114 int tempstartline,tempendline;
115
116 for (i = 0; i < 12; i++)
117 {
118 r[i] = 0;
119 a[i] = 0;
120 }
121
122 for (i = 0; i < 12; i++)
123 {
124 r[i] = new doublearray;
125 if (r[i] == 0)
126 {
127 r[i] = 0;
128 ErrorFlag = true;
129 goto j1;
130 }
131 }
132
133 for (i = 0; i < 12; i++)
134 {
135 a[i] = new doublearray;
136 if (a[i] == 0)
137 {
138 a[i] = 0;
139 ErrorFlag = true;
140 goto j1;
141 }
142 }
143
144 inxrow = new doublearray;
145 if (inxrow == 0)
146 {
147 inxrow = 0;
148 ErrorFlag = true;
149 goto j1;
150 }
151
152 inyrow = new doublearray;
153 if (inyrow == 0)
154 {
155 inyrow = 0;
156 ErrorFlag = true;
157 goto j1;
158 }
159
160 // remove invalid data
161 tempstartline = StartLine;
162 if (EndLine == -1)
163 tempendline = n;
164 else
165 {
166 if (EndLine <= n)
167 tempendline = EndLine;
168 else
169 tempendline = n;
170 }
171
172 normalcount = 0;
173 inxrow->add(0);
174 inyrow->add(0);
175 for (i = tempstartline; i <= tempendline; i++)
176 {
177 if (isnormalvalue((*inx)[i - 1]) && isnormalvalue((*iny)[i - 1]))
178 {
179 normalcount++;
180 inxrow->addat(normalcount, (*inx)[i - 1]);
181 inyrow->addat(normalcount, (*iny)[i - 1]);
182 }
183 }
184
185 if (normalcount == 0)
186 {
187 ErrorFlag = true;
188 goto j1;
189 }
190
191 n = normalcount;
192
193 for (i = 0; i < 12; i++)
194 r[i]->setarray(n + 1);
195
196 // generate input function
197 // in case of polynomial
198 if (mode == 1)
199 {
200 if (!Fixed)
201 {
202 for (i = 1; i <= n; i++)
203 (*r[0])[i] = 1.0;
204
205 for (j = 1;j <= m; j++)
206 {
207 for (i = 1; i <= n; i++)
208 {
209 dbtemp = (*r[j - 1])[i] * (*inxrow)[i];
210
211 if (isnormalvalue(dbtemp))
212 {
213 (*r[j])[i] = dbtemp;
214 }
215 else
216 {
217 ErrorFlag = true;
218 goto j1;
219 }
220 }
221 }
222 }
223 else if (Fixed)
224 {
225 if (m == 0)
226 {
227 ErrorFlag = true;
228 goto j1;
229 }
230
231 m--;
232 for (j = 0; j <= m; j++)
233 {
234 for (i = 1; i <= n; i++)
235 {
236 dbtemp = pow((*inxrow)[i], j +1 ) - pow(x0, j + 1);
237 if (isnormalvalue(dbtemp))
238 {
239 (*r[j])[i] = dbtemp;
240 }
241 else
242 {
243 ErrorFlag = true;
244 goto j1;
245 }
246 }
247 }
248 }
249 }
250
251 //in case of exponential (y=Aexp(Bx) -> lny=lnA+Bx)
252 else if (mode == 2)
253 {
254 if (!Fixed)
255 {
256 m = 1;
257 for (i = 1; i <= n; i++)
258 {
259 (*r[0])[i] = 1.0;
260 (*r[1])[i] = (*inxrow)[i];
261 }
262 }
263 else if (Fixed)
264 {
265 m = 0;
266 for (i = 1; i <= n; i++)
267 (*r[0])[i] = (*inxrow)[i] - x0;
268 }
269 }
270
271 //in case of logalithm (y=A + Bln x )
272 else if (mode == 3)
273 {
274 if (!Fixed)
275 {
276 m = 1;
277 for (i = 1;i <= n; i++)
278 {
279 if ((*inxrow)[i] <= 0.0)
280 {
281 ErrorFlag = true;
282 goto j1;
283 }
284
285 (*r[0])[i] = 1.0;
286 (*r[1])[i] = log((*inxrow)[i]);
287 }
288 }
289 else if (Fixed)
290 {
291 if (x0 <= 0.0)
292 {
293 ErrorFlag = true;
294 goto j1;
295 }
296
297 m = 0;
298 for (i = 1; i <= n; i++)
299 {
300 if ((*inxrow)[i] <= 0.0)
301 {
302 ErrorFlag = true;
303 goto j1;
304 }
305
306 (*r[0])[i] = log((*inxrow)[i]) - log(x0);
307 }
308 }
309 }
310 // in case of power1 (y = A*b^x -> ln y = ln A + ln B * x )
311 else if (mode == 4)
312 {
313 if (!Fixed)
314 {
315 m = 1;
316 for (i = 1; i <= n; i++)
317 {
318 (*r[0])[i] = 1.0;
319 (*r[1])[i] = (*inxrow)[i];
320 }
321 }
322 else if (Fixed)
323 {
324 m = 0;
325 for (i = 1; i <= n; i++)
326 (*r[0])[i] = (*inxrow)[i] - x0;
327 }
328 }
329 //in case of power2 (y = A*x^B -> ln y = ln A + B * ln x )
330 else if (mode == 5)
331 {
332 if (!Fixed)
333 {
334 m = 1;
335 for (i = 1; i <= n; i++)
336 {
337 (*r[0])[i] = 1.0;
338 (*r[1])[i] = log((*inxrow)[i]);
339 }
340 }
341 else if (Fixed)
342 {
343 if (x0 <= 0.0)
344 {
345 ErrorFlag = true;
346 goto j1;
347 }
348
349 m = 0;
350 for (i = 1; i <= n; i++)
351 {
352 if ((*inxrow)[i] <= 0.0)
353 {
354 ErrorFlag = true;
355 goto j1;
356 }
357
358 (*r[0])[i] = log((*inxrow)[i]) - log(x0);
359 }
360 }
361 }
362
363 for (i = 0; i < 12; i++)
364 a[i]->setarray(n);
365
366 for (i = 0; i <= n; i++)
367 {
368 for (j = i; j <= m; j++)
369 {
370 sum = 0.0;
371 for (k = 1; k <= n; k++)
372 {
373 dbtemp = sum + (*r[i])[k] * (*r[j])[k];
374 if (isnormalvalue(dbtemp))
375 {
376 sum = dbtemp;
377 }
378 else
379 {
380 ErrorFlag = true;
381 goto j1;
382 }
383 }
384
385 (*a[j])[i] = sum;
386 (*a[i])[j] = sum;
387 }
388 }
389
390 for (i = 0; i <= m; i++)
391 {
392 sum = 0.0;
393 for (k = 1; k <= n; k++)
394 {
395 if (!Fixed)
396 {
397 if (mode == 1 || mode == 3)
398 dbtemp = sum + (*r[i])[k] * (*inyrow)[k];
399 else if (mode == 2 || mode == 4 || mode == 5)
400 {
401 if ((*inyrow)[k] <= 0.0)
402 {
403 ErrorFlag = true;
404 goto j1;
405 }
406 dbtemp = sum + (*r[i])[k] * log((*inyrow)[k]);
407 }
408 }
409 else
410 {
411 if (mode == 1 || mode == 3)
412 dbtemp = sum + (*r[i])[k] * ((*inyrow)[k] - y0);
413 else if (mode == 2 || mode == 4 || mode == 5)
414 {
415 if ((*inyrow)[k] <= 0.0 || y0 <= 0)
416 {
417 ErrorFlag = true;
418 goto j1;
419 }
420 dbtemp = sum + (*r[i])[k] * (log((*inyrow)[k]) - log(y0));
421 }
422 }
423
424 if (isnormalvalue(dbtemp))
425 {
426 sum = dbtemp;
427 }
428 else
429 {
430 ErrorFlag = true;
431 goto j1;
432 }
433
434 }
435
436 (*a[m+1])[i] = sum;
437 }
438
439 for (k = 0; k <= m; k++)
440 {
441 pivot = (*a[k])[k];
442
443 if (pivot == 0.0)
444 {
445 ErrorFlag = true;
446 goto j1;
447 }
448
449 for (j = k; j <= m + 1; j++)
450 (*a[j])[k] /= pivot;
451
452 for (i = 0; i <= m; i++)
453 {
454 if (i != k)
455 {
456 t = (*a[k])[i];
457 for (j = k; j <= m + 1; j++)
458 {
459 dbtemp = (*a[j])[i] - t * (*a[j])[k];
460 if (isnormalvalue(dbtemp))
461 {
462 (*a[j])[i] = dbtemp;
463 }
464 else
465 {
466 ErrorFlag = true;
467 goto j1;
468 }
469 }
470 }
471 }
472 }
473
474 // calculate coefficients
475 for (i = 0; i <= m; i++)
476 factor[i] = (*a[m + 1])[i];
477
478 // polynomial
479 if (mode == 1)
480 {
481 if (Fixed)
482 {
483 m++;
484 for (i = m; i >= 1; i--)
485 factor[i] = factor[i - 1];
486
487 factor[0] = y0;
488 for (i = 1; i <= m; i++)
489 factor[0] -= factor[i] * pow(x0, i);
490 }
491 }
492 else if (mode == 2)
493 {
494 m = 1;
495 if (!Fixed)
496 factor[0] = exp(factor[0]);
497 else
498 {
499 factor[1] = factor[0];
500 factor[0] = log(y0) - factor[1] * x0; factor[0] = exp(factor[0]);
501 }
502
503 }
504 else if (mode == 3)
505 {
506 m = 1;
507 if (Fixed)
508 {
509 if (x0 <= 0.0)
510 {
511 ErrorFlag = true;
512 goto j1;
513 }
514
515 factor[1] = factor[0];
516 factor[0] = y0 - factor[1] * log(x0);
517 }
518 }
519 else if (mode == 4)
520 {
521 m = 1;
522 if (!Fixed)
523 {
524 factor[0] = exp(factor[0]);
525 factor[1] = exp(factor[1]);
526 }
527 else if (Fixed)
528 {
529 factor[1] = exp(factor[0]);
530 if (factor[1] <= 0.0)
531 {
532 ErrorFlag = true;
533 goto j1;
534 }
535 factor[0] = exp(log(y0) - log(factor[1]) * x0);
536 }
537 }
538 else if (mode == 5)
539 {
540 m = 1;
541 if (!Fixed)
542 factor[0] = exp(factor[0]);
543 else if (Fixed)
544 {
545 factor[1] = factor[0];
546 factor[0] = log(y0) - factor[1] * log(x0); factor[0] = exp(factor[0]);
547 }
548 }
549
550 // correlation coefficient
551 if (!Fixed)
552 {
553 for (j = 0; j < n; j++)
554 {
555 if (mode == 1)
556 {
557 p = factor[m];
558 for (i = m - 1; i >= 0; i--)
559 p = p * (*inxrow)[j + 1] + factor[i];
560 rss += (p - (*inyrow)[j + 1]) * (p - (*inyrow)[j + 1]);
561 ave += (*inyrow)[j + 1];
562 }
563 else if (mode == 2)
564 {
565 p = (*inxrow)[j + 1] * factor[1];
566 p += log(factor[0]);
567
568 rss += (p - log((*inyrow)[j + 1])) * (p - log((*inyrow)[j + 1]));
569 ave += log((*inyrow)[j + 1]);
570 }
571 else if (mode == 3)
572 {
573 p = log((*inxrow)[j + 1]);
574 p *= factor[1];
575 p += factor[0];
576
577 rss += (p - (*inyrow)[j + 1]) * (p - (*inyrow)[j + 1]);
578 ave += (*inyrow)[j + 1];
579 }
580 else if (mode == 4)
581 {
582 p = (*inxrow)[j + 1] * log(factor[1]);
583 p += log(factor[0]);
584
585 rss += (p - log((*inyrow)[j + 1])) * (p - log((*inyrow)[j + 1]));
586 ave += log((*inyrow)[j + 1]);
587 }
588 else if (mode == 5)
589 {
590 p = log((*inxrow)[j + 1]) * factor[1];
591 p += log(factor[0]);
592
593 rss += (p - log((*inyrow)[j + 1])) * (p - log((*inyrow)[j + 1]));
594 ave += log((*inyrow)[j + 1]);
595 }
596 }
597
598 ave /= n;
599 for (j = 0; j < n; j++)
600 {
601 if (mode == 1 || mode == 3)
602 {
603 if (!Fixed)
604 var += ((*inyrow)[j + 1] - ave) * ((*inyrow)[j + 1] - ave);
605 else
606 var += ((*inyrow)[j + 1] - y0 - ave) * ((*inyrow)[j + 1] -y0 - ave);
607 }
608 else if(mode == 2 || mode == 4 || mode == 5)
609 {
610 if (!Fixed)
611 var += (log((*inyrow)[j + 1]) - ave) * (log((*inyrow)[j + 1]) - ave);
612 else
613 var += (log((*inyrow)[j + 1]) - log(y0) - ave) * (log((*inyrow)[j + 1]) - log(y0) - ave);
614 }
615 }
616
617 *cor = sqrt(1 - rss / var);
618 }
619 j1:;
620
621 if (inxrow != 0)
622 delete inxrow;
623
624 if (inyrow != 0)
625 delete inyrow;
626
627 for (i = 0; i < 12; i++)
628 {
629 if (r[i] != 0)
630 delete r[i];
631 }
632
633 for (i = 0; i < 12; i++)
634 {
635 if (a[i] != 0)
636 delete a[i];
637 }
638
639 if (ErrorFlag)
640 return false;
641 else
642 return true;
643 }
644
ExecLinearSQ(int fittingCurve,bool Fixed,double fixed_x,double fixed_y,int startline,int endline,std::string * outfnstr,std::string * infostr)645 bool data::ExecLinearSQ(int fittingCurve, bool Fixed, double fixed_x, double fixed_y, int startline, int endline, std::string* outfnstr, std::string *infostr)
646 {
647 doublearray *inx = 0, *iny = 0;
648 int inputno, i;
649 bool ErrorFlag = false;
650 double factor[20];
651 double cor;
652 double x0;
653 double y0;
654 int StartLine;
655 int EndLine;
656 char stemp[256];
657 std::string temp_pre_formula, temp_pre_formula_x, temp_pre_formula_y;
658 std::string mes;
659 int NumPoint;
660 double fa, fb;
661 frame *fr = (frame *)parent;
662 std::string outfn;
663 int m = 0;
664 int mode = 1;
665
666 inx = new doublearray;
667 if (inx == 0)
668 {
669 inx = 0;
670 ErrorFlag = true;
671 goto j1;
672 }
673
674 iny = new doublearray;
675 if (iny == 0)
676 {
677 iny = 0;
678 ErrorFlag = true;
679 goto j1;
680 }
681
682
683 //calclation
684 NumPoint = x->getarraylength() - 1;
685 if (trans_x)
686 temp_pre_formula_x = *format_x;
687 else
688 temp_pre_formula_x = std::string("x");
689 if (trans_y)
690 temp_pre_formula_y = *format_y;
691 else
692 temp_pre_formula_y = std::string("y");
693
694 for (i = 0; i <= NumPoint; i++)
695 {
696 CalcFormulastring((char *)temp_pre_formula_x.c_str(), &fa, true ,this, 0, fr->getdataindex(id), 0);
697 CalcFormulastring((char *)temp_pre_formula_x.c_str(), &fa, false ,this, i, fr->getdataindex(id), 0);
698 CalcFormulastring((char *)temp_pre_formula_y.c_str(), &fb, true ,this, 0, fr->getdataindex(id), 0);
699 CalcFormulastring((char *)temp_pre_formula_y.c_str(), &fb, false ,this, i, fr->getdataindex(id), 0);
700 if (isnormalvalue(fa) && (*x)[i].status == 0 &&
701 isnormalvalue(fb) && (*y)[i].status == 0
702 )
703 {
704 if (!inx->add(fa))
705 goto j1;
706 if (!iny->add(fb))
707 goto j1;
708 }
709 }
710
711 // data no.
712 inputno = inx->getarraylength();
713
714 // order
715 switch (fittingCurve)
716 {
717 case 0:
718 case 1:
719 case 2:
720 case 3:
721 case 4:
722 case 5:
723 case 6:
724 case 7:
725 case 8:
726 case 9:
727 case 10:
728 m = fittingCurve;
729 break;
730 case 11:
731 case 12:
732 case 13:
733 case 14:
734 m = 1;
735 break;
736 }
737 //mode
738 switch (fittingCurve)
739 {
740 case 0:
741 case 1:
742 case 2:
743 case 3:
744 case 4:
745 case 5:
746 case 6:
747 case 7:
748 case 8:
749 case 9:
750 case 10:
751 mode = 1;
752 break;
753 case 11:
754 mode = 2;
755 break;
756 case 12:
757 mode = 3;
758 break;
759 case 13:
760 mode = 4;
761 break;
762 case 14:
763 mode = 5;
764 break;
765 }
766
767 //fixed point X
768 x0 = fixed_x;
769 //fixed point Y
770 y0 = fixed_y;
771
772 //startline
773 StartLine = startline;
774 //endline
775 EndLine = endline;
776
777 if (CalcLinearLeastSq(inx
778 ,iny
779 ,inputno
780 , m
781 , mode
782 , factor
783 , &cor
784 , x0
785 , y0
786 , Fixed
787 , StartLine
788 , EndLine
789 ))
790 {
791
792 // generate message
793 // source data
794 mes = std::string("--------- fitting results ---------\\n");
795 mes += std::string("\\n");
796 mes += std::string("Original data: ");
797 mes += std::string(filename);
798 mes += std::string("\\n");
799 mes += std::string("\\n");
800 outfn = std::string("");
801
802 switch(mode)
803 {
804 case 1:
805 mes += std::string("y = sigma Ai * x^i\\n");
806 mes += std::string("\\n");
807 for (i = 0; i <= m; i++)
808 {
809 sprintf(stemp, " A%d = %.15g\\n", i, factor[i]);
810 mes += std::string(stemp);
811 }
812 for (i = m; i >= 0; i--)
813 {
814 sprintf(stemp,"%+.15g*x^(%d)", factor[i], i);
815 outfn += std::string(stemp);
816 }
817 break;
818 case 2:
819 mes += std::string("y = A * exp( B * x)\\n");
820 mes += std::string("\\n");
821 sprintf(stemp, " A = %.15g\\n", factor[0]);
822 mes += std::string(stemp);
823 sprintf(stemp, " B = %.15g\\n", factor[1]);
824 mes += std::string(stemp);
825 sprintf(stemp,"%+.15g*exp(%.15g*x)", factor[0], factor[1]);
826 outfn += std::string(stemp);
827 break;
828 case 3:
829 mes += std::string("y = A + B * ln( x )\\n");
830 mes += std::string("\\n");
831 sprintf(stemp, " A = %.15g\\n", factor[0]);
832 mes += std::string(stemp);
833 sprintf(stemp, " B = %.15g\\n", factor[1]);
834 mes += std::string(stemp);
835 sprintf(stemp,"%+.15g%+.15g*ln(x)", factor[0], factor[1]);
836 outfn += std::string(stemp);
837 break;
838 case 4:
839 mes += std::string("y = A * x ^ B\\n");
840 mes += std::string("\\n");
841 sprintf(stemp, " A = %.15g\\n", factor[0]);
842 mes += std::string(stemp);
843 sprintf(stemp, " B = %.15g\\n", factor[1]);
844 mes += std::string(stemp);
845 sprintf(stemp,"%.15g*%15g^x", factor[0], factor[1]);
846 outfn += std::string(stemp);
847 break;
848 case 5:
849 mes += std::string("y = A * B ^ x\\n");
850 mes += std::string("\\n");
851 sprintf(stemp, " A = %.15g\\n", factor[0]);
852 mes += std::string(stemp);
853 sprintf(stemp, " B = %.15g\\n", factor[1]);
854 mes += std::string(stemp);
855 sprintf(stemp,"%.15g*x^(%.15g)", factor[0], factor[1]);
856 outfn += std::string(stemp);
857 break;
858 }
859 if (!Fixed)
860 {
861 mes += std::string("\\n");
862 sprintf(stemp, " Corr. Coeff. = %.15g\\n", cor);
863 mes += std::string(stemp);
864 }
865
866 mes += std::string("\\n");
867 mes += std::string("\\n");
868
869 sprintf(stemp, "Fitting line: %d -> %d\\n", StartLine, EndLine);
870 mes += std::string(stemp);
871
872 if (Fixed)
873 {
874 mes += std::string("Go through fixed point: true\\n");
875 sprintf(stemp, "Fixed point: (%.15g, %.15g)\\n", x0, y0);
876 mes += std::string(stemp);
877 }
878 mes += std::string("\\n");
879 mes += std::string("--------------- end ---------------\\n");
880
881 *outfnstr = outfn;
882 *infostr = mes;
883 }
884 else
885 {
886 ErrorFlag = true;
887 goto j1;
888 }
889 j1:;
890
891 if (inx != 0)
892 delete inx;
893 if (iny != 0)
894 delete iny;
895
896 if (ErrorFlag)
897 return false;
898 else
899 return true;
900 }
901
902