extfieldsofstudy: Mathematics Physics Biology

text:
Cryo electron microscopy: 3-D reconstruction of stochastic objects having statistics with geometric symmetry from 2-D noisy conditionally-linear imaging

The images obtained from cryo electron microscopy are the projection of many heterogeneous instances of the object under study (e.g., a virus). When the object is heterogeneous but has an overall geometric symmetry, it is natural to describe the object as stochastic with symmetrical statistics. A maximum likelihood reconstruction approach is described in which all instances of the object lack symmetry but the statistics of the ensemble of objects are constrained to have symmetry. This algorithm is demonstrated on bacteriophage HK97 and is contrasted with existing algorithms in which each object, while still heterogeneous, has the symmetry. Reconstruction results show that the proposed algorithm eliminates long-standing distortions in previous heterogeneity calculations associated with symmetry axes, and provides estimates that make more biological sense than the estimates of existing algorithms.

I. INTRODUCTION
T HE goal of structural biology is to determine the 3-D geometry of a biological macromolecular complex. Such complexes, for example, a virus, function as machines and are typically called "particles". The 3-D geometry, typically called the "structure", can be as detailed as the 3-D location of each atom in the particle or, if only lower resolution is achievable, the electron density as a function of position in 3-D (in xray crystallography) or the electron scattering intensity as a function of position in 3-D (in electron microscopy). Understanding the 3-D geometry of synthetic nanoparticles is also becoming important in materials science [1]- [3]. Structural biology has long been dominated by x-ray crystallography which has very limited capability for dealing with heterogeneity and dynamics among the particles. In particular, the particles making up the crystal must be sufficiently uniform to allow the formation of a crystal and then the additional non-biological interactions between neighboring particles in the crystal further restrict the individual particles. Finally the data is a function of the entire crystal, i.e., a function of all of the particles. Therefore x-ray crystallography provides a 3-D snapshot of a particular state of the particle rather than following a single instance of the particle through the 4-D trajectory that it would experience in its natural environment.
Single-particle cryo electron microscopy (cryo EM) is transforming structural biology [4] because, due to better image detectors [5], phase plates [6] for improved image contrast, and more sophisticated statistical computing [7]- [9], cryo EM is beginning to provide atomic-resolution structures similar to those of x-ray crystallography without the need to create a crystal. This is a very important advance because many particles simply do not form crystals for a variety of biophysical reasons. However, more important for understanding the functioning of the particles, cryo EM has the capability of displaying heterogeneity and dynamics among the particles.
In a cryo EM experiment, typically 10 4 -10 6 noisy 2-D projection images of the 3-D electron scattering intensity of the flash-frozen particles are recorded, one image for each instance of the particle. Typical particle weights are in the range of 10 5 -10 8 Daltons. The experimental design is dominated by the desire to avoid damage to the particle by the electron beam. In particular, damage minimization is achieved by low dose (leading to noisy images), by recording only one projection image per particle (unlike medical x-ray computed tomography), and by performing the experiment on flash-frozen particles at cryogenic temperatures. Because the orientation of the particle during the flash freezing is not controlled, the orientation of the projection image is random and because the image is noisy, the projection direction cannot be determined from an individual image alone. Standard computational methods assume that all instances of the particle are identical or that all instances of the particle belong to one of a few discrete classes ("discrete heterogeneity", the class label for each image is not known, typically there are less than ten classes) and that every instance within a class is identical. Then a structure is computed for each class, possibly at atomic resolution. More sophisticated approaches characterize the continuous heterogeneity of the particles within each class.
The geometric symmetry of the particle is an important aspect of structural studies of both biological and synthetic particles. Many so-called spherical viruses, which include polio virus and other human pathogens, have a so-called capsid that is constructed out of many copies of a peptide building block which are arranged in pattern that has icosahedral symmetry much like a human-made geodesic dome. Symmetry is described mathematically as a rotational point group. This paper describes how to impose symmetry on the statistics of the particle rather than on each realization of the particle (Sections II, III, IV, and VII) and how to integrate this mathematical model of the particle into a maximum likelihood reconstruction method (Section V). The result is demonstrated on images of a Virus Like Particle of bacteriophage HK97 (Section VI).
The approach of this paper is focused on the estimation of first and second order statistics of the electron scattering intensity of the particles in order to characterize the heterogeneity of the particles and related work is discussed in Section II. However there are other approaches to characterizing the differences among the particles visualized by cryo EM. Prominent approaches include direct reconstruction of principal components [10], low-dimensional manifold embedding of multiple reconstructions with interpolation to create trajectories [11], [12], and normal modes of a references model [13]. These papers do not concern symmetry (the numerical examples are all ribosome problems) and so are outside of the scope of this paper.
Notation: E denotes expectation, i.i.d. stands for "independent and identically distributed", x ∈ IR 3 , and x = x 2 . The quantity x/x is a unit vector and, since x ∈ IR 3 , is shorthand for the two angles of spherical coordinates. I n is the n × n identity matrix and 0 m,n is the m × n zero matrix. dΩ is integration over the surface of the sphere in 3-D which, in spherical coordinates, is π θ=0 2π φ=0 sin(θ)dθdφ. Define the Kronecker delta function δ m,n by δ m,n = 1 if m = n and = 0 otherwise where m and n are integers or vectors of integers.

II. SYMMETRY AND SYMMETRICAL STATISTICS
The symmetries that are relevant to the cryo EM application are the 3-D rotational point group symmetries which are the icosahedral I, octahedral O, tetrahedral T , cyclic C n (n ∈ {1, 2, . . . }), and dihedral D n (n ∈ {1, 2, . . . }) groups. Once a coordinate system is chosen, each of these groups can be described as a collection of matrices denoted by R β ∈ IR 3×3 (β ∈ {1, . . . , N g }) which are rotation matrices, i.e., R −1 β = R T β and det R β = +1. The value of N g is 60, 24, 12, n, and 2n for I, O, T , C n , and D n , respectively.
(2) Assume that each particle in a class is different but that each particle has the symmetry so that the electron scattering intensity of the ith particle, which is denoted by ρ i (x), satisfies , . . . , N g } and x ∈ IR 3 . This is available in some software systems by resampling [17]- [21] (implemented, for instance, in Ref. [7]), moment estimators [22]- [24], or maximum likelihood estimation of the parameters of a probabilistic model for the Fourier series coefficients of the electron scattering intensity [25]- [30]. It is known that the resulting spatial variance functions (i.e., the variance as a function of position) have anomalous peaks on and near symmetry axes of the particle [31, p. 173]. Since the symmetry axes can be the location of important biology, e.g., in Flock House Virus (FHV) [32], [33], the biology concerning the pentamer of γ peptides near the 5-fold symmetry axis of the icosahedrally symmetric particle [34], [35], so anomalous results near symmetry axes is not a satisfactory situation.
(3) Assume that each particle in a class is different, that no particle has the symmetry, but the statistics of the i.i.d. ensemble of particles has the symmetry as is described in the following paragraph. The authors are unaware of prior work of this type.

