Electronics > Projects, Designs, and Technical Stuff

4th order polynomial coefficients for pressure at temperature readings, Help!!!

<< < (9/12) > >>

itsbiodiversity:
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.

Nominal Animal:
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: ---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

--- End code ---
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: ---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

--- End code ---
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!)

itsbiodiversity:
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!

itsbiodiversity:
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.

Nominal Animal:
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: ----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

--- End code ---

such a program reports something like

--- Code: ---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

--- End code ---

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: ---// 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;
}

--- End code ---

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.

Navigation

[0] Message Index

[#] Next page

[*] Previous page

There was an error while thanking
Thanking...
Go to full version
Powered by SMFPacks Advanced Attachments Uploader Mod