spectrum.cpp 16.2 KB
Newer Older
1
2
#include <fstream>
#include <boost/format.hpp>
Carl Poelking's avatar
Carl Poelking committed
3
4
#include <boost/archive/binary_oarchive.hpp>
#include <boost/archive/binary_iarchive.hpp>
5

6
#include "soap/spectrum.hpp"
Carl Poelking's avatar
Setup.  
Carl Poelking committed
7
8
9

namespace soap {

Carl Poelking's avatar
Carl Poelking committed
10
Spectrum::Spectrum(Structure &structure, Options &options) :
11
    _log(NULL), _options(&options), _structure(&structure), _own_basis(true), _global_atomic(NULL) {
Carl Poelking's avatar
Setup.  
Carl Poelking committed
12
	GLOG() << "Configuring spectrum ..." << std::endl;
13
14
	// CREATE & CONFIGURE BASIS
	_basis = new Basis(&options);
Carl Poelking's avatar
Setup.  
Carl Poelking committed
15
16
}

Carl Poelking's avatar
Carl Poelking committed
17
Spectrum::Spectrum(Structure &structure, Options &options, Basis &basis) :
18
	_log(NULL), _options(&options), _structure(&structure), _basis(&basis), _own_basis(false), _global_atomic(NULL) {
Carl Poelking's avatar
Carl Poelking committed
19
20
21
	;
}

Carl Poelking's avatar
Carl Poelking committed
22
Spectrum::Spectrum(std::string archfile) :
23
	_log(NULL), _options(NULL), _structure(NULL), _basis(NULL), _own_basis(true), _global_atomic(NULL) {
Carl Poelking's avatar
Carl Poelking committed
24
25
26
	this->load(archfile);
}

27
28
29
30
31
Spectrum::Spectrum() :
	_log(NULL), _options(NULL), _structure(NULL), _basis(NULL), _own_basis(true), _global_atomic(NULL) { 
    ;
}

Carl Poelking's avatar
Setup.  
Carl Poelking committed
32
33
34
Spectrum::~Spectrum() {
	delete _log;
	_log = NULL;
Carl Poelking's avatar
Carl Poelking committed
35
36
37
38
	if (_own_basis) {
		delete _basis;
		_basis = NULL;
	}
Carl Poelking's avatar
Carl Poelking committed
39
40
41
42
43
44
45
	atomspec_array_t::iterator it;
	for (it = _atomspec_array.begin(); it != _atomspec_array.end(); ++it) {
		delete *it;
	}
	_atomspec_array.clear();
	// AtomicSpectra in type map already deleted above, clear only:
	_map_atomspec_array.clear();
46
47
48

	if (_global_atomic) delete _global_atomic;
	_global_atomic = NULL;
Carl Poelking's avatar
Setup.  
Carl Poelking committed
49
50
51
}

void Spectrum::compute() {
Carl Poelking's avatar
Carl Poelking committed
52
53
54
    this->compute(_structure->particles(), _structure->particles());
}

Carl Poelking's avatar
Carl Poelking committed
55
56
57
58
void Spectrum::compute(Segment *center) {
    this->compute(center->particles(), _structure->particles());
}

Carl Poelking's avatar
Carl Poelking committed
59
60
61
62
63
64
65
void Spectrum::compute(Segment *center, Segment *target) {
    this->compute(center->particles(), target->particles());
}

void Spectrum::compute(Structure::particle_array_t &centers, Structure::particle_array_t &targets) {
    GLOG() << "Compute spectrum "
        << "(centers " << centers.size() << ", targets " << targets.size() << ") ..." << std::endl;
Carl Poelking's avatar
Setup.  
Carl Poelking committed
66
    GLOG() << _options->summarizeOptions() << std::endl;
67
68
    GLOG() << "Using radial basis of type '" << _basis->getRadBasis()->identify() << "'" << std::endl;
    GLOG() << "Using angular basis of type '" << _basis->getAngBasis()->identify() << "'" << std::endl;
Carl Poelking's avatar
Carl Poelking committed
69
    GLOG() << "Using cutoff function of type '" << _basis->getCutoff()->identify() << "'" << std::endl;
70
71

    Structure::particle_it_t pit;
Carl Poelking's avatar
Carl Poelking committed
72
    for (pit = centers.begin(); pit != centers.end(); ++pit) {
73
        // Continue if exclusion defined ...
74
75
        if (_options->doExcludeCenter((*pit)->getType()) ||
            _options->doExcludeCenterId((*pit)->getId())) continue;
76
        // Compute ...
Carl Poelking's avatar
Carl Poelking committed
77
78
79
80
81
82
83
84
85
86
87
88
89
        AtomicSpectrum *atomic_spectrum = this->computeAtomic(*pit, targets);
        this->addAtomic(atomic_spectrum);
    }
}

void Spectrum::computePower() {
    atomspec_array_t::iterator it;
    for (it = _atomspec_array.begin(); it != _atomspec_array.end(); ++it) {
        (*it)->computePower();
    }
    return;
}

Carl Poelking's avatar
Carl Poelking committed
90
91
92
93
94
95
void Spectrum::computePowerGradients() {
    for (auto it = _atomspec_array.begin(); it != _atomspec_array.end(); ++it) {
        (*it)->computePowerGradients();
    }
}

Carl Poelking's avatar
Carl Poelking committed
96
97
98
99
100
101
102
103
AtomicSpectrum *Spectrum::computeAtomic(Particle *center) {
    return this->computeAtomic(center, _structure->particles());
}

AtomicSpectrum *Spectrum::computeAtomic(Particle *center, Structure::particle_array_t &targets) {
    GLOG() << "Compute atomic spectrum for particle " << center->getId()
        << " (type " << center->getType() << ", targets " << targets.size() << ") ..." << std::endl;

104
    // FIND IMAGE REPITIONS REQUIRED TO SATISFY CUTOFF
105
106
107
108
109
    vec box_a = _structure->getBoundary()->getBox().getCol(0);
    vec box_b = _structure->getBoundary()->getBox().getCol(1);
    vec box_c = _structure->getBoundary()->getBox().getCol(2);

    double rc = _basis->getCutoff()->getCutoff();
110
111
112
113
    std::vector<int> na_nb_nc = _structure->getBoundary()->calculateRepetitions(rc);
    int na_max = na_nb_nc[0];
    int nb_max = na_nb_nc[1];
    int nc_max = na_nb_nc[2];
114

115
116
117
    //GLOG() << box_a << " " << box_b << " " << box_c << std::endl;
    //GLOG() << rc << std::endl;
    //GLOG() << na_max << " " << nb_max << " " << nc_max << std::endl;
118

119
    // CREATE BLANK
Carl Poelking's avatar
Carl Poelking committed
120
121
122
    AtomicSpectrum *atomic_spectrum = new AtomicSpectrum(center, this->_basis);

    Structure::particle_it_t pit;
123
    for (pit = targets.begin(); pit != targets.end(); ++pit) { // TODO Consider images
Carl Poelking's avatar
Carl Poelking committed
124
125

        // CHECK FOR EXCLUSIONS
126
127
        if (_options->doExcludeTarget((*pit)->getType()) ||
            _options->doExcludeTargetId((*pit)->getId())) continue;
Carl Poelking's avatar
Carl Poelking committed
128

129
130
131
132
133
134
135
    for (int na=-na_max; na<na_max+1; ++na) {
    for (int nb=-nb_max; nb<nb_max+1; ++nb) {
    for (int nc=-nc_max; nc<nc_max+1; ++nc) {

        //GLOG() << na << " " << nb << " " << nc << std::endl;
        vec L = na*box_a + nb*box_b + nc*box_c;

Carl Poelking's avatar
Carl Poelking committed
136
        // FIND DISTANCE & DIRECTION, CHECK CUTOFF
137
        vec dr = _structure->connect(center->getPos(), (*pit)->getPos()) + L;  // TODO Consider images
Carl Poelking's avatar
Carl Poelking committed
138
        double r = soap::linalg::abs(dr);
Carl Poelking's avatar
Carl Poelking committed
139
        if (! this->_basis->getCutoff()->isWithinCutoff(r)) continue;
Carl Poelking's avatar
Carl Poelking committed
140
        vec d = (r > 0.) ? dr/r : vec(0.,0.,1.);
Carl Poelking's avatar
Carl Poelking committed
141

Carl Poelking's avatar
Carl Poelking committed
142
        // APPLY CUTOFF (= WEIGHT REDUCTION)
143
144
        bool is_image = (*pit == center);
        bool is_center = (*pit == center && na==0 && nb==0 && nc==0); // TODO Consider images
Carl Poelking's avatar
Carl Poelking committed
145
146
        double weight0 = (*pit)->getWeight();
        double weight_scale = _basis->getCutoff()->calculateWeight(r);
Carl Poelking's avatar
Carl Poelking committed
147
        if (is_center) {
Carl Poelking's avatar
Carl Poelking committed
148
            weight0 *= _basis->getCutoff()->getCenterWeight();
Carl Poelking's avatar
Carl Poelking committed
149
150
        }

Carl Poelking's avatar
Carl Poelking committed
151
        GLOG() << (*pit)->getType() << " X " << dr.getX() << " Y " << dr.getY() << " Z " << dr.getZ() << " W " << (*pit)->getWeight() << " S " << (*pit)->getSigma() << std::endl;
152

Carl Poelking's avatar
Carl Poelking committed
153
        // COMPUTE EXPANSION & ADD TO SPECTRUM
154
        bool gradients = (is_image) ? false : _options->get<bool>("spectrum.gradients");
Carl Poelking's avatar
Carl Poelking committed
155
156
        BasisExpansion *nb_expansion = new BasisExpansion(this->_basis); // <- kept by AtomicSpectrum
        nb_expansion->computeCoefficients(r, d, weight0, weight_scale, (*pit)->getSigma(), gradients);
157
        atomic_spectrum->addQnlmNeighbour(*pit, nb_expansion); // TODO Consider images
158
159
160

    }}} // Close loop over images
    } // Close loop over particles
Carl Poelking's avatar
Carl Poelking committed
161
162

    return atomic_spectrum;
163
164
}

165
166
167
168
AtomicSpectrum *Spectrum::computeGlobal() {
    if (_global_atomic) throw soap::base::APIError("<Spectrum::computeGlobal> Already initialised.");
    _global_atomic = new AtomicSpectrum(_basis);
    GLOG() << "Computing global spectrum ..." << std::endl;
Carl Poelking's avatar
Carl Poelking committed
169
    bool gradients = _options->get<bool>("spectrum.gradients");
170
171
172
    for (auto it = _atomspec_array.begin(); it != _atomspec_array.end(); ++it) {
        GLOG() << "  Adding center " << (*it)->getCenter()->getId()
            << " (type " << (*it)->getCenter()->getType() << ")" << std::endl;
Carl Poelking's avatar
Carl Poelking committed
173
        _global_atomic->mergeQnlm(*it, 0.5, gradients); // <- Scale factor 0.5 necessary in order not to overcount pairs
174
175
    }
    _global_atomic->computePower();
Carl Poelking's avatar
Carl Poelking committed
176
    if (gradients) _global_atomic->computePowerGradients();
177
178
179
    return _global_atomic;
}

Carl Poelking's avatar
Carl Poelking committed
180
AtomicSpectrum *Spectrum::getAtomic(int slot_idx, std::string center_type) {
Carl Poelking's avatar
Carl Poelking committed
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
	AtomicSpectrum *atomic_spectrum = NULL;
	// FIND SPECTRUM
	if (center_type == "") {
		// NO TYPE => ACCESS ARRAY
		// Slot index valid?
		if (slot_idx < _atomspec_array.size()) {
			atomic_spectrum = _atomspec_array[slot_idx];
		}
		else {
			throw soap::base::OutOfRange("Spectrum slot index");
		}
	}
	else {
		// TYPE => ACCESS TYPE MAP
		map_atomspec_array_t::iterator it = _map_atomspec_array.find(center_type);
		// Any such type?
		if (it == _map_atomspec_array.end()) {
			throw soap::base::OutOfRange("No spectrum of type '" + center_type + "'");
		}
		else {
			// Slot index valid?
			if (slot_idx < it->second.size()) {
				atomic_spectrum = it->second[slot_idx];
			}
			else {
				throw soap::base::OutOfRange("Spectrum slot index, type '" + center_type + "'");
			}
		}
	}
Carl Poelking's avatar
Carl Poelking committed
210
211
212
213
214
	return atomic_spectrum;
}

void Spectrum::writeDensityOnGrid(int slot_idx, std::string center_type, std::string density_type) {
	AtomicSpectrum *atomic_spectrum = this->getAtomic(slot_idx, center_type);
Carl Poelking's avatar
Carl Poelking committed
215
216
	// WRITE CUBE FILES
	if (atomic_spectrum) {
217
218
		//atomic_spectrum->getQnlm(density_type)->writeDensityOnGrid(
	    //	"density.expanded.cube", _options, _structure, atomic_spectrum->getCenter(), true);
Carl Poelking's avatar
Carl Poelking committed
219
		atomic_spectrum->getQnlm(density_type)->writeDensityOnGrid(
Carl Poelking's avatar
Carl Poelking committed
220
221
222
223
224
			"density.explicit.cube", _options, _structure, atomic_spectrum->getCenter(), false);
	}
	return;
}

225
226
227
228
229
230
231
232
233
234
235
236
void Spectrum::writeDensityCubeFile(int atom_idx, std::string density_type, std::string filename, bool from_expansion) {
    if (atom_idx < _atomspec_array.size()) {
        AtomicSpectrum *atomic_spectrum = this->_atomspec_array[atom_idx];
		atomic_spectrum->getQnlm(density_type)->writeDensityOnGrid(
			filename, _options, _structure, atomic_spectrum->getCenter(), from_expansion);
    }
    else {
        throw soap::base::OutOfRange("Spectrum::writeDensityCubeFile");
    }
    return;
}

Carl Poelking's avatar
Carl Poelking committed
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
void Spectrum::writeDensityOnGridInverse(int slot_idx, std::string center_type, std::string type1, std::string type2) {
	AtomicSpectrum *atomic_spectrum = this->getAtomic(slot_idx, center_type);
	AtomicSpectrum inverse_atomic_spectrum(atomic_spectrum->getCenter(), atomic_spectrum->getBasis());
	inverse_atomic_spectrum.invert(atomic_spectrum->getXnklMap(), atomic_spectrum->getXnklGenericCoherent(), type1, type2);
	inverse_atomic_spectrum.getQnlm("")->writeDensityOnGrid(
	    "density.inverse.cube", _options, _structure, inverse_atomic_spectrum.getCenter(), true);
	inverse_atomic_spectrum.getQnlm("")->writeDensity(
		"density.inverse.coeff", _options, _structure, inverse_atomic_spectrum.getCenter());
	return;
}

void Spectrum::writeDensity(int slot_idx, std::string center_type, std::string density_type) {
	AtomicSpectrum *atomic_spectrum = this->getAtomic(slot_idx, center_type);
	// WRITE COEFF FILE
	if (atomic_spectrum) {
		atomic_spectrum->getQnlm(density_type)->writeDensity(
			"density.expanded.coeff", _options, _structure, atomic_spectrum->getCenter());
	}
	return;
}

void Spectrum::writePowerDensity(int slot_idx, std::string center_type, std::string type1, std::string type2) {
	AtomicSpectrum *atomic_spectrum = this->getAtomic(slot_idx, center_type);
	if (atomic_spectrum) {
		AtomicSpectrum::type_pair_t types(type1, type2);
		atomic_spectrum->getXnkl(types)->writeDensity(
			"density.power.coeff", _options, _structure, atomic_spectrum->getCenter());
	}
	return;
}

void Spectrum::addAtomic(AtomicSpectrum *atomspec) {
Carl Poelking's avatar
Carl Poelking committed
269
270
271
272
273
274
275
276
277
278
	assert(atomspec->getBasis() == _basis &&
		"Should not append atomic spectrum linked against different basis.");
	std::string atomspec_type = atomspec->getCenterType();
	map_atomspec_array_t::iterator it = _map_atomspec_array.find(atomspec_type);
	if (it == _map_atomspec_array.end()) {
		_map_atomspec_array[atomspec_type] = atomspec_array_t();
		it = _map_atomspec_array.find(atomspec_type);
	}
	it->second.push_back(atomspec);
	_atomspec_array.push_back(atomspec);
279
	return;
Carl Poelking's avatar
Setup.  
Carl Poelking committed
280
281
}

Carl Poelking's avatar
Carl Poelking committed
282
283
284
285
286
287
288
void Spectrum::save(std::string archfile) {
	std::ofstream ofs(archfile.c_str());
	boost::archive::binary_oarchive arch(ofs);
	arch << (*this);
	return;
}

289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
std::string Spectrum::saves(bool prune) {
    if (prune) {
        // Delete all pid-resolved data from atomic spectra.
        // Note that this data is required to compute gradients.
        // Pruning will no longer allow computing gradients
        // for serialised objects.
        for (auto it = _atomspec_array.begin(); it != _atomspec_array.end(); ++it) {
            (*it)->prunePidData();
        }
        if (_global_atomic) {
            _global_atomic->prunePidData();
        }
    }
    std::stringstream bstream;
    boost::archive::binary_oarchive arch(bstream);
    arch << (*this);
    return bstream.str();
}

Carl Poelking's avatar
Carl Poelking committed
308
309
310
311
312
313
314
void Spectrum::load(std::string archfile) {
	std::ifstream ifs(archfile.c_str());
	boost::archive::binary_iarchive arch(ifs);
	arch >> (*this);
	return;
}

315
316
317
318
319
320
321
322
Spectrum &Spectrum::loads(std::string bstr) {
    std::stringstream bstream;
    bstream << bstr;
    boost::archive::binary_iarchive arch(bstream);
    arch >> (*this);
    return (*this);
}

Carl Poelking's avatar
Setup.  
Carl Poelking committed
323
324
void Spectrum::registerPython() {
    using namespace boost::python;
Carl Poelking's avatar
Carl Poelking committed
325
    void (Spectrum::*computeAll)() = &Spectrum::compute;
Carl Poelking's avatar
Carl Poelking committed
326
    void (Spectrum::*computeSeg)(Segment*) = &Spectrum::compute;
Carl Poelking's avatar
Carl Poelking committed
327
328
329
    void (Spectrum::*computeSegPair)(Segment*, Segment*) = &Spectrum::compute;
    void (Spectrum::*computeCentersTargets)(Structure::particle_array_t&, Structure::particle_array_t&) = &Spectrum::compute;

Carl Poelking's avatar
Setup.  
Carl Poelking committed
330
    class_<Spectrum>("Spectrum", init<Structure &, Options &>())
Carl Poelking's avatar
Carl Poelking committed
331
    	.def(init<Structure &, Options &, Basis &>())
Carl Poelking's avatar
Carl Poelking committed
332
    	.def(init<std::string>())
333
        .def(init<>())
334
    	.def("__iter__", range<return_value_policy<reference_existing_object> >(&Spectrum::beginAtomic, &Spectrum::endAtomic))
335
        .def("__len__", &Spectrum::length)
Carl Poelking's avatar
Carl Poelking committed
336
	    .def("compute", computeAll)
Carl Poelking's avatar
Carl Poelking committed
337
        .def("compute", computeSeg)
Carl Poelking's avatar
Carl Poelking committed
338
339
	    .def("compute", computeSegPair)
	    .def("compute", computeCentersTargets)
340
		.def("computePower", &Spectrum::computePower)
Carl Poelking's avatar
Carl Poelking committed
341
		.def("computePowerGradients", &Spectrum::computePowerGradients)
342
        .def("deleteGlobal", &Spectrum::deleteGlobal)
343
		.def("computeGlobal", &Spectrum::computeGlobal, return_value_policy<reference_existing_object>())
Carl Poelking's avatar
Carl Poelking committed
344
345
		.def("addAtomic", &Spectrum::addAtomic)
		.def("getAtomic", &Spectrum::getAtomic, return_value_policy<reference_existing_object>())
346
		.def("getGlobal", &Spectrum::getGlobal, return_value_policy<reference_existing_object>())
Carl Poelking's avatar
Carl Poelking committed
347
	    .def("saveAndClean", &Spectrum::saveAndClean)
Carl Poelking's avatar
Carl Poelking committed
348
349
		.def("save", &Spectrum::save)
		.def("load", &Spectrum::load)
350
351
        .def("saves", &Spectrum::saves)
        .def("loads", &Spectrum::loads, ref_existing())
Carl Poelking's avatar
Carl Poelking committed
352
        .def("writeDensityOnGrid", &Spectrum::writeDensityOnGrid)
353
        .def("writeDensityCubeFile", &Spectrum::writeDensityCubeFile)
Carl Poelking's avatar
Carl Poelking committed
354
355
356
357
358
359
		.def("writeDensity", &Spectrum::writeDensity)
		.def("writePowerDensity", &Spectrum::writePowerDensity)
		.def("writeDensityOnGridInverse", &Spectrum::writeDensityOnGridInverse)
		.add_property("options", make_function(&Spectrum::getOptions, ref_existing()))
		.add_property("basis", make_function(&Spectrum::getBasis, ref_existing()))
		.add_property("structure", make_function(&Spectrum::getStructure, ref_existing()));
360
361
    class_<atomspec_array_t>("AtomicSpectrumContainer")
           .def(vector_indexing_suite<atomspec_array_t>());
Carl Poelking's avatar
Setup.  
Carl Poelking committed
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
}

/* STORAGE, BASIS, COMPUTATION, PARALLELIZATION */
/*
 *
 * Spectrum, PowerSpectrum
 * SpectrumK, PowerSpectrumK
 * StructureK, ParticleK (with weight, type)
 *
 * Basis
 * BasisFactory
 *
 * RadialBasis
 * RadialBasisLegendre
 * RadialBasisGaussian
 * RadialBasisHermite
 *
 * AngularBasis
 * AngularBasisHarmonic
 */


/*
 * Parallelization
 *     based on centers
 *     based on wavevectors
 */

/*
 * Spectrum->Setup(system, options)
 *     -> Basis->Setup()
 * Compute(patterns)
 *     for pattern in patterns:
 *        single out neighbours that match pattern
 *        for each neighbour:
 *            -> Compute(pos, refpos)
 */

/*
 * PowerSpectrum->Setup(Expansion)
 *      [a][b][n][n'][l] = sum{m} [a][nlm]*[b][n'lm]
 *
 * Spectrum in linear storage:
 * Start n at offset_n = (n-1)*(L+1)^2
 *     Start l at offset_l = offset_n + l^2
 *         Start m at offset_m = offset_l + (m+l)
 *
 * PowerSpectrum in linear storage:
 * Start n at offset_n = (n-1)*N*(L+1)
 *     Start n' at offset_n' = offset_n + (n'-1)*(L+1)
 *         Start l at offset_l = offset_n'
 *
 * map< pattern1, spectrum  >
 * for (nit in spectrum)
 *     for (lit in nit)
 *         for (mit in lit)
 *             nit->n
 *             lit->l
 *             mit->m
 *
 *
 */

/*
 * Storage:
 *     -> Serialize()
 *     For each center:
 *         for each pattern:
 *             [n][l][m]
 *     Store [n] and [l][m] separately?
 *
 *     -> PowerSerialize()
 *     For each center:
 *         for each pair of patterns:
 *             [n][n'][l]
 */

/*
 *
 * Types of basis functions: Legendre, Hermite, Gaussian*
 * Make sure to include appropriate n-factors, e.g. sqrt(2*n+1) for Legendre
 * Normalized on which interval?
 *
 * e.g. RadialBasisLegendre : RadialBasis
 *
 * RadialBasis->Setup()
 *     -> Orthogonalize()
 * RadialBasis -> EvaluateAt(radius)
 *     -> AtomicBasis->EvaluateAt(radius)
 *     -> Transform coefficients
 *     -> Return coefficients (N-vector)
 * RadialBasis->Iterator() ?
 *
 */

/*
 * AngularBasis->Setup()
 * AngularBasis->EvaluateAt(direction)
 *     -> AtomicBasisL->EvaluateAt(radius)
 *          -> AtomicBasisLM -> EvaluateAt(radius)
 *     -> Return coefficients (L*M-vector)
 * AngularBasis->Iterator()
 */

}