III. BASIS FUNCTIONS
The particle is described as a linear combination of basis functions, specifically, where c ω are random variables and F ω (x) are known functions. The goal of the reconstruction calculation (Section V) is to estimate the joint probability distribution of the c ω . To achieve the desired constraints on the first and second order statistics of ρ(x) (Eqs. 1 and 2) requires constraints on the first and second order statistics of c ω . Determining these constraints is simplified if F ω (x) have specific properties under rotations from the symmetry point group of the particle. The numerical example of Section VI concerns a virus with icosahedral (I) symmetry. For I (and also for O and T ), the rotational symmetry axes are oriented in a variety of directions leading to particles that tend to fill their minimum bounding sphere. Therefore the choice of F ω (x) that is considered in detail uses spherical coordinates and is the product of a so-called angular function (which depends on the angles of spherical coordinates and is a linear combination of spherical harmonics 1 ) and a so-called radial function (which depends on the radius of spherical coordinates), i.e., F ω (x) = I ω (x/x)ψ ω (x) where x = x 2 and ω is an index. For C n and possibly D n , it is possible to use the same framework but, especially for C n , it might be more natural to use cylindrical coordinates and a product of a so-called axial function (in z), a polar function (in φ), and a radial function (in r). Because rotations of x do not effect x = x , in this section we focus on I ω (x/x).
Ref. [38] derives a basis I ω (x/x) that has precise properties under the operations of a rotational point group. In order to describe the properties of the basis, it is necessary to describe the the idea of representations. A representation (rep) of a rotational point group, which is a finite group, is a set of matrices which are homomorphic under matrix multiplication to the group elements [39, p. 61]. The irreducible representations (irreps) of a finite group are a set of unitary reps such that all unitary reps can be decomposed into a direct sum of the irreps [39, Theorem 3.3 p. 69]. For the icosahedral and octahedral cases, all irreps can be chosen to be real-valued orthonormal rather than complex-valued unitary and, because of the biological importance of the icosahedral case, the remainder of this section focuses on the real-valued orthonormal icosahedral case. Let N rep be the number of irreps of the N g -order group G. Let the set of real-valued matrices in the pth irrep be denoted by Γ p (g) ∈ R dp×dp for all g ∈ G where p ∈ {1, . . . , N rep }. For the icosahedral group, N rep = 5, N g = 60, and d p = 1, 3, 3, 4, 5 for p = 1, 2, 3, 4, 5, respectively. The basis functions I p,ω (x/x) were derived [38] by applying real-valued generalized projection operators constructed from real-valued irreps to real-valued spherical harmonics resulting in the index ω being the ordered pair (l, n) where l indexes the subspace defined by spherical harmonics of order l which contains the basis function and n indexes the different basis functions within the subspace. The basis functions I p;l,n (x/x) have the properties 1) Each I p;l,n is a d p -dimensional real-valued vector function, i.e., I p;l,n ∈ R dp . 2) The I p;l,n functions are orthonormal on the surface of the sphere, i.e., I p,l,n (x/x)I T p ,l ,n (x/x)dΩ = δ p,p δ l,l δ n,n I dp (4) 3) The subspace of square integrable functions on the surface of the sphere defined by spherical harmonics of index l, contains a set of I p;l,n functions with a total of 2l + 1 components. 4) Each I p;l,n function has a specific transformation property under rotations from the rotational point group [40, p. 20], in particular, where T is transpose not Hermitian transpose and R g is the 3×3 rotation matrix corresponding to the gth element of the group. Examples of I p,j;l,n functions for l = 15 are visualized in Figure 1. Note that p = 1 exhibits all of the symmetries of an icosahedron.
and 0 otherwise, where κ 1 and κ 2 are chosen so that κ 1 + κ 2 I p;l,n (x/x) varies between 0.5 and 1.0, are visualized by UCSF Chimera [41] where the color indicates the distance from the center of the object.

