correlated_fields.py 15.4 KB
Newer Older
Martin Reinecke's avatar
Martin Reinecke committed
1 2 3 4 5 6 7 8 9 10 11 12 13
# 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/>.
#
14
# Copyright(C) 2013-2019 Max-Planck-Society
Philipp Arras's avatar
Philipp Arras committed
15
# Authors: Philipp Frank, Philipp Arras
Martin Reinecke's avatar
Martin Reinecke committed
16
#
17
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
Philipp Arras's avatar
Philipp Arras committed
18

Philipp Arras's avatar
Philipp Arras committed
19
import numpy as np
20

Philipp Arras's avatar
Philipp Arras committed
21
from ..domain_tuple import DomainTuple
Philipp Arras's avatar
Philipp Arras committed
22 23
from ..domains.power_space import PowerSpace
from ..domains.unstructured_domain import UnstructuredDomain
24
from ..field import Field
Philipp Arras's avatar
Philipp Arras committed
25
from ..operators.adder import Adder
26
from ..operators.contraction_operator import ContractionOperator
Philipp Arras's avatar
Philipp Arras committed
27
from ..operators.distributors import PowerDistributor
Philipp Arras's avatar
Philipp Arras committed
28
from ..operators.endomorphic_operator import EndomorphicOperator
Martin Reinecke's avatar
Martin Reinecke committed
29
from ..operators.harmonic_operators import HarmonicTransformOperator
Philipp Arras's avatar
Philipp Arras committed
30
from ..operators.linear_operator import LinearOperator
Philipp Arras's avatar
Philipp Arras committed
31 32
from ..operators.operator import Operator
from ..operators.simple_linear_operators import VdotOperator, ducktape
Philipp Arras's avatar
Philipp Arras committed
33
from ..operators.value_inserter import ValueInserter
34
from ..sugar import from_global_data, full, makeDomain
35
from ..probing import StatCalculator
36

Philipp Arras's avatar
Philipp Arras committed
37

38 39
def _lognormal_moments(mean, sig):
    mean, sig = float(mean), float(sig)
Philipp Arras's avatar
Philipp Arras committed
40 41 42
    assert sig > 0
    logsig = np.sqrt(np.log((sig/mean)**2 + 1))
    logmean = np.log(mean) - logsig**2/2
43 44
    return logmean, logsig

Philipp Arras's avatar
Philipp Arras committed
45

46
class _lognormal_moment_matching(Operator):
Philipp Arras's avatar
Philipp Arras committed
47
    def __init__(self, mean, sig, key):
48 49 50 51 52 53 54 55
        key = str(key)
        logmean, logsig = _lognormal_moments(mean, sig)
        self._mean = mean
        self._sig = sig
        op = _normal(logmean, logsig, key).exp()
        self._domain = op.domain
        self._target = op.target
        self.apply = op.apply
56

57 58 59
    @property
    def mean(self):
        return self._mean
Philipp Arras's avatar
Philipp Arras committed
60

61 62 63
    @property
    def std(self):
        return self._sig
64

Philipp Arras's avatar
Philipp Arras committed
65

Philipp Arras's avatar
Philipp Arras committed
66 67 68 69 70
def _normal(mean, sig, key):
    return Adder(Field.scalar(mean)) @ (
        sig*ducktape(DomainTuple.scalar_domain(), None, key))


Philipp Arras's avatar
Philipp Arras committed
71 72 73 74
def _log_k_lengths(pspace):
    return np.log(pspace.k_lengths[1:])


Philipp Arras's avatar
Philipp Arras committed
75 76 77 78 79 80 81 82 83 84 85
def _logkl(power_space):
    power_space = DomainTuple.make(power_space)
    assert isinstance(power_space[0], PowerSpace)
    assert len(power_space) == 1
    logkl = _log_k_lengths(power_space[0])
    assert logkl.shape[0] == power_space[0].shape[0] - 1
    logkl -= logkl[0]
    logkl = np.insert(logkl, 0, 0)
    return logkl


Philipp Arras's avatar
Philipp Arras committed
86 87 88 89 90 91 92 93
def _log_vol(power_space):
    power_space = DomainTuple.make(power_space)
    assert isinstance(power_space[0], PowerSpace)
    assert len(power_space) == 1
    logk_lengths = _log_k_lengths(power_space[0])
    return logk_lengths[1:] - logk_lengths[:-1]


