[fityk-devel] new built-in function

[Scholz, Mirko]

Hello,

if you accept new built-in functions, here is code for log-normal shaped
peaks, which I use
to model the ground state absorption of some molecules.

Mirko Scholz

{{{

diff -p -r fityk-0.8.6/src/bfunc.cpp fityk-0.8.6-LogNormal/src/bfunc.cpp
*** fityk-0.8.6/src/bfunc.cpp 2007-04-08 01:16:36.000000000 +0200
--- fityk-0.8.6-LogNormal/src/bfunc.cpp 2009-03-21 15:23:07.000000000 +0100
*************** FUNC_CALCULATE_VALUE_DERIV_BEGIN(Pielasz
*** 927,930 ****
--- 927,1002 ----
dy_dx = dcenter;
FUNC_CALCULATE_VALUE_DERIV_END(height*t);

+ ///////////////////////////////////////////////////////////////////////
+
+ const char *FuncLogNormal::formula
+ = "LogNormal(height, center, bandwidth=fwhm, asymmetry = 0.1) = "
+ "height*exp(-ln(2)*(ln(2*asymmetry*(x-center)/bandwidth+1)/asymmetry)^2)";
+
+ void FuncLogNormal::more_precomputations()
+ {
+ if (vv.size() != 4)
+ vv.resize(4);
+ if (fabs(vv[2]) < epsilon)
+ vv[2] = epsilon;
+ if (fabs(vv[3]) < epsilon)
+ vv[3] = 0.001;
+ }
+
+ FUNC_CALCULATE_VALUE_BEGIN(LogNormal)
+ fp a = 2.0 * vv[3] * (x - vv[1]) / vv[2];
+ fp ex = 0.0;
+ if (a > -1.0) {
+ fp b = log(1 + a) / vv[3];
+ ex = vv[0] * exp_(-M_LN2 * b * b);
+ }
+ FUNC_CALCULATE_VALUE_END(ex)
+
+ FUNC_CALCULATE_VALUE_DERIV_BEGIN(LogNormal)
+ fp a = 2.0 * vv[3] * (x - vv[1]) / vv[2];
+ fp ex;
+ if (a > -1.0) {
+ fp b = log(1 + a) / vv[3];
+ ex = exp_(-M_LN2 * b * b);
+ dy_dv[0] = ex;
+ ex *= vv[0];
+ dy_dv[1] = 4.0*M_LN2/(vv[2]*(a+1))*ex*b;
+ dy_dv[2] = 4.0*(x-vv[1])*M_LN2/(vv[2]*vv[2]*(a+1))*ex*b;
+ dy_dv[3] = ex*(2.0*M_LN2*b*b/vv[3]
+ -4.0*(x-vv[1])*log(a+1)*M_LN2/(vv[2]*vv[3]*vv[3]*(a+1)));
+ dy_dx = -4.0*M_LN2/(vv[2]*(a+1))*ex*b;
+ }
+ else {
+ ex = 0.0;
+ dy_dv[0] = 0.0;
+ dy_dv[1] = 0.0;
+ dy_dv[2] = 0.0;
+ dy_dv[3] = 0.0;
+ dy_dx = 0.0;
+ }
+ FUNC_CALCULATE_VALUE_DERIV_END(ex)
+
+ bool FuncLogNormal::get_nonzero_range (fp level, fp &left, fp &right) const
+ { /* untested */
+ if (level == 0)
+ return false;
+ else if (fabs(level) >= fabs(vv[0]))
+ left = right = 0;
+ else {
+ fp w = sqrt (log (fabs(vv[0]/level)) / M_LN2) * vv[2];
+ fp w1 = (1-exp_(sqrt(log(fabs(vv[0]/level))/M_LN2)*vv[3]))*vv[2]
+ /2.0/vv[3]+vv[1];
+ fp w0 = (1-exp_(-sqrt(log(fabs(vv[0]/level))/M_LN2)*vv[3]))*vv[2]
+ /2.0/vv[3]+vv[1];
+ if (w1>w0) {
+ left = w0;
+ right = w1;
+ }
+ else {
+ left = w1;
+ right = w0;
+ }
+ }
+ return true;
+ }

diff -p -r fityk-0.8.6/src/bfunc.h fityk-0.8.6-LogNormal/src/bfunc.h
*** fityk-0.8.6/src/bfunc.h 2007-07-24 02:17:13.000000000 +0200
--- fityk-0.8.6-LogNormal/src/bfunc.h 2009-03-21 15:45:37.000000000 +0100
*************** class FuncValente : public Function
*** 212,216 ****
--- 212,231 ----
fp center() const { return vv[3]; }
};

+ class FuncLogNormal : public Function
+ {
+ DECLARE_FUNC_OBLIGATORY_METHODS(LogNormal)
+ void more_precomputations();
+ bool get_nonzero_range (fp level, fp &left, fp &right) const;
+ fp center() const { return vv[1]; }
+ bool has_height() const { return true; }
+ fp height() const { return vv[0]; }
+ bool has_fwhm() const { return true; }
+ fp fwhm() const { return fabs(vv[2]); }
+ bool has_area() const { return true; }
+ fp area() const { return vv[0]/sqrt(M_LN2/M_PI)
+ /(2.0/vv[2])/exp(-vv[3]*vv[3]/4.0/M_LN2); }
+ };
+

#endif
diff -p -r fityk-0.8.6/src/func.cpp fityk-0.8.6-LogNormal/src/func.cpp
*** fityk-0.8.6/src/func.cpp 2008-04-15 13:08:44.000000000 +0200
--- fityk-0.8.6-LogNormal/src/func.cpp 2009-03-21 15:37:01.000000000 +0100
*************** Function* Function::factory (Ftk const*
*** 123,128 ****
--- 123,129 ----
FACTORY_FUNC(EMG)
FACTORY_FUNC(DoniachSunjic)
FACTORY_FUNC(PielaszekCube)
+ FACTORY_FUNC(LogNormal)
else if (UdfContainer::is_defined(type_name)) {
UdfContainer::UDF const* udf = UdfContainer::get_udf(type_name);
if (udf->is_compound)
*************** const char* builtin_formulas[] = {
*** 153,159 ****
FuncVoigtA::formula,
FuncEMG::formula,
FuncDoniachSunjic::formula,
! FuncPielaszekCube::formula
};

vector<string> Function::get_all_types()
--- 154,161 ----
FuncVoigtA::formula,
FuncEMG::formula,
FuncDoniachSunjic::formula,
! FuncPielaszekCube::formula,
! FuncLogNormal::formula
};

vector<string> Function::get_all_types()
*************** void initialize_udfs()
*** 535,541 ****
"LorentzianA(area, center, hwhm) = Lorentzian(area/hwhm/pi, center,
hwhm)\n"
"Pearson7A(area, center, hwhm, shape=2) =
Pearson7(area/(hwhm*exp(lgamma(shape-0.5)-lgamma(shape))*sqrt(pi/(2^(1/shape)
-1))), center, hwhm, shape)\n"
"PseudoVoigtA(area, center, hwhm, shape=0.5) = GaussianA(area*(1-shape),
center, hwhm) + LorentzianA(area*shape, center, hwhm)\n"
! "ExpDecay(a=0, t=1) = a*exp(-x/t)",
"\n");
udfs.clear();
for (vector<string>::const_iterator i = formulae.begin();
--- 537,544 ----
"LorentzianA(area, center, hwhm) = Lorentzian(area/hwhm/pi, center,
hwhm)\n"
"Pearson7A(area, center, hwhm, shape=2) =
Pearson7(area/(hwhm*exp(lgamma(shape-0.5)-lgamma(shape))*sqrt(pi/(2^(1/shape)
-1))), center, hwhm, shape)\n"
"PseudoVoigtA(area, center, hwhm, shape=0.5) = GaussianA(area*(1-shape),
center, hwhm) + LorentzianA(area*shape, center, hwhm)\n"
! "ExpDecay(a=0, t=1) = a*exp(-x/t)\n"
! "LogNormalA(area, center, bandwidth=fwhm, asymmetry=0.1) =
LogNormal(sqrt(ln(2)/pi)*(2*area/bandwidth)*exp(-asymmetry^2/4/ln(2)),
center, bandwidth, asymmetry)",
"\n");
udfs.clear();
for (vector<string>::const_iterator i = formulae.begin();

}}}