IV. CONSTRAINTS ON THE MOMENTS OF THE WEIGHTS CAUSED BY SYMMETRIC STATISTICS
In this section we suppress the index η which describes the class of the particle. The electron scattering intensity, denoted by ρ(x), is described as a linear combination of basis functions where the basis functions are denoted by I p,ζ (x/x) (p indexes the irrep, ζ is a shorthand for all other indices, and x = x 2 ) and take values in IR dp (d p is the dimension of the pth irrep) and the weights in the linear combination are denoted by c p,ζ (x), also take values in IR dp , and depend on x = x 2 . In summary, The electron scattering intensity ρ(x) is real-valued. In this section we deal with the situation where there are a complete set of real-valued orthonormal irreps (including the important icosahedral case) so that I p,ζ (x/x) can be taken as real and therefore c p,ζ can also be taken as real. In Section VII the case where I p,ζ (x/x) and c p,ζ are complex-valued is described. The constraint on the mean function of and from Eqs. 6 and 5 (g ∈ {1, . . . , N g }) From Eq. 1 it follows that the left hand sides of Eqs. 7 and 8 are equal so that for all g ∈ {1, . . . , N g }, Multiplying by I T p ,ζ (x/x) and integrating over the surface of the sphere implies [via the orthogonality of the Taking transposes and renaming p to p and ζ to ζ gives that for all g ∈ {1, . . . , N g }, The identity irrep is the p = 1 case and for this case, d p=1 = 1 and Γ p=1 (g) = 1 for all g ∈ {1 . . . , N g } so thatc p=1,ζ (x) is unconstrained by Eq. 11. Eq. 11 implies that for all g ∈ {1, . . . , N g }, Consider p ∈ {2, . . . , N rep } such that d p > 1. Ifc p,ζ (x) = 0 dp,1 then there is at least a one-dimensional subspace which is invariant under the action of Γ p (g) for all g ∈ {1, . . . , N g } and therefore, by definition [39, p. 67], Γ p (g) is not an irrep. But this contradicts the assumption that Γ p (g) is an irrep. Therefore, the statementc p,ζ (x) = 0 dp,1 must be false. Finally, consider p ∈ {2, . . . , N rep } such that d p = 1. Ifc p,ζ (x) = 0 then Eq. 12 implies that Γ p (g) = 1 for all g ∈ {1, . . . , N g } but this is the identity irrep which is the p = 1 irrep. Therefore, c p,ζ (x) = 0. In summary, The constraint on the covariance function of c p,ζ (x) is derived in this paragraph acting as if c p,ζ (x), I p,ζ (x/x), and Γ p (g) are complex valued which simplifies the presentation in Section VII.
Eq. 29 implies structure for the complete covariance matrix V . Suppose that the basis functions are enumerated with indices changing in the order p (slowest), ζ, q, and j (fastest) where j is the index that enumerates the elements of the vector I p,ζ (x/x). Then the V matrix is a block matrix with N rep ×N rep blocks (determined by the values of the p indices) and only the diagonal blocks are nonzero (due to the "p 1 = p 2 " condition in Eq. 29). Furthermore, each block is constructed of N ζ,q × N ζ,q subblocks (where N ζ,q is the total number of (ζ, q) pairs) and each subblock is proportional to the d p × d p identity matrix with proportionality constant v p1 (ζ 1 , q 1 ; ζ 2 , q 2 ).

V. HETEROGENEOUS RECONSTRUCTION WITH SYMMETRICAL STATISTICS
This section describes the maximum likelihood estimator that computes the first and second order statistics of the weights in the orthonormal expansion of the electron scattering intensity ρ(x) from the image data under the constraint of symmetrical statistics for ρ(x) (Eqs. 1 and 2).
The standard imaging model is that the 2-D image is the projection along the axis of the microscope of the 3-D electron scattering intensity of the particle linearly modified by the so-called contrast transfer function of the microscope and observed with an unknown translation in additive noise. Let N v be the number of 2-D projection images. For the ith particle and ith image, let c i be the random variable coefficients in the orthonormal expansion for ρ i (x). Since ρ i (x) is a linear function of c i and the standard imaging model is a linear function of ρ i (x), there is a linear relationship between the ith image and c i , which can be written in the spatial frequency domain (so-called reciprocal space) by the projection slice theorem. Specifically [25], where y i is the ith 2-D projection image in reciprocal space, z i is a vector of nuisance parameters which include the projection orientation of the ith image and the projected location of the particle in the ith image, c i ∈ R Nc is the vector of weights of the ith image, L(z i ) is the transformation matrix from weights to the ith reciprocal-space image, and w i is the noise which corrupts the ith image. We assume that the noise w i and the weights c i are jointly Gaussian and independent with pdfs w i ∼ N (0, Q) and c i ∼ N (c, V ). A maximum likelihood estimator is used to estimatec and V from the image data, as is described in Ref. [25]. The maximum likelihood estimator is computed by an expectationmaximization algorithm with nuisance parameters {z i } Nv i=1 . The jointly Gaussian assumption on c i and w i is quite important because then the stochastic part of c i can be combined with w i , specifically, where from the set of nuisance parameters at the cost of having a noise w i with a structured covariance matrix. Alternating updates ofc and V are used. Additional algorithmic results are described in Ref. [25] including results for multiple discrete classes of particle where the a priori distribution is to be estimated.
Our previous calculations [25]- [27], [29] were for icosahedral symmetry and used only the p = 1 subset of the basis functions described in Section III. By direct calculation or based on the results in Section IV, no constraints onc or V are required in order to achieve Eqs. 1 and 2. The optimization of V uses Matlab's fmincon (option "trust-region-reflective") with symbolic cost, gradient of the cost, and Hessian of the cost. In order to incorporate the complete set of basis functions (i.e., all values of p) into this software [42], we think of the mean vector and covariance matrix as functions of a further level of parameterization. The further level of parameterization is quite simple. Specifically, the mean is either arbitrary or zero (Eq. 25) so the further parameters are the arbitrary values while the covariance has certain elements that much be equal but are otherwise arbitrary (Eq. 29) so the further parameters are the arbitrary values. Using the additional level of parameterization, the gradient and Hessian can be computed via the chain rule. These equations compute the necessary gradient and Hessian in terms of larger vectors and matrices which are then reduced in size. While this approach fits the software easily, more efficient approaches may be possible.

VI. RECONSTRUCTION RESULTS
Numerical results on both synthetic and experimental images of a Virus Like Particle (VLP) derived from bacteriophage HK97 [43] are described in this section. The VLP is essentially the bacteriophage minus the bacteriophage's portal, tail, and genome leaving only the icosahedrally symmetric capsid. HK97 has a complicated lifecycle wherein it first selfassembles in near equilibrium conditions and then undergoes a sequence of essentially irreversible maturation transformations which result in a robust particle capable of surviving outside of the host cell. One of the first steps of maturation is the digestion by a virally-encoded protease of the so-called δ domain of the 60×7 copies of the capsid peptide that together make up the capsid of the bacteriophage. The experimental images [44] come from the particle denoted by PhI Pro+ [29] which contains a protease that is defective so that the particle is trapped in the Prohead I step of maturation. The average outer radius of the capsid is 254Å.
Both the synthetic and experimental image calculations are performed in two steps. First a homogeneous reconstruction, i.e., a reconstruction in which all particles are assumed to be identical and have icosahedral symmetry, is performed. In terms of Section V, this is the case where p = 1, c i =c, and V = 0 Nc,Nc so this algorithm is denoted by Hetero (zero covar). Second, a heterogeneous reconstruction, i.e., a reconstruction in which no particles are the same, is performed using the homogeneous reconstructionc as an initial condition. Two types of heterogeneous reconstruction are computed and contrasted, one uses only the p = 1 basis functions (Hetero (SymPart)) so that each particle has the icosahedral symmetry while the other uses the full set of p ∈ {1, . . . , N rep = 5} basis functions with the constraints of Section IV (Hetero (SymStat)) so that no particle has the icosahedral symmetry but the first and second order statistics of the electron scattering intensity ρ(x) do have the icosahedral symmetry. The Hetero (zero covar) followed by Hetero (SymPart) reconstructions used basis functions with p = 1, l ≤ 55, and q ≤ 20 (1060 total basis functions) while the Hetero (zero covar) followed by Hetero [47, p. 879]) between two structures (the two estimates of mean electron scattering intensity) computed from disjoint sets of images. This measure is applied in all calculations, evaluated by the formulas of Ref. [48,: FSC(k) = S P1,P2 (k)/ S P1,P1 (k)S P2,P2 (k) where S P1,P2 (k) = P 1 (k)P * 2 (k)dΩ and the integral can be done symbolically [48,Eqs. 23].

A. Simulated images
Simulated 2-D projection images directly in reciprocal space were generated from Eq. 30 which is the mathematical model used in the reconstruction algorithm. There are multiple possible sources forc, including from atomic resolution coordinates from x-ray crystallography experiments, but few sources for V . The values ofc and V used in the simulation ("ground truth") come from the Hetero (SymStat) reconstruction results. Each of the 1200 images measures 100 × 100 pixels with a pixel sampling interval of 5.52Å and a SNR of 0.5 measured on variance not standard deviation. Examples of the simulated real-space 2-D images are shown in Fig. 2.
The sampling interval is 5.52Å so the Nyquist spatial frequency for the images is 0.5/5.52 = 0.0906Å −1 . Over this range the FSC curve based on reconstructions from the non-overlapping subsets of images is very close to 1 so that the structures are essentially identical. Due to computation limits imposed by our computer system, the number of basis functions that can be used is limited which restricts the spatial resolution that can be achieved. One measure of this restriction is that the spherically averaged energy (the denominator of FSC) is quite small by k = 0.0906Å −1 , approximately 10 −6 times it peak value.

B. Experimental images
Limitations of our computer system make it impossible to process the entire set of images in Ref. [44] so two nonoverlaping sets of 1200 images each showing one instance of the particle were randomly selected for processing. Each image measures 200 × 200 pixels with a pixel sampling interval of 2.76Å.
FSC was computed for the Hetero (SymPart) reconstructions from the non-overlapping subsets of images with the result that resolution is 0.061Å −1 (equivalently 16.4Å) based on the frequency at which the FSC curve crosses the 0.5 level. The same calculation for the Hetero (SymStat) reconstruction leads to an FSC curve that is almost 1 for the entire range of frequencies out to the Nyquist frequency for the images which is 0.5/2.76 = 0.181Å −1 since the pixel sampling interval is 2.76Å. Therefore these two structures are essentially identical.

Figs. 3(b) and (c)
show the standard deviation s ρ (x) computed by Hetero (SymPart) and Hetero (SymStat), respectively. In Fig. 4,ρ(x) computed by Hetero (SymPart) and the two heterogeneous solutions s ρ (x) computed by Hetero (SymPart) and Hetero (SymStat) are jointly visualized by computing a 3-D surface of constant value ofρ(x) and coloring the surface by s ρ (x). Theρ and s ρ functions in Figs. 3-4 have the same orientation as is shown in Fig. 3(a), wherein the x-, y-, and z-axes are all 2-fold rotational symmetry axes of the icosahedron which is one of two possible so-called 222 coordinate systems. Fig. 5 shows cross sections of the standard deviation function s ρ (x) with contours of the mean function ρ(x). Each cross section is either perpendicular to one of the three types of symmetry axis (5-, 3-, and 2-fold) or includes the symmetry axis as is described in the diagrams in the upper ofρ(x) colored by the standard deviation s ρ (x), which is visualized by UCSF Chimera [41]. All markings are scaled by 10 −3 .
row. Results are shown for both a Hetero (SymPart) and a Hetero (SymStat) reconstruction. The primary conclusion from Figs. 3-5 is that the assumption of symmetrical particles (i.e., Hetero (SymPart)) leads to unphysical peaks in the variance s ρ (x) near and on symmetry axes while these peaks are not present when using the assumption of symmetrical statistics (i.e., Hetero (SymStat)). Using a different approach, specifically, a resampling approach, but still assuming that each particle is symmetric, this same type of artifact is well-known to be a common result [31, p. 173] from standard software systems, e.g., EMAN2 [7]. The radially-directed artifacts in Figs. 3-5 are inconsistent with the annularly structured PhI Pro+ particle, which has, from outside to inside, annuli of well ordered capsid protein, poorly order capsid δ domain and protease protein, and poorly ordered DNA genome. The annular characteristic of the standard deviation function s ρ (x) from the Hetero (SymStat) reconstruction is particularly clear in Fig. 3(b) and even more so in Fig. 4(a). The annularly-directed artifacts are the origin of the light blue radial lines in Fig. 5(a) while the actual physical structure of the particle is the origin of the annular character of Fig. 5(b). In summary, the results based on assuming symmetric statistics (Hetero (SymStat)) rather than symmetric particles (Hetero (SymPart)) make more biological sense. Previous analyses based on Hetero (SymPart) primarily used the square root of the spherical average of the variance C(x, x) and therefore the peaks along radially-directed lines were averaged out leading to interpretable information. However, as we try to increase the spatial detail of our interpretation, such behavior is difficult to understand. for PhI Pro+ using both Hetero (SymPart) and Hetero (SymStat) and displayed with a common color map. The white lines depict contour levels of the meanρ(x). The geometrical orientation of the s ρ (x) function is shown in the top row and for each orientation, cross sections in the x-y and xz planes are displayed so that, for instance, a cross section perpendicular to a 5-fold axis and a cross section that contains a 5-fold axis are both shown. Panel (a) shows the reconstruction results based on symmetrical particles, i.e., Hetero (SymPart), whereas Panel (b) shows the reconstruction results based on symmetrical statistics, i.e., Hetero (SymStat).

VII. CASES WHERE REAL-VALUED IRREPS ARE NOT AVAILABLE AND NON-SPHERICAL COORDINATE SYSTEMS
Standard references for rotational point group irreps give the irrep matrices in complex-valued unitary form. In some instances, there is no similarity transformation to a realvalued orthonormal form [38]. While it is possible to consider representation theory over a real-valued vector space rather than a complex-valued vector space, the real-valued theory is more limited [39, p. 70], e.g., a matrix may not have any eigenvalues since the roots of the characteristic polynomial could potentially all be complex. In this section, the constraints on the moments of the weights caused by symmetric statistics are investigated when the irreps and therefore the basis functions and the weights are complex.
The results for the mean function of c p,ζ (x) are unaltered from Eqs. 7-13.
The results for the covariance function of c p,ζ (x) are changed from Eqs. 15-20. Because the variables are now complex, it is necessary to compute and to apply symmetry constraints to both C ρ,ρ (x 1 , x 2 ) and C ρ,ρ * (x 1 , x 2 ) which are defined by C ρ,ρ (x 1 , , respectively. The calculations for C ρ,ρ are the same as Eqs. 15-20. In a more detailed notation where C c,c p1,ζ1;p2,ζ2 (x 1 , It is not possible to continue on to Eq. 21 because Γ p (g) is no longer real and a new solution is described in the following paragraph. The calculations for C ρ,ρ * follow the same plan with the result that [C c,c * p1,ζ1;p2,ζ2 (x 1 , = c c,c * p1 (ζ 1 , x 1 ; ζ 2 , x 2 )I dp 1 , p 1 = p 2 0 dp 1 ,dp 2 , otherwise .
The calculations in this section allow ρ(x) to be complex. However, ρ(x) is real in the cryo EM application. A "penalty function" approach would be to create synthetic measurements of ρ(x) with value 0 and variance σ 2 and process the resulting complex-valued ρ(x), possibly allowing σ 2 to decrease as the iterative algorithm converges. A constraint approach would be to determine additional constraints on c p,ζ (x) such that ρ(x) is guaranteed to be real. We believe that such constraints involve the behavior of basis functions under conjugation (i.e., for Fourier series with period T in 1-D, ψ * n (t) = [exp(i2πnt/T )] * = exp(−i2πnt/T ) = ψ −n (t)) and are not closely related to the rotational point-group symmetries discussed in this paper.
Except in this paragraph, results in this paper are stated in the spherical coordinate system (e.g., Eq. 6). Spheres, which are contours of constant radius in spherical coordinates, enclose typical so-called spherical viruses reasonably closely and therefore lead to parsimonious parameterizations. This is typical for particles obeying the symmetries of I, O, and T because the symmetry axes have widely varying directions. However, for particles obeying the symmetries of C n (and possibly D n ), e.g., the tail of bacteriophage P22 [49], cylindrical coordinates (r, φ, z) can be more natural with the single symmetry axis of C n on the z axis of the coordinate system. For the case of C n , Eq. 6 is replaced by ρ(x) = Ng p=1 ζ c T p,ζ (r, z)I p,ζ (φ) and the constraint of symmetrical statistics specifies the mean and covariance functions of c p,ζ (r, z) rather than the mean (Eq. 13) and covariance (Eqs. 22, 34, and 36) functions of c p,ζ (x). The case of D n is somewhat more complicated because there are additional rotational symmetry axes beyond z axis. In either case, the effect of the symmetrical statistics constraint is a constraint on the covariance function of the coordinate-dependent weights where the weights depend on the coordinates that are not involved in the symmetry operation. Finally, if basis functions F p,ζ (x) are used which directly have symmetry properties that are the obvious generalization of Eq. 5, then the orthonormal expansion is ρ( with random-variable weights c p,ζ and there is no mean or covariance function intermediate on the path to mean vector and covariance matrix constraints analogous to Eqs. 25 and 29.

VIII. SYMMETRIZING NON-SYMMETRIC FUNCTIONS
Suppose s is a function from IR 3 to IR. Suppose R β (β ∈ {1, . . . , N g }) are the rotation matrices for a point-group symmetry. An important result in theory and software is that has the symmetry described by Eq. 1. The proof is via the Rearrangement Theorem [40, Theorem II, p. 24] which states that a sum over all the elements of a group, i.e., Ng β=1 f (R β ), is unaltered (although the terms in the sum occur in a permuted order) if R β is replaced by R β0 R β where β 0 is a fixed element of the group. The purpose of this section is to describe a similar symmetrization procedure for stochastic processes such that the mean (covariance) of the result has the symmetry described by Eq. 1 (Eq. 2).
Define ρ(·) by We will directly compute the first and second order statistics of ρ(·) and find that they satisfy Eq. 1 and Eq. 2, respectively, as is desired. For the first order statistics,ρ(x) . = E[ρ(x)] = Ng β=1s (R −1 β x) which is in the form of Eq. 37 so thatρ(·) has the symmetry described by Eq. 1. For the second order statistics, Independent of the properties of C = and C = , application of the Rearrangement Theorem [40, Theorem II, p. 24] demonstrates the C(x 1 , x 2 ) satisfies the symmetry of Eq. 2. For Gaussian stochastic processes, mean and covariance satisfying the constraints of Section IV or VII provides a simulation method based on the basis functions with special symmetry properties but for non-Gaussian processes, problems with many basis functions having symmetry properties, or problems with basis functions that lack symmetry properties, Eq. 38 provides an alternative method.

IX. CONCLUSION
In this paper we formulate and solve 3-D modeling and reconstruction problems for stochastic signals where the statistics of the stochastic signal are required to have a geometric symmetry. The work is motivated by cryo electron microscopy, which is of great importance in structural biology and of growing importance in the material science of nanoscale particles.
Previous algorithms assumed that the realizations of the stochastic signal have the geometric symmetry. While this implies that the statistics have the symmetry, it is a sufficient condition rather than a necessary condition and it has lead to inaccurate results for the second order statistics of the electron scattering intensity on and near symmetry axes [31, p. 173]. Since these are often locations of important biology, inaccurate results greatly reduce the value of the entire experiment.
Beyond correcting the problem in existing computations, this is the first publication of which the authors are aware in which a full 6-D covariance function is computed. Interpreting the 6-D information is challenging and we are working with our biological collaborators on this problem.


Energy Detection Based Spectrum Sensing Technique for Software Defined Radio

A lot of hardware components are used to make a conventional radio communication system, the complexity of these hardware can be reduced with help of Software Defined Radio (SDR) which is a form of wireless communication. Software Defined Radio implemented system allows us to make alteration in it or operate it with help of a software, thus making the system flexible and reconfigurable. The system changes from single mode or single band to multi-mode, multi-band as we can change the functionality and parameters of the system with the help of a software. The hardware components of a traditional radio communication system such as converters , modulator, demodulator, detector, filter, etc, make the platform cost very high and increases the complexity of a system. Software Defined Radio helps in excluding all the hardware from the system as it replaces them by pure software. In SDR modulation plays a very vital role as different modulation techniques are very important part of wireless technology and radio communication. This paper discusses about the amplitude modulation in brief and also discusses the optimized utilization of spectrum resources with help of software defined radio.

Introduction
Amplitude Modulation transmits information from one system to another through radio carrier wave and is widely used in electronic communication. In amplitude modulation, in line with variations in intensity of sound wave, the amplitude or signal strength of carrier wave is varied. In this way overall amplitude or envelope of carrier is modulated. A Canadian engineer named Reginald Fesseden was the first person who recorded AM in 1900. In the process he placed a carbon microphone in the antenna lead by which amplitude modulation was done through continuous spark transmission. Nowadays, wireless applications require spectrum resources. Each new service is provided with its own frequency block. This pre-decided spectrum allotment to users lead to wastage of spectrum, as huge portion of spectrum is not productively used; some of it is used to a limited extent, while some of it is left completely unused. The unlicensed users are not allowed to access the licensed spectrum as it is restricted by current government policies; as a result, they use interferenceprone and heavily populated frequency bands, which may be already in use. Thus in those bands there is huge spectrum scarcity problem. The major government agency such as Telecom Regulation Authority of India (TRAI) is responsible for maintaining the radio spectrum as it has been a limited resource. It is observed that some of the frequency bands in the spectrum are unutilized for some time, most of the frequency bands are partially occupied and some of the bands are heavily used. This observation is done by scanning the radio spectrum in which the revenue rich urban areas are included.
This paper discusses about the project work in which a Software Defined Radio (SDR) is developed which will help in improving the utilization of spectrum by allowing secondary user to share the spectrum with the primary users or if the radio spectrum is not in use by them then to borrow it from primary users. Software Defined radio is acquainted with the radio frequency environment, the communication parameters (such as modulation type, carrier frequency, transmission power and bandwidth) are chosen appropriately to make the optimized utilization of the spectrum and adapt transmission and reception accordingly. The SDR need to fill in the available spectrum holes by allotting these gaps to new upcoming users and hence can avoid the harmful interference to licensed user. For this, the Software Defined radio senses the spectrum continuously in order to detect the re-appearance of the primary user and once the primary user is detected, SDR should pull out itself from the spectrum gap it is using; this will help in avoiding the interference. While communication there may be variable propagation environment and also secondary users may generate interference. Thus accomplishment of this task becomes very difficult because of such environment conditions and also different primary users will be employed with no. of parameters such as different modulation techniques, varying data rates and transmission powers. Communication system is a broad area and is collection of many different processes including modulation, conversion, amplification etc. Modulation is a very important and basic process of communication, a type of modulation such as amplitude modulation has a great significance in communication system, hence this technique is also performed with help of software defined radio. 1) In the paper -SDR -Implementation of low frequency trans-receiver on FPGA‖ by Sunil G Yankanchiand Mrs. Kavitha .D (2015), it is discussed that FPGA is the best suited platform to achieve reprogrammed radio components. For complex design systems of the radio, SDR (Software Defined Radio) is the basic element to be considered as it can be programmed again and again to perform many functions making it configurable. In this Licensed Under Creative Commons Attribution CC BY project they have done the design implementation of low frequency transceiver using vitrex5 device. System level models and designs are carried out using system generator. 2) In paper -Investigating Latency in GNU Software Radio with USRP Embedded Series SDR Platform‖ by Nguyen B.Truong and Chansu Yu (2013), it is studied that besides the advantages like re-configurability and flexibility, a SDR system has to face a critical challenge on real-time processing due to high jitter and latency. Thus, understanding latency in SDR systems is very important. The paper focuses on investigating the latency of a popular SDR platform -GNU Software Radio with a new Universal Software Radio Peripherals (USRP) device -USRP Embedded Series.

Literature Review
3) In paper -SDR Implementation for DCF77‖ by Filip Zaplata and Miroslav Kasal (2013), The DCF77 signal is used to carry out the information of precise time for long distances. Standard DCF77 receivers consist of AM long wave receiver circuit, connected to a data processing unit. The idea is to remove the analog demodulator and process the amplified RF signal also on the processor. Specialized AM radio chip then would be replaced by a simple amplifier, e.g. operational amplifier circuit. 4) In the paper -Spectrum Sensing Techniques in cognitive

Concept of Software Defined Radio (SDR)
In this project a Software Defined radio is developed as a form of wireless communication where it will perform the amplitude modulation and a transceiver is made to detect which of the channels are being used and which are not in use for communication and while avoiding the channels which are occupied, access through the unused channels.

Spectrum Hole
Spectrum hole (SH) is the un-utilized slot of spectrum or that spectrum band which can be accessed and utilized by unlicensed users. The spectrum holes are detected by sensing whether a primary signal is present or absent in particular band and if the primary user is absent then try to access them. The primary users utilize the spectrum band either at different time slots or in different geographic regions.
Spectrum holes can be contemplated as multidimensional regions within frequency, time and space and because of them we can make non-interfering or safe use of spectrum as these vacant bands represent absence of primary users. The main task for software defined radio systems is to be able to strongly sense when they are within such a spectrum hole and once found , they can make use of it. The concept of Spectrum Hole is depicted in Figure (1).

Power Spectral Density
The energy of the frequencies, whether the variations (energy) are strong or weak is depicted by the Power Spectral Density of the signal. In other words PSD characterises the frequency content of a signal which is estimated by Fast Fourier Transform (FFT), it also represents the distribution of power of a signal over the different frequencies. The integral of power spectral density or the average power is fed as input to energy detector which in output would determine the total energy of the signal.

Energy detection
Since long time, people are recommending methods to identify spectrum detection as well as for presence of signal transmission. This is done with FFT (Fast Fourier Transform) which helps in indicating the presence of licensed signal. The technique which uses this method is energy detection technique. This basically determines the strength of Power Spectral Density (PSD, which is defined as average power of the signal) in each frequency of the signal.  The strength of PSD determines the presence or absence of primary user. Primary user is present if the average power is high and if it is low, it indicates absence of primary user. This phenomenon is depicted with the hypothesis. H0 is a null hypothesis, which determines that there is no licensed user signal i.e. primary user is absent and H1 states that signal is transmitted, which means primary user is present . A pre-decided threshold value is compared with the final output of the energy detector and on the basis of that we get aware of the signal status. The output signal V from the integrator is Comparing this output signal to threshold value we get the presence or absence of signal.

Project Work
For the analysis, a message frequency, sampling frequency and 5 Carrier Frequency Bands for users are initialized first . Different values are given to different frequencies .

In the first step
AM modulates the data in the first step over the respective frequency bands . After that , all the modulated signals are added in adder to produce a transmitting signal.

Results
After performing Amplitude Modulation , the project leads to second step i.e. calculation of Power Spectral Density (PSD). The results are compared with predefined threshold value. This would provide us with the status of presence of primary user signal. The simulation is done for five signals. These signals will be having sampling frequency of 12KHz and carrier frequencies will be 1KHz, 2KHz, 3KHz, 4KHz and 5KHz respectively.
After calculating absence or presence of primary user, in this analysis, it has been taken that primary users at 3 slots i.e. second , third and fourth are absent. And at two slots i.e. 1 st and 5 th the primary users are present.

Figure 5: Adder Output
The graph gives the adder output. Green arrow bands indicate the occupied or used bands (1st and 5th), red ones are the unused bands (2nd, 3rd, and 4th). The amplitude or the peak value is between (15 to 30) dB/Hz for the occupied spectral gaps or where the primary user is already present.After searching for the first available spectral gap (Spectrum Hole), the Software Defined Radio (SDR) system will automatically assign that hole to the secondary user (SU) in the spectrum. It is shown in the figure below that the first available spectral gap was occupied by the secondary user (SU) 1.

Figure 6:
Secondary User 1(SU1) occupies 1st spectral gap A greater dip in the peak indicated by the green arrow shows that secondary user has occupied that place.
As the first available spectrum hole was assigned to new user by system, in the same way the next available spectrum gaps would also be assigned one by one automatically to the upcoming users in spectrum. As soon as all the slots are allotted to the users, system will not entertain any other user and will be able to free up the spectral gap (slot) one by one. If the primary user reappears and the system needs an empty slot, it will eradicate the data in the first spectral gap by asking the secondary user to leave the hole and make it ready for the next upcoming assignment. Graph after all the spectrum holes are occupied one by one is as shown.  It can be observed through the graphs obtained that as the SNR increases, the disturbance in the spectrum decreases. This makes it obvious that the probability of the error in the received signal is increased due to noise. If we take the attenuation percentage values as 20% and 40% we get the following results. It is observed that as the attenuation increases, the signal peaks reduce proportionately. This reduces the signal power which results in weakening of the proper signal reception.

Conclusion
The idea in this project was to perform the amplitude modulation using SDR and take decisions on the basis of power spectral density (PSD) of the channel which can be used to search the available spectral gaps/holes .These available spectrum holes can be provided to new incoming users i.e. Secondary Users (SU) which will thus improve the entire channel's throughput as now the spectrum band is properly utilized. In this work the method of energy detection for spectrum sensing is done using FFT within the specified frequency band with help of MATLAB. It is successfully demonstrated in simulation results that Software Defined Radio can work in a very dynamic manner with changing frequency band from one to another. The values for the Signal to Noise ratio (SNR) are taken as 10dB, 15dB and Attenuation percentages 20% and 40% has been used to show the simulation results.


extfieldsofstudy:
Computer Science
Engineering

text:
A Whole-Body Software Abstraction layer for Control Design of free-floating Mechanical Systems

In this paper, we propose a software abstraction layer to simplify the design and synthesis of whole-body controllers without making any preliminary assumptions on the control law to be implemented. The main advantage of the proposed library is the decoupling of the control software from implementation details, which are related to the robotic platform. Furthermore, the resulting code is more clean and concise than ad-hoc code, as it focuses only on the implementation of the control law. In addition, we present a reference implementation of the abstraction layer together with a Simulink interface to provide support to Model-Driven based development. We also show the implementation of a simple proportional-derivative plus gravity compensation control together with a more complex momentum-based bipedal balance controller.

I. INTRODUCTION
Nowadays, robotics is moving from the original industrial context to more human-like environments. Foreseen applications involve robots with augmented autonomy and physical mobility. Within this novel context, physical interaction influences stability and balance. Consequently, the requirements and tasks that we expect from some platforms are changing as well. Instead of precise positioning tasks confined in cages in industrial assemblies, robots are foreseen to help in everydaylife tasks such as cleaning houses or elderly assistance.
The increase in complexity of robotic systems demands an increase in complexity of the corresponding control software. While ad-hoc solutions can be easy to implement, it is important to consider scalability, flexibility and portability of the developed software. The possibility to use the same software to control more than one platform can be of enormous importance in simplifying the testing, tuning, and deployment of the same controller on different robots.
Whole-body control has received an increased attention by the robotics community because of the possibility it offers to accomplish tasks coordination and to fully take advantage of the robots dynamics in presence of contacts. Indeed, the possibility to specify multiple objectives, even conflicting, at the same time opens the possibility to properly exploit robots for complex scenarios. In particular, citing the definition from the RAS Technical Committee [1] "Whole-Body Control aims to i) define a small set of simple, low-dimensional rules (e.g., equilibrium, self collision avoidance, etc.) ii) that are sufficient to guarantee the correct execution of any single task, whenever feasible [...], and of simultaneous multiple tasks [...], iii) exploiting the full capabilities of the entire body of redundant, floating-based robots in compliant multi-contact interaction with the environment".
In the context of whole-body control the Task function approach [2] has been successfully used. In this method, the control objectives are represented as n-dimensional continuous output functions, called tasks, to be regulated to zero. All the tasks, together with possible constraints are then transformed into a constrained optimization problem. From a software perspective, different implementations exist nowadays, among which the Stack of Task (SoT) [3], OpenSoT [4], Con-trolIt! [5] and the Instantaneous Task Specification using Constraints (iTaSC) [6]. The above softwares allow the user to specify the objectives and constraints but they solve the control problem internally. A disadvantage is that they force the user to choose a specific task-based approach to obtain the control solution thus denying the control designer the possibility to synthesize different control laws.
In this paper we propose a different approach for the wholebody control of mechanical systems. We deal with the control problem from a more general perspective, without limiting the user to the use of a task-based approach. Indeed, when we consider a generic control system, we usually identify three main building blocks: • Plant model. If we consider a model-based control system, in this block the information about the plant model given the current plant state are computed. • Feedback from the plant. This usually implies the possibility to obtain the current state of the controlled plant. • Actuation. The control system must interact with the plant. Any control-oriented software library must provide the above features to be of any use. Given the complexity of robotic systems it can be difficult, time consuming and error prone to write the controller directly in a low-level programming language such as C++. Nevertheless the control library must be efficient as it is usually required to have fast control loops. The aforementioned requirements serve as motivation for a model-based driven approach in such control libraries.
In this paper we propose a software abstraction layer which is responsible of decoupling the control software from i) the actual interface used to obtain the state feedback; ii) the actual interface used to command the actuation; iii) the dy-namic software library used to represent the robot dynamical model. Furthermore the proposed library is scalable and easily portable to other robots or different configurations. This paper is structured as follows. Section II introduces the mathematical formulation of the dynamics of mechanical systems and it shows an example of a simple classic controller. Section III describes the architecture of the proposed wholebody abstraction library and its key elements. A specific implementation is instead presented in Section IV. The controller mathematically introduced in Section II is implemented with the proposed library in Section V. Finally Section VI draws the conclusions.