Philipp Frank's avatar
Philipp Frank committed
94
class _SlopeRemover(EndomorphicOperator):
Philipp Arras's avatar
Philipp Arras committed
95
    def __init__(self, domain):
Philipp Frank's avatar
Philipp Frank committed
96
        self._domain = makeDomain(domain)
Philipp Arras's avatar
Philipp Arras committed
97 98
        assert len(self._domain) == 1
        assert isinstance(self._domain[0], PowerSpace)
Philipp Arras's avatar
Philipp Arras committed
99
        logkl = _logkl(self._domain)
100
        self._sc = logkl/float(logkl[-1])
Philipp Frank's avatar
Philipp Frank committed
101
        self._capability = self.TIMES | self.ADJOINT_TIMES
Philipp Arras's avatar
Philipp Arras committed
102

103 104
    def apply(self, x, mode):
        self._check_input(x, mode)
Philipp Frank's avatar
Philipp Frank committed
105 106
        x = x.to_global_data()
        if mode == self.TIMES:
107
            res = x - x[-1]*self._sc
Philipp Frank's avatar
Philipp Frank committed
108
        else:
109
            res = np.zeros(x.shape, dtype=x.dtype)
Philipp Frank's avatar
Philipp Frank committed
110 111
            res += x
            res[-1] -= (x*self._sc).sum()
112
        return from_global_data(self._tgt(mode), res)
Philipp Frank's avatar
Philipp Frank committed
113

114

Philipp Arras's avatar
Philipp Arras committed
115 116 117
class _TwoLogIntegrations(LinearOperator):
    def __init__(self, target):
        self._target = makeDomain(target)
Philipp Arras's avatar
Philipp Arras committed
118 119
        assert len(self._target) == 1
        assert isinstance(self._target[0], PowerSpace)
Philipp Arras's avatar
Philipp Arras committed
120 121 122 123 124
        self._domain = makeDomain(
            UnstructuredDomain((2, self.target.shape[0] - 2)))
        self._capability = self.TIMES | self.ADJOINT_TIMES
        if not isinstance(self._target[0], PowerSpace):
            raise TypeError
Philipp Arras's avatar
Philipp Arras committed
125
        self._log_vol = _log_vol(self._target[0])
Philipp Arras's avatar
Philipp Arras committed
126 127 128 129 130 131 132 133 134

    def apply(self, x, mode):
        self._check_input(x, mode)
        if mode == self.TIMES:
            x = x.to_global_data()
            res = np.empty(self._target.shape)
            res[0] = 0
            res[1] = 0
            res[2:] = np.cumsum(x[1])
Philipp Arras's avatar
Philipp Arras committed
135
            res[2:] = (res[2:] + res[1:-1])/2*self._log_vol + x[0]
Philipp Arras's avatar
Philipp Arras committed
136 137 138 139 140 141 142
            res[2:] = np.cumsum(res[2:])
            return from_global_data(self._target, res)
        else:
            x = x.to_global_data_rw()
            res = np.zeros(self._domain.shape)
            x[2:] = np.cumsum(x[2:][::-1])[::-1]
            res[0] += x[2:]
Philipp Arras's avatar
Philipp Arras committed
143
            x[2:] *= self._log_vol/2.
144 145
            x[1:-1] += x[2:]
            res[1] += np.cumsum(x[2:][::-1])[::-1]
Philipp Arras's avatar
Philipp Arras committed
146 147 148 149 150 151
            return from_global_data(self._domain, res)


class _Normalization(Operator):
    def __init__(self, domain):
        self._domain = self._target = makeDomain(domain)
Philipp Arras's avatar
Philipp Arras committed
152
        assert len(self._domain) == 1
Philipp Arras's avatar
Philipp Arras committed
153
        assert isinstance(self._domain[0], PowerSpace)
Philipp Arras's avatar
Philipp Arras committed
154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172
        hspace = self._domain[0].harmonic_partner
        pd = PowerDistributor(hspace, power_space=self._domain[0])
        cst = pd.adjoint(full(pd.target, 1.)).to_global_data_rw()
        cst[0] = 0
        self._cst = from_global_data(self._domain, cst)
        self._specsum = _SpecialSum(self._domain)

    def apply(self, x):
        self._check_input(x)
        amp = x.exp()
        spec = (2*x).exp()
        # FIXME This normalizes also the zeromode which is supposed to be left
        # untouched by this operator
        return self._specsum(self._cst*spec)**(-0.5)*amp