[Marcin Wojdyr]

On Mon, Mar 23, 2009 at 08:55:52AM +0100, Scholz, Mirko wrote:
>
> Hello,
>
> if you accept new built-in functions, here is code for log-normal shaped

Hi Mirko, thanks, but the patch is malformed and can be easily applied.
Probably your mail program added a few additional line breaks into the
patch. Could you resend it as an attachment, preferably in unified
format (diff -u)?
Marcin

[Mirko Scholz]

> > if you accept new built-in functions, here is code for log-normal shaped

> Could you resend it as an attachment, preferably in unified
> format (diff -u)?

Ok, here with a real mua.

Sorry for the inconvenience

Mirko

[Marcin Wojdyr]

On Mon, Mar 23, 2009 at 09:07:22PM +0100, Mirko Scholz wrote:
> > > if you accept new built-in functions, here is code for log-normal shaped
> > Could you resend it as an attachment, preferably in unified
> > format (diff -u)?
> Ok, here with a real mua.

Thanks, I applied the patch.
The only change I did was commenting out one line, because of the warning:
bfunc.cpp:986: warning: unused variable ‘w’

I didn't change the names of variables, but `bandwidth' seems to
refer to a specific data, and I'd prefer a more generic name (it can
be a single letter) and I'd like to change it later.