II. DYNAMICS OF A MECHANICAL SYSTEM
This section introduces the mathematical formulation commonly used in the robotics literature to describe the dynamics of mechanical systems, such as robots. Because a precise formulation of the mathematical problem is out of the scope of the present paper, we refer the interested reader to books on dynamics of mechanical systems [7], [8], [9] and control systems [10], [11] for further readings.

A. Notation
Throughout the section we will use the following definitions: • I denotes an inertial frame, with its z axis pointing against the gravity. • 1 n ∈ R n×n is the identity matrix of size n; 0 m×n ∈ R m×n is the zero matrix of size m × n and 0 n = 0 n×1 . • Given two orientation frames A and B, and vectors of coordinates expressed in these orientation frames, i.e. A p and B p, respectively, the rotation matrix A R B is such that • We denote with S(x) ∈ R 3×3 the skew-symmetric matrix such that S(x)y = x × y, where × denotes the cross product operator in R 3 .

B. System modelling
We assume that the mechanical model is composed of n+1 rigid bodies -called links -connected by n joints with one degree of freedom each. In addition, we also assume that the multi-body system is free floating, i.e. none of the links has an a priori constant pose with respect to the inertial frame. This implies that the multi-body system possesses n + 6 degrees of freedom. The configuration space of the multi-body system can then be characterized by the position and the orientation of a frame attached to a robot's link -called base frame B -and the joint configurations. More precisely, the robot configuration can be represented by the triplet where ( I p B , I R B ) denotes the origin and orientation of the base frame expressed in the inertial frame, and q j denotes the joint angles.
The velocity of the multi-body system can then be characterized by the triplet where I ω B is the angular velocity of the base frame expressed w.r.t. the inertial frame, i.e. IṘ B = S( I ω B ) I R B . We also assume that the robot is interacting with the environment through n c distinct contacts. The application of the Euler-Poincaré formalism [12,Ch. 13.5] to the multi-body system yields the following equations of motion: where M ∈ R n+6×n+6 is the mass matrix, C ∈ R n+6×n+6 is the Coriolis matrix and G ∈ R n+6 is the gravity term. τ are the internal actuation torques and B is a selector matrix which depends on the available actuation, e.g. in case all joints are actuated it is equal to respectively the force and corresponding moment of the force, denotes an external wrench applied by the environment on the link of the k-th contact. The Jacobian is the map between the robot velocity ν and the linear and angular velocity of the frame C k , i.e.

