| 117 | } // End empty namespace |
| 118 | |
| 119 | bool Interface2C(void) |
| 120 | { // This routine is intentionally coded as if it were a C routine |
| 121 | // except for the fact that it uses the predefined type bool. |
| 122 | bool ok = true; |
| 123 | |
| 124 | // declare variables |
| 125 | float x, a[6], y, dyda[6], tmp[6]; |
| 126 | size_t na, i; |
| 127 | |
| 128 | // number of parameters (3 for each Gaussian) |
| 129 | na = 6; |
| 130 | |
| 131 | // number of Gaussians: n = na / 3; |
| 132 | |
| 133 | // value of x |
| 134 | x = 1.; |
| 135 | |
| 136 | // value of the parameter vector a |
| 137 | for(i = 0; i < na; i++) |
| 138 | a[i] = (float) (i+1); |
| 139 | |
| 140 | // evaluate function and derivative |
| 141 | sumGauss(x, a, &y, dyda, na); |
| 142 | |
| 143 | // compare dyda to central difference approximation for deriative |
| 144 | for(i = 0; i < na; i++) |
| 145 | { // local variables |
| 146 | float eps, ai, yp, ym, dy_da; |
| 147 | |
| 148 | // We assume that the type float has at least 7 digits of |
| 149 | // precision, so we choose eps to be about pow(10., -7./2.). |
| 150 | eps = (float) 3e-4; |
| 151 | |
| 152 | // value of this component of a |
| 153 | ai = a[i]; |
| 154 | |
| 155 | // evaluate F( a + eps * ei ) |
| 156 | a[i] = ai + eps; |
| 157 | sumGauss(x, a, &yp, tmp, na); |
| 158 | |
| 159 | // evaluate F( a - eps * ei ) |
| 160 | a[i] = ai - eps; |
| 161 | sumGauss(x, a, &ym, tmp, na); |
| 162 | |
| 163 | // evaluate central difference approximates for partial |
| 164 | dy_da = (yp - ym) / (2 * eps); |
| 165 | |
| 166 | // restore this component of a |
| 167 | a[i] = ai; |
| 168 | |
| 169 | ok &= NearEqual(dyda[i], dy_da, eps, eps); |
| 170 | } |
| 171 | return ok; |
| 172 | } |
| 173 | // END C++ |