class _SpecialSum(EndomorphicOperator):
    def __init__(self, domain):
        self._domain = makeDomain(domain)
Philipp Arras's avatar
Philipp Arras committed
173
        assert len(self._domain) == 1
Philipp Arras's avatar
Philipp Arras committed
174 175 176 177 178 179 180
        self._capability = self.TIMES | self.ADJOINT_TIMES

    def apply(self, x, mode):
        self._check_input(x, mode)
        return full(self._tgt(mode), x.sum())


181 182 183
class _Amplitude(Operator):
    def __init__(self, target, fluctuations, flexibility, asperity,
                 loglogavgslope, key):
Philipp Arras's avatar
Philipp Arras committed
184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199
        """
        fluctuations > 0
        flexibility > 0
        asperity > 0
        loglogavgslope probably negative
        """
        assert isinstance(fluctuations, Operator)
        assert isinstance(flexibility, Operator)
        assert isinstance(asperity, Operator)
        assert isinstance(loglogavgslope, Operator)
        target = makeDomain(target)
        assert len(target) == 1
        assert isinstance(target[0], PowerSpace)

        twolog = _TwoLogIntegrations(target)
        sc = np.zeros(twolog.domain.shape)
Philipp Arras's avatar
Philipp Arras committed
200
        sc[0] = sc[1] = np.sqrt(_log_vol(target))
Philipp Arras's avatar
Philipp Arras committed
201 202
        sc = from_global_data(twolog.domain, sc)
        expander = VdotOperator(sc).adjoint
Philipp Arras's avatar
Philipp Arras committed
203
        sig_flex = expander @ flexibility
Philipp Arras's avatar
Philipp Arras committed
204

Philipp Arras's avatar
Philipp Arras committed
205 206
        dist = np.zeros(twolog.domain.shape, dtype=np.float64)
        dist[0] += 1
Philipp Arras's avatar
Philipp Arras committed
207 208 209 210 211 212 213 214
        sig_asp = VdotOperator(from_global_data(twolog.domain,
                                                dist)).adjoint @ asperity

        shift = np.ones(twolog.domain.shape)
        shift[0] = _log_vol(target)**2/12.
        shift = from_global_data(twolog.domain, shift)
        sigma = sig_flex*(Adder(shift) @ sig_asp).sqrt()
        smooth = twolog @ (sigma*ducktape(twolog.domain, None, key))
Philipp Arras's avatar
Philipp Arras committed
215

Philipp Arras's avatar
Philipp Arras committed
216
        tg = smooth.target
Philipp Arras's avatar
Philipp Arras committed
217 218
        noslope = _SlopeRemover(tg) @ smooth
        _t = VdotOperator(from_global_data(tg, _logkl(tg))).adjoint
Philipp Arras's avatar
Philipp Arras committed
219
        smoothslope = _t @ loglogavgslope + noslope
Philipp Arras's avatar
Philipp Arras committed
220 221 222 223 224 225 226 227 228

        normal_ampl = _Normalization(target) @ smoothslope
        vol = target[0].harmonic_partner.get_default_codomain().total_volume
        arr = np.zeros(target.shape)
        arr[1:] = vol
        expander = VdotOperator(from_global_data(target, arr)).adjoint
        mask = np.zeros(target.shape)
        mask[0] = vol
        adder = Adder(from_global_data(target, mask))
Philipp Arras's avatar
Philipp Arras committed
229 230
        op = adder @ ((expander @ fluctuations)*normal_ampl)
        self.apply = op.apply
231
        self.fluctuation_amplitude = fluctuations
Philipp Arras's avatar
Philipp Arras committed
232
        self._domain, self._target = op.domain, op.target
Philipp Arras's avatar
Philipp Arras committed
233

234 235 236 237

class CorrelatedFieldMaker:
    def __init__(self):
        self._a = []
238
        self._azm = None
239 240 241 242 243 244 245 246 247 248 249

    def add_fluctuations(self,
                         target,
                         fluctuations_mean,
                         fluctuations_stddev,
                         flexibility_mean,
                         flexibility_stddev,
                         asperity_mean,
                         asperity_stddev,
                         loglogavgslope_mean,
                         loglogavgslope_stddev,
250
                         prefix='',
Philipp Arras's avatar
Philipp Arras committed
251
                         index=None):
