Author Topic: 4th order polynomial coefficients for pressure at temperature readings, Help!!!  (Read 10589 times)

0 Members and 1 Guest are viewing this topic.

Offline itsbiodiversityTopic starter

  • Regular Contributor
  • *
  • Posts: 99
  • Country: us
This may be very dense and my apologies as I am new in this area of science (obviously), but in your calculations, what would the a0, a1, a2, a3, and a4 be?  I am trying to learn how to identify what those are.
 

Offline Nominal Animal

  • Super Contributor
  • ***
  • Posts: 7198
  • Country: fi
    • My home page and email address
The only real difficulty here is getting the picture right.  Hopefully, this one will get closer!

Did you happen to look at the spreadsheet I attached on the last post?
Yes, it opens in LibreOffice Calc just fine (don't have Windows), but it is difficult to understand what the fields mean.  Remember, we don't have the context you do.

If I've understood you correctly, you have a very high precision pressure sensor.  It has internal quartic polynomial compensation (via five coefficients) for the pressure sensor at a specific temperature, and you are looking for how to compute these coefficients for different temperatures.



Let's assume the following:

Column A, "Standard", is the known standard pressure being measured.
Column B, "Temp", is the temperature at which the measurements are made.
Column C, "Temp Counts", is the raw temperature ADC value obtained from the device.
Column D, "Reading", is the pressure reported by the device.

We can see that the temperature compensation at temperature 21 is correct, i.e. the reported pressure matches the pressure exactly.
At lower temperatures, the magnitude of the reading is too large.
At higher temperatures, the magnitude of the reading is too small.

To me, this indicates that the pressure sensor needs to be calibrated (using five coefficients) to work correctly at a given specific temperature, i.e. that the reported pressure is
$$P = P_0 + P_1 x + P_2 x^2 + P_3 x^3 + P_4 x^4$$
with the user able to supply \$P_0\$, \$P_1\$, \$P_2\$, \$P_3\$, and \$P_4\$ to the device.  It is these coefficients that depend on temperature; each set of five coefficients work at a given temperature.  For temperature 21, the coefficients are 0 1 0 0 0, i.e. no compensation: \$P = x\$.

Most importantly: We are supplying parameters to compensate the sensor behaviour at a specific temperature, not specifying the sensor behaviour at different temperatures.

Let's save the first four columns of each temperature to a text file that Gnuplot can eat. temperature-15.txt:
Code: [Select]
#Standard	Temp	Temp Counts	Reading
-15 15 5485 -15.017
-13 15 5485 -13.014
-11 15 5485 -11.011
-9 15 5485 -9.01
-7 15 5485 -7.008
-5 15 5485 -5.006
-3 15 5485 -3.004
-2 15 5485 -2.002
0 15 5485 0
2 15 5485 2.001
3 15 5485 3.002
5 15 5485 5.004
7 15 5485 7.006
9 15 5485 9.007
11 15 5485 11.008
13 15 5485 13.008
15 15 5485 15.008
The first and fourth columns are the same for the temperature 21 table, so we can skip that. (It'll just give us a linear function anyway, 0 1 0 0 0).
Note that Gnuplot ignores lines beginning with a # and those lines are just comments for us humans.

temperature-28.txt:
Code: [Select]
#Standard	Temp	Temp Counts	Reading
-15 28 9471 -14.99
-13 28 9471 -12.991
-11 28 9471 -10.993
-9 28 9471 -8.994
-7 28 9471 -6.996
-5 28 9471 -4.997
-3 28 9471 -2.998
-2 28 9471 -1.999
0 28 9471 0
2 28 9471 1.999
3 28 9471 2.998
5 28 9471 4.997
7 28 9471 6.996
9 28 9471 8.996
11 28 9471 10.996
13 28 9471 12.996
15 28 9471 14.996

What we want to do, is fit a fourth-degree polynomial on the fourth column, so that it yields the first column.  However, we can see that there is no zero drift: when the first column is zero, fourth column is zero also.  We can "enforce" that by using a polynomial without a constant coefficient.  In Gnuplot:
Code: [Select]
P(p) = P1*p + P2*p*p + P3*p**3 + P4*p**4
fit P(x) 'temperature-15.txt' using 4:1 via P1, P2, P3, P4
where we use fourth column as x , and result in first column.  This outputs
Code: [Select]
iter      chisq       delta/lim  lambda   P1            P2            P3            P4           
   0 7.4442087182e+09   0.00e+00  1.04e+04    1.000000e+00   1.000000e+00   1.000000e+00   1.000000e+00
   1 2.2294220617e+07  -3.33e+07  1.04e+03    9.985765e-01   9.947791e-01   7.430715e-01   9.651151e-03
   2 3.1706096563e+04  -7.02e+07  1.04e+02    9.945782e-01   9.802179e-01   2.086389e-02  -5.007133e-03
   3 2.5079458204e+03  -1.16e+06  1.04e+01    9.946186e-01   3.960710e-01   1.535787e-05  -2.025147e-03
   4 1.1420062378e-01  -2.20e+09  1.04e+00    9.975663e-01   2.682431e-03   8.654369e-06  -1.362051e-05
   5 4.7602533798e-06  -2.40e+09  1.04e-01    9.990046e-01   1.557076e-05   7.926488e-07   1.567490e-08
   6 4.7485651793e-06  -2.46e+02  1.04e-02    9.990118e-01   1.539077e-05   7.529785e-07   1.659336e-08
   7 4.7485651793e-06  -4.80e-07  1.04e-03    9.990118e-01   1.539077e-05   7.529765e-07   1.659336e-08
iter      chisq       delta/lim  lambda   P1            P2            P3            P4           

After 7 iterations the fit converged.
final sum of squares of residuals : 4.74857e-06
rel. change during last iteration : -4.80261e-12

degrees of freedom    (FIT_NDF)                        : 13
rms of residuals      (FIT_STDFIT) = sqrt(WSSR/ndf)    : 0.000604379
variance of residuals (reduced chisquare) = WSSR/ndf   : 3.65274e-07

Final set of parameters            Asymptotic Standard Error
=======================            ==========================
P1              = 0.999012         +/- 4.097e-05    (0.004101%)
P2              = 1.53908e-05      +/- 4.78e-06     (31.06%)
P3              = 7.52977e-07      +/- 2.475e-07    (32.87%)
P4              = 1.65934e-08      +/- 2.544e-08    (153.3%)

correlation matrix of the fit parameters:
                P1     P2     P3     P4     
P1              1.000
P2              0.001  1.000
P3             -0.917 -0.002  1.000
P4             -0.001 -0.961  0.002  1.000
which means that at temperature 15, the coefficients we should supply to the device are 0 0.999012 0.0000153908 0.000000752977 0.0000000165934

Running the same for temperature-28.txt outputs (shortened to include only the important bit)
Code: [Select]
Final set of parameters            Asymptotic Standard Error
=======================            ==========================
P1              = 1.0006           +/- 2.148e-05    (0.002147%)
P2              = -8.84714e-06     +/- 2.51e-06     (28.37%)
P3              = -5.92832e-07     +/- 1.301e-07    (21.95%)
P4              = -2.39269e-08     +/- 1.339e-08    (55.97%)
i.e. the device coefficients needed to get correct pressure readings at temperature 28 are 0 1.0006 -0.00000884714 -0.000000592832 -0.0000000239269.



The next question is, obviously, how to obtain the five coefficients for any specific temperature.  We have three sets of parameters, so we can interpolate between them (and extrapolate outside them) using quadratic interpolation:
$$\begin{aligned}
p_{15}(x) &=  0.999012 x + 0.0000153908 x^2 + 0.000000752977 x^3 + 0.0000000165934 x^4 \\
p_{21}(x) &= x \\
p_{28}(x) &= 1.0006 x - 0.00000884714 x^2 - 0.000000592832 x^3 - 0.0000000239269 x^4 \\
w_{15}(t) &= 98/13 - 49/78 t + 1/78 t^2 \\
w_{21}(t) &= -10 + 43/42 t - 1/42 t^2 \\
w_{28}(t) &= 45/13 - 36/91 t + 1/91 t^2 \\
p(x, t) &= \frac{ p_{15}(x) w_{15}(t) + p_{21}(x) w_{21}(t) + p_{28}(x) w_{28} t }{ w_{15}(t) + w_{21}(t) + w_{28}(t) }
\end{aligned}$$
where \$w_{15}(15) = 1\$, \$w_{15}(21) = 0\$, \$w_{15}(28) = 0\$; \$w_{21}(15) = 0\$, \$w_{21}(21) = 1\$, \$w_{21}(28) = 0\$; \$w_{28}(15) = 0\$, \$w_{28}(21) = 0\$, and \$w_{28}(28) = 1\$.  These fully define the nine coefficients.  It turns out that \$w_{15}(t) + w_{21}(t) + w_{28}(t) = 1\$ for all \$t\$, so we can omit the divisor.

When combining polynomials, evaluating the polynomials first and then scaling and summing the results is equivalent to scaling and summing the coefficients first, then evaluating the resulting polynomial.  So, with the above, we can immediately describe the quadratic interpolation/extrapolation of the five configuration values to the device as a function of ambient temperature \$t\$:
  • First coefficient is always \$0\$ (because the above fit observed there is no zero drift)
  • Second coefficient is \$0.9946289230769221 + 0.0003833040293040879 t - 0.000006073260073258604 t^2\$
  • Third coefficient is \$0.00008539823846153845 - 0.000006168612014652013 t + 0.00000010009663003663 t^2\$
  • Fourth coefficient is \$0.000003624177384615384 - 0.0000002384970677655678 t + 0.000000003138913919413919 t^2\$
  • Fifth coefficient is \$ 0.00000004226482307692308 - 0.0000000009584721611721613 t - 0.00000000005019706959706964 t^2\$
These I obtained by telling Maxima to expand the expression \$p(x, t)\$, and collecting the coefficients for each exponent of \$x\$ and \$t\$.

Note that if you had more samples, i.e. more temperatures (and perhaps at more than the 17 pressures), you could model the temperature dependence even better: this is the "best" approximation we can do, that should match the configuration data at the configuration temperatures and pressures, but makes no other assumptions.  (Even the quadratic weight polynomials \$w_{15}(t)\$, \$w_{21}(t)\$, and \$w_{28}(t)\$ are completely defined by having them be 1 at that specific temperature, and zero elsewhere.)
« Last Edit: August 19, 2020, 04:39:31 pm by Nominal Animal »
 
The following users thanked this post: itsbiodiversity

Offline Nominal Animal

  • Super Contributor
  • ***
  • Posts: 7198
  • Country: fi
    • My home page and email address
in your calculations, what would the a0, a1, a2, a3, and a4 be?
I assumed that you are in full control of the data conversion, obtaining raw pressure and temperature values.

In other words, I showed the model I believe the device should already have incorporated, with say 5×3 = 15 coefficients.

Instead, it looks like the device needs to be configured for a specific ambient temperature.  The pressure sensor is factory calibrated to give correct, linear results at 21°C, but the compensation needs to be adjusted for other temperatures.

My latter response shows how to calculate the five coefficients for a given temperature (in Celsius centigrade).  If you want to change that to the temperature reading (5485 at 15°C, 7386 at 25°C, and 9471 at 28°C), you need to change the weight functions.  In Maxima:
Code: [Select]
w15(t) := W00 + W01*t + W02*t^2 $
w21(t) := W10 + W11*t + W12*t^2 $
w28(t) := W20 + W21*t + W22*t^2 $
solve([ w15(5485)=1, w15(7386)=0, w15(9471)=0,
        w21(5485)=0, w21(7386)=1, w21(9471)=0,
        w28(5485)=0, w28(7386)=0, w28(9471)=1 ], [ W00,W01,W02, W10,W11,W12, W20,W21,W22 ]);
which yields
Code: [Select]
[[W00 = 34976403/3788693, W01 = -16857/7577386, W02 = 1/7577386,
  W10 = -3463229/264239, W11 = 14956/3963585, W12 = -1/3963585,
  W20 = 1350407/277027, W21 = -12871/8310810, W22 = 1/8310810]]
So, we construct those polynomials, and see what \$p(x, t)\$ looks like:
Code: [Select]
w15(t) := 34976403/3788693 - 16857/7577386*t + 1/7577386*t^2 $
w21(t) := -3463229/264239 + 14956/3963585*t - 1/3963585*t^2 $
w28(t) := 1350407/277027 -12871/8310810*t + 1/8310810*t^2 $
If we check expand(w15(t) + w21(t) + w28(t)); Maxima tells us the result is always 1, so we can omit the divisor.  Next, we expand \$p(x, t)\$:
Code: [Select]
p15(x) := 0.999012*x + 0.0000153908*x^2 + 0.000000752977*x^3 + 0.0000000165934*x^4 $
p21(x) := x $
p28(x) := 1.0006*x - 0.00000884714*x^2 - 0.000000592832*x^3 - 0.0000000239269*x^4 $
expand( w15(t)*p15(x) + w21(t)*p21(x) + w28(t)*p28(x) );
which yields
Code: [Select]
-6.891513009452426E-16*t^2*x^4
+1.412920651658525E-13*t*x^4
+3.655168497054565E-8*x^4
+2.803896602141317E-14*t^2*x^3
-7.569847447073744E-10*t*x^3
+4.061479723198378E-6*x^3
+9.666149779038884E-13*t^2*x^2
-2.053746129643471E-8*t*x^2
+9.895794611233592E-5*x^2
-5.8192843849845E-11*t^2*x
+1.268726552948896E-6*t*x
+0.9938037796576857*x
The first coefficient given to the device is always 0, \$C_0 = 0\$, because there is no constant term (there was no zero drift in the data).  Terms for \$x\$ give the first coefficient,
$$C_1 = 0.9938037796576857 +1.268726552948896 \cdot 10^{-6} t - 5.8192843849845\cdot 10^{-11} t^2$$
Terms for \$x^2\$ give the third coefficient, and so on, with terms for \$x^4\$ giving the fifth coefficient,
$$C_4 = 3.655168497054565\cdot 10^{-8} + 1.412920651658525 \cdot 10^{-13} t - 6.891513009452426 \cdot 10^{-16} t^2$$
These use the standard scientific notation, \$\text{1.2e-3} = 1.2 \cdot 10^{-3} = 0.0012\$.

Because the device does this compensation, we cannot do the bivariate polynomial.  (Why? Odd device.  It's not like 15 multiplications and 14 additions is that much more expensive, per pressure reading, than 5 multiplications and 4 additions.)

Instead, we must find the quartic polynomial the device can use to do the compensation itself, for each temperature separately.  Pretty annoying, but scientifically quite valid, and I kinda know the reasoning why (it keeps policy in the userland, so to speak).  Here, the cross terms appear in the interpolated form above (the last one we expanded - the one with terms including powers of \$x\$ and \$t\$), and are "hidden" in the changes of the five coefficients supplied to the device.  So, as long as the numerical computation is done with more digits than are correct in the result plus sufficient digits to cover the rounding errors (1 ULP per multiplication or addition – so 4 bits or 1.2 digits per 16 multiplications or additions), this sequential fitting model is just as precise as the single bivariate one would be.  So no worries!
« Last Edit: August 19, 2020, 05:12:48 pm by Nominal Animal »
 

Offline itsbiodiversityTopic starter

  • Regular Contributor
  • *
  • Posts: 99
  • Country: us
We are so so so so close.  I am going to use Ambient temperature and the calculations below.  On return to the office I will apply the coefficients and let you know the result.  Thank you!
 

Offline itsbiodiversityTopic starter

  • Regular Contributor
  • *
  • Posts: 99
  • Country: us
For all you've done to get me this close I truly appreciate it.  I am trying to calculate these formulas provided using t=21 (is that correct).  If so I get the following:

Coefficient #   Calculated (a)                                                  Temp         
0   0             
1   0.993830397   0.99380378   1.27E-06          5.82E-11   21    FORMULA : =0.99380378 + (1.27E-06*21)-(5.82E-11*21^2) = 0.993830397
2   9.85E-05           9.90E-05           -2.05E-08          9.67E-13    21    FORMULA : =9.90E-05 + (-2.05E-08*21)+(9.67E-13*21^2) = 9.85E-05
3   4.05E-06           4.06E-06           -7.57E-10          2.80E-14      21    FORMULA : =4.06E-06 + (-7.57E-10*21)+(2.8E-14*21^2) = 4.05E-06
4   3.6555E-08   3.66E-08            1.41E-13         -6.89E-16     21    FORMULA : =3.66E-08 + (1.41E-13*21)-(6.89E-16*21^2) = 3.6555E-08

I tried entering the seqence a0 = 0, a1 = 0.993830397, a2 = 9.85E-05, a3 = 4.05E-06, and a4 = 3.6555E-08 and the pressure went offscale.  I also tried using Counts rather than degrees C and it also went offscale.  Am I misunderstanding how to calculate those formulas above?  If you're sick of helping I don't blame you!!!  I will have another unit in house soon that was programmed by them prior to leaving the organization and am interested to see what is stored.  Default is 0 0 1 0 0, and the module meets specifications across temperatures but hoping to "dial it in".  I am kind of "pigeon holed" because the inventor/programmer did the design but all information regarding their application of the statistics is gone to the wind. 
 

Offline ahbushnell

  • Frequent Contributor
  • **
  • Posts: 751
  • Country: us
Thank you for such a detailed explanation, and I can clearly see how that technique would be superior.  Would employing that technique produce similarly effective temperature coefficients?  Is there any way to help me identify what those a0-a4 coefficients would be in both of these scenarios.  I can apply this to the unit and report back the effects.

Sounds like going back to the source of the parameters and getting clarification would be a good idea. 
 

Offline Nominal Animal

  • Super Contributor
  • ***
  • Posts: 7198
  • Country: fi
    • My home page and email address
For all you've done to get me this close I truly appreciate it.  I am trying to calculate these formulas provided using t=21 (is that correct).
No, not exactly.  Also check if your environment uses ** for exponentiation and ^ for exclusive-or.

For temperature \$t\$ in °C, we use the formulae
Code: [Select]
a0 = 0
a1 = 0.9946289230769221 + 0.0003833040293040879 * t - 0.000006073260073258604 * t * t
a2 = 0.00008539823846153845 - 0.000006168612014652013 * t + 0.00000010009663003663 * t * t
a3 = 0.000003624177384615384 - 0.0000002384970677655678 * t + 0.000000003138913919413919 * t * t
a4 = 0.00000004226482307692308 - 0.0000000009584721611721613 * t - 0.00000000005019706959706964 * t * t
so for \$t = 21\$ (°C), they are a0=0, a1=1, a2=0, a3=0, and a4=0.

As an example, for \$t = 24\$ (°C), they are 0 +1.000330021978023 -4.992790989010987E-6 -2.91737824175825E-7 -9.652020879120904E-9.

For \$t = 18\$(°C), they are 0 +0.9995606593406599 +6.794530329670339E-6 +3.482382747252735E-7 +8.748473626373614E-9.

The coefficients you specified, best match temperatures \$t \approx -2\$ (°C) or \$t \approx 64\$ (°C), so it is not a surprise the sensor reported pressure out of range, unless the ambient temperature was that extreme (compared to the calibrated range, +15°C to +28°C).



If you want to use the temperature counts (\$c\$) instead, then use
Code: [Select]
a0 = 0
a1 = 0.9938037796576857 + 1.268726552948896E-6 * c - 5.8192843849845E-11 * c * c
a2 = 9.895794611233592E-5 - 2.053746129643471E-8 * c + 9.666149779038884E-13 * c * c
a3 = 4.061479723198378E-6 - 7.569847447073744E-10 * c + 2.803896602141317E-14 * c * c
a4 = 3.655168497054565E-8 + 1.412920651658525E-13 * c - 6.891513009452426E-16 * c * c
At c = 7386 (which corresponds to \$t = 21°C\$), we have a0 = 0, a1 = 0.999999999999997 ≃ 1.000000, a2 = 7.453889935837843e-20 ≃ 0.000000, a3 = -2.117582368135751e-22 ≃ -0.000000, and a4 = 1.985233470127266E-23 ≃ 0.000000 .



Let's check the math anyways, that's always the best option.  I'll use Maxima syntax, but you can adjust it trivially to SageMath, Maple, Mathematica etc.

Let's use the temperature counts.  The weight functions, the fitted pressure compensation functions, and the combined pressure estimate are
Code: [Select]
w15(t) := 34976403/3788693 - 16857/7577386*t + 1/7577386*t^2 $
w21(t) := -3463229/264239 + 14956/3963585*t - 1/3963585*t^2 $
w28(t) := 1350407/277027 -12871/8310810*t + 1/8310810*t^2 $
p15(x) := 0.999012*x + 0.0000153908*x^2 + 0.000000752977*x^3 + 0.0000000165934*x^4 $
p21(x) := x $
p28(x) := 1.0006*x - 0.00000884714*x^2 - 0.000000592832*x^3 - 0.0000000239269*x^4 $
p(x,t) := w15(t)*p15(x) + w21(t)*p21(x) + w28(t)*p28(x) $
The temperature count at \$21°C\$ is \$7386\$, so let's look what the pressure compensation function looks like at that temperature:
Code: [Select]
expand(p(x,7386));
It yields simply x , like it should: no compensation.

The temperature count at \$15°C\$ is \$5485\$, and expand(p(x,5485)); yields 1.65934E-8*x^4+7.52977E-7*x^3+1.53908E-5*x^2+0.999012*x , which is the same as p15(x) :
Code: [Select]
expand(p(x,5485) - p15(x));
yields 0.0.  Similarly for \$28°C\$ at \$9471\$: expand(p(x,9471)-p28(x)) yields 0.

So, we know the weighting function works, we haven't goofed there.

Let's check some of the calibration values, using p(x,t); with x from the fourth column in the tabulated data, t from the third column, and the expected result from the first column. In tabulated form, the results are
Code: [Select]
   x     │   t  │   p(x,t)   │ Expectd │ Error
─────────┼──────┼────────────┼─────────┼─────────────────
 -15.017 │ 5485 │ -15.000398 │  -15    │ -0.000398
  -9.01  │ 5485 │  -9.000290 │   -9    │ -0.000290
  -3.004 │ 5485 │  -3.000912 │   -3    │ -0.000912
   0     │ 5485 │   0        │    0    │  0
   3.002 │ 5485 │   2.999194 │    3    │  0.000806
   9.007 │ 5485 │   9.000009 │    9    │  0.000009
  15.008 │ 5485 │  15.000026 │   25    │  0.000026
─────────┼──────┼────────────┼─────────┼─────────────────
 -15     │ 7386 │ -15        │  -15    │  0
  -9     │ 7386 │  -9        │   -9    │  0
  -3     │ 7386 │  -3        │   -3    │  0
   3     │ 7386 │   3        │    3    │  0
   9     │ 7386 │   9        │    9    │  0
  15     │ 7386 │  15        │   15    │  0
─────────┼──────┼────────────┼─────────┼─────────────────
 -14.99  │ 9471 │ -15.000193 │  -15    │ -0.000193
  -8.994 │ 9471 │  -8.999837 │   -9    │  0.000163
  -2.998 │ 9471 │  -2.999864 │   -3    │  0.000136
   0     │ 9471 │   0        │    0    │  0
   2.998 │ 9471 │   2.999701 │    3    │ -0.000299
   8.996 │ 9471 │   9.000093 │    9    │  0.000093
  14.996 │ 9471 │  14.999799 │   15    │ -0.000201
From this, we can see that the absolute errors are less than 0.001.  The spreadsheet says tolerance is 0.035% of full range, which corresponds to 30×0.035/100 = 0.0105.  Using the function shown above, the results are correct (compared to calibration values) with less than 0.001/30×100% = 0.0033% of (absolute) error.



Let's verify the coefficient functions I gave earlier, for the temperature counts:
Code: [Select]
A4(t) := -6.891513009452426E-16*t^2
         +1.412920651658525E-13*t
         +3.655168497054565E-8      $
A3(t) := +2.803896602141317E-14*t^2
         -7.569847447073744E-10*t
         +4.061479723198378E-6      $
A2(t) := +9.666149779038884E-13*t^2
         -2.053746129643471E-8*t
         +9.895794611233592E-5      $
A1(t) := -5.8192843849845E-11*t^2
         +1.268726552948896E-6*t
         +0.9938037796576857        $
A0(t) := 0 $
q(x, t) := A0(t) + A1(t)*x + A2(t)*x^2 + A3(t)*x^3 + A4(t)*x^4 $
and repeating the table, but using q(x, t) instead of p(x, t), we see that the results are the same.  In fact, expand(p(x,t) - q(x,t)); gives a polynomial with coefficients of the order 10-21 to 10-29, which is close enough to zero for us to say the two polynomials are the same to within the precision we represent them.



If the correctly calculated coefficients – please don't round them, all digits matter – still produce erroneous pressures, then a core assumption I made is incorrect.

I assume the formula the device uses to report the pressure is \$P(p) = C_0 + C_1 p + C_2 p^2 + C_3 p^3 + C_4 p^4\$ (with \$p\$ the internal factory-calibrated sensor value, and \$P(p)\$ the pressure it provides to the user).  (Mainly, because the defaults being 0 1 0 0 0 , it would make sense for it to correspond to \$P(p) = p\$.  However, if OP is not getting sane results, this assumption is wrong.)

If we look at the spreadsheet, it has a separate cell for "constant" (zero – equal to my assumption about no zero drift!), and four X coefficients, but it is the second X coefficient (corresponding to \$C_2\$).  Could it be that the coefficients the device wants are for the temperature and not the pressure?  That it provides \$P(p, t)\$ to the user, internally calculating the pressure as \$P(p, t) = p + T_0 + T_1 t + T_2 t^2 + T_3 t^3 + T_4 t^4\$, with \$T_0\$ to \$T_4\$ being the four coefficients we can specify?  Well, no. You see, if that was the case, we could not get \$P(p_{-}, t) = -15\$, \$P(0, t) = 0\$, and \$P(p_{+},t) = +15\$ for a constant \$t\$ but \$p_{-} \ne -p_{+}\$.

In fact, I cannot imagine right now what the exact formula for \$P(p, t)\$ the device uses, if not \$C_0 + C_1 p + C_2 p^2 + C_3 p^3 + C_4 p^4\$.  This does not mean it cannot be anything else; just that it is something not immediately obvious.  Perhaps \$C_0 + C_1 p + C_2 p^2 + C_3 t + C_4 t^2\$ or something?  I don't know, and I don't see anything specific in the spreadsheet that would help us find what it is.

I am kind of "pigeon holed" because the inventor/programmer did the design but all information regarding their application of the statistics is gone to the wind.
I know how that feels.  I'd like to help one team to measure pressure in a fusor, and temperature of the vacuum cask, but I seem to be the only one willing to use non-off-the-shelf parts and non-Raspberry Pi SBCs, which makes everything harder.  (Although I do hear that the "I shall randomly drop USB packets because I'm overburdened right now" bugs have mostly worked around nowadays.)
Modular hardware you can build and test in units (with calibration voltage sources, compare to known good measurements) and in increasingly complex combinations, but if you use COTS hardware with LabView-only interfaces or arse-backwards 'Pi libraries that do not even check for short reads and writes, you can only hope stuff works for the duration of the experiment.  (On the other hand, when something goes wrong, you get to blame other people, and not tarnish your own image.)
Me, I like to make robust, reliable stuff on the cheap.
« Last Edit: August 19, 2020, 09:05:24 pm by Nominal Animal »
 

Offline IanB

  • Super Contributor
  • ***
  • Posts: 12539
  • Country: us
This may be very dense and my apologies as I am new in this area of science (obviously), but in your calculations, what would the a0, a1, a2, a3, and a4 be?  I am trying to learn how to identify what those are.

I'm having trouble following all the details Nominal Animal is posting, but let's simplify it a bit.

There is some difference between the actual pressure reading (call it P) and the assumed true pressure (call it P*). So there is an offset between P and P* (call it deltaP).

So we have:

   deltaP = P - P*

Looking at the data table, it is clear that deltaP depends both on T and P.

So we need to fit a correlating function of the form:

   deltaP = f(P, T)

We might assume, for example:

  deltaP ~ (T-Tref) * (a0 + a1 * P + a2 * P * (T-Tref) + a3 * P^2 + a4 * P^2 * (T-Tref))  (say)

You have to know the actual polynomial assumed by the designer of the equipment. What I wrote here is just for instance.

Then your job is to find values of a0..a4 that best fit the data in the table. You can use Excel for that, or any other tool of your choice.

To simplify matters I chose to use T-Tref because the offset appears to be zero at Tref = 21°C. But really, Tref is just another coefficient in the polynomial.

At the end of the day, the only way to know what a0..a4 are is to know the form of the interpolating polynomial, and the designer has to tell you that.
« Last Edit: August 19, 2020, 10:14:36 pm by IanB »
 
The following users thanked this post: itsbiodiversity

Offline Nominal Animal

  • Super Contributor
  • ***
  • Posts: 7198
  • Country: fi
    • My home page and email address
Yeah, sorry about the excess verbosity, everyone.

I basically just explained, step by step, how to obtain \$P(p) = P_1(t) p + P_2(t) p^2 + P_3(t) p^3 + P_4(t) p^4\$, where \$p\$ is the device-internal pressure calibrated to be correct at \$t = 21°C\$ (which corresponds to \$c = 7386\$), and how to derive \$P_1(t)\$, \$P_2(t)\$, \$P_3(t)\$, and \$P_4(t)\$ from the samples at \$t = 15°C\$, \$t = 21°C\$, and \$t = 28°C\$.  (And the same for using temperature counts \$c\$ instead of temperature \$t\$.)

My core assumption is/was that the five (\$P_0(t) = P_0(c) = 0\$, as there is no zero drift in the calibration data) are coefficients for powers of \$p\$, but that may not be true.  It could be \$P(p, t) = P_0 + P_1 p + P_2 p^2 + P_3 t + P_4 t^2\$, or \$P(p, c) = P_0 + P_1 p + P_2 p^2 + P_3 c + P_3 c^2\$, or something else entirely.

The spreadsheet sorta-kinda suggests the last one (\$P(p,c)\$) – it has no formulae for the derivation of the five coefficients used, just using the column labels as hints –, but \$P(p)\$ would be the most logical one, considering how the device is used.



Actually, if we look at the spreadsheet, a possible form is \$P(p, c) = A_0 p^2 + A_1 p + A_2 c^2 + A_3 c + A_4\$.  If so, the coefficients fitted to the calibration data are
      0.000000669808265654805  0.999876931661194  -0.0000000000881848278027671  0.000000710797927540773  0
If these produce sane/good results, I'll be happy to show how to obtain these using the dataset and Gnuplot, although it is almost trivial.
« Last Edit: August 20, 2020, 12:28:44 pm by Nominal Animal »
 
The following users thanked this post: itsbiodiversity

Offline itsbiodiversityTopic starter

  • Regular Contributor
  • *
  • Posts: 99
  • Country: us
HOLY CRAP MAN!!!  That was close.  It actually did not cause the positive pressure linearity to vary from what I set as the "ambient" linearity; however, it caused the formerly near perfect "ambient" negative pressure values to change near identical to the un-adjust reading in Column D (-15.017 psi at a standard of -15 psi).  Basically it went "backwards" to the unadjusted cold temperature value even though I am operating at ambient.  Does this explanation help guide to a portion of the formula???  I am very interested how you got to this point!  Of course the negative pressure was the last data point I took I thought you had it.
 

Offline Jay_Diddy_B

  • Super Contributor
  • ***
  • Posts: 2766
  • Country: ca
Hi,

It make intuitive sense.

You have a two second order polynomials, one for temperature and one for pressure.

The constant term A4 is the sum of the constant terms in each of the second order polynomials.

On the surface plot that i shared no correction is needed at Temp = 21 (all pressures) or pressure = 0 (all temperatures)

Jay_Diddy_B
 
The following users thanked this post: itsbiodiversity

Offline Nominal Animal

  • Super Contributor
  • ***
  • Posts: 7198
  • Country: fi
    • My home page and email address
Let's recap.  We have five coefficients that can be defined in the device.  Let's call these a0 a1 a2 a3 a4 .
We do not know what these coefficients are, exactly, but we want to find out.  Because of the default, we can be almost certain that a1 is the first-degree coefficient for the pressure.

Let's assume the internal compensation function is P(p, c) = P2 p2 + P1 p + C2 c2 + C1 c + Oo.  Then, we can assume that a1 = P1, i.e. P1 is always the second one.  The simplest option, then, would be to try the most likely possibilities, and see if they make sense.  Here's the ones I would try:
  • 0.000000669808265654805  0.999876931661194  -0.0000000000881848278027671  0.000000710797927540773  0
    This corresponds to P2 P1 C2 C1 Oo, but we already know this didn't work.
  • 0  0.999876931661194  0.000000669808265654805  0.000000710797927540773  -0.0000000000881848278027671
    This corresponds to Oo P1 P2 C1 C2 .
  • 0  0.999876931661194  0.000000710797927540773  0.000000669808265654805  -0.0000000000881848278027671
    This corresponds to Oo P1 C1 P2 C2 .
  • 0.000000669808265654805  0.999876931661194  0  -0.0000000000881848278027671  0.000000710797927540773
    This corresponds to P2 P1 Oo C2 C1 .
  • 0.000000669808265654805  0.999876931661194  0  0.000000710797927540773  -0.0000000000881848278027671
    This corresponds to P2 P1 Oo C1 C2 .
  • 0.000000710797927540773  0.999876931661194  -0.0000000000881848278027671  0.000000669808265654805  0
    This corresponds to C1 P1 C2 P2 O0 .
There are a total of 24 possible combinations with P1 in the second position, but I think there has to be some kind of logic in them, not just a random jumble, and the above ones are the ones I would test.

If anyone has additional suggestions I didn't immediately see, please pipe up!  Me err often, me be uncle bumblespork.
(I don't just like peer review; my output requires it.)



Another option is to do ten sets of samples at different temperatures and pressures, each set with a specific coefficient set:
Code: [Select]
0 1 0 0 0
+0.00000001  0  0  0  0
-0.00000001  0  0  0  0
0 0 +0.00000001 0 0
0 0 -0.00000001 0 0
0 0 0 +0.00000001 0
0 0 0 -0.00000001 0
0 0 0 0 +0.00000001
0 0 0 0 -0.00000001
Instead of ±0.00000001, it would be better to use as large value in magnitude as possible, without error in the reported pressure.

Note that each of the ten sets are independent: they don't need to have samples at the same temperature and pressure, as long as there are at least five measurements, at at least three different pressures and three different temperatures.  You can either do a sweep of temperature or pressure, or measure a specific pressure and temperature using each of the ten sets of parameters.

The idea is that when we have a table with eight columns – the five coefficients used, the temperature (count) and pressure received, and the actual pressure being measured –, we can deduce the exact algebraic form of the compensation function the device uses.  For each set of five coefficients, we do need samples at at least three different temperatures and three different pressures, so that we can differentiate between linear terms and quadratic terms.  (If we suspect there are third degree terms, or higher, we need correspondingly more samples to deduce the exponent/power.)
« Last Edit: August 20, 2020, 06:59:42 pm by Nominal Animal »
 

Offline IanB

  • Super Contributor
  • ***
  • Posts: 12539
  • Country: us
Let's recap.  We have five coefficients that can be defined in the device.  Let's call these a0 a1 a2 a3 a4 .
We do not know what these coefficients are, exactly, but we want to find out.  Because of the default, we can be almost certain that a1 is the first-degree coefficient for the pressure.

Unfortunately, it doesn't seem like the meaning of a0..a4 should be unknown. If we look at the OP's earlier posts, there is apparently access to the inventor of the device, and there is access to a manual. This should not be a mystery. Either the inventor or the manual should be able to say, "This is the form of the correction polynomial programmed into the device."

I am perplexed that not once has the OP provided this information in plain text. It would be simple to answer if only that information were given.

I am using a previous manual that used the term "correction factor" - I spoke to the inventor yesterday

The UUT counts is used to calibrate to engineering units of psi.  The Counts of PSI wasn't used in the author's formula, just the counts for temperature at varying pressures.
 

Offline itsbiodiversityTopic starter

  • Regular Contributor
  • *
  • Posts: 99
  • Country: us
My apologies having issues quoting.  But to answer a few questions.  First off the person I spoke to initially isn't responding anymore.  Here is verbiage concerning the coefficients from the admin manual.

Enter Temperature Coefficients
This command provides for the entry of temperature coefficients that will compensate the sensor for ambient temperature conditions. The
coefficients are determined during the factory calibration process.
tcomp <a0> <a1> <a2> <a3> <a4>
The operator specifies five coefficients, which are used in a fourth order polynomial that corrects temperature readings for the ambient temperature at the sensor.

From the cal procedure:
"Take your recorded readings from the 3 temperatures and enter them into the polynomial linearization program to get your coefficients."

The "linearization program" I am told is a Quattro spreadsheet that I have access to.  I am not able to upload the spreadsheet here or I would.  I am attaching that pasted in excel to this comment.  In Quattro I can see the following is used, and to great effect actually, to normalize the pressure readings on the UUT to match the standard.
+$A:$N$11*C7+$A:$O$11*D7+$A:$P$11*E7+$A:$Q$11*F7

where N11 is an x coefficient, C7 is the temperature counts, O11 is another x coefficient, D7 is the UUT PSI reading, P11 is the third x coefficient, E7 is Temperature Counts (C7) x UUT PSI reading (D7), Q11 is the fourth x coefficient, and F7 is temperature counts squared (C7 x C7) multiplied by UUT PSI reading (D7).

The x coefficients referenced above are the four x coefficients in ZERO regression output left to right.

This is all of the information I have - as it appears this team abandoned the project at some point and I'm working off of what they did.  My apologies for not having further information - it has been frustrating trying to get information from them.  You've all helped more than you know already.

In the attached file I went ahead and wrote the equation for correction to column G labeled CORRECTED so that one can see what is going on.  The confusion I am having is what within the simple formula is considered a0 through a4.  Does this help at all?

« Last Edit: August 20, 2020, 09:05:17 pm by itsbiodiversity »
 

Offline itsbiodiversityTopic starter

  • Regular Contributor
  • *
  • Posts: 99
  • Country: us
I just received this info and I hope this helps.

Here is how the unit uses the coefficients:

"/* 3.2 - Function to Apply Temperature Compensation */
float ApplyCompensation (float rdg, int tempAdc)
    {
    float tmp = (float)tempAdc;
    float tmp2 = tmp * tmp;
    float limit, result;

    /* Limit the result to comform to the linearization table */
    limit = tblCnt ? fabs(calTbl[tblCnt-1].value) * 1.5 : 17.0;
   
    /* Adjust the pressure reading using temperture polynomial */
    result = coefTbl[0]
            + (coefTbl[1] * tmp)
            + (rdg * (coefTbl[2] + (coefTbl[3] * tmp) + (coefTbl[4] * tmp2)));

    /* Don't allow the result to stray too far from the table's limits */
    if (result < -limit)
        result = -limit;
    else if (result > limit)
        result = limit;

    /* Provide the temperature compensated pressure reading */
    return (result);
    }

You may not be a connoisseur of the C language, so the distilled version is:

result = a0 + (a1 * tmp) + X * (a2 + (a3 * tmp) + (a4 * tmp2 ))

   where: tmp is temperature ADC and X is the pressure reading


I seem to recall that A0 was usually zero". 
 
The following users thanked this post: Nominal Animal

Offline itsbiodiversityTopic starter

  • Regular Contributor
  • *
  • Posts: 99
  • Country: us
I have just successfully programmed the first of these units thanks to everyone's help on this thread.

In the end, for a similar unit to the one data was provided for the coefficients were:
a0 = 0
a1 = 6.18E-10
a2 = 0.9938
a3 = 1.16E-06
a4 = -5.43749E-11

a2 has a large effect on the unit.  Once these were established I reprogrammed certain lines of the calibration table and the unit successfully met the 0.011 psi requirement.  I am going to be diving more in to all of the information provided thus far to gain a better understanding of this type of statistics.  Again - thanks.  And if anyone has any input on these coefficients I am interested in what your thoughts are.
 

Offline Nominal Animal

  • Super Contributor
  • ***
  • Posts: 7198
  • Country: fi
    • My home page and email address
Yes!  This is the key: result = a0 + (a1 * tmp) + X * (a2 + (a3 * tmp) + (a4 * tmp2 )) where tmp is the temperature count and tmp2 is that squared, and X is the uncorrected pressure reading.

Using \$c\$ for temperature counts (\$c = 7386\$ corresponding to 21°C), and \$p\$ for the uncorrected pressure reading, the correction function is
$$P(p, c) = a_0 + a_1 c + a_2 p + a_3 c p + a_4 c^2 p = a_0 + a_1 c + p (a_2 + a_3 c + a_4 c^2 )$$



Using Gnuplot with the four-column calibration data (all three sets) in data.txt,
Code: [Select]
P(p, c) = a0 + a1*c + p*(a2 + a3*c + a4*c*c)
fit P(x,y) 'data.txt' using 4:3:1 via a0, a1, a2, a3, a4
print a0, a1, a2, a3, a4
gives us a nice set of coefficients:
0.00504498959142752 -6.45913664596841e-07 0.994401298866226 1.15958714191414e-06 -5.43693236648548e-11

However, keeping a0 zero since there is no zero drift makes still sense, so
Code: [Select]
P(p, c) = a1*c + p*(a2 + a3*c + a4*c*c)
fit P(x,y) 'data.txt' using 4:3:1 via a1, a2, a3, a4
print 0, a1, a2, a3, a4
gives us the set
0 6.18461239919483e-10 0.994400954117506 1.15967494730689e-06 -5.4374892229142e-11

Could you try this one, itsbiodiversity?

(You can use print P(-15.017,5485) (fourth column, third column, compare result to first column) to see how well the compensation matches the calibration dataset.  I'll omit the comparisons for now, to keep my posts more readable!)
« Last Edit: August 20, 2020, 10:01:56 pm by Nominal Animal »
 
The following users thanked this post: itsbiodiversity

Offline itsbiodiversityTopic starter

  • Regular Contributor
  • *
  • Posts: 99
  • Country: us
I just disconnected after battling with this one for days!  I have four more to program tomorrow and will let you know how those go using the coefficients you provided.  Thanks again!
 

Offline itsbiodiversityTopic starter

  • Regular Contributor
  • *
  • Posts: 99
  • Country: us
Those are the coefficients NOMINAL ANIMAL. Exactly. Now I need to figure out how to calculate those myself based on any dataset using excel.  I tried to download Gnuplot but had some issues installing.
 

Offline Nominal Animal

  • Super Contributor
  • ***
  • Posts: 7198
  • Country: fi
    • My home page and email address
Well, instead of a linear regression, you could do a brute-force (kinda-sorta gradient descent) numerical search that minimizes the largest error in magnitude at calibration samples.

Given a text file with the 51 samples (each line has one sample: actual_pressure  temperature_count  measured_pressure obtained using default compensation coefficients), i.e.
Code: [Select]
-15 5485 -15.017
-13 5485 -13.014
-11 5485 -11.011
-9 5485 -9.01
-7 5485 -7.008
-5 5485 -5.006
-3 5485 -3.004
-2 5485 -2.002
0 5485 0
2 5485 2.001
3 5485 3.002
5 5485 5.004
7 5485 7.006
9 5485 9.007
11 5485 11.008
13 5485 13.008
15 5485 15.008
-15 7386 -15
-13 7386 -13
-11 7386 -11
-9 7386 -9
-7 7386 -7
-5 7386 -5
-3 7386 -3
-2 7386 -2
0 7386 0
2 7386 2
3 7386 3
5 7386 5
7 7386 7
9 7386 9
11 7386 11
13 7386 13
15 7386 15
-15 9471 -14.99
-13 9471 -12.991
-11 9471 -10.993
-9 9471 -8.994
-7 9471 -6.996
-5 9471 -4.997
-3 9471 -2.998
-2 9471 -1.999
0 9471 0
2 9471 1.999
3 9471 2.998
5 9471 4.997
7 9471 6.996
9 9471 8.996
11 9471 10.996
13 9471 12.996
15 9471 14.996

such a program reports something like
Code: [Select]
Maximum error minimized to 0.002420993, using
0.007697997673805127 -9.893590127282349e-07 0.9937816855224063 1.353571043339379e-06 -6.784307221818112e-11

  p     |    c |     P(p, c)   |  P  |  Absolute error
--------+------+---------------+-----+-----------------
-15.017 | 5485 | -15.002188670 | -15 | -0.002188669678
-13.014 | 5485 | -13.000861278 | -13 | -0.000861277627
-11.011 | 5485 | -10.999533886 | -11 | +0.000466114424
 -9.010 | 5485 |  -9.000204823 |  -9 | -0.000204823422
 -7.008 | 5485 |  -6.999876596 |  -7 | +0.000123403681
 -5.006 | 5485 |  -4.999548369 |  -5 | +0.000451630784
 -3.004 | 5485 |  -2.999220142 |  -3 | +0.000779857886
 -2.002 | 5485 |  -1.998056864 |  -2 | +0.001943136386
  0.000 | 5485 |   0.002271363 |   0 | +0.002271363489
  2.001 | 5485 |   2.001600426 |   2 | +0.001600425643
  3.002 | 5485 |   3.001764539 |   3 | +0.001764539194
  5.004 | 5485 |   5.002092766 |   5 | +0.002092766297
  7.006 | 5485 |   7.002420993 |   7 | +0.002420993400
  9.007 | 5485 |   9.001750056 |   9 | +0.001750055554
 11.008 | 5485 |  11.001079118 |  11 | +0.001079117708
 13.008 | 5485 |  12.999409015 |  13 | -0.000590985086
 15.008 | 5485 |  14.997738912 |  15 | -0.002261087881
-15.000 | 7386 | -15.000781184 | -15 | -0.000781184012
-13.000 | 7386 | -13.000624947 | -13 | -0.000624947209
-11.000 | 7386 | -11.000468710 | -11 | -0.000468710407
 -9.000 | 7386 |  -9.000312474 |  -9 | -0.000312473605
 -7.000 | 7386 |  -7.000156237 |  -7 | -0.000156236802
 -5.000 | 7386 |  -5.000000000 |  -5 | -0.000000000000
 -3.000 | 7386 |  -2.999843763 |  -3 | +0.000156236802
 -2.000 | 7386 |  -1.999765645 |  -2 | +0.000234355203
  0.000 | 7386 |   0.000390592 |   0 | +0.000390592006
  2.000 | 7386 |   2.000546829 |   2 | +0.000546828808
  3.000 | 7386 |   3.000624947 |   3 | +0.000624947209
  5.000 | 7386 |   5.000781184 |   5 | +0.000781184012
  7.000 | 7386 |   7.000937421 |   7 | +0.000937420814
  9.000 | 7386 |   9.001093658 |   9 | +0.001093657616
 11.000 | 7386 |  11.001249894 |  11 | +0.001249894419
 13.000 | 7386 |  13.001406131 |  13 | +0.001406131221
 15.000 | 7386 |  15.001562368 |  15 | +0.001562368023
-14.990 | 9471 | -14.999404724 | -15 | +0.000595275661
-12.991 | 9471 | -12.999373552 | -13 | +0.000626447983
-10.993 | 9471 | -11.000342896 | -11 | -0.000342895540
 -8.994 | 9471 |  -9.000311723 |  -9 | -0.000311723218
 -6.996 | 9471 |  -7.001281067 |  -7 | -0.001281066740
 -4.997 | 9471 |  -5.001249894 |  -5 | -0.001249894418
 -2.998 | 9471 |  -3.001218722 |  -3 | -0.001218722096
 -1.999 | 9471 |  -2.001703394 |  -2 | -0.001703393858
  0.000 | 9471 |  -0.001672222 |   0 | -0.001672221536
  1.999 | 9471 |   1.998358951 |   2 | -0.001641049214
  2.998 | 9471 |   2.997874279 |   3 | -0.002125720975
  4.997 | 9471 |   4.997905451 |   5 | -0.002094548653
  6.996 | 9471 |   6.997936624 |   7 | -0.002063376331
  8.996 | 9471 |   8.998968312 |   9 | -0.001031688166
 10.996 | 9471 |  11.000000000 |  11 | +0.000000000000
 12.996 | 9471 |  13.001031688 |  13 | +0.001031688166
 14.996 | 9471 |  15.002063376 |  15 | +0.002063376332

i.e., if you use
    0.007697997673805127  -9.893590127282349e-07  0.9937816855224063  1.353571043339379e-06  -6.784307221818112e-11
the compensated pressure reading is within ±0.002421 of actual pressure at every calibration point.  This is just 0.008033% of the full pressure range.  In other words, the absolute error is less than 1/12000th of the full range.

I wrote such a program in C99 (or later) – edited to include a progress estimate:
Code: [Select]
// SPDX-License-Identifier: CC0-1.0
// This file is in Public Domain.
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <stdio.h>
#include <errno.h>
#include <math.h>

/* Each configuration is an array of five doubles. */
typedef struct {
    double  arg[5];
} config;

#define  ROW_FORMAT "%14.11f %18.15f %12.9f %17.14f %22.19f"

#define  ERROR_FACTOR  1.50

/* Calibration data entry */
struct sample {
    double  actual_pressure;
    double  temperature_count;
    double  measured_pressure;
};

static size_t       max_samples = 0;
static size_t       num_samples = 0;
static struct sample   *samples = NULL;

static inline void add_sample(const double actual_pressure,
                              const double temperature_count,
                              const double measured_pressure)
{
    /* Check if already known. */
    size_t  i = num_samples;
    while (i-->0)
        if (samples[i].actual_pressure == actual_pressure &&
            samples[i].temperature_count == temperature_count &&
            samples[i].measured_pressure == measured_pressure)
            return;

    /* Grow sample buffer whenever necessary. */
    if (num_samples >= max_samples) {
        const size_t  new_max = (num_samples | 1023) + 1025 - 8;
        void         *new_ptr;

        new_ptr = realloc(samples, new_max * sizeof samples[0]);
        if (!new_ptr) {
            fprintf(stderr, "Warning: Out of memory (%zu calibration samples)\n", num_samples);
            return;
        }

        max_samples = new_max;
        samples = new_ptr;
    }

    /* Add sample. */
    samples[num_samples].actual_pressure = actual_pressure;
    samples[num_samples].temperature_count = temperature_count;
    samples[num_samples].measured_pressure = measured_pressure;
    ++num_samples;
}

/* Returns the magnitude of the maximum absolute error for a given config, and the direction of the steepest descent. */
static double  error_magnitude_and_descent(const config given, config *const descent)
{
    config  gradient = {{ 0.0, 0.0, 0.0, 0.0, 0.0 }};
    double  result = 0.0;
    size_t  i;

    for (i = 0; i < num_samples; i++) {
        const double  p        = samples[i].measured_pressure;
        const double  c        = samples[i].temperature_count;
        const double  actual   = samples[i].actual_pressure;
        const double  reported = given.arg[0] + given.arg[1]*c
                               + p*(given.arg[2] + c*(given.arg[3] + c*given.arg[4]));
        const double  err      = fabs(reported - actual);

        if (result < err) {
            result = err;
            gradient.arg[0] = 1.0;
            gradient.arg[1] = c;
            gradient.arg[2] = p;
            gradient.arg[3] = p*c;
            gradient.arg[4] = p*c*c;
        }
    }

    if (descent) {
        const double  scale = -result / ( gradient.arg[0]*gradient.arg[0]
                                        + gradient.arg[1]*gradient.arg[1]
                                        + gradient.arg[2]*gradient.arg[2]
                                        + gradient.arg[3]*gradient.arg[3]
                                        + gradient.arg[4]*gradient.arg[4] );
        descent->arg[0] = gradient.arg[0] * scale;
        descent->arg[1] = gradient.arg[1] * scale;
        descent->arg[2] = gradient.arg[2] * scale;
        descent->arg[3] = gradient.arg[3] * scale;
        descent->arg[4] = gradient.arg[4] * scale;
    }

    return result;
}

/* Returns the magnitude of the maximum absolute error for given coefficients. */
static double  error_magnitude(const config  given)
{
    double  result = 0.0;
    size_t  i;

    for (i = 0; i < num_samples; i++) {
        const double  p        = samples[i].measured_pressure;
        const double  c        = samples[i].temperature_count;
        const double  actual   = samples[i].actual_pressure;
        const double  reported = given.arg[0] + given.arg[1]*c
                               + p*(given.arg[2] + c*(given.arg[3] + c*given.arg[4]));
        const double  err      = fabs(reported - actual);
        if (result < err)
            result = err;
        else
        if (result < -err)
            result = -err;
    }

    return result;
}

/* Point cache for dense minima. */
static size_t    max_points = 0;
static size_t    num_points = 0;
static config   *points = NULL;

static inline void push_point(const config point)
{
    if (num_points >= max_points) {
        const size_t  new_max = (num_points | 4095) + 4097 - 8;
        void         *new_ptr;

        new_ptr = realloc(points, new_max * sizeof points[0]);
        if (!new_ptr) {
            fprintf(stderr, "Warning: Out of memory (%zu points).\n", num_points);
            fflush(stderr);
            return;
        }

        max_points = new_max;
        points     = new_ptr;
    }

    points[num_points] = point;
    ++num_points;
}

static inline int pop_point(config *const to)
{
    if (num_points < 1)
        return 0;

    --num_points;

    if (to) *to = points[num_points];

    return 1;
}

/* Known best configuration and error magnitude. Updated after each walk. */
static config  best_point = {{ 0, 0, 0, 0, 0 }};
static double  best_error = -1.0;

/* Find local minima. */
static void  walk(config const from)
{
    config  min_point, max_point, the_point, local_point, descent;
    double  min_error, max_error, the_error, local_error, req_error;
    size_t  n, nmax;

    min_point = from;
    min_error = error_magnitude(min_point);
    push_point(min_point);

    local_point = min_point;
    local_error = min_error;

    nmax = 1;
    while (pop_point(&min_point)) {

        min_error = error_magnitude_and_descent(min_point, &descent);
        req_error = min_error;

        max_point.arg[0] = min_point.arg[0] + descent.arg[0];
        max_point.arg[1] = min_point.arg[1] + descent.arg[1];
        max_point.arg[2] = min_point.arg[2] + descent.arg[2];
        max_point.arg[3] = min_point.arg[3] + descent.arg[3];
        max_point.arg[4] = min_point.arg[4] + descent.arg[4];
        max_error = error_magnitude(max_point);

        while (1) {
            the_point.arg[0] = max_point.arg[0] + descent.arg[0];
            the_point.arg[1] = max_point.arg[1] + descent.arg[1];
            the_point.arg[2] = max_point.arg[2] + descent.arg[2];
            the_point.arg[3] = max_point.arg[3] + descent.arg[3];
            the_point.arg[4] = max_point.arg[4] + descent.arg[4];
            the_error = error_magnitude(the_point);
            if (the_error >= max_error)
                break;

            max_point = the_point;
            max_error = the_error;
        }

        if (min_error > max_error) {
            the_point = max_point;  the_error = max_error;
            max_point = min_point;  max_error = min_error;
            min_point = the_point;  min_error = the_error;
        }

        if (local_error > min_error) {
            local_error = min_error;
            local_point = min_point;
        }

        /* Binary search for minimum between curr_point and curr_point+descent. */
        n = 64;  /* Number of iterations */
        while (n-->0) {
            the_point.arg[0] = (double)(0.5*min_point.arg[0]) + (double)(0.5*max_point.arg[0]);
            the_point.arg[1] = (double)(0.5*min_point.arg[1]) + (double)(0.5*max_point.arg[1]);
            the_point.arg[2] = (double)(0.5*min_point.arg[2]) + (double)(0.5*max_point.arg[2]);
            the_point.arg[3] = (double)(0.5*min_point.arg[3]) + (double)(0.5*max_point.arg[3]);
            the_point.arg[4] = (double)(0.5*min_point.arg[4]) + (double)(0.5*max_point.arg[4]);
            the_error = error_magnitude(the_point);

            if (min_error >= the_error) {
                min_error = the_error;
                min_point = the_point;
            } else
            if (max_error >= the_error) {
                max_error = the_error;
                max_point = the_point;
            } else
            if (max_error < req_error) {
                push_point(max_point);
                max_error = the_error;
                max_point = the_point;
            }
        }

        if (nmax < num_points)
            nmax = num_points;
    }

    fprintf(stderr, "\r" ROW_FORMAT " %.9f [%zu] ",
                    local_point.arg[0], local_point.arg[1],
                    local_point.arg[2], local_point.arg[3],
                    local_point.arg[4], local_error, nmax);
    if (best_error > local_error || best_error < 0.0) {
        best_error = local_error;
        best_point = local_point;
        fprintf(stderr, "\n");
    }
    fflush(stderr);
}

static inline size_t unique5(const double v0, const double v1, const double v2, const double v3, const double v4)
{
    size_t  result = 5;

    if (v1 == v0) result--;
    if (v2 == v0 || v2 == v1) result--;
    if (v3 == v0 || v3 == v1 || v3 == v2) result--;
    if (v4 == v0 || v4 == v1 || v4 == v2 || v4 == v3) result--;

    return result;
}

/* This solves the system of equations
        a0 = X0 + X1*c0 + p0*(X2 + c0*(X3 + c0*X4))
        a1 = X1 + X1*c1 + p1*(X2 + c1*(X3 + c1*X4))
        a2 = X2 + X1*c2 + p2*(X2 + c2*(X3 + c2*X4))
        a3 = X3 + X1*c3 + p3*(X2 + c3*(X3 + c3*X4))
        a4 = X4 + X1*c4 + p4*(X2 + c4*(X3 + c4*X4))
   If a solution exists, the function returns 1, with X0,X1,X2,X3,X4 in *to.
   Otherwise, the function returns 0.
*/
static int solve(const double p0, const double p1, const double p2, const double p3, const double p4,
                 const double c0, const double c1, const double c2, const double c3, const double c4,
                 const double a0, const double a1, const double a2, const double a3, const double a4,
                 config *const to)
{
    config  result;

    const double c00 = c0 * c0, c01 = c0 * c1, c02 = c0 * c2, c03 = c0 * c3, c04 = c0 * c4,
                 c11 = c1 * c1, c12 = c1 * c2, c13 = c1 * c3, c14 = c1 * c4,
                 c22 = c2 * c2, c23 = c2 * c3, c24 = c2 * c4,
                 c33 = c3 * c3, c34 = c3 * c4,
                 c44 = c4 * c4;

    const double c001 = c00 * c1,
                 c002 = c00 * c2,
                 c003 = c00 * c3,
                 c004 = c00 * c4,
                 c011 = c11 * c0,
                 c022 = c22 * c0,
                 c033 = c33 * c0,
                 c044 = c44 * c0,
                 c112 = c11 * c2,
                 c113 = c11 * c3,
                 c114 = c11 * c4,
                 c122 = c22 * c1,
                 c133 = c33 * c1,
                 c144 = c44 * c1,
                 c223 = c22 * c3,
                 c224 = c22 * c4,
                 c233 = c33 * c2,
                 c244 = c44 * c2,
                 c334 = c33 * c4,
                 c344 = c44 * c3;
    const double c0012 = c00 * c12,
                 c0013 = c00 * c13,
                 c0014 = c00 * c14,
                 c0023 = c00 * c23,
                 c0024 = c00 * c24,
                 c0034 = c00 * c34,
                 c0112 = c11 * c02,
                 c0113 = c11 * c03,
                 c0114 = c11 * c04,
                 c0122 = c01 * c22,
                 c0133 = c01 * c33,
                 c0144 = c01 * c44,
                 c0223 = c03 * c22,
                 c0224 = c04 * c22,
                 c0233 = c02 * c33,
                 c0244 = c02 * c44,
                 c0334 = c04 * c33,
                 c0344 = c03 * c44,
                 c1123 = c11 * c23,
                 c1124 = c11 * c24,
                 c1134 = c11 * c34,
                 c1223 = c13 * c22,
                 c1224 = c14 * c22,
                 c1233 = c12 * c33,
                 c1244 = c12 * c44,
                 c1334 = c14 * c33,
                 c1344 = c13 * c44,
                 c2234 = c22 * c34,
                 c2334 = c24 * c33,
                 c2344 = c23 * c44;

    const double p01 = p0 * p1, p02 = p0 * p2, p03 = p0 * p3, p04 = p0 * p4,
                 p12 = p1 * p2, p13 = p1 * p3, p14 = p1 * p4,
                 p23 = p2 * p3, p24 = p2 * p4,
                 p34 = p3 * p4;
    const double p012 = p01 * p2,
                 p013 = p01 * p3,
                 p014 = p01 * p3,
                 p023 = p02 * p3,
                 p024 = p02 * p4,
                 p034 = p03 * p4,
                 p123 = p12 * p3,
                 p124 = p12 * p4,
                 p134 = p13 * p4,
                 p234 = p23 * p4;

    const double divisor = (c1224 - c0224 - c1124 + c0024 + c0114 - c0014 - c1223 + c0223 + c1123 - c0023 - c0113 + c0013) * p012
                         - (c0012 + c1334 - c0334 - c1134 + c0034 + c0114 - c0014 - c1233 + c0233 + c1123 - c0023 - c0112) * p013
                         + (c0012 + c1344 - c0344 - c1244 + c0244 - c1134 + c0034 + c1124 - c0024 + c0113 - c0013 - c0112) * p014
                         + (c0012 + c2334 - c0334 - c2234 + c0034 + c0224 - c0024 - c1233 + c0133 + c1223 - c0013 - c0122) * p023
                         - (c0012 + c2344 - c0344 - c1244 + c0144 - c2234 + c0034 + c1224 - c0014 + c0223 - c0023 - c0122) * p024
                         + (c2344 - c1344 - c0244 + c0144 - c2334 + c1334 + c0024 - c0014 + c0233 - c0133 - c0023 + c0013) * p034
                         - (c2334 - c1334 - c2234 + c1134 + c1224 - c1124 - c0233 + c0133 + c0223 - c0113 - c0122 + c0112) * p123
                         + (c2344 - c1344 - c0244 + c0144 - c2234 + c1134 + c0224 - c0114 + c1223 - c1123 - c0122 + c0112) * p124
                         - (c2344 - c0344 - c1244 + c0144 - c2334 + c0334 + c1124 - c0114 + c1233 - c0133 - c1123 + c0113) * p134
                         + (c1344 - c0344 - c1244 + c0244 - c1334 + c0334 + c1224 - c0224 + c1233 - c0233 - c1223 + c0223) * p234;

    if (divisor == 0.0)
        return 0;

    result.arg[0] = ( (  a4*c0013 - a3*c0014 - a4*c0023 + a3*c0024 - a4*c0113 + a3*c0114 + a4*c0223 - a3*c0224 + a4*c1123 - a3*c1124 - a4*c1223 + a3*c1224) * p012
                    + ( -a4*c0012 + a2*c0014 + a4*c0023 - a2*c0034 + a4*c0112 - a2*c0114 - a4*c0233 + a2*c0334 - a4*c1123 + a2*c1134 + a4*c1233 - a2*c1334) * p013
                    + (  a3*c0012 - a2*c0013 - a3*c0024 + a2*c0034 - a3*c0112 + a2*c0113 + a3*c0244 - a2*c0344 + a3*c1124 - a2*c1134 - a3*c1244 + a2*c1344) * p014
                    + (  a4*c0012 - a4*c0013 - a1*c0024 + a1*c0034 - a4*c0122 + a4*c0133 + a1*c0224 - a1*c0334 + a4*c1223 - a4*c1233 - a1*c2234 + a1*c2334) * p023
                    + ( -a3*c0012 + a3*c0014 + a1*c0023 - a1*c0034 + a3*c0122 - a3*c0144 - a1*c0223 + a1*c0344 - a3*c1224 + a3*c1244 + a1*c2234 - a1*c2344) * p024
                    + (  a2*c0013 - a2*c0014 - a1*c0023 + a1*c0024 - a2*c0133 + a2*c0144 + a1*c0233 - a1*c0244 + a2*c1334 - a2*c1344 - a1*c2334 + a1*c2344) * p034
                    + ( -a4*c0112 + a4*c0113 + a4*c0122 - a4*c0133 - a4*c0223 + a4*c0233 + a0*c1124 - a0*c1134 - a0*c1224 + a0*c1334 + a0*c2234 - a0*c2334) * p123
                    + (  a3*c0112 - a3*c0114 - a3*c0122 + a3*c0144 + a3*c0224 - a3*c0244 - a0*c1123 + a0*c1134 + a0*c1223 - a0*c1344 - a0*c2234 + a0*c2344) * p124
                    + ( -a2*c0113 + a2*c0114 + a2*c0133 - a2*c0144 - a2*c0334 + a2*c0344 + a0*c1123 - a0*c1124 - a0*c1233 + a0*c1244 + a0*c2334 - a0*c2344) * p134
                    + (  a1*c0223 - a1*c0224 - a1*c0233 + a1*c0244 + a1*c0334 - a1*c0344 - a0*c1223 + a0*c1224 + a0*c1233 - a0*c1244 - a0*c1334 + a0*c1344) * p234
                    ) / (divisor);

    result.arg[1] = ( ( a3*c001 - a4*c001 + a4*c002 - a3*c002 + a4*c011 - a3*c011 + a3*c022 - a4*c022 + a3*c112 - a4*c112 + a4*c122 - a3*c122) * p012
                    + ( a4*c001 - a2*c001 + a2*c003 - a4*c003 + a2*c011 - a4*c011 + a4*c033 - a2*c033 + a4*c113 - a2*c113 + a2*c133 - a4*c133) * p013
                    + ( a2*c001 - a3*c001 + a3*c004 - a2*c004 + a3*c011 - a2*c011 + a2*c044 - a3*c044 + a2*c114 - a3*c114 + a3*c144 - a2*c144) * p014
                    + ( a1*c002 - a4*c002 + a4*c003 - a1*c003 + a4*c022 - a1*c022 + a1*c033 - a4*c033 + a1*c223 - a4*c223 + a4*c233 - a1*c233) * p023
                    + ( a3*c002 - a1*c002 + a1*c004 - a3*c004 + a1*c022 - a3*c022 + a3*c044 - a1*c044 + a3*c224 - a1*c224 + a1*c244 - a3*c244) * p024
                    + ( a1*c003 - a2*c003 + a2*c004 - a1*c004 + a2*c033 - a1*c033 + a1*c044 - a2*c044 + a1*c334 - a2*c334 + a2*c344 - a1*c344) * p034
                    + ( a4*c112 - a0*c112 + a0*c113 - a4*c113 + a0*c122 - a4*c122 + a4*c133 - a0*c133 + a4*c223 - a0*c223 + a0*c233 - a4*c233) * p123
                    + ( a0*c112 - a3*c112 + a3*c114 - a0*c114 + a3*c122 - a0*c122 + a0*c144 - a3*c144 + a0*c224 - a3*c224 + a3*c244 - a0*c244) * p124
                    + ( a2*c113 - a0*c113 + a0*c114 - a2*c114 + a0*c133 - a2*c133 + a2*c144 - a0*c144 + a2*c334 - a0*c334 + a0*c344 - a2*c344) * p134
                    + ( a0*c223 - a1*c223 + a1*c224 - a0*c224 + a1*c233 - a0*c233 + a0*c244 - a1*c244 + a0*c334 - a1*c334 - a0*c344 + a1*c344) * p234
                    ) / (divisor);

    result.arg[2] = ( ( a4*c0012 - a3*c0012 + a2*c0013 - a4*c0013 + a3*c0014 - a2*c0014 + a3*c0112 - a4*c0112 + a4*c0113 - a2*c0113 + a2*c0114 - a3*c0114) * p01
                    + ( a3*c0012 - a4*c0012 + a4*c0023 - a1*c0023 + a1*c0024 - a3*c0024 + a4*c0122 - a3*c0122 + a1*c0223 - a4*c0223 + a3*c0224 - a1*c0224) * p02
                    + ( a4*c0013 - a2*c0013 + a1*c0023 - a4*c0023 + a2*c0034 - a1*c0034 + a2*c0133 - a4*c0133 + a4*c0233 - a1*c0233 + a1*c0334 - a2*c0334) * p03
                    + ( a2*c0014 - a3*c0014 + a3*c0024 - a1*c0024 + a1*c0034 - a2*c0034 + a3*c0144 - a2*c0144 + a1*c0244 - a3*c0244 + a2*c0344 - a1*c0344) * p04
                    + ( a4*c0112 - a3*c0112 + a3*c0122 - a4*c0122 + a0*c1123 - a4*c1123 + a3*c1124 - a0*c1124 + a4*c1223 - a0*c1223 + a0*c1224 - a3*c1224) * p12
                    + ( a2*c0113 - a4*c0113 + a4*c0133 - a2*c0133 + a4*c1123 - a0*c1123 + a0*c1134 - a2*c1134 + a0*c1233 - a4*c1233 + a2*c1334 - a0*c1334) * p13
                    + ( a3*c0114 - a2*c0114 + a2*c0144 - a3*c0144 + a0*c1124 - a3*c1124 + a2*c1134 - a0*c1134 + a3*c1244 - a0*c1244 + a0*c1344 - a2*c1344) * p14
                    + ( a4*c0223 - a1*c0223 + a1*c0233 - a4*c0233 + a0*c1223 - a4*c1223 + a4*c1233 - a0*c1233 + a1*c2234 - a0*c2234 + a0*c2334 - a1*c2334) * p23
                    + ( a1*c0224 - a3*c0224 + a3*c0244 - a1*c0244 + a3*c1224 - a0*c1224 + a0*c1244 - a3*c1244 + a0*c2234 - a1*c2234 + a1*c2344 - a0*c2344) * p24
                    + ( a2*c0334 - a1*c0334 + a1*c0344 - a2*c0344 + a0*c1334 - a2*c1334 + a2*c1344 - a0*c1344 + a1*c2334 - a0*c2334 + a0*c2344 - a1*c2344) * p34
                    ) / (divisor);

    result.arg[3] = ( ( a3*c002 - a4*c002 + a4*c003 - a2*c003 + a2*c004 - a3*c004 + a4*c112 - a3*c112 + a2*c113 - a4*c113 + a3*c114 - a2*c114) * p01
                    + ( a4*c001 - a3*c001 + a1*c003 - a4*c003 + a3*c004 - a1*c004 + a3*c122 - a4*c122 + a4*c223 - a1*c223 + a1*c224 - a3*c224) * p02
                    + ( a2*c001 - a4*c001 + a4*c002 - a1*c002 + a1*c004 - a2*c004 + a4*c133 - a2*c133 + a1*c233 - a4*c233 + a2*c334 - a1*c334) * p03
                    + ( a3*c001 - a2*c001 + a1*c002 - a3*c002 + a2*c003 - a1*c003 + a2*c144 - a3*c144 + a3*c244 - a1*c244 + a1*c344 - a2*c344) * p04
                    + ( a3*c011 - a4*c011 + a4*c022 - a3*c022 + a4*c113 - a0*c113 + a0*c114 - a3*c114 + a0*c223 - a4*c223 + a3*c224 - a0*c224) * p12
                    + ( a4*c011 - a2*c011 + a2*c033 - a4*c033 + a0*c112 - a4*c112 + a2*c114 - a0*c114 + a4*c233 - a0*c233 + a0*c334 - a2*c334) * p13
                    + ( a2*c011 - a3*c011 + a3*c044 - a2*c044 + a3*c112 - a0*c112 + a0*c113 - a2*c113 + a0*c244 - a3*c244 + a2*c344 - a0*c344) * p14
                    + ( a1*c022 - a4*c022 + a4*c033 - a1*c033 + a4*c122 - a0*c122 + a0*c133 - a4*c133 + a0*c224 - a1*c224 + a1*c334 - a0*c334) * p23
                    + ( a3*c022 - a1*c022 + a1*c044 - a3*c044 + a0*c122 - a3*c122 + a3*c144 - a0*c144 + a1*c223 - a0*c223 + a0*c344 - a1*c344) * p24
                    + ( a1*c033 - a2*c033 + a2*c044 - a1*c044 + a2*c133 - a0*c133 + a0*c144 - a2*c144 + a0*c233 - a1*c233 - a0*c244 + a1*c244) * p34
                    ) / (divisor);

    result.arg[4] = ( ( a4*c02 - a3*c02 + a2*c03 - a4*c03 + a3*c04 - a2*c04 + a3*c12 - a4*c12 + a4*c13 - a2*c13 + a2*c14 - a3*c14) * p01
                    + ( a3*c01 - a4*c01 + a4*c03 - a1*c03 + a1*c04 - a3*c04 + a4*c12 - a3*c12 + a1*c23 - a4*c23 + a3*c24 - a1*c24) * p02
                    + ( a4*c01 - a2*c01 + a1*c02 - a4*c02 + a2*c04 - a1*c04 + a2*c13 - a4*c13 + a4*c23 - a1*c23 + a1*c34 - a2*c34) * p03
                    + ( a2*c01 - a3*c01 + a3*c02 - a1*c02 + a1*c03 - a2*c03 + a3*c14 - a2*c14 + a1*c24 - a3*c24 + a2*c34 - a1*c34) * p04
                    + ( a4*c01 - a3*c01 + a3*c02 - a4*c02 + a0*c13 - a4*c13 + a3*c14 - a0*c14 + a4*c23 - a0*c23 + a0*c24 - a3*c24) * p12
                    + ( a2*c01 - a4*c01 + a4*c03 - a2*c03 + a4*c12 - a0*c12 + a0*c14 - a2*c14 + a0*c23 - a4*c23 + a2*c34 - a0*c34) * p13
                    + ( a3*c01 - a2*c01 + a2*c04 - a3*c04 + a0*c12 - a3*c12 + a2*c13 - a0*c13 + a3*c24 - a0*c24 + a0*c34 - a2*c34) * p14
                    + ( a4*c02 - a1*c02 + a1*c03 - a4*c03 + a0*c12 - a4*c12 + a4*c13 - a0*c13 + a1*c24 - a0*c24 + a0*c34 - a1*c34) * p23
                    + ( a1*c02 - a3*c02 + a3*c04 - a1*c04 + a3*c12 - a0*c12 + a0*c14 - a3*c14 + a0*c23 - a1*c23 + a1*c34 - a0*c34) * p24
                    + ( a2*c03 - a1*c03 + a1*c04 - a2*c04 + a0*c13 - a2*c13 + a2*c14 - a0*c14 + a1*c23 - a0*c23 + a0*c24 - a1*c24) * p34
                    ) / (divisor);

    if (to) *to = result;
    return 1;
}

static size_t parse_fields(const char *line, double *values, size_t count)
{
    const char *ends;
    double      value;
    size_t      n = 0;

    while (1) {

        ends = line;
        errno = 0;
        value = strtod(line, (char **)(&ends));
        if (errno || ends == line)
            return n;

        if (*ends == '\t' || *ends == '\n' || *ends == '\v' ||
            *ends == '\f' || *ends == '\r' || *ends == ' ') {
            line = ends + 1;
            if (n < count)
                values[n] = value;
            n++;

        } else
            return n;
    }
}


int main(int argc, char *argv[])
{
    double   nmax;
    uint64_t n;
    size_t   i;
    int      arg;

    if (argc < 2 || !strcmp(argv[1], "-h") || !strcmp(argv[1], "--help") || !strcmp(argv[1], "/?")) {
        const char *argv0 = (argc > 0 && argv && argv[0]) ? argv[0] : "(this)";
        fprintf(stderr, "\n");
        fprintf(stderr, "Usage: %s [ -h | --help | /? ]\n", argv0);
        fprintf(stderr, "       %s INPUT-FILE [ INPUT-FILE ]\n", argv0);
        fprintf(stderr, "\n");
        fprintf(stderr, "Each input file should contain rows of form\n");
        fprintf(stderr, "    ACTUAL-PRESSURE  TEMPERATURE-COUNT  MEASURED-PRESSURE\n");
        fprintf(stderr, "and nothing else.\n");
        fprintf(stderr, "\n");
        fprintf(stderr, "There must be at least five unique rows,\n");
        fprintf(stderr, "with at least three unique actual pressures,\n");
        fprintf(stderr, "and at least three unique temperature counts.\n");
        fprintf(stderr, "The more samples, the better.\n");
        fprintf(stderr, "\n");
        return EXIT_SUCCESS;
    }

    for (arg = 1; arg < argc; arg++) {
        double field[3];
        char   buffer[1024], *line;
       
        FILE  *in;

        in = fopen(argv[arg], "r");
        if (!in) {
            fprintf(stderr, "%s: %s.\n", argv[arg], strerror(errno));
            return EXIT_FAILURE;
        }

        while (1) {
            line = fgets(buffer, sizeof buffer, in);
            if (!line)
                break;

            /* Skip leading whitespace. */
            line += strspn(line, "\t\n\v\f\r ");

            /* Skip lines beginning with # or ;. */
            if (*line == '#' || *line == ';')
                continue;

            /* Parse three numeric fields. */
            if (parse_fields(line, field, 3) != 3) {
                /* Remove trailing whitespace. */
                line[strcspn(line, "\t\n\v\f\r ")] = '\0';
                fprintf(stderr, "%s: Invalid input: \"%s\".\n", argv[arg], line);
                return EXIT_FAILURE;
            }

            /* Add as a new sample. */
            add_sample(field[0], field[1], field[2]);
        }

        if (!feof(in) || ferror(in)) {
            fclose(in);
            fprintf(stderr, "%s: Read error.\n", argv[arg]);
            return EXIT_FAILURE;
        } else
        if (fclose(in)) {
            fprintf(stderr, "%s: Read error.\n", argv[arg]);
            return EXIT_FAILURE;
        }
    }

    if (num_samples < 5) {
        fprintf(stderr, "Not enough calibration samples.\n");
        return EXIT_FAILURE;
    }

    fprintf(stderr, "Read %zu calibration samples.\n", num_samples);

    nmax = (double)(num_samples    )
         * (double)(num_samples - 1)
         * (double)(num_samples - 2)
         * (double)(num_samples - 3)
         * (double)(num_samples - 4)
         * 0.66; /* Approximate, depends on sample distribution -- good'enuf */
    n = 0;

    {
        config point;
        size_t i0, i1, i2, i3, i4;

        for (i0 = 0; i0 < num_samples; i0++) {
            const double p0 = samples[i0].measured_pressure;
            const double c0 = samples[i0].temperature_count;
            const double a0 = samples[i0].actual_pressure;
            for (i1 = 0; i1 < num_samples; i1++) {
                const double p1 = samples[i1].measured_pressure;
                const double c1 = samples[i1].temperature_count;
                const double a1 = samples[i1].actual_pressure;
                if (i1 == i0)
                    continue;
                for (i2 = 0; i2 < num_samples; i2++) {
                    const double p2 = samples[i2].measured_pressure;
                    const double c2 = samples[i2].temperature_count;
                    const double a2 = samples[i2].actual_pressure;
                    if (i2 == i0 || i2 == i1)
                        continue;
                    for (i3 = 0; i3 < num_samples; i3++) {
                        const double p3 = samples[i3].measured_pressure;
                        const double c3 = samples[i3].temperature_count;
                        const double a3 = samples[i3].actual_pressure;
                        if (i3 == i0 || i3 == i1 || i3 == i2)
                            continue;
                        for (i4 = 0; i4 < num_samples; i4++) {
                            const double p4 = samples[i4].measured_pressure;
                            const double c4 = samples[i4].temperature_count;
                            const double a4 = samples[i4].actual_pressure;
                            if (i4 == i0 || i4 == i1 || i4 == i2 || i4 == i3)
                                continue;
                            if (unique5(a0, a1, a2, a3, a4) < 3)
                                continue;
                            if (unique5(c0, c1, c2, c3, c4) < 3)
                                continue;
                            if (unique5(p0, p1, p2, p3, p4) < 3)
                                continue;
                            n++;
                            if (solve(p0, p1, p2, p3, p4, c0, c1, c2, c3, c4, a0, a1, a2, a3, a4, &point)) {
                                if (best_error < 0 || error_magnitude(point) < ERROR_FACTOR * best_error)
                                   walk(point);
                            }
                        }
                    }
                }
            }
            fprintf(stderr, "\r[ Estimating %.3f%% complete ]\n", (double)n * 100.0 / nmax);
            fflush(stderr);
        }
    }
    fprintf(stderr, "\n");
    fflush(stderr);

    printf("  p     |    c |     P(p, c)   |  P  |  Absolute error\n");
    printf("--------+------+---------------+-----+-----------------\n");
    for (i = 0; i < num_samples; i++) {
        const double  p = samples[i].measured_pressure;
        const double  c = samples[i].temperature_count;
        const double  actual = samples[i].actual_pressure;
        const double  r = best_point.arg[0] + best_point.arg[1]*c
                        + p*(best_point.arg[2] + c*(best_point.arg[3] + best_point.arg[4]*c));
        printf("%7.3f | %4.0f | %13.9f | %3.0f | %+15.12f\n", p, c, r, actual, r - actual);
    }
    printf("\n");
    printf("Maximum error minimized to %.9f, using\n", best_error);
    printf("%.16g %.16g %.16g %.16g %.16g\n",
            best_point.arg[0], best_point.arg[1], best_point.arg[2],
            best_point.arg[3], best_point.arg[4]);

    return EXIT_SUCCESS;
}

Note that this does not quarantee that the absolute error is that or smaller over the entire measurement range; this minimizes the absolute error (magnitude) at calibration sample points only.  It's up to you to provide enough samples to cover the device operational range.
« Last Edit: August 21, 2020, 02:56:44 pm by Nominal Animal »
 
The following users thanked this post: sorin, itsbiodiversity

Offline itsbiodiversityTopic starter

  • Regular Contributor
  • *
  • Posts: 99
  • Country: us
you are a MAD PERSON!  Wow.  Thanks so much for sticking with me on this.
 

Offline IanB

  • Super Contributor
  • ***
  • Posts: 12539
  • Country: us
Those are the coefficients NOMINAL ANIMAL. Exactly. Now I need to figure out how to calculate those myself based on any dataset using excel.  I tried to download Gnuplot but had some issues installing.

As I mentioned in an earlier post, this is easily done using the Excel solver add-in (it is included with every version of Excel, but it does not show up on ribbon by default--you have to add it).

With the solver add-in you can say "vary these cells to minimize the value of that cell". It becomes quite easy to do all sorts of regressions that way.
 

Offline Nominal Animal

  • Super Contributor
  • ***
  • Posts: 7198
  • Country: fi
    • My home page and email address
If you want to explore the gradient descent min-max optimizer, you can contact me via email, or continue discussion here.
(I'm only interested in this as an optimization problem, not looking to get anything.)
I noticed that I can get well over 50% acceleration using explicit AVX2 intrinsics; also a thread pool approach should work really well here.
My limitation is that I don't do Windows-only code.

In particular, if you have a calibration setup where you have more than one known good setup measuring the same thing as the device being calibrated, you could probably make the entire calibration process continuous.
 

Offline b_force

  • Super Contributor
  • ***
  • Posts: 1381
  • Country: 00
    • One World Concepts
Noble amount of labor to figure all these things out  :-+

But I don't want to be much of a party pooper, but there are quite some programs available that can do very complicated curve fitting.

Offline Nominal Animal

  • Super Contributor
  • ***
  • Posts: 7198
  • Country: fi
    • My home page and email address
Note that the gradient descent optimization is not curve fitting, it is actually an optimizer that works on the full parameter space, finding the local optimum (minimizing maximum error to calibration samples) from a given starting point.

Mathematically, we have a function \$f(a_0, a_1, a_2, a_3, a_4)\$ that provides the maximum error among the calibration points when the compensation function for measured pressure \$p\$ and temperature count \$c\$ is \$P(p, c) = a_0 + a_1 c + a_2 p + a_3 p c + a_4 p c^2\$.  The program assumes that only one calibration point causes the maximum error, and uses the gradient of that point – i.e., \$\left(\frac{\partial f}{\partial a_0}, \frac{\partial f}{\partial a_1}, \frac{\partial f}{\partial a_2}, \frac{\partial f}{\partial a_3}, \frac{\partial f}{\partial a_4} \right)\$ – to determine the direction in which to move in the parameter space.  The length of the step is first maximized (to the maximum length that provides a smaller maximum error than the starting point), then a binary search along this hyperline segment is used to find the \$(a_0, a_1, a_2, a_3, a_4)\$ that yields the smallest maximum error along this hyperline segment.

The program does a very, very large number of these descents.  It uses each unique quintet among the calibration samples (each unique 5-tuple of calibration points) to use the \$\left(\begin{matrix}N\\5\end{matrix}\right) = N(N-1)(N-2)(N-3)(N-4)/(1\cdot 2 \cdot 3 \cdot 4 \cdot 5)\$ starting points for the descents.  (Because of numerical instability of the naive implementation of the 5-point solver, it actually uses each of the 120 variants.)

In essence, the curve \$P(p, c) = a_0 + a_1 c + a_2 p + a_3 p c + a_4 p c^2\$ is only used with sets of 5 calibration samples, and then only to get useful starting points for the gradient descent.  There is no fitting at all, because the set of 5 samples define the five unknowns exactly.  (The code implements this algebraic solver itself.)

To get similar results from curve fitting, you use the maximum error in magnitude, \$\max\left(\left\lvert P(p_i,c_i) - A_i\right\rvert\right)\$ across the \$i\$ calibration samples (\$A_i\$ being the actual pressure being measured, \$p_i\$ the corresponding uncompensated measured pressure, and \$c_i\$ the temperature count).  But, since \$\max\left(\left\lvert P(p_i,c_i) - A_i\right\rvert\right)\$ is nonlinear (piecewise quadratic with a cross term), with a non-contiguous second derivative (because the \$i\$ that causes the maximum error varies), with LOTS of local minima (thousands if not more), normal methods (including anything based on linear algebra) will not produce as good results as this gradient descent does.

No need to believe me, though.  Just use the curve fitting package you think does better or even just as well, and see if it discovers the minimum I reported a couple of messages back, that minimizes the maximum error among the calibration samples to 0.002421 pressure units.  In comparison, least squares fitting only reaches a minimum error a magnitude larger (around 0.02, depending on details).

(Actually, if you can, please do; I would appreciate the results.  The ones I have access to – gnuplot, maxima, sagemath, maple, and mathematica – all use least squares fitting, which produces the same minimum error, about 0.02; or an order of magnitude worse maximum error than the gradient descent does.)
« Last Edit: August 26, 2020, 02:15:52 pm by Nominal Animal »
 


Share me

Digg  Facebook  SlashDot  Delicious  Technorati  Twitter  Google  Yahoo
Smf