1 /***************************************************************************
2 * fit.cpp is part of Math Graphic Library
3 * Copyright (C) 2007-2016 Alexey Balakin <mathgl.abalakin@gmail.ru> *
4 * *
5 * This program is free software; you can redistribute it and/or modify *
6 * it under the terms of the GNU Lesser General Public License as *
7 * published by the Free Software Foundation; either version 3 of the *
8 * License, or (at your option) any later version. *
9 * *
10 * This program is distributed in the hope that it will be useful, *
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
13 * GNU General Public License for more details. *
14 * *
15 * You should have received a copy of the GNU Lesser General Public *
16 * License along with this program; if not, write to the *
17 * Free Software Foundation, Inc., *
18 * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
19 ***************************************************************************/
20 #include <ctype.h>
21 #include "mgl2/fit.h"
22 #include "mgl2/prim.h"
23 #include "mgl2/eval.h"
24 #include "mgl2/data.h"
25 #include "mgl2/base.h"
26
27 #if MGL_HAVE_GSL
28 #include <gsl/gsl_multifit_nlin.h>
29 #include <gsl/gsl_blas.h>
30 #endif
31 //-----------------------------------------------------------------------------
32 int mglFitPnts=100; ///< Number of output points in fitting
33 char mglFitRes[1024]; ///< Last fitted formula
34 mreal mglFitChi=NAN; ///< Chi value for last fitted formula
35 mglData mglFitCovar; ///< Covar matrix for lat fitted formula
36 //-----------------------------------------------------------------------------
mgl_get_fit_chi()37 mreal MGL_EXPORT_PURE mgl_get_fit_chi() { return mglFitChi; }
mgl_get_fit_chi_()38 mreal MGL_EXPORT_PURE mgl_get_fit_chi_() { return mglFitChi; }
39 //-----------------------------------------------------------------------------
mgl_get_fit_covar()40 HCDT MGL_EXPORT_CONST mgl_get_fit_covar() { return &mglFitCovar; }
mgl_get_fit_covar_()41 uintptr_t MGL_EXPORT_CONST mgl_get_fit_covar_() { return (uintptr_t)&mglFitCovar; }
42 //-----------------------------------------------------------------------------
mgl_puts_fit(HMGL gr,double x,double y,double z,const char * pre,const char * font,double size)43 void MGL_EXPORT mgl_puts_fit(HMGL gr, double x, double y, double z, const char *pre, const char *font, double size)
44 {
45 long n = strlen(mglFitRes)+(pre?strlen(pre):0)+1;
46 char *buf = new char[n];
47 if(pre) snprintf(buf,n,"%s%s",pre,mglFitRes);
48 else mgl_strncpy(buf,mglFitRes,n);
49 buf[n-1]=0; mgl_puts(gr,x,y,z,buf,font,size);
50 delete []buf;
51 }
mgl_puts_fit_(uintptr_t * gr,mreal * x,mreal * y,mreal * z,const char * prefix,const char * font,mreal * size,int l,int n)52 void MGL_EXPORT mgl_puts_fit_(uintptr_t* gr, mreal *x, mreal *y, mreal *z, const char *prefix, const char *font, mreal *size, int l, int n)
53 {
54 char *s=new char[l+1]; memcpy(s,prefix,l); s[l]=0;
55 char *d=new char[n+1]; memcpy(d,font,n); d[n]=0;
56 mgl_puts_fit(_GR_, *x,*y,*z, s, d, *size);
57 delete []s; delete []d;
58 }
59 //-----------------------------------------------------------------------------
60 /// Structure for keeping data and precompiled fitted formula
61 struct mglFitData
62 {
63 long n; ///< number of points
64 mglDataA *x,*y,*z; ///< x, y, z values
65 mreal *a; ///< function values
66 mreal *s; ///< value dispersions (sigma)
67 const char *eq; ///< approximation formula
68 int m; ///< number of variables
69 const char *var; ///< variables for fitting
70 };
71 //-----------------------------------------------------------------------------
72 #if MGL_HAVE_GSL
mgl_fit__f(const gsl_vector * x,void * data,gsl_vector * f)73 int mgl_fit__f (const gsl_vector *x, void *data, gsl_vector *f)
74 {
75 mglFitData *fd = (mglFitData *)data;
76 mglDataV *var = new mglDataV[fd->m];
77 std::vector<mglDataA*> list;
78 for(long i=0;i<fd->m;i++)
79 { var[i].s = fd->var[i]; var[i].Fill(gsl_vector_get(x,i)); list.push_back(var+i); }
80 if(fd->x) list.push_back(fd->x);
81 if(fd->y) list.push_back(fd->y);
82 if(fd->z) list.push_back(fd->z);
83 HMDT res = mglFormulaCalc(fd->eq, list);
84 #pragma omp parallel for
85 for(long i=0;i<fd->n;i++)
86 {
87 mreal aa = fd->a[i], ss = fd->s[i];
88 if(mgl_isnum(aa) && ss==ss && ss!=0)
89 gsl_vector_set (f, i, (res->a[i] - aa)/ss);
90 else gsl_vector_set (f, i, 0);
91 }
92 delete []var; mgl_delete_data(res);
93 return GSL_SUCCESS;
94 }
95 //-----------------------------------------------------------------------------
mgl_fit__df(const gsl_vector * x,void * data,gsl_matrix * J)96 int static mgl_fit__df (const gsl_vector * x, void *data, gsl_matrix * J)
97 {
98 mglFitData *fd = (mglFitData *)data;
99 mglDataV *var = new mglDataV[fd->m];
100 std::vector<mglDataA*> list;
101 for(long i=0;i<fd->m;i++)
102 { var[i].s = fd->var[i]; var[i].Fill(gsl_vector_get(x,i)); list.push_back(var+i); }
103 if(fd->x) list.push_back(fd->x);
104 if(fd->y) list.push_back(fd->y);
105 if(fd->z) list.push_back(fd->z);
106 HMDT res = mglFormulaCalc(fd->eq, list);
107 const mreal eps = 1e-5;
108 for(long j=0;j<fd->m;j++)
109 {
110 var[j].Fill(gsl_vector_get(x,j)+eps);
111 HMDT dif = mglFormulaCalc(fd->eq, list);
112 var[j].Fill(gsl_vector_get(x,j));
113 #pragma omp parallel for
114 for(long i=0;i<fd->n;i++)
115 {
116 mreal aa = fd->a[i], ss = fd->s[i];
117 if(mgl_isnum(aa) && ss==ss && ss!=0)
118 gsl_matrix_set (J, i, j, (dif->a[i]-res->a[i])/(eps*ss));
119 else gsl_matrix_set (J, i, j, 0);
120 }
121 mgl_delete_data(dif);
122 }
123 delete []var; mgl_delete_data(res);
124 return GSL_SUCCESS;
125 }
126 //-----------------------------------------------------------------------------
mgl_fit__fdf(const gsl_vector * x,void * data,gsl_vector * f,gsl_matrix * J)127 int static mgl_fit__fdf (const gsl_vector * x, void *data, gsl_vector * f, gsl_matrix * J)
128 {
129 mglFitData *fd = (mglFitData *)data;
130 mglDataV *var = new mglDataV[fd->m];
131 std::vector<mglDataA*> list;
132 for(long i=0;i<fd->m;i++)
133 { var[i].s = fd->var[i]; var[i].Fill(gsl_vector_get(x,i)); list.push_back(var+i); }
134 if(fd->x) list.push_back(fd->x);
135 if(fd->y) list.push_back(fd->y);
136 if(fd->z) list.push_back(fd->z);
137 HMDT res = mglFormulaCalc(fd->eq, list);
138 #pragma omp parallel for
139 for(long i=0;i<fd->n;i++)
140 {
141 mreal aa = fd->a[i], ss = fd->s[i];
142 if(mgl_isnum(aa) && ss==ss && ss!=0)
143 gsl_vector_set (f, i, (res->a[i] - aa)/ss);
144 else gsl_vector_set (f, i, 0);
145 }
146 const mreal eps = 1e-5;
147 for(long j=0;j<fd->m;j++)
148 {
149 var[j].Fill(gsl_vector_get(x,j)+eps);
150 HMDT dif = mglFormulaCalc(fd->eq, list);
151 var[j].Fill(gsl_vector_get(x,j));
152 #pragma omp parallel for
153 for(long i=0;i<fd->n;i++)
154 {
155 mreal aa = fd->a[i], ss = fd->s[i];
156 if(mgl_isnum(aa) && ss==ss && ss!=0)
157 gsl_matrix_set (J, i, j, (dif->a[i]-res->a[i])/(eps*ss));
158 else gsl_matrix_set (J, i, j, 0);
159 }
160 mgl_delete_data(dif);
161 }
162 delete []var; mgl_delete_data(res);
163 return GSL_SUCCESS;
164 }
165 #endif
166 //-----------------------------------------------------------------------------
167 /// GSL based fitting procedure for formula/arguments specified by string
mgl_fit_base(mglFitData & fd,mreal * ini)168 mreal static mgl_fit_base(mglFitData &fd, mreal *ini)
169 {
170 #if MGL_HAVE_GSL
171 long m=fd.m,n=fd.n,iter=0;
172 if(n<1 || ini==0) return -1;
173 // setup data
174 double *x_init = new double[fd.m];
175 for(long i=0;i<m;i++) x_init[i] = ini[i];
176 // setup fitting
177 gsl_vector_view vx = gsl_vector_view_array(x_init, m);
178 const gsl_multifit_fdfsolver_type *T = gsl_multifit_fdfsolver_lmsder;
179 gsl_multifit_fdfsolver *s = gsl_multifit_fdfsolver_alloc(T, n, m);
180 gsl_multifit_function_fdf f;
181 f.f = mgl_fit__f; f.df = mgl_fit__df;
182 f.fdf = mgl_fit__fdf; f.n = n; f.p = m;
183 f.params = &fd;
184 gsl_multifit_fdfsolver_set(s, &f, &vx.vector);
185 int status; // start fitting
186 do
187 {
188 iter++;
189 status = gsl_multifit_fdfsolver_iterate(s);
190 if ( status ) break;
191 status = gsl_multifit_test_delta (s->dx, s->x, 1e-4, 1e-4 );
192 }
193 while ( status == GSL_CONTINUE && iter < 500 );
194
195 gsl_matrix *covar = gsl_matrix_alloc(m, m);
196 #ifdef MGL_HAVE_GSL2
197 gsl_matrix *J = gsl_matrix_alloc(s->fdf->n, s->fdf->p);
198 gsl_multifit_fdfsolver_jac(s, J);
199 gsl_multifit_covar (J, 0.0, covar);
200 gsl_matrix_free (J);
201 #else
202 gsl_multifit_covar(s->J, 0.0, covar);
203 #endif
204 mglFitCovar.Set(covar);
205 gsl_matrix_free(covar);
206
207 mreal res = gsl_blas_dnrm2(s->f);
208 for(long i=0;i<m;i++) ini[i] = gsl_vector_get(s->x, i);
209 // free memory
210 gsl_multifit_fdfsolver_free(s);
211 delete []x_init;
212 return res;
213 #else
214 return 0.0;
215 #endif
216 }
217 //-----------------------------------------------------------------------------
mglPrepareFitEq(mglBase * gr,mreal chi,const char * eq,const char * var,mreal * par)218 void mglPrepareFitEq(mglBase *gr,mreal chi, const char *eq, const char *var, mreal *par)
219 {
220 char buf[32]="";
221 mglFitChi = chi;
222 snprintf(mglFitRes,1024,"chi=%g",chi); mglFitRes[1023]=0;
223 size_t i,k,len=strlen(var);
224 for(i=0;i<len;i++)
225 {
226 snprintf(buf,32,", %c=%g",var[i],par[i]);
227 buf[31]=0; strcat(mglFitRes,buf);
228 }
229 gr->SetWarn(-1,mglFitRes);
230
231 memset(mglFitRes, 0, 1024); //mglFitRes[0] = 0;
232 len=strlen(eq);
233 for(i=k=0;i<len;i++)
234 {
235 const char *c = strchr(var,eq[i]);
236 if(c && (i==0 || !isalnum(eq[i-1])) && (i==len-1 || !isalnum(eq[i+1])))
237 {
238 snprintf(buf,32,"%g",par[c-var]);
239 buf[31]=0; strcat(mglFitRes+k, buf); k+=strlen(buf);
240 }
241 else { mglFitRes[k] = eq[i]; k++; }
242 }
243 mglFitRes[k]=0;
244 }
245 //-----------------------------------------------------------------------------
mgl_fit_1(HMGL gr,HCDT y,const char * eq,const char * var,HMDT ini,const char * opt)246 HMDT MGL_EXPORT mgl_fit_1(HMGL gr, HCDT y, const char *eq, const char *var, HMDT ini, const char *opt)
247 {
248 gr->SaveState(opt);
249 mglData x(y->GetNx()); x.Fill(gr->Min.x, gr->Max.x);
250 mglData s(y); s.Fill(1,1);
251 return mgl_fit_xys(gr,&x,y,&s,eq,var,ini,0);
252 }
253 //-----------------------------------------------------------------------------
mgl_fit_2(HMGL gr,HCDT z,const char * eq,const char * var,HMDT ini,const char * opt)254 HMDT MGL_EXPORT mgl_fit_2(HMGL gr, HCDT z, const char *eq, const char *var, HMDT ini, const char *opt)
255 {
256 gr->SaveState(opt);
257 mglData x(z->GetNx()); x.Fill(gr->Min.x, gr->Max.x);
258 mglData y(z->GetNy()); y.Fill(gr->Min.y, gr->Max.y);
259 mglData s(z); s.Fill(1,1);
260 return mgl_fit_xyzs(gr,&x,&y,z,&s,eq,var,ini,0);
261 }
262 //-----------------------------------------------------------------------------
mgl_fit_3(HMGL gr,HCDT a,const char * eq,const char * var,HMDT ini,const char * opt)263 HMDT MGL_EXPORT mgl_fit_3(HMGL gr, HCDT a, const char *eq, const char *var, HMDT ini, const char *opt)
264 {
265 gr->SaveState(opt);
266 mglData x(a->GetNx()); x.Fill(gr->Min.x, gr->Max.x);
267 mglData y(a->GetNy()); y.Fill(gr->Min.y, gr->Max.y);
268 mglData z(a->GetNz()); z.Fill(gr->Min.z, gr->Max.z);
269 mglData s(a); s.Fill(1,1);
270 return mgl_fit_xyzas(gr,&x,&y,&z,a,&s,eq,var,ini,0);
271 }
272 //-----------------------------------------------------------------------------
mgl_fit_xy(HMGL gr,HCDT x,HCDT y,const char * eq,const char * var,HMDT ini,const char * opt)273 HMDT MGL_EXPORT mgl_fit_xy(HMGL gr, HCDT x, HCDT y, const char *eq, const char *var, HMDT ini, const char *opt)
274 {
275 mglData s(y); s.Fill(1,1);
276 return mgl_fit_xys(gr,x,y,&s,eq,var,ini,opt);
277 }
278 //-----------------------------------------------------------------------------
mgl_fit_xyz(HMGL gr,HCDT x,HCDT y,HCDT z,const char * eq,const char * var,HMDT ini,const char * opt)279 HMDT MGL_EXPORT mgl_fit_xyz(HMGL gr, HCDT x, HCDT y, HCDT z, const char *eq, const char *var, HMDT ini, const char *opt)
280 {
281 mglData s(z); s.Fill(1,1);
282 return mgl_fit_xyzs(gr,x,y,z,&s,eq,var,ini,opt);
283 }
284 //-----------------------------------------------------------------------------
mgl_fit_xyza(HMGL gr,HCDT x,HCDT y,HCDT z,HCDT a,const char * eq,const char * var,HMDT ini,const char * opt)285 HMDT MGL_EXPORT mgl_fit_xyza(HMGL gr, HCDT x, HCDT y, HCDT z, HCDT a, const char *eq, const char *var, HMDT ini, const char *opt)
286 {
287 mglData s(a); s.Fill(1,1);
288 return mgl_fit_xyzas(gr,x,y,z,a,&s,eq,var,ini,opt);
289 }
290 //-----------------------------------------------------------------------------
mgl_fit_ys(HMGL gr,HCDT y,HCDT s,const char * eq,const char * var,HMDT ini,const char * opt)291 HMDT MGL_EXPORT mgl_fit_ys(HMGL gr, HCDT y, HCDT s, const char *eq, const char *var, HMDT ini, const char *opt)
292 {
293 gr->SaveState(opt);
294 mglData x(y->GetNx()); x.Fill(gr->Min.x, gr->Max.x);
295 return mgl_fit_xys(gr,&x,y,s,eq,var,ini,0);
296 }
297 //-----------------------------------------------------------------------------
mgl_fill_fit(HMGL gr,mglData & fit,mglData & in,mglFitData & fd,const char * var,long nx,long ny,long nz,long k)298 void static mgl_fill_fit(HMGL gr, mglData &fit, mglData &in, mglFitData &fd, const char *var, long nx, long ny, long nz, long k)
299 {
300 mglDataV *vv = new mglDataV[fd.m];
301 std::vector<mglDataA*> list;
302 for(long i=0;i<fd.m;i++)
303 { vv[i].s = var[i]; vv[i].Fill(in.a[i]); list.push_back(vv+i); }
304 mglDataV x(nx,ny,nz, gr->Min.x,gr->Max.x,'x'); x.Name(L"x"); list.push_back(&x);
305 mglDataV y(nx,ny,nz, gr->Min.y,gr->Max.y,'y'); y.Name(L"y"); list.push_back(&y);
306 mglDataV z(nx,ny,nz, gr->Min.z,gr->Max.z,'z'); z.Name(L"z"); list.push_back(&z);
307 HMDT res = mglFormulaCalc(fd.eq, list);
308 long nn = nx*ny*nz;
309 memcpy(fit.a+k*nn,res->a,nn*sizeof(mreal));
310 delete []vv; mgl_delete_data(res);
311 }
312 //-----------------------------------------------------------------------------
mgl_fit_xys(HMGL gr,HCDT xx,HCDT yy,HCDT ss,const char * eq,const char * var,HMDT ini,const char * opt)313 HMDT MGL_EXPORT mgl_fit_xys(HMGL gr, HCDT xx, HCDT yy, HCDT ss, const char *eq, const char *var, HMDT ini, const char *opt)
314 {
315 long m = yy->GetNx();
316 mreal rr = gr->SaveState(opt);
317 long nn = (mgl_isnan(rr) || rr<=0) ? mglFitPnts:long(rr+0.5);
318 if(xx->GetNx()!=m)
319 { gr->SetWarn(mglWarnDim,"Fit[S]"); return 0; }
320 if(m<2)
321 { gr->SetWarn(mglWarnLow,"Fit[S]"); return 0; }
322 if(ss->GetNN() != yy->GetNN())
323 { gr->SetWarn(mglWarnDim,"Fit[S]"); return 0; }
324 if(!var || *var==0)
325 { gr->SetWarn(mglWarnNull,"Fit[S]"); return 0; }
326
327 mglData x(xx), y(yy), s(ss); x.Name(L"x");
328 long mm = yy->GetNy()*yy->GetNz();
329 #pragma omp parallel for
330 for(long i=0;i<m;i++) if(mgl_isnan(x.a[i]))
331 for(long j=0;j<mm;j++) y.a[i+m*j] = NAN;
332 mglFitData fd;
333 fd.n = m; fd.x = &x; fd.y = 0;
334 fd.z = 0; fd.a = y.a; fd.s = s.a;
335 fd.eq = eq; fd.var = var; fd.m = strlen(var);
336 mglData in(fd.m), *fit=new mglData(nn, yy->GetNy(), yy->GetNz());
337 mreal res=-1;
338 mglDataR xc(x);
339 for(long i=0;i<yy->GetNy()*yy->GetNz();i++)
340 {
341 if(ini && ini->nx>=fd.m) in.Set(ini->a,fd.m);
342 else in.Fill(0.,0);
343 xc.SetInd(i%x.ny, L"x");
344 fd.a = y.a+i*m; fd.x = &xc; //x.a+(i%x.ny)*m;
345 fd.s = s.a+i*m;
346 res = mgl_fit_base(fd,in.a);
347 mgl_fill_fit(gr,*fit,in,fd,var,nn,1,1,i);
348 if(ini && ini->nx>=fd.m) memcpy(ini->a,in.a,fd.m*sizeof(mreal));
349 }
350 mglPrepareFitEq(gr,res,eq,var,in.a);
351 gr->LoadState(); return fit;
352 }
353 //-----------------------------------------------------------------------------
mgl_fit_xyzs(HMGL gr,HCDT xx,HCDT yy,HCDT zz,HCDT ss,const char * eq,const char * var,HMDT ini,const char * opt)354 HMDT MGL_EXPORT mgl_fit_xyzs(HMGL gr, HCDT xx, HCDT yy, HCDT zz, HCDT ss, const char *eq, const char *var, HMDT ini, const char *opt)
355 {
356 long m=zz->GetNx(),n=zz->GetNy();
357 mreal rr = gr->SaveState(opt);
358 long nn = (mgl_isnan(rr) || rr<=0) ? mglFitPnts:long(rr+0.5);
359 if(xx->GetNx()!=m)
360 { gr->SetWarn(mglWarnDim,"Fit[S]"); return 0; }
361 if(ss->GetNN() != zz->GetNN())
362 { gr->SetWarn(mglWarnDim,"Fit[S]"); return 0; }
363 if(yy->GetNx()!=n && (xx->GetNy()!=n || yy->GetNx()!=m || yy->GetNy()!=n))
364 { gr->SetWarn(mglWarnDim,"Fit[S]"); return 0; }
365 if(m<2|| n<2)
366 { gr->SetWarn(mglWarnLow,"Fit[S]"); return 0; }
367 if(!var || *var==0)
368 { gr->SetWarn(mglWarnNull,"Fit[S]"); return 0; }
369
370 mglData x(m, n), y(m, n), z(zz), s(ss); x.Name(L"x"); y.Name(L"y");
371 long nz = zz->GetNz(), mm = n*m;
372 #pragma omp parallel for collapse(2)
373 for(long j=0;j<n;j++) for(long i=0;i<m;i++)
374 {
375 long i0 = i+m*j;
376 x.a[i0] = GetX(xx,i,j,0).x;
377 y.a[i0] = GetY(yy,i,j,0).x;
378 if(mgl_isnan(x.a[i0]) || mgl_isnan(y.a[i0]))
379 for(long k=0;k<nz;k++) z.a[i0+mm*k] = NAN;
380 }
381 mglFitData fd;
382 fd.n = m*n; fd.x = &x; fd.y = &y;
383 fd.z = 0; fd.a = z.a; fd.s = s.a;
384 fd.eq = eq; fd.var=var; fd.m = strlen(var);
385
386 mglData in(fd.m), *fit=new mglData(nn, nn, zz->GetNz());
387 mreal res = -1;
388 for(long i=0;i<nz;i++)
389 {
390 if(ini && ini->nx>=fd.m) in.Set(ini->a,fd.m);
391 else in.Fill(0.,0);
392 fd.a = z.a+i*m*n; fd.s = s.a+i*m*n;
393 res = mgl_fit_base(fd,in.a);
394 mgl_fill_fit(gr,*fit,in,fd,var,nn,nn,1,i);
395 if(ini && ini->nx>=fd.m) memcpy(ini->a,in.a,fd.m*sizeof(mreal));
396 }
397 mglPrepareFitEq(gr,res, eq,var,in.a);
398 gr->LoadState(); return fit;
399 }
400 //-----------------------------------------------------------------------------
mgl_fit_xyzas(HMGL gr,HCDT xx,HCDT yy,HCDT zz,HCDT aa,HCDT ss,const char * eq,const char * var,HMDT ini,const char * opt)401 HMDT MGL_EXPORT mgl_fit_xyzas(HMGL gr, HCDT xx, HCDT yy, HCDT zz, HCDT aa, HCDT ss, const char *eq, const char *var, HMDT ini, const char *opt)
402 {
403 long m=aa->GetNx(), n=aa->GetNy(), l=aa->GetNz(), i = n*m*l;
404 mreal rr = gr->SaveState(opt);
405 long nn = (mgl_isnan(rr) || rr<=0) ? mglFitPnts:long(rr+0.5);
406 if(m<2 || n<2 || l<2)
407 { gr->SetWarn(mglWarnLow,"Fit[S]"); return 0; }
408 if(ss->GetNN() != i)
409 { gr->SetWarn(mglWarnDim,"Fit[S]"); return 0; }
410 bool both = xx->GetNN()==i && yy->GetNN()==i && zz->GetNN()==i;
411 if(!(both || (xx->GetNx()==m && yy->GetNx()==n && zz->GetNx()==l)))
412 { gr->SetWarn(mglWarnDim,"Fit[S]"); return 0; }
413 if(!var || *var==0)
414 { gr->SetWarn(mglWarnNull,"Fit[S]"); return 0; }
415
416 mglData x(m,n,l), y(m,n,l), z(m,n,l), a(aa), s(ss);
417 x.Name(L"x"); y.Name(L"y"); z.Name(L"z");
418 #pragma omp parallel for collapse(3)
419 for(long k=0;k<l;k++) for(long j=0;j<n;j++) for(long i=0;i<m;i++)
420 {
421 long i0 = i+m*(j+n*k);
422 x.a[i0] = GetX(xx,i,j,k).x;
423 y.a[i0] = GetY(yy,i,j,k).x;
424 z.a[i0] = GetZ(zz,i,j,k).x;
425 if(mgl_isnan(x.a[i0]) || mgl_isnan(y.a[i0]) || mgl_isnan(z.a[i0])) a.a[i0] = NAN;
426 }
427 mglFitData fd;
428 fd.n = m*n*l; fd.x = &x; fd.y = &y;
429 fd.z = &z; fd.a = a.a; fd.s = s.a;
430 fd.eq = eq; fd.var=var; fd.m = strlen(var);
431 mglData in(fd.m), *fit=new mglData(nn, nn, nn);
432 mreal res = -1;
433
434 if(ini && ini->nx>=fd.m) in.Set(ini->a,fd.m);
435 else in.Fill(0.,0);
436 res = mgl_fit_base(fd,in.a);
437 mgl_fill_fit(gr,*fit,in,fd,var,nn,nn,nn,0);
438 if(ini && ini->nx>=fd.m) memcpy(ini->a,in.a,fd.m*sizeof(mreal));
439
440 mglPrepareFitEq(gr,res, eq,var,in.a);
441 gr->LoadState(); return fit;
442 }
443 //-----------------------------------------------------------------------------
mgl_hist_x(HMGL gr,HCDT x,HCDT a,const char * opt)444 HMDT MGL_EXPORT mgl_hist_x(HMGL gr, HCDT x, HCDT a, const char *opt)
445 {
446 long nn=a->GetNN();
447 if(nn!=x->GetNN())
448 { gr->SetWarn(mglWarnDim,"Hist"); return 0; }
449 mreal rr = gr->SaveState(opt);
450 long n = (mgl_isnan(rr) || rr<=0) ? mglFitPnts:long(rr+0.5);
451 mglData *res = new mglData(n);
452
453 mreal vx = n/(gr->Max.x-gr->Min.x);
454 for(long i=0;i<nn;i++)
455 {
456 long j1 = long((x->vthr(i)-gr->Min.x)*vx);
457 if(j1>=0 && j1<n) res->a[j1] += a->vthr(i);
458 }
459 gr->LoadState(); return res;
460 }
461 //-----------------------------------------------------------------------------
mgl_hist_xy(HMGL gr,HCDT x,HCDT y,HCDT a,const char * opt)462 HMDT MGL_EXPORT mgl_hist_xy(HMGL gr, HCDT x, HCDT y, HCDT a, const char *opt)
463 {
464 long nn=a->GetNN();
465 if(nn!=x->GetNN() || nn!=y->GetNN())
466 { gr->SetWarn(mglWarnDim,"Hist"); return 0; }
467 mreal rr = gr->SaveState(opt);
468 long n = (mgl_isnan(rr) || rr<=0) ? mglFitPnts:long(rr+0.5);
469 mglData *res = new mglData(n, n);
470 mreal vx = n/(gr->Max.x-gr->Min.x);
471 mreal vy = n/(gr->Max.y-gr->Min.y);
472 for(long i=0;i<nn;i++)
473 {
474 long j1 = long((x->vthr(i)-gr->Min.x)*vx);
475 long j2 = long((y->vthr(i)-gr->Min.y)*vy);
476 if(j1>=0 && j1<n && j2>=0 && j2<n) res->a[j1+n*j2] += a->vthr(i);
477 }
478 gr->LoadState(); return res;
479 }
480 //-----------------------------------------------------------------------------
mgl_hist_xyz(HMGL gr,HCDT x,HCDT y,HCDT z,HCDT a,const char * opt)481 HMDT MGL_EXPORT mgl_hist_xyz(HMGL gr, HCDT x, HCDT y, HCDT z, HCDT a, const char *opt)
482 {
483 long nn=a->GetNN();
484 if(nn!=x->GetNN() || nn!=y->GetNN() || nn!=z->GetNN())
485 { gr->SetWarn(mglWarnDim,"Hist"); return 0; }
486 mreal rr = gr->SaveState(opt);
487 long n = (mgl_isnan(rr) || rr<=0) ? mglFitPnts:long(rr+0.5);
488 mglData *res = new mglData(n, n, n);
489 mreal vx = n/(gr->Max.x-gr->Min.x), vy = n/(gr->Max.y-gr->Min.y), vz = n/(gr->Max.z-gr->Min.z);
490 for(long i=0;i<nn;i++)
491 {
492 long j1 = long((x->vthr(i)-gr->Min.x)*vx);
493 long j2 = long((y->vthr(i)-gr->Min.y)*vy);
494 long j3 = long((z->vthr(i)-gr->Min.z)*vz);
495 if(j1>=0 && j1<n && j2>=0 && j2<n && j3>=0 && j3<n)
496 res->a[j1+n*(j2+n*j3)] += a->vthr(i);
497 }
498 gr->LoadState(); return res;
499 }
500 //-----------------------------------------------------------------------------
mgl_hist_x_(uintptr_t * gr,uintptr_t * x,uintptr_t * a,const char * opt,int lo)501 uintptr_t MGL_EXPORT mgl_hist_x_(uintptr_t* gr, uintptr_t* x, uintptr_t* a, const char *opt, int lo)
502 { char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
503 uintptr_t r = (uintptr_t)mgl_hist_x(_GR_, _DA_(x), _DA_(a), o);
504 delete []o; return r; }
mgl_hist_xy_(uintptr_t * gr,uintptr_t * x,uintptr_t * y,uintptr_t * a,const char * opt,int lo)505 uintptr_t MGL_EXPORT mgl_hist_xy_(uintptr_t* gr, uintptr_t* x, uintptr_t* y, uintptr_t* a, const char *opt, int lo)
506 { char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
507 uintptr_t r = (uintptr_t)mgl_hist_xy(_GR_, _DA_(x), _DA_(y), _DA_(a), o);
508 delete []o; return r; }
mgl_hist_xyz_(uintptr_t * gr,uintptr_t * x,uintptr_t * y,uintptr_t * z,uintptr_t * a,const char * opt,int lo)509 uintptr_t MGL_EXPORT mgl_hist_xyz_(uintptr_t* gr, uintptr_t* x, uintptr_t* y, uintptr_t* z, uintptr_t* a, const char *opt, int lo)
510 { char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
511 uintptr_t r = (uintptr_t)mgl_hist_xyz(_GR_, _DA_(x), _DA_(y), _DA_(z), _DA_(a), o);
512 delete []o; return r; }
513 //-----------------------------------------------------------------------------
mgl_get_fit(HMGL)514 MGL_EXPORT_CONST const char *mgl_get_fit(HMGL ) { return mglFitRes; }
mgl_get_fit_(uintptr_t * gr,char * out,int len)515 int MGL_EXPORT mgl_get_fit_(uintptr_t *gr, char *out, int len)
516 {
517 const char *res = mgl_get_fit(_GR_);
518 if(out) mgl_strncpy(out,res,len);
519 return strlen(res);
520 }
521 //-----------------------------------------------------------------------------
mgl_fit_1_(uintptr_t * gr,uintptr_t * y,const char * eq,const char * var,uintptr_t * ini,const char * opt,int l,int n,int lo)522 uintptr_t MGL_EXPORT mgl_fit_1_(uintptr_t* gr, uintptr_t* y, const char *eq, const char *var, uintptr_t *ini, const char *opt, int l, int n, int lo)
523 {
524 char *s=new char[l+1]; memcpy(s,eq,l); s[l]=0;
525 char *d=new char[n+1]; memcpy(d,var,n); d[n]=0;
526 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
527 uintptr_t r = (uintptr_t)mgl_fit_1(_GR_, _DA_(y), s, d, _DM_(ini), o);
528 delete []o; delete []s; delete []d; return r;
529 }
mgl_fit_2_(uintptr_t * gr,uintptr_t * z,const char * eq,const char * var,uintptr_t * ini,const char * opt,int l,int n,int lo)530 uintptr_t MGL_EXPORT mgl_fit_2_(uintptr_t* gr, uintptr_t* z, const char *eq, const char *var, uintptr_t *ini, const char *opt, int l, int n, int lo)
531 {
532 char *s=new char[l+1]; memcpy(s,eq,l); s[l]=0;
533 char *d=new char[n+1]; memcpy(d,var,n); d[n]=0;
534 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
535 uintptr_t r = (uintptr_t)mgl_fit_2(_GR_, _DA_(z), s, d, _DM_(ini), o);
536 delete []o; delete []s; delete []d; return r;
537 }
mgl_fit_3_(uintptr_t * gr,uintptr_t * a,const char * eq,const char * var,uintptr_t * ini,const char * opt,int l,int n,int lo)538 uintptr_t MGL_EXPORT mgl_fit_3_(uintptr_t* gr, uintptr_t* a, const char *eq, const char *var, uintptr_t *ini, const char *opt, int l, int n, int lo)
539 {
540 char *s=new char[l+1]; memcpy(s,eq,l); s[l]=0;
541 char *d=new char[n+1]; memcpy(d,var,n); d[n]=0;
542 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
543 uintptr_t r = (uintptr_t)mgl_fit_3(_GR_, _DA_(a), s, d, _DM_(ini), o);
544 delete []o; delete []s; delete []d; return r;
545 }
mgl_fit_xy_(uintptr_t * gr,uintptr_t * x,uintptr_t * y,const char * eq,const char * var,uintptr_t * ini,const char * opt,int l,int n,int lo)546 uintptr_t MGL_EXPORT mgl_fit_xy_(uintptr_t* gr, uintptr_t* x, uintptr_t* y, const char *eq, const char *var, uintptr_t *ini, const char *opt, int l, int n, int lo)
547 {
548 char *s=new char[l+1]; memcpy(s,eq,l); s[l]=0;
549 char *d=new char[n+1]; memcpy(d,var,n); d[n]=0;
550 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
551 uintptr_t r = (uintptr_t)mgl_fit_xy(_GR_, _DA_(x), _DA_(y), s, d, _DM_(ini), o);
552 delete []o; delete []s; delete []d; return r;
553 }
mgl_fit_xyz_(uintptr_t * gr,uintptr_t * x,uintptr_t * y,uintptr_t * z,const char * eq,const char * var,uintptr_t * ini,const char * opt,int l,int n,int lo)554 uintptr_t MGL_EXPORT mgl_fit_xyz_(uintptr_t* gr, uintptr_t* x, uintptr_t* y, uintptr_t* z, const char *eq, const char *var, uintptr_t *ini, const char *opt, int l, int n, int lo)
555 {
556 char *s=new char[l+1]; memcpy(s,eq,l); s[l]=0;
557 char *d=new char[n+1]; memcpy(d,var,n); d[n]=0;
558 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
559 uintptr_t r = (uintptr_t)mgl_fit_xyz(_GR_, _DA_(x), _DA_(y), _DA_(z), s, d, _DM_(ini), o);
560 delete []o; delete []s; delete []d; return r;
561 }
mgl_fit_xyza_(uintptr_t * gr,uintptr_t * x,uintptr_t * y,uintptr_t * z,uintptr_t * a,const char * eq,const char * var,uintptr_t * ini,const char * opt,int l,int n,int lo)562 uintptr_t MGL_EXPORT mgl_fit_xyza_(uintptr_t* gr, uintptr_t* x, uintptr_t* y, uintptr_t* z, uintptr_t* a, const char *eq, const char *var, uintptr_t *ini, const char *opt, int l, int n, int lo)
563 {
564 char *s=new char[l+1]; memcpy(s,eq,l); s[l]=0;
565 char *d=new char[n+1]; memcpy(d,var,n); d[n]=0;
566 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
567 uintptr_t r = (uintptr_t)mgl_fit_xyza(_GR_, _DA_(x), _DA_(y), _DA_(z), _DA_(a), s, d, _DM_(ini), o);
568 delete []o; delete []s; delete []d; return r;
569 }
mgl_fit_ys_(uintptr_t * gr,uintptr_t * y,uintptr_t * ss,const char * eq,const char * var,uintptr_t * ini,const char * opt,int l,int n,int lo)570 uintptr_t MGL_EXPORT mgl_fit_ys_(uintptr_t* gr, uintptr_t* y, uintptr_t* ss, const char *eq, const char *var, uintptr_t *ini, const char *opt, int l, int n, int lo)
571 {
572 char *s=new char[l+1]; memcpy(s,eq,l); s[l]=0;
573 char *d=new char[n+1]; memcpy(d,var,n); d[n]=0;
574 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
575 uintptr_t r = (uintptr_t)mgl_fit_ys(_GR_, _DA_(y), _DA_(ss), s, d, _DM_(ini), o);
576 delete []o; delete []s; delete []d; return r;
577 }
mgl_fit_xys_(uintptr_t * gr,uintptr_t * x,uintptr_t * y,uintptr_t * ss,const char * eq,const char * var,uintptr_t * ini,const char * opt,int l,int n,int lo)578 uintptr_t MGL_EXPORT mgl_fit_xys_(uintptr_t* gr, uintptr_t* x, uintptr_t* y, uintptr_t* ss, const char *eq, const char *var, uintptr_t *ini, const char *opt, int l, int n, int lo)
579 {
580 char *s=new char[l+1]; memcpy(s,eq,l); s[l]=0;
581 char *d=new char[n+1]; memcpy(d,var,n); d[n]=0;
582 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
583 uintptr_t r = (uintptr_t)mgl_fit_xys(_GR_, _DA_(x), _DA_(y), _DA_(ss), s, d, _DM_(ini), o);
584 delete []o; delete []s; delete []d; return r;
585 }
mgl_fit_xyzs_(uintptr_t * gr,uintptr_t * x,uintptr_t * y,uintptr_t * z,uintptr_t * ss,const char * eq,const char * var,uintptr_t * ini,const char * opt,int l,int n,int lo)586 uintptr_t MGL_EXPORT mgl_fit_xyzs_(uintptr_t* gr, uintptr_t* x, uintptr_t* y, uintptr_t* z, uintptr_t* ss, const char *eq, const char *var, uintptr_t *ini, const char *opt, int l, int n, int lo)
587 {
588 char *s=new char[l+1]; memcpy(s,eq,l); s[l]=0;
589 char *d=new char[n+1]; memcpy(d,var,n); d[n]=0;
590 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
591 uintptr_t r = (uintptr_t)mgl_fit_xyzs(_GR_, _DA_(x), _DA_(y), _DA_(z), _DA_(ss), s, d, _DM_(ini), o);
592 delete []o; delete []s; delete []d; return r;
593 }
mgl_fit_xyzas_(uintptr_t * gr,uintptr_t * x,uintptr_t * y,uintptr_t * z,uintptr_t * a,uintptr_t * ss,const char * eq,const char * var,uintptr_t * ini,const char * opt,int l,int n,int lo)594 uintptr_t MGL_EXPORT mgl_fit_xyzas_(uintptr_t* gr, uintptr_t* x, uintptr_t* y, uintptr_t* z, uintptr_t* a, uintptr_t* ss, const char *eq, const char *var, uintptr_t *ini, const char *opt, int l, int n, int lo)
595 {
596 char *s=new char[l+1]; memcpy(s,eq,l); s[l]=0;
597 char *d=new char[n+1]; memcpy(d,var,n); d[n]=0;
598 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
599 uintptr_t r = (uintptr_t)mgl_fit_xyzas(_GR_, _DA_(x), _DA_(y), _DA_(z), _DA_(a), _DA_(ss), s, d, _DM_(ini), o);
600 delete []o; delete []s; delete []d; return r;
601 }
602 //-----------------------------------------------------------------------------
603