Philipp Arras's avatar
Philipp Arras committed
252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275
        fluctuations_mean = float(fluctuations_mean)
        fluctuations_stddev = float(fluctuations_stddev)
        flexibility_mean = float(flexibility_mean)
        flexibility_stddev = float(flexibility_stddev)
        asperity_mean = float(asperity_mean)
        asperity_stddev = float(asperity_stddev)
        loglogavgslope_mean = float(loglogavgslope_mean)
        loglogavgslope_stddev = float(loglogavgslope_stddev)
        prefix = str(prefix)
        assert fluctuations_stddev > 0
        assert fluctuations_mean > 0
        assert flexibility_stddev > 0
        assert flexibility_mean > 0
        assert asperity_stddev > 0
        assert asperity_mean > 0
        assert loglogavgslope_stddev > 0

        fluct = _lognormal_moment_matching(fluctuations_mean,
                                           fluctuations_stddev,
                                           prefix + 'fluctuations')
        flex = _lognormal_moment_matching(flexibility_mean, flexibility_stddev,
                                          prefix + 'flexibility')
        asp = _lognormal_moment_matching(asperity_mean, asperity_stddev,
                                         prefix + 'asperity')
276
        avgsl = _normal(loglogavgslope_mean, loglogavgslope_stddev,
Philipp Arras's avatar
Philipp Arras committed
277
                        prefix + 'loglogavgslope')
278 279 280 281 282
        amp = _Amplitude(target, fluct, flex, asp, avgsl, prefix + 'spectrum')
        if index is not None:
            self._a.insert(index, amp)
        else:
            self._a.append(amp)
283 284 285

    def finalize_from_op(self, zeromode, prefix=''):
        assert isinstance(zeromode, Operator)
286
        self._azm = zeromode
287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308
        hspace = makeDomain([dd.target[0].harmonic_partner for dd in self._a])
        foo = np.ones(hspace.shape)
        zeroind = len(hspace.shape)*(0,)
        foo[zeroind] = 0
        azm = Adder(from_global_data(hspace, foo)) @ ValueInserter(
            hspace, zeroind) @ zeromode

        n_amplitudes = len(self._a)
        ht = HarmonicTransformOperator(hspace, space=0)
        for i in range(1, n_amplitudes):
            ht = HarmonicTransformOperator(ht.target, space=i) @ ht

        pd = PowerDistributor(hspace, self._a[0].target[0], 0)
        for i in range(1, n_amplitudes):
            foo = PowerDistributor(pd.domain, self._a[i].target[0], space=i)
            pd = pd @ foo

        spaces = tuple(range(n_amplitudes))
        a = ContractionOperator(pd.domain, spaces[1:]).adjoint @ self._a[0]
        for i in range(1, n_amplitudes):
            co = ContractionOperator(pd.domain, spaces[:i] + spaces[(i + 1):])
            a = a*(co.adjoint @ self._a[i])
Philipp Arras's avatar
Philipp Arras committed
309

310
        return ht(azm*(pd @ a)*ducktape(hspace, None, prefix + 'xi'))
Philipp Arras's avatar
Philipp Arras committed
311 312 313 314

    def finalize(self,
                 offset_amplitude_mean,
                 offset_amplitude_stddev,
315
                 prefix='',
Philipp Arras's avatar
Philipp Arras committed
316 317 318 319 320 321 322 323 324
                 offset=None):
        """
        offset vs zeromode: volume factor
        """
        offset_amplitude_stddev = float(offset_amplitude_stddev)
        offset_amplitude_mean = float(offset_amplitude_mean)
        assert offset_amplitude_stddev > 0
        assert offset_amplitude_mean > 0
        if offset is not None:
325
            raise NotImplementedError
Philipp Arras's avatar
Philipp Arras committed
326 327 328 329
            offset = float(offset)
        azm = _lognormal_moment_matching(offset_amplitude_mean,
                                         offset_amplitude_stddev,
                                         prefix + 'zeromode')
330
        return self.finalize_from_op(azm, prefix)
Philipp Arras's avatar
Philipp Arras committed
331 332 333

    @property
    def amplitudes(self):
334
        return self._a
335

336 337 338 339 340 341 342
    @property
    def amplitude_total_offset(self):
        return self._azm

    @property
    def total_fluctuation(self):
        if len(self._a) == 0:
Philipp Arras's avatar
Philipp Arras committed
343
            raise NotImplementedError
344 345 346 347 348
        if len(self._a) == 1:
            return self._a[0].fluctuation_amplitude
        q = 1.
        for a in self._a:
            fl = a.fluctuation_amplitude
