From 10d3e1b75433740fc07e15c787ab19c1e9469dba Mon Sep 17 00:00:00 2001 From: Giles Atkinson <“gatk555@gmail.com”> Date: Fri, 15 Jul 2022 19:17:49 +0100 Subject: [PATCH] Extend the output and partial derivatives of the output from the 2/3-D table lookup models (XSPICE) continuously outside the defined rectangle or cube, as it may help convergence and avoid ambiguity caused by rounding at the boundaries. See bug #591 "Wrong 2D table model output". The value is that of the nearest point in the region while the reported derivative is ramped to zero outside the boundaries. Also fixes incorrect output at the upper limit in each dimension. --- src/xspice/icm/table/table2D/cfunc.mod | 68 +++++++++++++----- src/xspice/icm/table/table3D/cfunc.mod | 95 ++++++++++++++++++++------ 2 files changed, 124 insertions(+), 39 deletions(-) diff --git a/src/xspice/icm/table/table2D/cfunc.mod b/src/xspice/icm/table/table2D/cfunc.mod index cc771692a..a5375c3e7 100644 --- a/src/xspice/icm/table/table2D/cfunc.mod +++ b/src/xspice/icm/table/table2D/cfunc.mod @@ -245,14 +245,19 @@ x0yiy-1 x1yiy-1 x2yiy-1 ... xix-1yiy-1 ==============================================================================*/ +/* How quickly partial derivative ramps to zero at boundary. */ + +#define RAMP_WIDTH 0.125 /*=== CM_table2D ROUTINE ===*/ +static const char exceed_fmt[] = "%c value %g exceeds table limits,\n" + " please enlarge range of your table."; void cm_table2D(ARGS) /* structure holding parms, inputs, outputs, etc. */ { int size, xind, yind; - double xval, yval, xoff, yoff, xdiff, ydiff; + double xval, yval, xoff, yoff, xdiff, ydiff, xramp, yramp; double derivval[2], outval; Table2_Data_t *loc; /* Pointer to local static data, not to be included @@ -280,23 +285,50 @@ void cm_table2D(ARGS) /* structure holding parms, inputs, outputs, etc. */ /* check table ranges */ if (xval < loc->xcol[0] || xval > loc->xcol[loc->ix - 1]) { - if (PARAM(verbose) > 0) { - cm_message_printf("x value %g exceeds table limits,\n" - " please enlarge range of your table", - xval); + /* x input out of region: use nearest point in region, ramping + * partial derivative to zero at edge. + */ + + if (PARAM(verbose) > 0) + cm_message_printf(exceed_fmt, 'x', xval); + if (xval < loc->xcol[0]) { + xramp = 1 - ((loc->xcol[0] - xval) / + (RAMP_WIDTH * (loc->xcol[1] - loc->xcol[0]))); + if (xramp < 0.0) + xramp = 0.0; + xval = loc->xcol[0]; + } else { + xramp = 1 - ((xval - loc->xcol[loc->ix - 1]) / + (RAMP_WIDTH * (loc->xcol[loc->ix - 1] - + loc->xcol[loc->ix - 2]))); + if (xramp < 0.0) + xramp = 0.0; + xval = loc->xcol[loc->ix - 1]; } - return; + } else { + xramp = 1.0; } if (yval < loc->ycol[0] || yval > loc->ycol[loc->iy - 1]) { - if (PARAM(verbose) > 0) { - cm_message_printf("y value %g exceeds table limits,\n" - " please enlarge range of your table", - yval); + if (PARAM(verbose) > 0) + cm_message_printf(exceed_fmt, 'y', yval); + if (yval < loc->ycol[0]) { + yramp = 1 - ((loc->ycol[0] - yval) / + (RAMP_WIDTH * (loc->ycol[1] - loc->ycol[0]))); + if (yramp < 0.0) + yramp = 0.0; + yval = loc->ycol[0]; + } else { + yramp = 1 - ((yval - loc->ycol[loc->iy - 1]) / + (RAMP_WIDTH * (loc->ycol[loc->iy - 1] - + loc->ycol[loc->iy - 2]))); + if (yramp < 0.0) + yramp = 0.0; + yval = loc->ycol[loc->iy - 1]; } - return; + } else { + yramp = 1.0; } - /*** find indices where interpolation will be done ***/ /* something like binary search to get the index */ xind = findCrossOver(loc->xcol, loc->ix, xval); @@ -319,12 +351,14 @@ void cm_table2D(ARGS) /* structure holding parms, inputs, outputs, etc. */ DER /* what to compute [FUNC, DER, BOTH] */ ); - /* xind and yind may become too large */ + /* xind may become too large when xval == xcol[loc->ix - 1] */ if (xind == loc->ix - 1) { --xind; + xoff += loc->xcol[xind + 1] - loc->xcol[xind]; } if (yind == loc->iy - 1) { --yind; + yoff += loc->ycol[yind + 1] - loc->ycol[yind]; } /* Overwrite outval from sf_eno2_apply by bilinear interpolation */ @@ -337,9 +371,9 @@ void cm_table2D(ARGS) /* structure holding parms, inputs, outputs, etc. */ double xderiv, yderiv, outv; outv = PARAM(offset) + PARAM(gain) * outval; OUTPUT(out) = outv; - xderiv = PARAM(gain) * derivval[0] / xdiff; + xderiv = xramp * PARAM(gain) * derivval[0] / xdiff; PARTIAL(out, inx) = xderiv; - yderiv = PARAM(gain) * derivval[1] / ydiff; + yderiv = yramp * PARAM(gain) * derivval[1] / ydiff; PARTIAL(out, iny) = yderiv; if (PARAM(verbose) > 1) { @@ -350,10 +384,10 @@ void cm_table2D(ARGS) /* structure holding parms, inputs, outputs, etc. */ } else { Mif_Complex_t ac_gain; - ac_gain.real = PARAM(gain) * derivval[0] / xdiff; + ac_gain.real = xramp * PARAM(gain) * derivval[0] / xdiff; ac_gain.imag= 0.0; AC_GAIN(out, inx) = ac_gain; - ac_gain.real = PARAM(gain) * derivval[1] / ydiff; + ac_gain.real = yramp * PARAM(gain) * derivval[1] / ydiff; ac_gain.imag= 0.0; AC_GAIN(out, iny) = ac_gain; } diff --git a/src/xspice/icm/table/table3D/cfunc.mod b/src/xspice/icm/table/table3D/cfunc.mod index 751637622..2792cfdcc 100644 --- a/src/xspice/icm/table/table3D/cfunc.mod +++ b/src/xspice/icm/table/table3D/cfunc.mod @@ -247,10 +247,18 @@ x0yiy-1 x1yiy-1 x2yiy-1 ... xix-1yiy-1 /*=== CM_table3D ROUTINE ===*/ +/* How quickly partial derivative ramps to zero at boundary. */ + +#define RAMP_WIDTH 0.125 + +static const char exceed_fmt[] = "%c value %g exceeds table limits,\n" + " please enlarge range of your table."; + void cm_table3D(ARGS) /* structure holding parms, inputs, outputs, etc. */ { int size, xind, yind, zind; double xval, yval, zval, xoff, yoff, zoff, xdiff, ydiff, zdiff; + double xramp, yramp, zramp; double derivval[3], outval; Table3_Data_t *loc; /* Pointer to local static data, not to be included @@ -279,28 +287,68 @@ void cm_table3D(ARGS) /* structure holding parms, inputs, outputs, etc. */ /* check table ranges */ if (xval < loc->xcol[0] || xval > loc->xcol[loc->ix - 1]) { - if (PARAM(verbose) > 0) { - cm_message_printf("x value %g exceeds table limits,\n" - " please enlarge range of your table", - xval); + /* x input out of region: use nearest point in region, ramping + * partial derivative to zero at edge. + */ + + if (PARAM(verbose) > 0) + cm_message_printf(exceed_fmt, 'x', xval); + if (xval < loc->xcol[0]) { + xramp = 1 - ((loc->xcol[0] - xval) / + (RAMP_WIDTH * (loc->xcol[1] - loc->xcol[0]))); + if (xramp < 0.0) + xramp = 0.0; + xval = loc->xcol[0]; + } else { + xramp = 1 - ((xval - loc->xcol[loc->ix - 1]) / + (RAMP_WIDTH * (loc->xcol[loc->ix - 1] - + loc->xcol[loc->ix - 2]))); + if (xramp < 0.0) + xramp = 0.0; + xval = loc->xcol[loc->ix - 1]; } - return; + } else { + xramp = 1.0; } if (yval < loc->ycol[0] || yval > loc->ycol[loc->iy - 1]) { - if (PARAM(verbose) > 0) { - cm_message_printf("y value %g exceeds table limits,\n" - " please enlarge range of your table", - yval); + if (PARAM(verbose) > 0) + cm_message_printf(exceed_fmt, 'y', yval); + if (yval < loc->ycol[0]) { + yramp = 1 - ((loc->ycol[0] - yval) / + (RAMP_WIDTH * (loc->ycol[1] - loc->ycol[0]))); + if (yramp < 0.0) + yramp = 0.0; + yval = loc->ycol[0]; + } else { + yramp = 1 - ((yval - loc->ycol[loc->iy - 1]) / + (RAMP_WIDTH * (loc->ycol[loc->iy - 1] - + loc->ycol[loc->iy - 2]))); + if (yramp < 0.0) + yramp = 0.0; + yval = loc->ycol[loc->iy - 1]; } - return; + } else { + yramp = 1.0; } if (zval < loc->zcol[0] || zval > loc->zcol[loc->iz - 1]) { - if (PARAM(verbose) > 0) { - cm_message_printf("z value %g exceeds table limits,\n" - " please enlarge range of your table", - zval); + if (PARAM(verbose) > 0) + cm_message_printf(exceed_fmt, 'z', zval); + if (zval < loc->zcol[0]) { + zramp = 1 - ((loc->zcol[0] - zval) / + (RAMP_WIDTH * (loc->zcol[1] - loc->zcol[0]))); + if (zramp < 0.0) + zramp = 0.0; + zval = loc->zcol[0]; + } else { + zramp = 1 - ((zval - loc->zcol[loc->iz - 1]) / + (RAMP_WIDTH * (loc->zcol[loc->iz - 1] - + loc->zcol[loc->iz - 2]))); + if (zramp < 0.0) + zramp = 0.0; + zval = loc->zcol[loc->iz - 1]; } - return; + } else { + zramp = 1.0; } @@ -329,15 +377,18 @@ void cm_table3D(ARGS) /* structure holding parms, inputs, outputs, etc. */ DER /* what to compute [FUNC, DER, BOTH] */ ); - /* xind and yind may become too large */ + /* xind may become too large when xval == xcol[loc->ix - 1] */ if (xind == loc->ix - 1) { --xind; + xoff += loc->xcol[xind + 1] - loc->xcol[xind]; } if (yind == loc->iy - 1) { --yind; + yoff += loc->ycol[yind + 1] - loc->ycol[yind]; } if (zind == loc->iz - 1) { --zind; + zoff += loc->zcol[zind + 1] - loc->zcol[zind]; } /* overwrite outval from sf_eno3_apply by trilinear interpolation */ @@ -351,11 +402,11 @@ void cm_table3D(ARGS) /* structure holding parms, inputs, outputs, etc. */ double xderiv, yderiv, zderiv, outv; outv = PARAM(offset) + PARAM(gain) * outval; OUTPUT(out) = outv; - xderiv = PARAM(gain) * derivval[0] / xdiff; + xderiv = xramp * PARAM(gain) * derivval[0] / xdiff; PARTIAL(out, inx) = xderiv; - yderiv = PARAM(gain) * derivval[1] / ydiff; + yderiv = yramp * PARAM(gain) * derivval[1] / ydiff; PARTIAL(out, iny) = yderiv; - zderiv = PARAM(gain) * derivval[2] / zdiff; + zderiv = zramp * PARAM(gain) * derivval[2] / zdiff; PARTIAL(out, inz) = zderiv; if (PARAM(verbose) > 1) { @@ -366,13 +417,13 @@ void cm_table3D(ARGS) /* structure holding parms, inputs, outputs, etc. */ } else { Mif_Complex_t ac_gain; - ac_gain.real = PARAM(gain) * derivval[0] / xdiff; + ac_gain.real = xramp * PARAM(gain) * derivval[0] / xdiff; ac_gain.imag= 0.0; AC_GAIN(out, inx) = ac_gain; - ac_gain.real = PARAM(gain) * derivval[1] / ydiff; + ac_gain.real = yramp * PARAM(gain) * derivval[1] / ydiff; ac_gain.imag= 0.0; AC_GAIN(out, iny) = ac_gain; - ac_gain.real = PARAM(gain) * derivval[2] / zdiff; + ac_gain.real = zramp * PARAM(gain) * derivval[2] / zdiff; ac_gain.imag= 0.0; AC_GAIN(out, iny) = ac_gain; }