C. Control Example
To illustrate the use of the dynamical model presented in Section II-B we present the classic Proportional Derivative (PD) plus Gravity compensation controller as example.
This kind of controller has been usually applied to fullyactuated fixed-base robots. Considering the model presented in Section II-B this means that the base frame position and orientation are constant and known a-priori and thus they are not part of the robot state.
The control objective is the asymptotical stabilization of a desired constant joint configuration q d j or equivalently the asymptotical stabilization to zero of the error The choice of the following control action where K p , K d ∈ R n×n are the positive definite proportional and derivative gain matrices and G j (q) = [0 n×6 1 n ] G(q), satisfy the control objective, i.e. the stabilization to zero ofq j , and it can be proved by Lyapunov arguments [7, Sec. 6.5.1].

III. SOFTWARE ARCHITECTURE
We propose a software abstraction layer to simplify creating whole-body controllers for highly redundant mechanical systems. Given the requirements introduced in the previous sections we highlight four main elements that must be present in the library, i.e. Actuators, Sensors, State and Model. The abstraction offered by the library allows one to also easily implement higher-level interfaces such as a Simulink R interface. Figure 1 summarizes the whole software architecture.
Note that the paradigm described by this library does not assume any particular robot operating system or underlining software as this is left to the actual library implementation.
A crucial feature of the proposed abstraction layer is related to the ordering of the information provided from and to the robot. In fact, the elements that directly interface the hardware, i.e. the Actuators, Sensors and State have to represent the information in a robot-dependent suitable way. On the other hand, the Model element usually interfaces with libraries that represent the information with the formalism of Eq. (1). To further complicate the problem, the control software may want to access only a subset of the degrees of freedom modelled by the dynamics library, or provided by the robot. The wholebody abstraction library must thus orchestrate all the various elements to provide a unified interface to the control software.
We now describe in detail the role that each element has in the proposed library.