Philipp Arras's avatar
Philipp Arras committed
349 350
            q = q*(Adder(full(fl.target, 1.)) @ fl**2)
        return (Adder(full(q.target, -1.)) @ q).sqrt()
351

Philipp Arras's avatar
Philipp Arras committed
352
    def slice_fluctuation(self, space):
353
        if len(self._a) == 0:
Philipp Arras's avatar
Philipp Arras committed
354
            raise NotImplementedError
355 356 357 358 359 360 361
        assert space < len(self._a)
        if len(self._a) == 1:
            return self._a[0].fluctuation_amplitude
        q = 1.
        for j in range(len(self._a)):
            fl = self._a[j].fluctuation_amplitude
            if j == space:
Philipp Arras's avatar
Philipp Arras committed
362
                q = q*fl**2
363
            else:
Philipp Arras's avatar
Philipp Arras committed
364
                q = q*(Adder(full(fl.target, 1.)) @ fl**2)
365
        return q.sqrt()
Philipp Arras's avatar
Philipp Arras committed
366 367

    def average_fluctuation(self, space):
368
        if len(self._a) == 0:
Philipp Arras's avatar
Philipp Arras committed
369
            raise NotImplementedError
370 371 372 373 374
        assert space < len(self._a)
        if len(self._a) == 1:
            return self._a[0].fluctuation_amplitude
        return self._a[space].fluctuation_amplitude

Philipp Arras's avatar
Philipp Arras committed
375
    def offset_amplitude_realized(self, samples):
376 377 378 379
        res = 0.
        for s in samples:
            res += s.mean()**2
        return np.sqrt(res/len(samples))
Philipp Arras's avatar
Philipp Arras committed
380 381

    def total_fluctuation_realized(self, samples):
382 383
        res = 0.
        for s in samples:
Philipp Arras's avatar
Philipp Arras committed
384
            res = res + (s - s.mean())**2
385 386
        res = res/len(samples)
        return np.sqrt(res.mean())
Philipp Arras's avatar
Philipp Arras committed
387 388

    def average_fluctuation_realized(self, samples, space):
389 390 391 392
        ldom = len(samples[0].domain)
        assert space < ldom
        if ldom == 1:
            return self.total_fluctuation_realized(samples)
Philipp Arras's avatar
Philipp Arras committed
393
        spaces = ()
394 395 396 397 398 399
        for i in range(ldom):
            if i != space:
                spaces += (i,)
        res = 0.
        for s in samples:
            r = s.mean(spaces)
Philipp Arras's avatar
Philipp Arras committed
400
            res = res + (r - r.mean())**2
401 402
        res = res/len(samples)
        return np.sqrt(res.mean())
Philipp Arras's avatar
Philipp Arras committed
403 404

    def slice_fluctuation_realized(self, samples, space):
405 406 407 408 409 410 411 412 413 414 415
        ldom = len(samples[0].domain)
        assert space < ldom
        if ldom == 1:
            return self.total_fluctuation_realized(samples)
        res1 = 0.
        res2 = 0.
        for s in samples:
            res1 = res1 + s**2
            res2 = res2 + s.mean(space)**2
        res1 = res1/len(samples)
        res2 = res2/len(samples)
Philipp Arras's avatar
Philipp Arras committed
416
        res = res1.mean() - res2.mean()
417 418
        return np.sqrt(res)

Philipp Arras's avatar
Philipp Arras committed
419
    def stats(self, op, samples):
420 421 422
        sc = StatCalculator()
        for s in samples:
            sc.add(op(s.extract(op.domain)))
423
        return sc.mean.to_global_data(), sc.var.sqrt().to_global_data()
Philipp Arras's avatar
Philipp Arras committed
424 425

    def moment_slice_to_average(self, fluctuations_slice_mean, nsamples=1000):
426 427 428 429 430
        fluctuations_slice_mean = float(fluctuations_slice_mean)
        assert fluctuations_slice_mean > 0
        scm = 1.
        for a in self._a:
            m, std = a.fluctuation_amplitude.mean, a.fluctuation_amplitude.std
Philipp Arras's avatar
Philipp Arras committed
431 432
            mu, sig = _lognormal_moments(m, std)
            flm = np.exp(mu + sig*np.random.normal(size=nsamples))
433 434
            scm *= flm**2 + 1.
        scm = np.mean(np.sqrt(scm))
Philipp Arras's avatar
Philipp Arras committed
435
        return fluctuations_slice_mean/scm