polynomial_fit.py 4.36 KB
Newer Older
Martin Reinecke's avatar
Martin Reinecke committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2019 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.

18
19
20
import matplotlib.pyplot as plt
import numpy as np

Martin Reinecke's avatar
5->6  
Martin Reinecke committed
21
import nifty6 as ift
22
23
24
25
26
27
28
29
30
31
np.random.seed(12)


def polynomial(coefficients, sampling_points):
    """Computes values of polynomial whose coefficients are stored in
    coefficients at sampling points. This is a quick version of the
    PolynomialResponse.

    Parameters
    ----------
32
    coefficients: Field
33
34
35
    sampling_points: Numpy array
    """

36
    if not (isinstance(coefficients, ift.Field)
37
38
            and isinstance(sampling_points, np.ndarray)):
        raise TypeError
Martin Reinecke's avatar
Martin Reinecke committed
39
    params = coefficients.val
40
41
42
43
44
45
46
    out = np.zeros_like(sampling_points)
    for ii in range(len(params)):
        out += params[ii] * sampling_points**ii
    return out


class PolynomialResponse(ift.LinearOperator):
Martin Reinecke's avatar
Martin Reinecke committed
47
48
    """Calculates values of a polynomial parameterized by input at sampling
    points.
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64

    Parameters
    ----------
    domain: UnstructuredDomain
        The domain on which the coefficients of the polynomial are defined.
    sampling_points: Numpy array
        x-values of the sampling points.
    """

    def __init__(self, domain, sampling_points):
        if not (isinstance(domain, ift.UnstructuredDomain)
                and isinstance(x, np.ndarray)):
            raise TypeError
        self._domain = ift.DomainTuple.make(domain)
        tgt = ift.UnstructuredDomain(sampling_points.shape)
        self._target = ift.DomainTuple.make(tgt)
Martin Reinecke's avatar
Martin Reinecke committed
65
        self._capability = self.TIMES | self.ADJOINT_TIMES
66
67
68
69
70
71
72
73

        sh = (self.target.size, domain.size)
        self._mat = np.empty(sh)
        for d in range(domain.size):
            self._mat.T[d] = sampling_points**d

    def apply(self, x, mode):
        self._check_input(x, mode)
Martin Reinecke's avatar
Martin Reinecke committed
74
        val = x.val_rw()
75
76
77
78
79
80
        if mode == self.TIMES:
            # FIXME Use polynomial() here
            out = self._mat.dot(val)
        else:
            # FIXME Can this be optimized?
            out = self._mat.conj().T.dot(val)
Martin Reinecke's avatar
Martin Reinecke committed
81
        return ift.makeField(self._tgt(mode), out)
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96


# Generate some mock data
N_params = 10
N_samples = 100
size = (12,)
x = np.random.random(size) * 10
y = np.sin(x**2) * x**3
var = np.full_like(y, y.var() / 10)
var[-2] *= 4
var[5] /= 2
y[5] -= 0

# Set up minimization problem
p_space = ift.UnstructuredDomain(N_params)
97
params = ift.full(p_space, 0.)
98
99
100
101
R = PolynomialResponse(p_space, x)
ift.extra.consistency_check(R)

d_space = R.target
Martin Reinecke's avatar
Martin Reinecke committed
102
103
d = ift.makeField(d_space, y)
N = ift.DiagonalOperator(ift.makeField(d_space, var))
104

Martin Reinecke's avatar
Martin Reinecke committed
105
IC = ift.DeltaEnergyController(tol_rel_deltaE=1e-12, iteration_limit=200)
Martin Reinecke's avatar
Martin Reinecke committed
106
likelihood = ift.GaussianEnergy(d, N)(R)
107
Ham = ift.StandardHamiltonian(likelihood, IC)
Martin Reinecke's avatar
fixes  
Martin Reinecke committed
108
H = ift.EnergyAdapter(params, Ham, want_metric=True)
109
110

# Minimize
111
minimizer = ift.NewtonCG(IC)
112
113
114
H, _ = minimizer(H)

# Draw posterior samples
Martin Reinecke's avatar
fixes  
Martin Reinecke committed
115
metric = Ham(ift.Linearization.make_var(H.position, want_metric=True)).metric
Martin Reinecke's avatar
fix  
Martin Reinecke committed
116
samples = [metric.draw_sample(from_inverse=True) + H.position
117
118
119
120
121
122
123
124
125
           for _ in range(N_samples)]

# Plotting
plt.errorbar(x, y, np.sqrt(var), fmt='ko', label='Data with error bars')
xmin, xmax = x.min(), x.max()
xs = np.linspace(xmin, xmax, 100)

sc = ift.StatCalculator()
for ii in range(len(samples)):
126
127
    sc.add(samples[ii])
    ys = polynomial(samples[ii], xs)
128
129
130
131
    if ii == 0:
        plt.plot(xs, ys, 'k', alpha=.05, label='Posterior samples')
        continue
    plt.plot(xs, ys, 'k', alpha=.05)
132
ys = polynomial(H.position, xs)
133
134
135
136
137
138
plt.plot(xs, ys, 'r', linewidth=2., label='Interpolation')
plt.legend()
plt.savefig('fit.png')
plt.close()

# Print parameters
Martin Reinecke's avatar
Martin Reinecke committed
139
140
141
142
mean = sc.mean.val
sigma = np.sqrt(sc.var.val)
for ii in range(len(mean)):
    print('Coefficient x**{}: {:.2E} +/- {:.2E}'.format(ii, mean[ii],
143
                                                            sigma[ii]))