A. Actuators
The actuators element abstracts the actual control of the robot motors. In particular it exposes the possible motors controllable mode, e.g. position control, velocity control and torque control just to cite the most common. Of course, it also provides the possibility to specify the references for the low level controllers.

B. Sensors
The sensors element is the counterpart of the actuators element. In fact, it abstracts all the sensors available on the robot, usually the readings from encoders, force/torque sensors or accelerometers and it is responsible for providing access to the latest sensor measurements.

C. State
The state element represents all the possible information which can be measured or estimated on the robot. This implies that state encompasses the information provided by the sensor element. Furthermore, it provides additional information which can come from estimation or filtering of the data. For example, if the robot provides only joint position measurements, e.g. coming from the joint encoders, a first and second derivative filter can provide velocity and acceleration measurements. In case this information is provided by the robot itself, no additional processing is required from the interface. It is important to notice that in both cases the control software using the abstraction library will remain exactly the same.

D. Model
The last element is the model element. It abstracts the kinematic-and dynamic-related information that a controller needs while computing the control law. In general data are represented with the formalism of Eq. (1). Note that a common requirement for a control library is to control only a subset of the degrees of freedom of a robot, e.g. control only the lower body of a legged robot while walking. For this reason, the library must correctly compute the kinematics and dynamics of the whole system, while considering the possibility to expose only a subset of the quantities as requested by the control software.

