Least-squares minimisation.
1894 {
1895
1896
1897
1898
1899
1900
1901
1902
1903
1904
1905
1906
1907
1908
1909
1910
1911 const unsigned int n = par.size();
1912 const unsigned int m = x.size();
1913
1914 if (n > m) {
1915 std::cerr << "LeastSquaresFit: Number of parameters to be varied\n"
1916 << " exceeds the number of data points.\n";
1917 return false;
1918 }
1919
1920
1921 if (*std::min_element(ey.cbegin(), ey.cend()) <= 0.) {
1922 std::cerr << "LeastSquaresFit: Not all errors are > 0; no fit done.\n";
1923 return false;
1924 }
1925 chi2 = 0.;
1926
1927 double diffc = -1.;
1928
1929 std::vector<double> r(m, 0.);
1930 for (unsigned int i = 0; i < m; ++i) {
1931
1932 r[i] = (y[i] - f(x[i], par)) / ey[i];
1933
1934 diffc = std::max(std::abs(r[i]), diffc);
1935
1936 chi2 += r[i] * r[i];
1937 }
1938 if (debug) {
1939 std::cout << " Input data points:\n"
1940 << " X Y Y - F(X)\n";
1941 for (unsigned int i = 0; i < m; ++i) {
1942 std::printf(" %9u %15.8e %15.8e %15.8e\n", i, x[i], y[i], r[i]);
1943 }
1944 std::cout << " Initial values of the fit parameters:\n"
1945 << " Parameter Value\n";
1946 for (unsigned int i = 0; i < n; ++i) {
1947 std::printf(" %9u %15.8e\n", i, par[i]);
1948 }
1949 }
1950 if (verbose) {
1951 std::cout << " MINIMISATION SUMMARY\n"
1952 << " Initial situation:\n";
1953 std::printf(" Largest difference between fit and target: %15.8e\n",
1954 diffc);
1955 std::printf(" Sum of squares of these differences: %15.8e\n",
1956 chi2);
1957 }
1958
1959 bool converged = false;
1960 double chi2L = 0.;
1961 for (unsigned int iter = 1; iter <= nMaxIter; ++iter) {
1962
1963 if ((diffc < diff) || (iter > 1 && std::abs(chi2L - chi2) < eps * chi2)) {
1964 if (debug || verbose) {
1965 if (diffc < diff) {
1966 std::cout << " The maximum difference stopping criterion "
1967 << "is satisfied.\n";
1968 } else {
1969 std::cout << " The relative change in chi-squared has dropped "
1970 << "below the threshold.\n";
1971 }
1972 }
1973 converged = true;
1974 break;
1975 }
1976
1977 std::vector<std::vector<double> > d(n, std::vector<double>(m, 0.));
1978 for (unsigned int i = 0; i < n; ++i) {
1979 const double epsdif = eps * (1. + std::abs(par[i]));
1980 par[i] += 0.5 * epsdif;
1981 for (unsigned int j = 0; j < m; ++j) d[i][j] = f(x[j], par);
1982 par[i] -= epsdif;
1983 for (unsigned int j = 0; j < m; ++j) {
1984 d[i][j] = (d[i][j] - f(x[j], par)) / (epsdif * ey[j]);
1985 }
1986 par[i] += 0.5 * epsdif;
1987 }
1988
1989 std::vector<double> colsum(n, 0.);
1990 std::vector<int> pivot(n, 0);
1991 for (unsigned int i = 0; i < n; ++i) {
1992 colsum[i] = std::inner_product(d[i].cbegin(), d[i].cend(),
1993 d[i].cbegin(), 0.);
1994 pivot[i] = i;
1995 }
1996
1997 std::vector<double>
alpha(n, 0.);
1998 bool singular = false;
1999 for (unsigned int k = 0; k < n; ++k) {
2000 double sigma = colsum[k];
2001 unsigned int jbar = k;
2002 for (unsigned int j = k + 1; j < n; ++j) {
2003 if (sigma < colsum[j]) {
2004 sigma = colsum[j];
2005 jbar = j;
2006 }
2007 }
2008 if (jbar != k) {
2009
2010 std::swap(pivot[k], pivot[jbar]);
2011 std::swap(colsum[k], colsum[jbar]);
2012 std::swap(d[k], d[jbar]);
2013 }
2014 sigma = 0.;
2015 for (unsigned int i = k; i < m; ++i) sigma += d[k][i] * d[k][i];
2016 if (sigma == 0. ||
sqrt(sigma) < 1.e-8 * std::abs(d[k][k])) {
2017 singular = true;
2018 break;
2019 }
2021 const double beta = 1. / (sigma - d[k][k] *
alpha[k]);
2022 d[k][k] -=
alpha[k];
2023 std::vector<double> b(n, 0.);
2024 for (unsigned int j = k + 1; j < n; ++j) {
2025 for (unsigned int i = k; i < n; ++i) b[j] += d[k][i] * d[j][i];
2026 b[j] *= beta;
2027 }
2028 for (unsigned int j = k + 1; j < n; ++j) {
2029 for (unsigned int i = k; i < m; ++i) {
2030 d[j][i] -= d[k][i] * b[j];
2031 colsum[j] -= d[j][k] * d[j][k];
2032 }
2033 }
2034 }
2035 if (singular) {
2036 std::cerr << "LeastSquaresFit: Householder matrix (nearly) singular;\n"
2037 << " no further optimisation.\n"
2038 << " Ensure the function depends on the parameters\n"
2039 << " and try to supply reasonable starting values.\n";
2040 break;
2041 }
2042
2043 for (unsigned int j = 0; j < n; ++j) {
2044 double gamma = 0.;
2045 for (unsigned int i = j; i < m; ++i) gamma += d[j][i] * r[i];
2046 gamma *= 1. / (
alpha[j] * d[j][j]);
2047 for (unsigned int i = j; i < m; ++i) r[i] += gamma * d[j][i];
2048 }
2049 std::vector<double>
z(n, 0.);
2050 z[n - 1] = r[n - 1] /
alpha[n - 1];
2051 for (int i = n - 1; i >= 1; --i) {
2052 double sum = 0.;
2053 for (unsigned int j = i + 1; j <= n; ++j) {
2054 sum += d[j - 1][i - 1] *
z[j - 1];
2055 }
2056 z[i - 1] = (r[i - 1] - sum) / alpha[i - 1];
2057 }
2058
2059 std::vector<double> s(n, 0.);
2060 for (unsigned int i = 0; i < n; ++i) s[pivot[i]] = z[i];
2061
2062 if (debug) {
2063 std::cout << " Correction vector in iteration " << iter << ":\n";
2064 for (unsigned int i = 0; i < n; ++i) {
2065 std::printf(" %5u %15.8e\n", i, s[i]);
2066 }
2067 }
2068
2069 chi2L = chi2;
2070 chi2 *= 2;
2071 for (unsigned int i = 0; i < n; ++i) par[i] += s[i] * 2;
2072 for (unsigned int i = 0; i <= 10; ++i) {
2073 if (chi2 <= chi2L) break;
2074 if (std::abs(chi2L - chi2) < eps * chi2) {
2075 if (debug) {
2076 std::cout << " Too little improvement; reduction loop halted.\n";
2077 }
2078 break;
2079 }
2080 chi2 = 0.;
2081 const double scale = 1. /
pow(2, i);
2082 for (unsigned int j = 0; j < n; ++j) par[j] -= s[j] * scale;
2083 for (unsigned int j = 0; j < m; ++j) {
2084 r[j] = (
y[j] - f(x[j], par)) / ey[j];
2085 chi2 += r[j] * r[j];
2086 }
2087 if (debug) {
2088 std::printf(" Reduction loop %3i: chi2 = %15.8e\n", i, chi2);
2089 }
2090 }
2091
2092 diffc = std::abs(r[0]);
2093 for (unsigned int i = 1; i < m; ++i) {
2094 diffc = std::max(std::abs(r[i]), diffc);
2095 }
2096
2097 if (debug) {
2098 std::cout << " Values of the fit parameters after iteration "
2099 << iter << "\n Parameter Value\n";
2100 for (unsigned int i = 0; i < n; ++i) {
2101 std::printf(" %9u %15.8e\n", i, par[i]);
2102 }
2103 std::printf(" for which chi2 = %15.8e and diff = %15.8e\n",
2104 chi2, diffc);
2105 } else if (verbose) {
2106 std::printf(" Step %3u: largest deviation = %15.8e, chi2 = %15.8e\n",
2107 iter, diffc, chi2);
2108 }
2109 }
2110
2111 if (!converged) {
2112 std::cerr << "LeastSquaresFit: Maximum number of iterations reached.\n";
2113 }
2114
2115 std::vector<std::vector<double> > d(n, std::vector<double>(m, 0.));
2116 for (unsigned int i = 0; i < n; ++i) {
2117 const double epsdif = eps * (1. + std::abs(par[i]));
2118 par[i] += 0.5 * epsdif;
2119 for (unsigned int j = 0; j < m; ++j) d[i][j] = f(x[j], par);
2120 par[i] -= epsdif;
2121 for (unsigned int j = 0; j < m; ++j) {
2122 d[i][j] = (d[i][j] - f(x[j], par)) / (epsdif * ey[j]);
2123 }
2124 par[i] += 0.5 * epsdif;
2125 }
2126
2127 std::vector<std::vector<double> > cov(n, std::vector<double>(n, 0.));
2128 for (unsigned int i = 0; i < n; ++i) {
2129 for (unsigned int j = 0; j < n; ++j) {
2130 cov[i][j] = std::inner_product(d[i].cbegin(), d[i].cend(),
2131 d[j].cbegin(), 0.);
2132 }
2133 }
2134
2135 double scale = m > n ? chi2 / (m - n) : 1.;
2136
2137 epar.assign(n, 0.);
2139 std::cerr << "LeastSquaresFit: Singular covariance matrix; "
2140 << "no error calculation.\n";
2141 } else {
2142 for (unsigned int i = 0; i < n; ++i) {
2143 for (unsigned int j = 0; j < n; ++j) cov[i][j] *= scale;
2144 epar[i] =
sqrt(std::max(0., cov[i][i]));
2145 }
2146 }
2147
2148 if (debug) {
2149 std::cout << " Comparison between input and fit:\n"
2150 << " X Y F(X)\n";
2151 for (unsigned int i = 0; i < m; ++i) {
2152 const double fit = f(x[i], par);
2153 std::printf(" %5u %15.8e %15.8e %15.8e\n", i, x[i], y[i], fit);
2154 }
2155 }
2156 if (verbose) {
2157 std::cout << " Final values of the fit parameters:\n"
2158 << " Parameter Value Error\n";
2159 for (unsigned int i = 0; i < n; ++i) {
2160 std::printf(" %9u %15.8e %15.8e\n", i, par[i], epar[i]);
2161 }
2162 std::cout << " The errors have been scaled by a factor of "
2163 <<
sqrt(scale) <<
".\n";
2164 std::cout << " Covariance matrix:\n";
2165 for (unsigned int i = 0; i < n; ++i) {
2166 for (unsigned int j = 0; j < n; ++j) {
2167 std::printf(" %15.8e", cov[i][j]);
2168 }
2169 std::cout << "\n";
2170 }
2171 std::cout << " Correlation matrix:\n";
2172 for (unsigned int i = 0; i < n; ++i) {
2173 for (unsigned int j = 0; j < n; ++j) {
2174 double cor = 0.;
2175 if (cov[i][i] > 0. && cov[j][j] > 0.) {
2176 cor = cov[i][j] /
sqrt(cov[i][i] * cov[j][j]);
2177 }
2178 std::printf(" %15.8e", cor);
2179 }
2180 std::cout << "\n";
2181 }
2182 std::cout << " Minimisation finished.\n";
2183 }
2184 return converged;
2185}
int dinv(const int n, std::vector< std::vector< double > > &a)
Replace square matrix A by its inverse.
DoubleAc pow(const DoubleAc &f, double p)
DoubleAc sqrt(const DoubleAc &f)
1.13.2