1 /*
2 * $Id: ofit.c,v 1.6 2002/07/06 08:51:42 isizaka Exp isizaka $
3 *
4 * This file is part of "Ngraph for X11".
5 *
6 * Copyright (C) 2002, Satoshi ISHIZAKA. isizaka@msa.biglobe.ne.jp
7 *
8 * "Ngraph for X11" is free software; you can redistribute it and/or
9 * modify it under the terms of the GNU General Public License
10 * as published by the Free Software Foundation; either version 2
11 * of the License, or (at your option) any later version.
12 *
13 * "Ngraph for X11" 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 /**
25 *
26 * $Log: ofit.c,v $
27 * Revision 1.6 2002/07/06 08:51:42 isizaka
28 * change to GPL.
29 *
30 * Revision 1.5 2001/03/23 12:15:31 isizaka
31 * for 6.3.13
32 *
33 * Revision 1.4 1999/04/15 12:15:27 isizaka
34 * for release 6.03.01
35 *
36 * Revision 1.3 1999/04/11 06:08:10 isizaka
37 * *** empty log message ***
38 *
39 * Revision 1.2 1999/03/20 12:30:47 isizaka
40 * add 'profile' field
41 *
42 * Revision 1.1 1999/03/17 13:46:09 isizaka
43 * Initial revision
44 *
45 *
46 **/
47
48 #include <stdio.h>
49 #include <stdlib.h>
50 #include <string.h>
51 #include <ctype.h>
52 #include <math.h>
53 #include "ngraph.h"
54 #include "object.h"
55 #include "mathcode.h"
56 #include "mathfn.h"
57 #include "oroot.h"
58
59 #define NAME "fit"
60 #define PARENT "object"
61 #define VERSION "1.00.01"
62 #define TRUE 1
63 #define FALSE 0
64
65 #define ERRNUM 14
66
67 #define ERRSYNTAX 100
68 #define ERRILLEGAL 101
69 #define ERRNEST 102
70 #define ERRMANYPARAM 103
71 #define ERRLN 104
72 #define ERRSMLDATA 105
73 #define ERRPOINT 106
74 #define ERRMATH 107
75 #define ERRSINGULAR 108
76 #define ERRNOEQS 109
77 #define ERRTHROUGH 110
78 #define ERRRANGE 111
79 #define ERRNEGATIVEWEIGHT 112
80 #define ERRCONVERGE 113
81
82 char *fiterrorlist[ERRNUM]={
83 "syntax error.",
84 "not allowd function.",
85 "sum() or dif(): deep nest.",
86 "too many parameters.",
87 "illegal data -> ignored.",
88 "too small number of data.",
89 "illegal point",
90 "math error -> ignored.",
91 "singular matrix.",
92 "no fitting equation.",
93 "`through_point' for type `user' is not supported.",
94 "math range check.",
95 "negative or zero weight -> ignored.",
96 "convergence error.",
97 };
98
99 char *fittypechar[6]={
100 "poly",
101 "pow",
102 "exp",
103 "log",
104 "user",
105 NULL
106 };
107
108 struct fitlocal {
109 char *codef;
110 struct narray *needdata;
111 char *codedf[10];
112 int dim;
113 double coe[11];
114 int num;
115 double derror,correlation;
116 char *equation;
117 };
118
fitinit(struct objlist * obj,char * inst,char * rval,int argc,char ** argv)119 int fitinit(struct objlist *obj,char *inst,char *rval,int argc,char **argv)
120 {
121 int div,dimension;
122 double converge;
123 struct fitlocal *fitlocal;
124 int i,disp;
125
126 if (_exeparent(obj,(char *)argv[1],inst,rval,argc,argv)) return 1;
127 div=500;
128 dimension=1;
129 converge=1;
130 disp=TRUE;
131 if (_putobj(obj,"div",inst,&div)) return 1;
132 if (_putobj(obj,"poly_dimension",inst,&dimension)) return 1;
133 if (_putobj(obj,"converge",inst,&converge)) return 1;
134 if (_putobj(obj,"display",inst,&disp)) return 1;
135
136 if ((fitlocal=memalloc(sizeof(struct fitlocal)))==NULL) return 1;
137 fitlocal->codef=NULL;
138 for (i=0;i<10;i++) fitlocal->codedf[i]=NULL;
139 fitlocal->needdata=NULL;
140 fitlocal->equation=NULL;
141 if (_putobj(obj,"_local",inst,fitlocal)) {
142 memfree(fitlocal);
143 return 1;
144 }
145 return 0;
146 }
147
fitdone(struct objlist * obj,char * inst,char * rval,int argc,char ** argv)148 int fitdone(struct objlist *obj,char *inst,char *rval,int argc,char **argv)
149 {
150 struct fitlocal *fitlocal;
151 int i;
152
153 if (_exeparent(obj,(char *)argv[1],inst,rval,argc,argv)) return 1;
154 _getobj(obj,"_local",inst,&fitlocal);
155 memfree(fitlocal->codef);
156 memfree(fitlocal->equation);
157 arrayfree(fitlocal->needdata);
158 for (i=0;i<10;i++) memfree(fitlocal->codedf[i]);
159 return 0;
160 }
161
fitput(struct objlist * obj,char * inst,char * rval,int argc,char ** argv)162 int fitput(struct objlist *obj,char *inst,char *rval,
163 int argc,char **argv)
164 {
165 char *field;
166 int rcode,ecode,maxdim,need2pass;
167 struct narray *needdata;
168 char *math,*code;
169 struct fitlocal *fitlocal;
170 char *equation;
171
172 _getobj(obj,"_local",inst,&fitlocal);
173 field=argv[1];
174 if (strcmp(field,"poly_dimension")==0) {
175 if (*(int *)(argv[2])<1) *(int *)(argv[2])=1;
176 else if (*(int *)(argv[2])>9) *(int *)(argv[2])=9;
177 } else if ((strcmp(field,"user_func")==0)
178 || ((strncmp(field,"derivative",10)==0) && (field[10]!='\0'))) {
179 math=(char *)(argv[2]);
180 if (math!=NULL) {
181 needdata=arraynew(sizeof(int));
182 rcode=mathcode(math,&code,needdata,NULL,&maxdim,&need2pass,
183 TRUE,FALSE,FALSE,TRUE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE);
184 if (field[0]!='u') {
185 arrayfree(needdata);
186 }
187 if (rcode!=MCNOERR) {
188 if (rcode==MCSYNTAX) ecode=ERRSYNTAX;
189 else if (rcode==MCILLEGAL) ecode=ERRILLEGAL;
190 else if (rcode==MCNEST) ecode=ERRNEST;
191 error(obj,ecode);
192 return 1;
193 }
194 if (maxdim>9) {
195 error(obj,ERRMANYPARAM);
196 return 1;
197 }
198 } else {
199 code=NULL;
200 needdata=NULL;
201 }
202 if (field[0]=='u') {
203 memfree(fitlocal->codef);
204 fitlocal->codef=code;
205 arrayfree(fitlocal->needdata);
206 fitlocal->needdata=needdata;
207 } else {
208 memfree(fitlocal->codedf[field[10]-'0']);
209 fitlocal->codedf[field[10]-'0']=code;
210 }
211 }
212 _getobj(obj,"equation",inst,&equation);
213 if (_putobj(obj,"equation",inst,NULL)) return 1;
214 memfree(equation);
215 return 0;
216 }
217
fitpoly(struct fitlocal * fitlocal,int type,int dimension,int through,double x0,double y0,double * data,int num,int disp,int weight,double * wdata)218 int fitpoly(struct fitlocal *fitlocal,
219 int type,int dimension,int through,double x0,double y0,
220 double *data,int num,int disp,int weight,double *wdata)
221 /*
222 return
223 -5: range check
224 -2: singular matrix
225 -1: too small data
226 0: no error
227 1: fatail error
228 */
229 {
230 int i,j,k,dim;
231 double yy,y1,y2,derror,sy,correlation,wt,sum;
232 vector b,x1,x2,coe;
233 matrix m;
234 char *equation;
235 char buf[1024];
236
237 if (type==0) dim=dimension+1;
238 else dim=2;
239 if (!through && (num<dim)) return -1;
240 if (through && (num<dim+1)) return -1;
241 yy=0;
242 for (i=0;i<dim+1;i++) {
243 b[i]=0;
244 for (j=0;j<dim+1;j++) m[i][j]=0;
245 }
246 sum=0;
247 for (k=0;k<num;k++) {
248 if (weight) wt=wdata[k];
249 else wt=1;
250 sum+=wt;
251 y1=data[k*2+1];
252 x1[0]=1;
253 for (i=1;i<dim;i++) {
254 if ((fabs(x1[i-1])>1e100) || (fabs(data[k*2])>1e100)) return -5;
255 x1[i]=x1[i-1]*data[k*2];
256 if (fabs(x1[i])>1e100) return -5;
257 }
258 yy+=wt*y1*y1;
259 if (fabs(yy)>1e100) return -5;
260 for (i=0;i<dim;i++) {
261 if (fabs(y1)>1e100) return -5;
262 b[i]+=wt*y1*x1[i];
263 if (fabs(b[i])>1e100) return -5;
264 }
265 for (i=0;i<dim;i++)
266 for (j=0;j<dim;j++) {
267 m[i][j]=m[i][j]+wt*x1[i]*x1[j];
268 if (fabs(m[i][j])>1e100) return -5;
269 }
270 }
271 if (through) {
272 y2=y0;
273 x2[0]=1;
274 for (i=1;i<dim;i++) x2[i]=x2[i-1]*x0;
275 for (i=0;i<dim;i++) {
276 m[i][dim]=x2[i];
277 m[dim][i]=x2[i];
278 }
279 b[dim]=y2;
280 dim++;
281 }
282 if (matsolv(dim,m,b,coe)) return -2;
283 derror=yy;
284 for (i=0;i<dim;i++) derror-=coe[i]*b[i];
285 derror=fabs(derror)/sum;
286 sy=yy/sum-(b[0]/sum)*(b[0]/sum);
287 if ((sy==0) || (1-derror/sy<0)) correlation=-1;
288 else correlation=sqrt(1-derror/sy);
289 derror=sqrt(derror);
290 if (through) dim--;
291 for (i=0;i<dim;i++) fitlocal->coe[i]=coe[i];
292 for (;i<10;i++) fitlocal->coe[i]=0;
293 fitlocal->dim=dim;
294 fitlocal->derror=derror;
295 fitlocal->correlation=correlation;
296 fitlocal->num=num;
297 if ((equation=memalloc(512))==NULL) return 1;
298 equation[0]='\0';
299 j=0;
300 if (type==0) {
301 for (i=dim-1;i>1;i--) j+=sprintf(equation+j,"%.15e*X^%d+",coe[i],i);
302 j+=sprintf(equation+j,"%.15e*X+%.15e",coe[1],coe[0]);
303 } else if (type==1) sprintf(equation,"exp(%.15e)*X^%.15e",coe[0],coe[1]);
304 else if (type==2) sprintf(equation,"exp(%.15e*X+%.15e)",coe[1],coe[0]);
305 else if (type==3) sprintf(equation,"%.15e*Ln(X)+%.15e",coe[1],coe[0]);
306 fitlocal->equation=equation;
307
308 if (disp) {
309 i=0;
310 if (type==0) i+=sprintf(buf+i,"Eq: %%0i*X^i (i=0-%d)\n\n",dim-1);
311 else if (type==1) i+=sprintf(buf+i,"Eq: exp(%%00)*X^%%01\n\n");
312 else if (type==2) i+=sprintf(buf+i,"Eq: exp(%%01*X+%%00)\n\n");
313 else if (type==3) i+=sprintf(buf+i,"Eq: %%01*Ln(X)+%%00\n\n");
314 for (j=0;j<dim;j++)
315 i+=sprintf(buf+i," %%0%d = %+.7e\n",j,coe[j]);
316 i+=sprintf(buf+i,"\n");
317 i+=sprintf(buf+i," points = %d\n",num);
318 i+=sprintf(buf+i," <DY^2> = %+.7e\n",derror);
319 if (correlation>=0)
320 i+=sprintf(buf+i,"|r| or |R| = %+.7e\n",correlation);
321 else
322 i+=sprintf(buf+i,"|r| or |R| = -------------\n");
323 ndisplaydialog(buf);
324 }
325
326 return 0;
327 }
328
fituser(struct objlist * obj,struct fitlocal * fitlocal,char * func,int deriv,double converge,double * data,int num,int disp,int weight,double * wdata)329 int fituser(struct objlist *obj,struct fitlocal *fitlocal,char *func,
330 int deriv,double converge,double *data,int num,int disp,
331 int weight,double *wdata)
332 /*
333 return
334 -7: interrupt
335 -6: convergence
336 -5: range check
337 -4: syntax error
338 -3: eqaution is not specified
339 -2: singular matrix
340 -1: too small data
341 0: no error
342 1: fatal error
343
344 */
345 {
346 int ecode;
347 int *needdata;
348 int tbl[10],dim,n,count,rcode,err,err2,err3;
349 double yy,y,y1,y2,y3,sy,spx,spy,dxx,dxxc,xx,derror,correlation;
350 double b[10],x2[10],par[10],par2[10],parerr[10];
351 char parstat[10];
352 matrix m;
353 int i,j,k,pnum;
354 char *equation;
355 char buf[1024];
356 double wt,sum;
357 /*
358 matrix m2;
359 int count2;
360 double sum2;
361 double lambda,s0;
362 */
363 if (num<1) return -1;
364 if (fitlocal->codef==NULL) return -3;
365 if (fitlocal->needdata==NULL) return -3;
366 dim=arraynum(fitlocal->needdata);
367 needdata=arraydata(fitlocal->needdata);
368 if (deriv) {
369 for (i=0;i<dim;i++) {
370 if (fitlocal->codedf[needdata[i]]==NULL) return -3;
371 }
372 }
373 for (i=0;i<dim;i++) tbl[i]=needdata[i];
374 for (i=0;i<10;i++) {
375 parstat[i]=MNOERR;
376 par[i]=fitlocal->coe[i];
377 }
378
379 ecode=0;
380 /*
381 err2=FALSE;
382 n=0;
383 sum=0;
384 yy=0;
385 lambda=0.001;
386 for (k=0;k<num;k++) {
387 if (weight) wt=wdata[k];
388 else wt=1;
389 sum+=wt;
390 spx=data[k*2];
391 spy=data[k*2+1];
392 err=FALSE;
393 rcode=calculate(fitlocal->codef,1,spx,MNOERR,0,0,0,0,
394 0,0,0,0,0,0,0,0,0,0,0,0,
395 par,parstat,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,
396 NULL,NULL,NULL,0,NULL,NULL,NULL,0,&y1);
397 if (rcode==MSERR) return -4;
398 else if (rcode!=MNOERR) err=TRUE;
399 if (!err) {
400 if (fabs(yy)>1e100) return -5;
401 y2=spy-y1;
402 if ((fabs(y2)>1e50) || (fabs(spy)>1e50)) err=TRUE;
403 }
404 if (!err) {
405 n++;
406 yy+=wt*y2*y2;
407 } else err2=TRUE;
408 }
409 if (err2) error(obj,ERRMATH);
410 if (n<1) {
411 ecode=-1;
412 goto errexit;
413 }
414 s0=yy/sum;
415 */
416
417 count=0;
418 err3=FALSE;
419 do {
420 yy=0;
421 y=0;
422 y3=0;
423 for (i=0;i<dim;i++) {
424 b[i]=0;
425 for (j=0;j<dim;j++) m[j][i]=0;
426 }
427 err2=FALSE;
428 n=0;
429 sum=0;
430 for (k=0;k<num;k++) {
431 if (weight) wt=wdata[k];
432 else wt=1;
433 sum+=wt;
434 spx=data[k*2];
435 spy=data[k*2+1];
436 err=FALSE;
437 rcode=calculate(fitlocal->codef,1,spx,MNOERR,0,0,0,0,
438 0,0,0,0,0,0,0,0,0,0,0,0,0,0,
439 par,parstat,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,
440 NULL,NULL,NULL,0,NULL,NULL,NULL,0,&y1);
441 if (rcode==MSERR) return -4;
442 else if (rcode!=MNOERR) err=TRUE;
443 if (deriv) {
444 for (j=0;j<dim;j++) {
445 rcode=calculate(fitlocal->codedf[tbl[j]],1,spx,MNOERR,0,0,0,0,
446 0,0,0,0,0,0,0,0,0,0,0,0,0,0,
447 par,parstat,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,
448 NULL,NULL,NULL,0,NULL,NULL,NULL,0,&(x2[j]));
449 if (rcode==MSERR) return -4;
450 else if (rcode!=MNOERR) err=TRUE;
451 }
452 } else {
453 for (j=0;j<dim;j++) {
454 for (i=0;i<10;i++) par2[i]=par[i];
455 dxx=par2[j]*converge*1e-6;
456 if (dxx==0) dxx=1e-6;
457 par2[tbl[j]]+=dxx;
458 rcode=calculate(fitlocal->codef,1,spx,MNOERR,0,0,0,0,
459 0,0,0,0,0,0,0,0,0,0,0,0,0,0,
460 par2,parstat,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,
461 NULL,NULL,NULL,0,NULL,NULL,NULL,0,&(x2[j]));
462 if (rcode==MSERR) return -4;
463 else if (rcode!=MNOERR) err=TRUE;
464 x2[j]=(x2[j]-y1)/dxx;
465 }
466 }
467 if (!err) {
468 if ((fabs(yy)>1e100) || (fabs(y3)>1e100)) return -5;
469 y2=spy-y1;
470 if ((fabs(y2)>1e50) || (fabs(spy)>1e50)) err=TRUE;
471 }
472 if (!err) {
473 n++;
474 yy+=wt*y2*y2;
475 y+=wt*spy;
476 y3+=wt*spy*spy;
477 for (j=0;j<dim;j++) b[j]+=wt*x2[j]*y2;
478 for (j=0;j<dim;j++)
479 for (i=0;i<dim;i++) m[j][i]=m[j][i]+wt*x2[j]*x2[i];
480 } else err2=TRUE;
481 }
482 if (!err3 && err2) {
483 error(obj,ERRMATH);
484 err3=TRUE;
485 }
486 if (n<1) {
487 ecode=-1;
488 goto errexit;
489 }
490
491 count++;
492
493 derror=fabs(yy)/sum;
494 sy=y3/sum-(y/sum)*(y/sum);
495 if ((sy==0) || (1-derror/sy<0)) correlation=-1;
496 else correlation=sqrt(1-derror/sy);
497 if (correlation>1) correlation=-1;
498 derror=sqrt(derror);
499
500 if (disp) {
501 i=0;
502 i+=sprintf(buf+i,"Eq: User defined\n\n");
503 for (j=0;j<dim;j++)
504 i+=sprintf(buf+i," %%0%d = %+.7e\n",tbl[j],par[tbl[j]]);
505 i+=sprintf(buf+i,"\n");
506 i+=sprintf(buf+i," points = %d\n",n);
507 if (count==1) i+=sprintf(buf+i," delta = \n");
508 else i+=sprintf(buf+i," delta = %+.7e\n",dxxc);
509 i+=sprintf(buf+i," <DY^2> = %+.7e\n",derror);
510 if (correlation>=0)
511 i+=sprintf(buf+i,"|r| or |R| = %+.7e\n",correlation);
512 else
513 i+=sprintf(buf+i,"|r| or |R| = -------------\n");
514 i+=sprintf(buf+i," Iteration = %d\n",count);
515 ndisplaydialog(buf);
516 }
517 sprintf(buf,"Iteration = %d",count);
518 ndisplaystatus(buf);
519
520 if (ninterrupt()) {
521 ecode=-7;
522 goto errexit;
523 }
524
525 /*
526 count2=0;
527 sum2=sum;
528 while (TRUE) {
529 count2++;
530 sprintf(buf,"Iteration = %d:%d",count,count2);
531 ndisplaystatus(buf);
532 for (j=0;j<dim;j++)
533 for (i=0;i<dim;i++) m2[j][i]=m[j][i];
534 for (i=0;i<dim;i++) m2[i][i]=m[i][i]+sum2*lambda;
535 if (matsolv(dim,m2,b,parerr)) goto repeat;
536 for (i=0;i<dim;i++) par2[tbl[i]]=par[tbl[i]]+parerr[i];
537 n=0;
538 err2=FALSE;
539 sum=0;
540 yy=0;
541 for (k=0;k<num;k++) {
542 if (weight) wt=wdata[k];
543 else wt=1;
544 sum+=wt;
545 spx=data[k*2];
546 spy=data[k*2+1];
547 err=FALSE;
548 rcode=calculate(fitlocal->codef,1,spx,MNOERR,0,0,0,0,
549 0,0,0,0,0,0,0,0,0,0,0,0,
550 par2,parstat,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,
551 NULL,NULL,NULL,0,NULL,NULL,NULL,0,&y1);
552 if (rcode==MSERR) return -4;
553 else if (rcode!=MNOERR) err=TRUE;
554 if (!err) {
555 if (fabs(yy)>1e100) goto repeat;
556 y2=spy-y1;
557 if ((fabs(y2)>1e50) || (fabs(spy)>1e50)) err=TRUE;
558 }
559 if (!err) {
560 n++;
561 yy+=wt*y2*y2;
562 }
563 }
564 if (n<1) goto repeat;
565 if (yy/sum<=s0) {
566 s0=yy/sum;
567 break;
568 }
569 repeat:
570 lambda*=10;
571 if (lambda>1e100) return -6;
572 }
573
574 lambda/=10;
575 */
576 if (matsolv(dim,m,b,parerr)) return -2;
577
578 dxxc=0;
579 xx=0;
580 for (i=0;i<dim;i++) {
581 dxxc+=parerr[i]*parerr[i];
582 par[tbl[i]]+=parerr[i];
583 xx+=par[tbl[i]]*par[tbl[i]];
584 }
585 dxxc=sqrt(dxxc);
586 xx=sqrt(xx);
587
588
589 } while ((dxxc>xx*converge/100) && ((xx>1e-6) || (dxxc>1e-6*converge/100)));
590
591 if (disp) {
592 i=0;
593 i+=sprintf(buf+i,"Eq: User defined\n\n");
594 for (j=0;j<dim;j++)
595 i+=sprintf(buf+i," %%0%d = %+.7e\n",tbl[j],par[tbl[j]]);
596 i+=sprintf(buf+i,"\n");
597 i+=sprintf(buf+i," points = %d\n",n);
598 i+=sprintf(buf+i," delta = %+.7e\n",dxxc);
599 i+=sprintf(buf+i," <DY^2> = %+.7e\n",derror);
600 if (correlation>=0)
601 i+=sprintf(buf+i,"|r| or |R| = %+.7e\n",correlation);
602 else
603 i+=sprintf(buf+i,"|r| or |R| = -------------\n");
604 ndisplaydialog(buf);
605 }
606
607 errexit:
608 if ((ecode==0) || (ecode==-5)) {
609 for (i=0;i<10;i++) fitlocal->coe[i]=par[i];
610 fitlocal->dim=dim;
611 fitlocal->derror=derror;
612 fitlocal->correlation=correlation;
613 fitlocal->num=n;
614 pnum=0;
615 for (i=0;func[i]!='\0';i++)
616 if (func[i]=='%') {
617 pnum++;
618 i+=2;
619 }
620 if ((equation=memalloc(strlen(func)+25*pnum+1))==NULL) return 1;
621 j=0;
622 for (i=0;func[i]!='\0';i++)
623 if (func[i]!='%') {
624 equation[j]=func[i];
625 j++;
626 } else {
627 if (isdigit(func[i+1]) && isdigit(func[i+2]) && isdigit(func[i+3])) {
628 pnum=(func[i+1]-'0')*100+(func[i+2]-'0')*10+(func[i+3]-'0');
629 i+=3;
630 } else if (isdigit(func[i+1]) && isdigit(func[i+2])) {
631 pnum=(func[i+1]-'0')*10+(func[i+2]-'0');
632 i+=2;
633 } else {
634 pnum=(func[i+1]-'0');
635 i+=1;
636 }
637 j+=sprintf(equation+j,"%.15e",par[pnum]);
638 }
639 equation[j]='\0';
640 fitlocal->equation=equation;
641 }
642 return ecode;
643 }
644
fitfit(struct objlist * obj,char * inst,char * rval,int argc,char ** argv)645 int fitfit(struct objlist *obj,char *inst,char *rval,int argc,char **argv)
646 {
647 struct fitlocal *fitlocal;
648 int i,type,through,dimension,deriv,disp;
649 double x,y,x0,y0,converge,wt;
650 struct narray *darray;
651 double *data,*wdata;
652 char *equation,*func;
653 int dnum,num,err,err2,err3,rcode;
654 double derror,correlation,pp;
655 int ecode;
656 int weight,anum;
657
658 if (_exeparent(obj,(char *)argv[1],inst,rval,argc,argv)) return 1;
659 _getobj(obj,"equation",inst,&equation);
660 memfree(equation);
661 if (_putobj(obj,"equation",inst,NULL)) return 1;
662 if ((equation=memalloc(6))==NULL) return 1;
663 strcpy(equation,"undef");
664 if (_putobj(obj,"equation",inst,equation)) {
665 memfree(equation);
666 return 1;
667 }
668 num=0;
669 derror=0;
670 correlation=0;
671 pp=0;
672 if (_putobj(obj,"number",inst,&num)) return 1;
673 if (_putobj(obj,"error",inst,&derror)) return 1;
674 if (_putobj(obj,"correlation",inst,&correlation)) return 1;
675 if (_putobj(obj,"%00",inst,&pp)) return 1;
676 if (_putobj(obj,"%01",inst,&pp)) return 1;
677 if (_putobj(obj,"%02",inst,&pp)) return 1;
678 if (_putobj(obj,"%03",inst,&pp)) return 1;
679 if (_putobj(obj,"%04",inst,&pp)) return 1;
680 if (_putobj(obj,"%05",inst,&pp)) return 1;
681 if (_putobj(obj,"%06",inst,&pp)) return 1;
682 if (_putobj(obj,"%07",inst,&pp)) return 1;
683 if (_putobj(obj,"%08",inst,&pp)) return 1;
684 if (_putobj(obj,"%09",inst,&pp)) return 1;
685
686 _getobj(obj,"_local",inst,&fitlocal);
687 _getobj(obj,"type",inst,&type);
688 _getobj(obj,"through_point",inst,&through);
689 _getobj(obj,"point_x",inst,&x0);
690 _getobj(obj,"point_y",inst,&y0);
691 _getobj(obj,"poly_dimension",inst,&dimension);
692
693 _getobj(obj,"user_func",inst,&func);
694 _getobj(obj,"derivative",inst,&deriv);
695 _getobj(obj,"converge",inst,&converge);
696 _getobj(obj,"parameter0",inst,&(fitlocal->coe[0]));
697 _getobj(obj,"parameter1",inst,&(fitlocal->coe[1]));
698 _getobj(obj,"parameter2",inst,&(fitlocal->coe[2]));
699 _getobj(obj,"parameter3",inst,&(fitlocal->coe[3]));
700 _getobj(obj,"parameter4",inst,&(fitlocal->coe[4]));
701 _getobj(obj,"parameter5",inst,&(fitlocal->coe[5]));
702 _getobj(obj,"parameter6",inst,&(fitlocal->coe[6]));
703 _getobj(obj,"parameter7",inst,&(fitlocal->coe[7]));
704 _getobj(obj,"parameter8",inst,&(fitlocal->coe[8]));
705 _getobj(obj,"parameter9",inst,&(fitlocal->coe[9]));
706 _getobj(obj,"display",inst,&disp);
707
708 if (through && (type==4)) {
709 error(obj,ERRTHROUGH);
710 return 1;
711 }
712 darray=(struct narray *)(argv[2]);
713 if (arraynum(darray)<1) return -1;
714 data=arraydata(darray);
715 anum=arraynum(darray)-1;
716 dnum=nround(data[0]);
717 data+=1;
718 if (dnum==(anum/2)) weight=FALSE;
719 else if (dnum==(anum/3)) {
720 weight=TRUE;
721 wdata=data+2*dnum;
722 } else return 1;
723 num=0;
724 err2=err3=FALSE;
725 for (i=0;i<dnum;i++) {
726 x=data[i*2];
727 y=data[i*2+1];
728 if (weight) wt=wdata[i];
729 err=FALSE;
730 if (type==1) {
731 if (y<=0) err=TRUE;
732 else y=log(y);
733 if (x<=0) err=TRUE;
734 else x=log(x);
735 } else if (type==2) {
736 if (y<=0) err=TRUE;
737 else y=log(y);
738 } else if (type==3) {
739 if (x<=0) err=TRUE;
740 else x=log(x);
741 }
742 if (err) err2=TRUE;
743 else if (weight && (wt<=0)) {
744 err=TRUE;
745 err3=TRUE;
746 }
747 if (!err) {
748 data[num*2]=x;
749 data[num*2+1]=y;
750 if (weight) wdata[num]=wt;
751 num++;
752 }
753 }
754 if (err2) error(obj,ERRLN);
755 if (err3) error(obj,ERRNEGATIVEWEIGHT);
756 if (through) {
757 err=FALSE;
758 if (type==1) {
759 if (y0<=0) err=TRUE;
760 else y0=log(y0);
761 if (x0<=0) err=TRUE;
762 else x0=log(x0);
763 } else if (type==2) {
764 if (y0<=0) err=TRUE;
765 else y0=log(y0);
766 } else if (type==3) {
767 if (x0<=0) err=TRUE;
768 else x0=log(x0);
769 }
770 if (err) {
771 error(obj,ERRPOINT);
772 return 1;
773 }
774 }
775 if (type!=4)
776 rcode=fitpoly(fitlocal,type,dimension,through,x0,y0,data,num,disp,weight,wdata);
777 else
778 rcode=fituser(obj,fitlocal,func,deriv,converge,data,num,disp,weight,wdata);
779 if (rcode==1) return 1;
780 else if (rcode==-1) ecode=ERRSMLDATA;
781 else if (rcode==-2) ecode=ERRSINGULAR;
782 else if (rcode==-3) ecode=ERRNOEQS;
783 else if (rcode==-4) ecode=ERRSYNTAX;
784 else if (rcode==-5) ecode=ERRRANGE;
785 else if (rcode==-6) ecode=ERRCONVERGE;
786 else ecode=0;
787 if (ecode!=0) {
788 error(obj,ecode);
789 return 1;
790 }
791 if (_putobj(obj,"number",inst,&(fitlocal->num))) return 1;
792 if (_putobj(obj,"error",inst,&(fitlocal->derror))) return 1;
793 if (_putobj(obj,"correlation",inst,&(fitlocal->correlation))) return 1;
794 if (_putobj(obj,"%00",inst,&(fitlocal->coe[0]))) return 1;
795 if (_putobj(obj,"%01",inst,&(fitlocal->coe[1]))) return 1;
796 if (_putobj(obj,"%02",inst,&(fitlocal->coe[2]))) return 1;
797 if (_putobj(obj,"%03",inst,&(fitlocal->coe[3]))) return 1;
798 if (_putobj(obj,"%04",inst,&(fitlocal->coe[4]))) return 1;
799 if (_putobj(obj,"%05",inst,&(fitlocal->coe[5]))) return 1;
800 if (_putobj(obj,"%06",inst,&(fitlocal->coe[6]))) return 1;
801 if (_putobj(obj,"%07",inst,&(fitlocal->coe[7]))) return 1;
802 if (_putobj(obj,"%08",inst,&(fitlocal->coe[8]))) return 1;
803 if (_putobj(obj,"%09",inst,&(fitlocal->coe[9]))) return 1;
804 _getobj(obj,"equation",inst,&equation);
805 if (_putobj(obj,"equation",inst,fitlocal->equation)) return 1;
806 memfree(equation);
807 fitlocal->equation=NULL;
808 return 0;
809 }
810
811
812
813 #define TBLNUM 54
814
815 struct objtable fit[TBLNUM] = {
816 {"init",NVFUNC,NEXEC,fitinit,NULL,0},
817 {"done",NVFUNC,NEXEC,fitdone,NULL,0},
818 {"next",NPOINTER,0,NULL,NULL,0},
819
820 {"profile",NSTR,NREAD|NWRITE,NULL,NULL,0},
821
822 {"type",NENUM,NREAD|NWRITE,fitput,fittypechar,0},
823 {"min",NDOUBLE,NREAD|NWRITE,NULL,NULL,0},
824 {"max",NDOUBLE,NREAD|NWRITE,NULL,NULL,0},
825 {"div",NINT,NREAD|NWRITE,oputge1,NULL,0},
826 {"interpolation",NBOOL,NREAD|NWRITE,NULL,NULL,0},
827 {"through_point",NBOOL,NREAD|NWRITE,fitput,NULL,0},
828 {"point_x",NDOUBLE,NREAD|NWRITE,fitput,NULL,0},
829 {"point_y",NDOUBLE,NREAD|NWRITE,fitput,NULL,0},
830 {"equation",NSTR,NREAD|NWRITE,NULL,NULL,0},
831
832 {"poly_dimension",NINT,NREAD|NWRITE,fitput,NULL,0},
833
834 {"weight_func",NSTR,NREAD|NWRITE,NULL,NULL,0},
835 {"user_func",NSTR,NREAD|NWRITE,fitput,NULL,0},
836 {"derivative",NBOOL,NREAD|NWRITE,fitput,NULL,0},
837 {"derivative0",NSTR,NREAD|NWRITE,fitput,NULL,0},
838 {"derivative1",NSTR,NREAD|NWRITE,fitput,NULL,0},
839 {"derivative2",NSTR,NREAD|NWRITE,fitput,NULL,0},
840 {"derivative3",NSTR,NREAD|NWRITE,fitput,NULL,0},
841 {"derivative4",NSTR,NREAD|NWRITE,fitput,NULL,0},
842 {"derivative5",NSTR,NREAD|NWRITE,fitput,NULL,0},
843 {"derivative6",NSTR,NREAD|NWRITE,fitput,NULL,0},
844 {"derivative7",NSTR,NREAD|NWRITE,fitput,NULL,0},
845 {"derivative8",NSTR,NREAD|NWRITE,fitput,NULL,0},
846 {"derivative9",NSTR,NREAD|NWRITE,fitput,NULL,0},
847 {"converge",NDOUBLE,NREAD|NWRITE,fitput,NULL,0},
848 {"parameter0",NDOUBLE,NREAD|NWRITE,fitput,NULL,0},
849 {"parameter1",NDOUBLE,NREAD|NWRITE,fitput,NULL,0},
850 {"parameter2",NDOUBLE,NREAD|NWRITE,fitput,NULL,0},
851 {"parameter3",NDOUBLE,NREAD|NWRITE,fitput,NULL,0},
852 {"parameter4",NDOUBLE,NREAD|NWRITE,fitput,NULL,0},
853 {"parameter5",NDOUBLE,NREAD|NWRITE,fitput,NULL,0},
854 {"parameter6",NDOUBLE,NREAD|NWRITE,fitput,NULL,0},
855 {"parameter7",NDOUBLE,NREAD|NWRITE,fitput,NULL,0},
856 {"parameter8",NDOUBLE,NREAD|NWRITE,fitput,NULL,0},
857 {"parameter9",NDOUBLE,NREAD|NWRITE,fitput,NULL,0},
858 {"%00",NDOUBLE,NREAD,NULL,NULL,0},
859 {"%01",NDOUBLE,NREAD,NULL,NULL,0},
860 {"%02",NDOUBLE,NREAD,NULL,NULL,0},
861 {"%03",NDOUBLE,NREAD,NULL,NULL,0},
862 {"%04",NDOUBLE,NREAD,NULL,NULL,0},
863 {"%05",NDOUBLE,NREAD,NULL,NULL,0},
864 {"%06",NDOUBLE,NREAD,NULL,NULL,0},
865 {"%07",NDOUBLE,NREAD,NULL,NULL,0},
866 {"%08",NDOUBLE,NREAD,NULL,NULL,0},
867 {"%09",NDOUBLE,NREAD,NULL,NULL,0},
868 {"number",NINT,NREAD,NULL,NULL,0},
869 {"error",NDOUBLE,NREAD,NULL,NULL,0},
870 {"correlation",NDOUBLE,NREAD,NULL,NULL,0},
871 {"display",NBOOL,NREAD|NWRITE,NULL,NULL,0},
872
873 {"fit",NVFUNC,NREAD|NEXEC,fitfit,"da",0},
874 {"_local",NPOINTER,0,NULL,NULL,0},
875 };
876
addfit()877 void *addfit()
878 /* addfit() returns NULL on error */
879 {
880 return addobject(NAME,NULL,PARENT,VERSION,TBLNUM,fit,ERRNUM,fiterrorlist,NULL,NULL);
881 }
882