IV. IMPLEMENTATION
This section describes the current implementation of the whole-body abstraction library conceptually described in Section III. The code has been implemented in C++ because of its diffusion and computational performance while remaining a high-level programming language. The implementation has been divided in two libraries: the wholeBodyInterface [13] and yarpWholeBodyInterface libraries [14].

A. wholeBodyInterface
The wholeBodyInterface is the direct transposition in C++ of the abstract concepts described in Section III. The four elements, i.e. actuators, sensors, state and model, are represented as abstract classes. Additionally, the library provides utilities to identify the various degrees of freedom.
As it represents the coded counterpart of the abstraction library, wholeBodyInterface does not make any assumption on the underlining robot framework, or how data is organized, using only native C++ types.

B. yarpWholeBodyInterface
The yarpWholeBodyInterface is the actual implementation of wholeBodyInterface specifically considering YARPpowered mechanical systems [15]. Regarding the model implementation we choose as kinematic and dynamic library the iDynTree library [16] and information about the kinematic and dynamic model can be loaded from different sources, e.g. a URDF representation.
The actuators and sensors elements directly interact with YARP control boards. Because a robot possesses in general multiple control boards, these two elements are also responsible for mapping the information coming from the control boards to the degrees of freedom selected by the library user.
Note that, because of the dependency on the YARP library, in the current implementation the state element uses YARP data structures, e.g. vectors and matrices, but this dependency can be easily dropped in future implementations.

C. Simulink Interface for Model-Driven Engineering
C++ applications can leverage the advantages of the proposed abstraction layer while keeping full control of the performance of the control software by directly using the provided C++ implementations. On the other hand, coding and testing a complex control system directly in C++ can be prohibitive. For example, even the simple task of monitoring a signal over time can be complex and requires the use of a dedicated library. The use of software to design and simulate dynamical system models greatly helps the design and synthesis of control systems. Domain-specific software for dynamical systems is a specific case of model-driven engineering [17].
We currently implemented the Simulink interface to our proposed whole-body abstraction library, which can be found in [18]. Most of the features accessible in C++ are also accessible to Simulink models. Furthermore, because the connection with the robot or the simulator is handled by the underlining C++ library, the Simulink interface does not require any particular toolbox to command the robot, e.g. Simulink Real Time R .
A further advantage of using Simulink R with respect to the C++ code consists in the possibility to exploit the abundance of toolboxes and Matlab native functions out of the box.

V. EXPERIMENTS
This section presents the implementation of the PD plus gravity compensation controllers briefly described in Section II-C. We also discuss the results of a more complex controller, namely a momentum-based balancing controller which has been implemented with the Simulink interface described in Section IV-C.

A. PD plus Gravity Compensation
This section reports the code for the example presented in Section II-C, i.e. the code for the PD plus gravity compensation controller.
Because it is a simple example we show both the C++ code (see Code 1 and 2) and the Simulink model diagram (see Figure 2). Note that, while the Simulink diagram completely represents the controller, the C++ code snippet has been extracted from the main loop function, i.e. the function which runs at every iteration. How the control thread is created and managed depends on the particular system and it is outside the scope of the present paper.
The snippet of code in Code 1 shows how the specific YARP-based implementation is instantiated. In particular, the current implementation needs information about the URDF model representing the kinematic and dynamic information of the robot and the mapping between the model joints and the YARP control boards. This is provided by the object created at line 4 and passed to the interface constructor at line 7. Additionally, the list of controlled joints are passed to the interface at line 19, just before the interface initialization routine is called.
Reading the code in Code 2 it is possible to observe how all the details regarding the specific robot platform are hidden by the library. The object robot, in fact, is accessed through its abstract type, as it can be also seen during its instantiation, i.e. in line 7 of Code 1. In lines 4 − 7 the state of the robot, i.e. (q j ,q j ), is read. The feedforward term, corresponding to G(q) is computed at lines 10 − 14 where the last parameter is the resulting gravity compensation term. Finally the error and the feedback term necessary to implement Eq. (2) is computed in lines 17 − 22. Because we did not use any specific mathematical library we explicitly computed the term K pqj + K dqj in the for loop. Finally, at line 25 we send the torque command to the robot, which we previously setup to be controlled in torque mode. robot->setControlReference(torques); Figure 2 shows the same code implemented directly in Simulink. It is evident how the block-based diagram is clearer with respect to its C++ counterpart. Furthermore, the possibility to add scopes, or dump signal variables directly into Matlab workspace greatly increases its advantages with respect to directly coding in C++.

B. Momentum-based Balance Control
To show the power of the proposed architecture we present here a second example, i.e. we show the results of a momentum-based balancing controller which has been synthesized directly by using the Simulink interface. Given the complexity of the control problem we do not report here screenshots or code snippets of the Simulink model, but the model can be examined in [19], while the mathematical formulation can be found in [20].
The YouTube c video [21] shows the robot performing complex movements by using the controller implemented and running as a Simulink Model. By using the yarpWholeBody-Interface implementation we also leverage the capabilities of the YARP middleware to seamlessly connect to the real or simulated system. In particular the test platform is the iCub humanoid robot [22], endowed with 53 degrees of freedom, 6-axis force/torque sensors and distributed tactile skin. The robot is simulated on the Gazebo simulator [23] by means of Gazebo-YARP plugins [24]. The same demo has also been implemented on a different configuration of the iCub platform [25]. Note that the two robots have a different set of degrees of freedom. Thanks to the flexibility of the library, the controller code remains the same in both scenarios.
We encourage the interested reader to test the controller on the Gazebo Simulator. Instructions on how to run the controller can be found directly in the model repository readme [19].

VI. CONCLUSIONS
In this paper we presented a software abstraction layer to simplify the development of whole-body controllers. While there are already some whole-body control software libraries, they already define the controller structure and leave to the user only the possibility to specify objectives and constraints.
On the other hand the proposed library leaves complete freedom to the control designer by exposing all the information needed. It does not make any assumptions on the controller structure. The whole-body abstraction library presents also the following advantages: • it decouples the writing of the controller from a particular robot implementation • it decouples the writing of the controller from a specific dynamic library implementation • it allows more concise and clear code as it represents uniquely the code needed to implement the mathematical formulation of the controller. All the implementation details are left to the library • it allows to benchmark the controller on different platforms or with different implementations. Furthermore, the possibility to expose the functionality at an higher level than C++ facilitates the writing of controllers as the results on the iCub robot clearly prove.
We voluntarily did not consider some aspects as they are out of the scope of the present contribution. Nevertheless they must be taken into account when a controller is implemented and used on the real system. In particular the following details should be considered: • how are controllers run on the platform? Do they run as threads? • how are controllers configured and initialized? • how is communication with other software performed?
For example, how are desired values provided to the controller, coming from a planner or higher-level control loop?
By not considering these details in the abstraction library, we render the library portable to different systems. Indeed, the actual control law is not concerned by the previously listed implementation details.
While the more complex demos have been achieved by directly executing the Simulink model connected to the robot, we recognize the need to automatically generate self-contained C++ code. The advantage is twofold. On one side the autogenerated code is in general more optimized than the code directly executed in Simulink, even if less optimized than ad-hoc C++ code. On the other side, this would remove the requirement of having a Simulink installation on the computers controlling the robot.


Towards Reduced Instruction Sets for Synchronization

Contrary to common belief, a recent work by Ellen, Gelashvili, Shavit, and Zhu has shown that computability does not require multicore architectures to support"strong"synchronization instructions like compare-and-swap, as opposed to combinations of"weaker"instructions like decrement and multiply. However, this is the status quo, and in turn, most efficient concurrent data-structures heavily rely on compare-and-swap (e.g. for swinging pointers and in general, conflict resolution). We show that this need not be the case, by designing and implementing a concurrent linearizable Log data-structure (also known as a History object), supporting two operations: append(item), which appends the item to the log, and get-log(), which returns the appended items so far, in order. Readers are wait-free and writers are lock-free, and this data-structure can be used in a lock-free universal construction to implement any concurrent object with a given sequential specification. Our implementation uses atomic read, xor, decrement, and fetch-and-increment instructions supported on X86 architectures, and provides similar performance to a compare-and-swap-based solution on today's hardware. This raises a fundamental question about minimal set of synchronization instructions that the architectures have to support.

