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.
-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
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-11the 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:
// 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.