And out of curiosity: Do you know if the LogNormal formula in fityk is
the same as in other programs? E.g. as this formula I just found:
http://www.systat.co.kr/products/TableCurve2D/help/1184.html

Marcin

--
Marcin Wojdyr | http://www.unipress.waw.pl/~wojdyr/

[Marcin Wojdyr]

On Wed, Mar 25, 2009 at 18:34, Mirko Scholz
<mirko....@stud.uni-goettingen.de> wrote:

> On Tue, Mar 24, 2009 at 12:53:26PM -0500, Marcin Wojdyr wrote:
>> The only change I did was commenting out one line, because of the warning:
>> bfunc.cpp:986: warning: unused variable ‘w’

> A right, but get_nonzero_range is never called???

it's called when option cut-function-level is set, it can make
calculations much faster in some cases, but I guess not many people
are using.

>> I didn't change the names of variables, but `bandwidth' seems to
>> refer to a specific data, and I'd prefer a more generic name (it can
>> be a single letter) and I'd like to change it later.

> As you please.  The name is related to electronic spectroscopy.
> FWHM describes it, too.

>
>> And out of curiosity: Do you know if the LogNormal formula in fityk is
>> the same as in other programs? E.g. as this formula I just found:
>> http://www.systat.co.kr/products/TableCurve2D/help/1184.html

> I used that formula of the article
>
> Bingemann, D.; Ernsting, N. P. J. Chem. Phys. 1995, 102, 2691–2700.
>
> because my work is somewhat related to theirs.  I can imagine
> that they used the one available in Origin, although І cannot
> tell you since I don't use Origin.
>
> Thanks for the link.  Having a bit thought about it, I think
> the two formulae differ substancially. 1st) The integral.
> 2nd) I have no idea how to convert 1/ln(d)**2 into 1/asymmetry.
>
> Parameter d is constrained to values between ]0;1[ or ]1;infty[
> while parameter asymmetry is constrained to values != 0
>
> On the other hand the functions simply look to close to have
> totally unrelated parameters.
>
> Hope it helps
>
> Mirko
>

I changed names bandwidth to width and asymmetry to asym, to make it shorter.
Thanks again for the patch,