Introduction
In order to develop efficient concurrent algorithms and data-structures in multiprocessor systems, processes that take steps asynchronously need to coordinate their actions. In shared memory systems, this is accomplished by applying hardware-supported low-level atomic instructions to memory locations. An atomic instruction takes effect as a single indivisible step. The most natural and universally supported instructions are read and write, as these are useful even in uniprocessors to store and load data from memory.
A concurrent implementation is wait-free, if any process that takes infinitely many steps completes infinitely many operation invocations. An implementation is lock-free if in any infinite execution infinitely many operations are completed. The celebrated FLP impossibility result [FLP85] implies that in a system equipped with only read and write instructions, there is no deterministic algorithm to solve binary lock-free/wait-free consensus among n ≥ 2 processes. Binary consensus is a synchronization task where processes start with input bits, and must agree on an output bit that was an input to one of the processes. For one-shot tasks like consensus, wait-freedom and lock-freedom are equivalent.
Herlihy's Consensus Hierarchy [Her91] takes the FLP result further. It assigns a consensus number to each object, namely, the number of processes for which there is a wait-free binary consensus algorithm using only instances of this object and read -write registers. An object with a higher consensus number is hence a more powerful tool for synchronization. Moreover, Herlihy showed that consensus is a fundamental synchronization task, by developing a universal construction which allows n processes to wait-free implement any object with a sequential specification, provided that they can solve consensus among themselves.
Herlihy's hierarchy is simple, elegant and, for many years, has been our best explanation of synchronization power. It provides an intuitive explanation as to why, for instance, the compare-and-swap instuction can be viewed "stronger" than fetch-and-increment, as the consensus number of a Compare-and-Swap object is n, while the consensus number of Fetch-and-Increment is 2.
However, key to this hierarchy is treating synchronization instructions as distinct objects, an approach that is far from the real-world, where multiprocessors do let processes apply supported atomic instructions to arbitrary memory locations. In fact, a recent work by Ellen et al. [EGSZ16] has shown that a combination of instructions like decrement and multiply-by-n, whose corresponding objects have consensus number 1 in Herlihy's hierarchy, when applied to the same memory location, allows solving wait-free consensus for n processes. Thus, in terms of computability, a combination of instructions traditionally viewed as "weak" can be as powerful as a compare-and-swap instruction, for instance.
The practical question is whether we can really replace a compare-and-swap instruction in concurrent algorithms and data-structures with a combination of weaker instructions. This might seem improbable for two reasons. First, compare-and-swap is ubiquitous in practice and used heavily for various tasks like swinging a pointer. Second, the protocol given by Ellen et al. solves only binary n-process consensus. It is not clear how to use it for implementing complex concurrent objects, as utilizing Herlihy's universal construction is not a practical solution. On the optimistic side, there exists a concurrent queue implementation based on fetch-and-add that outperforms compare-and-swap-based alternatives [MA13]. Both a Queue and a Fetch-and-Add object have consensus number 2, and this construction does not "circumvent" Herlihy's hierarchy by applying different nontrivial synchronization instructions to the same location. Indeed, we are not aware of any practical construction that relies on applying different instructions to the same location.
As a proof of concept, we develop a lock-free universal construction using only read , xor , decrement, and fetch-and-increment instructions. The construction could be made waitfree by standard helping techniques. In particular, we implement a Log object [BMW + 13] (also known as a History object [Dav04]), which supports high-level operations get-log() and append (item), and is linearizable [HW90] to the sequential specification that get-log() returns all previously appended items in order. This interface can be used to agree on a simulated object state, and thus, provides the universal construction [Her91]. In practice, we require a get-log() for each thread to return a suffix of items after the last get-log() by this thread. We design a lock-free Log with wait-free readers, which performs as well as a compare-and-swap-based solution on modern hardware.
In our construction, we could replace both fetch-and-increment and decrement with the atomic fetch-and-add instruction, reducing the instruction set size even further.

Algorithm
We work in the bounded concurrency model where at most n processes will ever access the Log implementation. The object is implemented by a single fetch-and-increment-based counter C, and an array A of b-bit integers on which the hardware supports atomic xor and decrement instructions. We assume that A is unbounded. Otherwise, processes can use A to agree on the next array A to continue the construction. C and the elements of A are initialized by 0. We call an array location invalid if it contains a negative value, i.e., if its most significant bit is 1, empty if it contains value 0, and valid otherwise. The least significant m = log 2 (n + 1) bits are contention bits and have a special importance to the algorithm. The remaining b − m − 1 bits are used to store items. See Figure 1 for illustration. For every array location, at most one process will ever attempt to record a (b − m − 1)-bit item, and at most n − 1 processes will attempt to invalidate this location. No process will try to record to or invalidate the same location twice. In order to record item x, a process invokes xor (x ), where x is x shifted by m bits to the left, plus 2 m − 1 ≥ n, i.e., the contention bits set to 1. To invalidate a location, a process calls a decrement. The following properties hold: 1. After a xor or decrement is performed on a location, no read on it ever returns 0.
2. If a xor is performed first, no later read returns an invalid value. Ignoring the most significant bit, the next most significant b − m − 1 bits contain the item recorded by xor .
3. If a decrement is performed first, then all values returned by later read s are invalid.
A xor instruction fails to record an item if it is performed after a decrement.
To implement a get-log() operation, process p starts at index i = 0, and keeps reading the values of A[i] and incrementing i until it encounters an empty location A[i] = 0. By the above properties, from every valid location A[j], it can extract the item x j recorded by a xor , and it returns an ordered list of all such items (x i 1 , x i 2 , . . . , x i k ). In practice, we require p to return only a suffix of items appended after the last get-log() invocation by p. This can be accomplished by keeping i in static memory instead of initializing it to 0 in every invocation. To make get-log wait-free, p first performs l = C.read (). Then, if i becomes equal to l during the traversal, it stops and returns the items extracted so far.
To implement append (x), process p starts by = C.fetch-and-increment(). Then it attempts to record item x in A[ ] using an atomic xor instruction. If it fails to record an item, the process does another fetch-and-increment and attempts xor at that location, and so on, until it is able to successfully record x. Suppose this location is A[ ]. Then p iterates from j = − 1 down to j = 0, reading each A[j], and if A[j] is empty, performing a decrement on it. Afterwards, process p can safely return.
fetch-and-increment guarantees that each location is xor ed at most once, and it can be decremented at most n − 1 times, once by each process that did not xor . As a practical optimization, each process can store the maximum from its previous append operations and only iterate down to in the next invocation (all locations with lower indices will be non-empty). Our implementation of append is lock-free, because if an operation takes steps and does not terminate it must be repeatedly failing to record items in locations. This only happens if other xor operations successfully record their items and invalidate these locations.
At any time t during the execution, let us denote by f (t) as the maximum index such that, A[f (t)] is valid and A[j] is non-empty for all j ≤ f (t). By the first property f (t) is non-decreasing, i.e., for t > t we have f (t ) ≥ f (t). We linearize an append Consider a get-log() operation with the latest returned item x extracted from A[ ]. We show by contradiction that the execution interval of this get-log() must contain time t such that f (t) = . We then linearize get-log() at the smallest such t. It is an easy exercise to deal with the case when multiple operations are linearized at exactly the same point by slightly perturbing linerization points to enforce the correct ordering. Suppose the get-log() operation extracts x from A[ ] at time T . f (T ) ≥ as get-log() stops at an empty index, and by the contradiction assumption we must have = f (T ) > . get-log() then reaches valid location A[ ] and extracts an item x from it, contradicting the definition of x .
We implemented the algorithm on X86 processor and with 32 threads. It gave the same performance as an implementation that used compare-and-swap for recording items and invalidating locations. It turns out that in today's architectures, the cost of supporting compare-and-swap is not significantly higher than that of supporting xor or decrement. This may or may not be the case in future Processing-in-Memory architectures [PAC + 97].
Finding a compact set of synchronization instructions that, when supported, is equally powerful as the set of instructions used today is an important question to establish in future research.