text:
## help pls!! How many different numbers between 1/1000 and 1000 can be written either as a power of 2 or as a power of 3, where the exp

Question

help pls!!
How many different numbers between 1/1000 and 1000 can be written either as a power of 2 or as a power of 3, where the exponent is an integer? Remember that we’re counting all possible numbers, not just integers.

in progress 0
2 months 2021-07-31T07:36:55+00:00 1 Answers 0 views 0

1. How many different numbers between 1/1000 and 1000 can be written either as a power of 2 or as a power of 3, where the exponent is an integer?

2^9 = 512    So that is  9+9+1 = 19 powers of 2

3^6=729     so that is 6+6+1 = 13 powers of 3

Now you just have to work out how many are powers of 2 and 3 because they have  been counted twice.

If they are powers of 2 and 3 that means they are powers of 6

6^3 = 216 which is the biggest so that is  3+1+3 = 9

So I think the answer is      19+13-9 = 23

mark me as a brill please

format: plain


text:
# AIRS SRF Study

Figure Viewer

## 1 Purpose and Scope

To investigate the sensitivity of the AIRS radiances to changes in the AIRS SPectral Response Functions (SRF). To answer the question whether in the long wave region (AIRS detector modules 11 and 12) the knowledge of the SRFs is in error and whether they can be improved using independent measurements and intercomparison with CrIS and or IASI.

It has been observed in the AIRS:CrIS SNO mean bias, that there is a repeatable and significant variation in the region from about 645 to 750 cm-1, that is for detectors of modules 11 and 12. This variation is attributed to the AIRS observed radiances. But to-date this variation has not been accounted for nor has it been determined whether it can be corrected.

The scope of this technical note is to record the methods and results.

## 2 Methods

Working directory: /home/chepplew/projects/airs/srf_study/
Main routine therein: airs_srf_analysis.m
Sources: (from strow) /home/chepplew/gitLib/srf_model/ and /home/chepplew/gitLib/l1c_freq_adjust/
Intermedite data files: /home/chepplew/projects/airs/srf_study/outd/
Figures: /home/chepplew/projects/airs/srf_study/figs/


text:
# Beamer : Overlay-aware environments using \newtcbtheorem

I recently tried \newtcbtheorem to define new environments, however, these are not overlay-aware. Any idea why this is happening ?

Here is a MWE :

\documentclass{beamer}

\usepackage[theorems]{tcolorbox}

\newtcbtheorem[number within=section]{theo_beamer}{Theorem}{colframe=blue}{thm}
\newtcbtheorem[number within=section,use counter from=theo_beamer]{prop_beamer}{Proposition}{colframe=red}{prop}

\resetcounteronoverlays{tcb@cnt@theo_beamer}

\begin{document}
\begin{frame}

\begin{theo_beamer}{Title}{label1}
main result
\end{theo_beamer}

\onslide<2->{
\begin{prop_beamer}{Title}{label2}
Let's call Theorem \ref{label1}
\end{prop_beamer}
}
\end{frame}
\end{document}

• With \ref{label1} replaced by \ref{thm:label1}, compiling your example, I get a two-page pdf. First page shows theo_beamer environment only, and second page shows two environments. Hence I cannot reproduce your problem. Jul 7, 2020 at 23:42


text:
# B Complex Neuro Factors with Alpha Lipoic Acid and CoQ10 - 60 Veg Capsules

0.0
Quantity per Container: 60 Capsules
 Item Number: PT1225 Serving Size: 2 Capsules Servings per Container: 30 Price per Serving:  \$0.80 Price per Unit*: \$0.40 per Capsule * Price Per Unit may reflect product's price per capsule, tablet, ounce, fluid ounce, etc.
Out of Stock
\$24.00

As co-factors in numerous biochemical reactions‚ B Complex vitamins are critical for energy production and the proper metabolism of carbohydrates‚ fats‚ and proteins. B Vitamins are also important for the maintenance of normal nervous system function‚ blood cell production‚ and countless other metabolic processes. B Complex Neuro Factors has a balanced combination of B Vitamins in their activated CoEnzyme states. Acid-resistant capsules help to preserve ingredient integrity through the harsh acidic environment of the stomach for intact delivery to intestinal absorption sites.

These statements have not been evaluated by the Food and Drug Administration (FDA). These products are not meant to diagnose‚ treat or cure any disease or medical condition. Please consult your doctor before starting any exercise or nutritional supplement program or before using these or any product during pregnancy or if you have a serious medical condition.

• MFG # P0406
• UPC # 707359104068
Official Website: http://www.protocolforlife.com

Serving Size: 2 Veg Capsules
Servings Per Container: 30

Amount Per Serving / % DV

Vitamin C (from Magnesium Ascorbate) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60mg / 67%
Thiamin (Vitamin B-1) (from Thiamin HCI and Cocarboxylase) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50mg / 4167%
Riboflavin (As Riboflavin & from Riboflavin-5-Phosphate Sodium) . . . . . . . . . . . . . . . . . . . . . . . . . . 50mg / 3846%
Niacin (Vitamin B-3) (from Inositol Hexanicotinate and NAD Trihydrate) . . . . . . . . . . . . . . . . . . . . . 100mg / 625%
Vitamin B-6 (from Pyridoxine HCI and P-5-P Monohydrate) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50mg / 2941%
Folate [from Quatrefolic® (6S)-5-MTHF** Glucosamine Salt] [400mcg (6S)-5-MTHF**] . . . . . 667mcg DFE / 167%
Vitamin B-12 (as Methylcobalamin and Dibencozide) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 250mcg / 10‚417%
Biotin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1‚000mcg / 3333%
Pantothenic Acid (Vitamin B-5) (from Calcium Pantothenate and Pantethine) . . . . . . . . . . . . . . . . 50mg / 1000%
Choline (from Choline Bitartrate) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .25mg / 5%
Betaine Anhydrous . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50mg / †
Alpha Lipoic Acid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50mg / †
PABA (Para-Aminobenzoic Acid) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25mg / †
Inositol (from Inositol Hexanicotinate‚ Inositol Hexaphosphate and Inositol) . . . . . . . . . . . . . . . . . . 30mg / †
CoQ10 (Coenzyme Q10) 10mg / †

†Daily Value not established.
**MTHF = Methyltetrahydrofolate.

Other Ingredients: Cellulose (capsule) (acid-resistant)‚ cellulose powder‚ stearic acid (vegetable source) and silica.

Not manufactured with wheat‚ gluten‚ soy‚ milk‚ egg‚ fish‚ shellfish or tree nut ingredients. Produced in a GMP facility that processes other ingredients containing these allergens.

Vegan / Vegetarian.

Take 2 capsules daily with food‚ or as directed by your healthcare practitioner.

If you are pregnant/nursing‚ consult your healthcare practitioner before using this product.

Do not eat freshness packet. Keep in bottle.

Notice: Natural color variation may occur in this product.

Allergen Caution: Not manufactured with wheat‚ gluten‚ soy‚ milk‚ egg‚ fish‚ shellfish or tree nut ingredients. Produced in a GMP facility that processes other ingredients containing these allergens.

• Egg Free
• Fish Free
• Gluten Free
• GMP (Good Manufacturing Practices) Certification
• Soy Free
• Vegan
• Vegetarian
• Wheat Free
• Crustacean Shellfish Free

## Ratings and Reviews

Be the first to review this item.

Your question(s) and our corresponding answer(s) will not be confidential and will be posted publically on this specific product page.


text:
0
Review Article

# Geometric Algorithms for Robot Dynamics: A Tutorial ReviewPUBLIC ACCESS

[+] Author and Article Information
Frank C. Park

Department of Mechanical and
Aerospace Engineering,
Seoul National University,
Seoul 08826, Korea
e-mail: fcp@snu.ac.kr

Beobkyoon Kim, Cheongjae Jang, Jisoo Hong

Department of Mechanical and
Aerospace Engineering,
Seoul National University,
Seoul 08826, Korea

Manuscript received April 25, 2017; final manuscript received December 12, 2017; published online February 7, 2018. Editor: Harry Dankowicz.

Appl. Mech. Rev 70(1), 010803 (Feb 07, 2018) (18 pages) Paper No: AMR-17-1029; doi: 10.1115/1.4039078 History: Received April 25, 2017; Revised December 12, 2017

## Abstract

We provide a tutorial and review of the state-of-the-art in robot dynamics algorithms that rely on methods from differential geometry, particularly the theory of Lie groups. After reviewing the underlying Lie group structure of the rigid-body motions and the geometric formulation of the equations of motion for a single rigid body, we show how classical screw-theoretic concepts can be expressed in a reference frame-invariant way using Lie-theoretic concepts and derive recursive algorithms for the forward and inverse dynamics and their differentiation. These algorithms are extended to robots subject to closed-loop and other constraints, joints driven by variable stiffness actuators, and also to the modeling of contact between rigid bodies. We conclude with a demonstration of how the geometric formulations and algorithms can be effectively used for robot motion optimization.

<>

## Introduction

Robotics is defined by motion—controlled programmable motion—and that mechanics, particularly multibody dynamics, playing a fundamental role in robotics should come as no surprise. A typical robot is constructed by connecting rigid bodies, called links, together with joints. Actuators such as electric motors then deliver forces and torques to the joints, enabling the robot to move. Dynamic models of the robot—by this we mean the equations of motion relating the applied forces and torques to the motion of the robot—are needed for a wide range of applications, from the design of model-based control laws and trajectory generation algorithms to the simulation, design optimization, and motion planning of robots.

Because rigid-body dynamics is a classical and well-understood subject, there is a tendency to regard robot dynamics as a solved problem. Although at some level this may be true, what makes the subject of robot dynamics challenging and still an open area of research is that robots increasingly have high-dimensional, very complex structures. A typical humanoid robot, for example, can possess up to 60 degrees-of-freedom (DOFs) and intermittently form closed chains and be subject to contact and other constraints during its motion. Some of the robot's links and joints may purposefully contain elastic (“soft”) elements like variable stiffness springs and dampers, which may either be passive or actively controlled. The actuators themselves may possess complicated dynamics involving gears and various transmission devices. Model-based control laws often require real-time calculation of the dynamics. The high-dimensional and nonlinear nature of the equations of motion, together with the need to account for contact and other physical constraints, all are factors that make the formulation and computation of the robot dynamics less than straightforward.

This paper attempts to provide a review of the state-of-the-art in robot dynamics algorithms. The literature on robot dynamics is vast, and a review of the earlier foundational works can be found in a number of sources, e.g., see Ref. [1]. Our focus in this paper will be limited to those formulations and algorithms that rely on methods from differential geometry, particularly the theory of Lie groups.

This focus on geometry is motivated by the observation that often the most salient physical features of a robot are best captured by a geometric description. The advantages of the geometric approach have been recognized for quite some time by practitioners of classical screw theory. Featherstone (see Ref. [2] and the earlier references cited therein) showed how the dynamics of rigid multibody systems can be formulated by using and extending constructs from classical screw theory. The articulated body inertia algorithm by Featherstone [3] is one of the remarkable achievements of this machinery.

What has made these tools largely inaccessible to practitioners is that they require an entirely new language of notations and constructs (screws, twists, wrenches, reciprocity, transversality, conjugacy, etc.) and their often obscure rules for manipulation and transformation. On the other hand, the mostly algebraic alternatives to screw theory often mean the practitioner ends up buried in the details of calculation, losing the simple and intuitive interpretation that lies at the heart of what they are calculating.

The geometric framework described in this paper has its roots in the Lie bracket-based formulation of the attitude dynamics equation for a rotating rigid body [4] and allows one to express screw-theoretic concepts as well as Featherstone's constructs using the language of Lie groups and Lie algebras. By exploiting the Lie group structure of the rigid-body motions, it becomes possible to express the kinematics of an open chain robot as a product of matrix exponentials—this remarkable result was first derived by Brockett [5]; some of its many powerful applications are described in Ref. [6]—allowing one to reinvent screw theory simply by appealing to basic notions from linear algebra and linear differential equations.

In fact, the geometric framework allows one to go much further. The resulting algorithms, together with the formulations upon which they are based, can be expressed in a way that is invariant with respect to the choice of link reference frames (coordinate invariance). Moreover, by drawing upon more general results from Lie theory, new algorithms, such as a recursive formulation for gradients and Hessians of the dynamics, which are essential to motion optimization, can be straightforwardly obtained.

Rather than attempting a comprehensive review of geometric dynamics algorithms for the entire spectrum of robot designs, our objective will be to describe in some detail the underlying geometric concepts for a single rigid body and then to derive the basic forward and inverse dynamics algorithms for open chains. Summary descriptions and references are then provided for the corresponding formulations and algorithms of more complex robots, including those subject to closed-loop and other constraints and robots containing soft elements. We then show how these geometric algorithms can be used effectively for robot motion optimization.

We will not explore those aspects of geometric mechanics that do not relate in an explicit way to the algorithmic and computational aspects of robot dynamics; we instead refer the reader to many excellent references on geometric mechanics and control such as Refs. [4] and [79]. Papers that explore in more detail the geometric aspects of the algorithmic constructs in this paper will be pointed out throughout.

This paper is organized as follows: After reviewing the Lie group structure of the rotation and Euclidean groups and their connections to velocities and forces, we then review in detail the geometric formulation of the equations of motion for a single rigid body; these equations form the basis for the ensuing formulations and algorithms for the forward and inverse dynamics of open chain structures and also for their derivatives. We also examine from the geometric perspective the dynamics of rigid bodies in contact with each other. The corresponding algorithms for robots subject to constraints, such as closed chains, and robots containing soft elements, are then reviewed. The remainder of this paper is then devoted to robot motion optimization, which employs in a fundamental way the geometric algorithms for the dynamics and also their differentiation. We conclude with a brief discussion of some of the open research issues in robot dynamics.

## Single Rigid-Body Dynamics

We begin with a motivating example that illustrates the advantages of the geometric approach to robot dynamics. Consider a single rigid body of mass $m$ with a body-fixed frame x-y-z attached to the body's center of mass; unless stated otherwise, all vector and matrix quantities mentioned later are expressed in terms of this body frame. Let $I∈ℝ3×3$ be the 3 × 3 inertia matrix of the body, and let $f∈ℝ3$ and $m∈ℝ3$, respectively, denote the external force and moment applied to the body (here m is applied about the center of mass). The equations of motion are then given by the familiar pair of equations

$f=mam=Iα+ω×Iω$

where $v∈ℝ3$ and $a∈ℝ3$ are, respectively, the velocity and acceleration of the center of mass, and $ω∈ℝ3$ and $α∈ℝ3$ are, respectively, the angular velocity and angular acceleration.

The previously mentioned dynamic equations are expressed in terms of the body frame attached to the center of mass. Let us now see how the equations would appear if a different body frame was chosen. Affix the $x′ -y′ -z′$ frame to another point on the body; to simplify calculations, we will even assume that the orientation of this frame is identical to that of the x-y-z frame. The external force and moment expressed in the new frame are $f′$ and $m′$. The velocities and accelerations of this new frame are now represented by primes: $ω′, v′, α′$, and $a′$, all of which are vectors expressed in the new frame. Finally, let $p′$ denote the vector from the new frame origin to the center of mass, also expressed in the new frame. The equations of motion now become

$f′=−m(p′×α′)+ma′+ω′×{mv′−m(p′×ω′)}m′=Iα′−m(p′×(p′×α′))+m(p′×a)+ω×{Iω′−m(p′×(p′×ω′))+m(p′×v′)}+v′×{−m(p′×α′)+ma′}$

The derivation of these equations is quite standard and can be found in standard dynamics texts, e.g., see Ref. [10].

What is immediately apparent is that the equations become considerably more complex with respect to this new reference frame. This feature is precisely what is behind the complexity of multibody robot dynamics: the forces and moments acting on a specific link are typically expressed in different reference frames, and these must be expressed in terms of a common frame before they can be summed.

What we show later in this section is that by applying the basic geometric Lie group constructs to the special Euclidean group of rigid-body motions, the dynamics for a single rigid body can, regardless of the choice of reference frame, always be expressed in the formDisplay Formula

(1)$GV˙=adVT(GV)+F$

where $V∈ℝ6$ is a twist representing the linear and angular velocity of the rigid body and $V˙∈ℝ6$ its time derivative, $F∈ℝ6$ is a wrench representing the force and moment applied to the rigid body, $G∈ℝ6×6$ is a spatial inertia matrix, and $adV∈ℝ6×6$ is a linear twist transformation matrix and $adVT$ its transpose. There exist well-defined transformation rules for twists, wrenches, and spatial inertias under a change of reference frame. Equation (1) forms the basis for the geometrical formulations and algorithms that follow. We now review the fundamental concepts and constructs for deriving Eq. (1).

###### Matrix Lie Group Fundamentals.

This section covers the fundamentals of matrix Lie groups at a level targeted to the practitioner; for a more rigorous discussion and review of Lie groups and their properties, see, e.g., Refs. [6] and [1114]. Let GL(n), referred to as the general linear group, denote the set of real n × n nonsingular matrices. A subset G of GL(n) is a matrix Lie group if it has the structure of both a differentiable manifold (imagine a smooth multidimensional surface embedded in a higher-dimensional Euclidean space) and an algebraic group under matrix multiplication

1. (1)The identity matrix $1$ is an element of G.
2. (2)For every XG, its inverse X−1 exists and is also an element of G.

GL(n) itself is a matrix Lie group: it is a differentiable manifold of dimension n2 and also satisfies the properties of an algebraic group under matrix multiplication. The set of rotation matrices and also the set of rigid-body motions (or homogeneous transformations) are examples of matrix Lie groups that we will treat in detail later. The special linear group SL(n), consisting of real n × n matrices of unit determinant, is also a matrix Lie group.

A Lie algebra is a vector space $V$ along with a mapping called the Lie bracket, $[·,·]:V×V→V$, that satisfies for all scalars α, β and all elements A, B, C in $V$

1. (1)[αA + βB, C] = α[A, C] + β[B, C] and [C, αA + βB] = α[C, A] + β[C, B];
2. (2)[A, B] = –[B, A];
3. (3)[A, [B, C]] + [C, [A, B]] + [B, [C, A]] = 0.

As an example of a Lie algebra, the set of real n × n matrices, denoted gl(n), is a matrix Lie algebra under the Lie bracket defined by the matrix commutatorDisplay Formula

(2)$[A,B]=AB−BA$
for A, B ∈ gl(n).

Any matrix Lie group gives rise to a matrix Lie algebra, and conversely, to any matrix Lie algebra, there is a unique corresponding matrix Lie group. The matrix exponential makes explicit this connection between matrix Lie groups and Lie algebras. For any real n × n matrix A ∈ gl(n), its matrix exponential eA is defined by the seriesDisplay Formula

(3)$eA= 1+A+12!A2+13!A3+⋯$

The series is always guaranteed to converge for finite A. The matrix exponential satisfies the following properties:

1. (1)eA is always nonsingular with inverse given by eA. The exponential map thus maps A ∈ gl(n) into eA ∈ GL(n).
2. (2)$ePAP−1=PeAP−1$ for any nonsingular $P∈ℝn×n$.
3. (3)$(d/dt)eAt=AeAt=eAtA$ and $det(eA)=etr(A)$.
4. (4)For any A, B ∈ gl(n), eAeB = eC, where C ∈ gl(n) is of the formDisplay Formula
(4)$C=A+B+12[A,B]+112{[[B,A],A]+[[B,A],B]}+⋯$

If G is a matrix Lie group, its corresponding matrix Lie algebra, denoted g, is defined to be the set of all n × n real matrices whose exponential is an element of G. More preciselyDisplay Formula

(5)$g={A∈gl(n) | eAt∈G ∀t∈ℝ}$
with the Lie bracket taken to be the matrix Lie commutator. Equivalently, g is the tangent space to G at the identity element. That is, let X(t) be a differentiable curve on G such that X(0) = $I$. Then $X˙(0)$, the derivative of X(t) at t = 0, is a tangent vector to G at $1$. The collection of all such vectors $X˙(0)$, obtained as the derivative of all possible differentiable curves X(t) on G such that X(0) = $I$, then constitutes the Lie algebra g. g is a vector space—a vector subspace of gl(n), to be precise—and the exponential map $exp :g→G$ is defined from the Lie algebra to the Lie group. As an example, the Lie algebra of GL(n) is gl(n). The Lie algebra of the special linear group SL(n) is the set of n × n real matrices with trace zero, denoted sl(n).

Having defined the Lie algebra of a Lie group, we now define a pair of mappings that play a central role in describing how various physical quantities that arise in multibody dynamics transform under a change of reference frames. Given an element Ag of the Lie algebra, the small adjoint map adA: gg is simply the Lie bracketDisplay Formula

(6)$adA(B)=[A,B]$

Given some element XG of the Lie group, the large adjoint map AdX: gg is defined asDisplay Formula

(7)$AdX(A)=XAX−1$

If X ∈ GL(n) and A ∈ gl(n) satisfy eA = X, then A is said to be a matrix logarithm of X. The matrix logarithm, or matrix log for short, satisfies $log eA=A$ and $e log X=X$. The matrix log need not be unique, but if X is restricted to some sufficiently small neighborhood of the identity, then the matrix log is indeed unique and defined by the convergent series Display Formula

(8)$log X=(X− 1)−12(X− 1)2+13(X− 1)3−⋯$

If X is now restricted to be an element of some matrix Lie group G, then over some neighborhood $N$ of the identity element of G, the matrix log is uniquely defined and acts as a set of local coordinates for $N⊆G$. That is, let {Ω1,…, ΩN} be a basis for g. The map $exp :ℝN→N$ that takesDisplay Formula

(9)$(q1,…,qN)↦eΩ1q1+⋯+ΩNqN$
then defines a set of local coordinates q for $N$, with the matrix log as its inverse. Local coordinates constructed in this fashion are referred to as canonical coordinates of the first kind [14]. Another way to construct local coordinates for $N$ is via the mapDisplay Formula
(10)$(q1,…,qN)↦eΩ1q1⋯eΩNqN$

Such coordinates are called canonical coordinates of the second kind. It is not always necessary for the {Ω1,…, ΩN} to form a basis in the case of canonical coordinates of the second kind.

For any group element XG, the left translation map LX: GG is defined asDisplay Formula

(11)$LX(Y)=XY$

Since $LX−1X=X−1X= 1$, the derivative of $LX−1$ at X is a linear mapping from the tangent space TXG to $T 1G=g$ (see Fig. 1). Given a differentiable curve X(t) on G with velocity $X˙(t)$, the differential of $LX−1$, denoted $dLX−1(X):TXG→g$ and also denoted as the left differential at X, maps the tangent vector $X˙$ to an element of g as follows:Display Formula

(12)$dLX−1(X)(X˙)=X−1X˙$

The notion of the right-translation map and right differential can also be defined analogously: RX(Y) = Y X, andDisplay Formula

(13)$dRX−1(X)X˙=X˙X−1$

The left- and right-translation maps provide two natural ways of mapping any tangent vector $X˙∈TXG$ to the Lie algebra g: via the differential of the left translation map as $X−1X˙$ or the differential of the right-translation map as $X˙X−1$. If A(t) ∈ g is differentiable, and X(t) = eA(t)G, then it can be shown that

$X−1X˙=∫01e−A(t)sA˙(t)eA(t)s dsX˙X−1=∫01eA(t)sA˙(t)e−A(t)s ds$

The left and right differential maps offer a natural way of constructing a Riemannian metric on G—a Riemannian metric on a manifold $M$ is a smoothly varying family of inner products defined at each tangent space to $M$—simply by constructing an inner product on the Lie algebra g. Let X(t) denote a smooth curve on G with velocity $X˙(t)$. Then both $X−1X˙$ and $X˙X−1$ are elements of g, which recall also is a vector space. By choosing an inner product $〈·,·〉 1$ on g, an inner product $〈·,·〉X$ on the tangent space TXG can be constructed asDisplay Formula

(14)$〈X˙,X˙〉X=〈X−1X˙,X−1X˙〉 1$

The above is a left-invariant Riemannian metric: it is left-invariant in the sense that if the original curve X(t) is now multiplied from the left by some constant TG, to Y (t) = TX(t), thenDisplay Formula

(15)$〈Y˙,Y˙〉Y=〈X˙,X˙〉X$

That is, the resulting inner product is invariant with respect to left translations by T. The right-invariant Riemannian metric is, following an analogous construction but this time replacing $X−1X˙$ by $X˙X−1$, invariant with respect to right translations by T.

For the left- and right-invariant Riemannian metrics defined earlier, the minimal geodesics, i.e., minimum-length paths between two arbitrary endpoints on G, can be determined by solving an appropriate optimal control problem. To obtain minimal geodesics with respect to the left-invariant metric, given endpoints X(0) and X(1) on G, we seek the trajectory U(t) ∈ g, t ∈ [0, 1] that minimizesDisplay Formula

(16)$J(U)=∫01〈U(t),U(t)〉 dt$

subject to the matrix differential equationDisplay Formula

(17)$X˙(t)=X(t)U(t)$

where $〈·,·〉$ denotes an inner product on g. (If minimal geodesics with respect to the right-invariant Riemannian metric are desired, the corresponding differential equation would be expressed as $X˙=UX$.) Choosing an inner product on g of the formDisplay Formula

(18)$〈A,B〉=tr(ATQB)$
for some symmetric positive-definite $Q∈ℝn×n$, the optimal U(t) and X(t) must satisfy the differential equationsDisplay Formula
(19)$X˙=XU$
Display Formula
(20)$U˙=Q−1UTQU−UUT$
with boundary conditions X(0) and X(1) given, leading to a classical two-point boundary value problem. In the event that $Q= 1$, the second equation simplifies to $U˙=[UT,U]$.

The solution X(t) to Eqs. (19) and (20) represents a minimum velocity curve on SE(3). Minimum acceleration curves can be found by minimizing the objective functionDisplay Formula

(21)$J(U)=∫01〈U˙(t),U˙(t)〉 dt$
subject to the matrix differential equation (17) with boundary conditions X(0), $X˙(0)$, X(1), and $X˙(1)$ given. For the choice $Q= 1$ for the metric on g, it can be shown that the optimal U(t) and X(t) must then satisfy the following differential equations:Display Formula
(22)$X˙=XU$
Display Formula
(23)$U⃛=[UT, U¨]$

###### The Rotation Group.

The set of 3 × 3 rotation matrices forms the special orthogonal group SO(3)Display Formula

(24)$SO(3)={R∈ℝ3×3 | RTR= 1, det R=1}$

SO(3) is a three-dimensional Lie group under matrix multiplication and describes all the possible orientations of a rigid body in physical space. The Lie algebra of SO(3), denoted so(3), consists of the real 3 × 3 skew-symmetric matricesDisplay Formula

(25)$so(3)={Ω∈ℝ3×3| Ω+ΩT=0}$

Since any element of so(3) has precisely three independent components, so(3) can be identified with $ℝ3$. The following notation is used: given $ω=(ω1,ω2,ω3)T∈ℝ3$, its skew-symmetric representation will be denoted asDisplay Formula

(26)$[ω]=[0−ω3ω2ω30−ω1−ω2ω10]$

For any [ω], [v] ∈ so(3) (or $ω,v∈ℝ3$), and R ∈ SO(3), the Lie bracket [⋅,⋅], small adjoint adω: so(3) → so(3), and large adjoint AdR: so(3) → so(3) areDisplay Formula

(27)$[[ω],[v]]=[ω][v]−[v][ω]=[ω×v]$
Display Formula
(28)$adω(v)=ω×v=[ω]v$
Display Formula
(29)$AdR([ω])=R[ω]RT=[Rω]$

Equation (29) can also be written in the vector form AdR(ω) = .

As is well-known from classical kinematics, any orientation can be achieved by rotating a frame about some axis $ω̂$ by some specified angle θ. The axis-angle formula is an explicit representation of this property that can be expressed through the matrix exponential operator. Given $[ω̂]θ∈so(3)$ such that $θ∈ℝ$ and $‖ω̂‖=1$, the exponential map $exp :so(3)→SO(3)$ can be expressed in closed-form as follows:Display Formula

(30)$e[ω̂]θ= 1+sin θ [ω̂]+(1−cos θ) [ω̂]2$

The above formula is also known as the Rodrigues rotation formula.

A closed-form formula for the matrix log on SO(3) is given byDisplay Formula

(31)$[ω̂]=12 sin θ(R−RT)$
Display Formula
(32)$θ=cos−1(tr(R)−12)$
for any R such that tr(R) ≠ −1. The formula in the general case is covered in, e.g., Ref. [15]. Other mathematical representation for rotations, e.g., Euler angles (which correspond to canonical coordinates of the second kind on SO(3)), Cayley–Rodrigues parameters, unit quaternions, etc., are discussed in Refs. [6] and [1518].

We now examine the angular velocity of a rigid body from the Lie group perspective. Suppose a rigid body is rotating in space with angular velocity w. Let {s} denote the fixed (or space) frame, and {b} denote the moving (or body) frame attached to the rotating rigid body. The orientation of {b} relative to {s} is then described by the rotation R ∈ SO(3). Let $ωs∈ℝ3$ be the column vector representation of w in {s} frame coordinates. ωs is then related to R as follows:Display Formula

(33)$[ωs]=R˙R−1$
which corresponds to the differential of the right-translation map. If $ωb∈ℝ3$ is the column vector representation of w in {b} frame coordinates, then from Eq. (29) it follows that Display Formula
(34)$[ωb]=R−1R˙$
which corresponds to the differential of the left translation map.

###### The Euclidean Group.

Let R ∈ SO(3) and $p∈ℝ3$ be the orientation and position of the rigid body relative to the fixed frame. The special Euclidean group SE(3), which we shall refer to more simply as the Euclidean group, consists of 4 × 4 real matrices of the formDisplay Formula

(35)$[Rp01]$

where R ∈ SO(3) and $p∈ℝ3$ (the 0 in the last row denotes a three-dimensional row vector of zeros). For convenience, we shall at times also write (R, p) ∈ SE(3). SE(3) is a six-dimensional Lie group under matrix multiplication and describes the set of all possible relative displacements between two reference frames situated in physical space.

The Lie algebra of SE(3), denoted se(3), consists of 4 × 4 matrices of the formDisplay Formula

(36)$[S]=[[ω]v00]$

where [ω] ∈ so(3) and $v∈ℝ3$. $S$ will also be expressed in six-dimensional column vector form as $S=(ωT,vT)T$, or more conveniently, $S=(ω,v)$.

The Lie bracket of two elements $S1=(ω1,v1), S2=(ω2,v2)$ can be evaluated as follows:Display Formula

(37)$[[S1],[S2]]=[S1][S2]−[S2][S1]=[[ω1×ω2][ω1]v2−[ω2]v100]$

The above can be written more compactly using the small adjoint map $adS1(S2)$Display Formula

(38)$adS1(S2)=[[ω1]0[v1][ω1]][ω2v2]$

or $adS1(S2)=(ω1×ω2,ω1×v2−ω2×v1)$.

Given T = (R, p) ∈ SE(3) and $S=(ω,v)∈se(3)$, the large adjoint map AdT: se(3) → se(3) isDisplay Formula

(39)$AdT([S])=T[S]T−1=[[Rω]Rv+[p]Rω00]$

In vector form, the above can also be written asDisplay Formula

(40)$AdT(S)=[R0[p]RR][ωv]$

or $AdT(S)=(Rω,Rv+p×Rω)$. The notation [AdT] will also be used to denote the 6 × 6 matrixDisplay Formula

(41)$[AdT]=[R0[p]RR]$

In what follows, $[ω̂]∈so(3)$ denotes a vector of unit length, i.e., $‖ω̂‖=1$. Given $Sθ=(ω̂θ,vθ)∈se(3)$, $θ∈ℝ$, the closed-form formula for the exponential map $exp :se(3)→SE(3)$ is given byDisplay Formula

(42)$e[S]θ=[e[ω̂]θG(θ)v01]$

where $e[ω̂]θ$ is given by Eq. (30) andDisplay Formula

(43)$G(θ)= 1θ+(1−cos θ)[ω̂]+(θ−sin θ)[ω̂]2$

The matrix log is uniquely defined over some neighborhood $N⊂SE(3)$ of the identityDisplay Formula

(44)$log([Rp01])=[[ω̂]v00]θ$


text:
# How can I find the sum of an infinite series of products?

## Background

Today in my macroeconomics class my teacher taught us three concepts.

The first is very simple: consumption $c$ is a linear function of national income $y$. Mathematically, $$c = My + b$$

We will call $M$ the marginal propensity to consume, and refer to it accordingly as MPC. MPC can also be thought of as the fraction of income that is spent rather than saved. So MPC is going to be some number less or equal to one but greater than or equal to zero ($0 \leq M \leq 1$).

The second concept we learned is also quite simple. As income $y$ increases, MPC decreases.

The third concept is a little bit more complicated. Let Person 1 has $5$ dollars. Let MPC be universally constant and equal to $\frac{3}{4}$. Now, since MPC $= \frac{3}{4}$, Person 1 spends $\frac{3}{4}$ of his or her income on Person 2. So Person 2 has an income of $\frac{3}{4} * 5 = 3.75$ dollars. Person 2 will now spend $\frac{3}{4}$ of his or her income on Person 3. And so on to infinity. I instantly recognized this as the sum of a geometric series: $$5*\sum_{n=0}^\infty\left(\frac{3}{4}\right)^{n} = 20$$ This makes sense, but only when disregarding the second concept (MPC decreases as $y$ increases). When I asked about it my teacher said that since calculus isn't a prerequisite for the class we won't go any further in depth; we'll just let MPC be a constant to simplify things in an introductory-level macroeconomics class. So I thought about the problem for a while today and realized that if we did account for MPC decreasing as y increases, or in this case, MPC increasing as $n$ increases, then we could write the total money spent as an infinite series of kind of infinite product (??). My model can be found below.

## The Main Question

I just finished Calculus III last fall so I've never formally learned anything about infinite products. I'm not even sure if these are infinite products since many of them are finite. Anyways, can the following infinite series that I came up with be solved for convergence? If so, how can I figure out the number it converges too? Long explanations will be appreciated. Keep in mind I have a very coarse knowledge of mathematics beyond Calculus III.

$$\sum_{i=0}^{\infty} \left(\prod_{n=0}^{i} \frac{n + 2}{\sqrt{n^2 + 5n + 7}}\right)$$

Perhaps it would be better if I knew how to solve the infinite product first.

Also, I came up with the function $\frac{n + 2}{\sqrt{n^2 + 5n + 7}}$ as a model of MPC because I figured MPC should start around $\frac{3}{4}$ and increase asymptotically towards $1$ as $n$ approaches $\infty$.

• @Shahar Woops! I made a mistake. I've fixed it now. – user168210 Jan 22 '15 at 2:52
• The n+2 term gives (i+2)!. The denominator is harder (duh!). You can ignore the square root and then take the square root of the product, if you can find it. Somehow gamma or theta or q-theta functions come to mind. – marty cohen Jan 22 '15 at 2:59
• I think it may not converge because the asymptotic behaviour of $\sum_{k=0}^n(\log(1+2/k)-\log(1+5/k+7/k^2)/2)$ seems like the harmonic series, which has $-\log(n)$ behaviour, and then $\sum e^{-\log(n)}$ gives a harmonic series again. – user1537366 Jan 22 '15 at 3:13

Not really an answer to your question, but since you seem interested in modelling, I wanted to put some ideas to you. When trying to understand an idea, its good to start simple.

Perhaps you may want to start with a simpler model of a dynamic $M$ as a function of $y$: $M(y):=\min\{\frac{a}{y},1\}\; 0<a$, which implies that people spend roughly $\text{\$} a$regardless of income, unless they can't afford to (hence decreasing propensity). Then, if person$1$has$X$dollars, they spend either$a$or$X$dollars it on person$2$, who subsequently pays person$3$either$a$or$X$dollars. This is a very simple model, and you end up getting that everyone, except person 1, is living "hand-to-mouth"...not fun for the people...but man, what a great multiplier! (i.e.,$\infty\$) ;-)...good example of why maximizing GDP doesn't necessarily maximize psychological well being of a people.

Now, this model is obviously deficient, but in what way? Is the constant propensity model you were presented in class really all that unrealistic for relevant salary ranges (certainly many in the middle class would say its realistic). As an extreme example, consider philanthropy...its done by those with ridiculous incomes and they can spend ridiculous sums on their pet projects. What fraction would we be talking about? Given the small number of such people, are they relevant to economics? Or do we try to model the middle?

Anyway, its great that you are pushing beyond what you've been taught formally...that's the kind of curiosity that leads to new insights! Good luck with your studies.


text:
# Math Help - Linear Regression

1. ## Linear Regression

Hello guys,

I have a small question please:

An experiment was made to find correlation between X and Y.
2 Regression lines has been build:

1. Y=a+bx

2. X=n+mY

It is known that: m=-0.4 and b=-0.625

Find the correlation r.

cheers !

2. using the standard formula for a OLS slope coefficient:

$-0.4 = \frac{Cov(X,Y)}{Var(Y)}$
$-0.625 = \frac{Cov(X,Y)}{Var(X)}$

multiply together

$-0.4 * -0.625 = \frac{Cov(X,Y)^2}{Var(Y) Var(X)}$
$-0.4 * -0.625 = r^2$
r=0.5 or -0.5

interesting problem. I agree with you that you should choose the negative solution since both coefficients are negative.


text:
source("http://macosa.dima.unige.it/r.R")    # If I have not already loaded the library
---------- ---------- ---------- ---------- ---------- ---------- ---------- ----------
S 08 Integrals and Derivatives (and Limits). Taylor polynomial. Volumes.

# deriv(f,VariableName), deriv2, …, deriv6

#  Note:  deriv(f,"a")  equals  D(body(f),"a")

G = function(u) sin(u); deriv(G,"u")
# cos(u)
f = function(x) 3*x^3-5*x^2+x+3
graphF(f, -2,4, "blue")
deriv(f,"x")
#  3 * (3 * x^2) - 5 * (2 * x) + 1
# To use the expression produced by "deriv" you need to evaluate it with eval:
df = function(x) eval(deriv(f,"x"))
graph(df,-3,5, "red")
# I could also define  df  by writing  3*(3*x^2)-5*(2*x)+1   [ or 9*x^2-10*x+1 ]
deriv2(f,"x")
#  3 * (3 * (2 * x)) - 5 * 2
d2f = function(x) eval(deriv2(f,"x")); graph(d2f,-3,5, "violet")
deriv3(f,"x")
#  3 * (3 * 2)
d3f = function(x) eval(deriv3(f,"x"))+x-x      # I put the "x" with a trick
graph(d3f,-3,5, "green")

# To evaluate the derivative of log10(x) and log2(x) you must use log(x)/log(10) and
# log(x)/log(2), where log(x) is the inverse of exp(x)
f=function(x) log(x); g=function(x) log(x)/10; h=function(x) log(x)/2
deriv(f,"x"); deriv(g,"x"); deriv(h,"x")
# 1/x  1/x/10  1/x/2

# To evaluate the derivative of a function where abs(u) appears, it needs to be
# replaced with sqrt(u^2)
#
# If f is a piecewise-defined function we must derive every sub-function (see).

#
# Note. df calculated by the software can have a wider domain than f. Example:

f = function(x) log(1+tan(x/2))
BF=5; HF=3; a=-2*pi; b=2*pi; Plane(a,b, -5,5)
graph2(f, a,b, "brown"); df = function(x) eval( deriv(f,"x") ); deriv(f,"x")
# 1/2/cos(x/2)^2/(1 + tan(x/2))
graph2(df, a,b, "seagreen")
FDOM=f; graph(DOM, a,b, "red")                # DOM is defined 0 in the domain of FDOM
pointO(pi*c(-1,-1/2,1,2-1/2),rep(0,4),"blue")
text(-5,-1/2,"domain",font=2,col="red"); text(-1/2,-1.5,"f",font=2,col="brown")
text(-1/2,2.5,"df",font=2,col="seagreen"); text(-2.5,-1.5,"?",font=2,col="seagreen")
#
# I can calculate the integral of a piecewise linear function (with coordinates in x
# and y), using the areaPol command, if I define:
INTEGRAL = function(x,y) areaPol( c(x[length(x)],x[1],x), c(0,0,y) )
# Example:

x = c(10,13, 17, 19); y = c(7,1,9,8)
Plane(min(x),max(x), 0,max(y)); polyl(x,y,"blue")
INTEGRAL = function(x,y) areaPol( c(x[length(x)],x[1],x), c(0,0,y) )
INTEGRAL(x,y)
# 49

# We can use areaPol also for spline functions (see):
x = c(-0.2, 1.3, 1.7,  3.2,  4.1,  4.6)
y = c(9.1, 17.5, 20.7, 25.4, 35.7, 51.2)
Plane(-1,6,-10,70);  POINT(x,y,1); L=spline(x,y,1000); polyline(L$x,L$y,"red")
INTEGRAL(L$x,L$y)   # INTEGRAL defined above
# 108.6348

# But I can use the integral(f, a,b) command, if I define such function explicitly:

f = function(x) ifelse(x<13, 7-6/3*(x-10), ifelse(x<17, 1+8/4*(x-13), 9-1/2*(x-17)))
Plane(10,19, 0,9); graph(f, 10,19, "blue")

integral(f, 10,19)
# 49

# A function that can be described more easily:

g = function(x) sin(x)
integral(sin, 0,pi/2); integral(sin, 0,pi); integral(sin, 0,2*pi)
#           1                2                2.032977e-16  ("pratically" 0)

# Here, how to calculate the effective force on the mast of a racing sailbot
#       (from "Numerical Methods", Chapra & Canale, McGraw-Hill)

# See here for integrals of F from A to B when A or B are infinite or F is undefined
# at A or B .

# Definite integrals of many funcions can be calculated with the areaF command: see.
# Only when we have focused the theorem fundamental theorem of calculus (see here or
# here) we can interpret areas as defined integrals.

#
# Gintegral(f,A,B,col)  traces the graph of x -> ∫[A,x]f from A to B (A and B must be
# finite);  Gintegra, Gintegr plot slimmer graphics.

Plane(-pi,2*pi, -1,1); graph1(cos, -pi,2*pi, "blue")
Gintegr(cos, 0,2*pi, "red"); Gintegr(cos, 0,-pi, "violet")
POINT(0,0, "red")

# The graph on the right:

Plane(-pi,2*pi, -1,1); graph1(cos, -pi,2*pi, "blue")
POINT(2.5,0.5, "red")
GintegrK(cos, 2.5,2*pi, 0.5, "brown"); GintegrK(cos, 2.5,-pi, 0.5, "seagreen")

# With GintegralK(f,A,B, k, col) (and GintegraK, GintegrK) I can trace the graph
# of x -> k + ∫[A,x]f

# In very special cases Gintegral (Gintegra, Gintegr) may not work; in such cases you
# can try using GintegralA or, if it is not working, GintegralB or … GintegralD, or
# Gintegral0  (or GintegraA, GintegrA, … Gintegra0, Gintegr0, for slimmer graphics):

f = function(x) ifelse(x<13, 7-6/3*(x-10), ifelse(x<17, 1+8/4*(x-13), 9-1/2*(x-17)))
Plane(10,19, 0,50); Gintegra(f, 10,19, "blue")
# Error in integrate
Plane(10,19, 0,50); GintegraA(f, 10,19, "brown")

# For the same reasons if integral does not work you can use integralA, …, intergralD
# or intergral0.
#
# In some case, some manipulations are necessary to achieve a result: a simple example.
#
# With R we can compute definite integrals. We can also calculate the indefinite
# integrals (ie antiderivatives) of polynomial functions.
Q = c(-2,-3,2.4,-1/5)                       # the polynomial -2-3*x+2.4*x^2-1/5*x^3
q = function(x) Q[1]+x*Q[2]+x^2*Q[3]+x^3*Q[4]
# the integral:
I = function(p) {I=NULL; for (i in 1:length(p)) I = c(I,p[i]/i); I}
I(Q); fraction(I(Q))
#  -2.00 -1.50  0.80 -0.05
#     -2  -3/2   4/5 -1/20                   -2*x-3/2*x^2+4/5*x^3-1/20*x^4 + constant
Iq = function(x) x*I(Q)[1]+x^2*I(Q)[2]+x^3*I(Q)[3]+x^4*I(Q)[4]
# ∫[-1,2] q  using I(Q)
Iq(2)-Iq(-1)
#  -4.05
# I can calculate the defined integer in each case, with:
integral(q, -1,2)
#  -4.05

# For other functions I can use WolframAlpha. Let's see how to use it in this case:
#  integrate -2-3*x+2.4*x^2-1/5*x^3 dx
#  -0.05 x^4 + 0.8 x^3 - 1.5 x^2 - 2 x + constant
#  integrate -2-3*x+2.4*x^2-1/5*x^3 dx from x=-1 to 2
#  -0.45
#
# The elementary functions (functions of one variable which are the composition of
# arithmetic operations, powers, exponentials, logarithms, trigonometric functions,
# inverse trigonometric functions) are closed under differentiation but not under
# integration. Some examples of elementary functions that have not an elementary
# antiderivative:  sqrt(1+x^3)  sqrt(1-x^4)  1/sqrt(1+x^4)  1/log(x)  log(log(x))
# exp(x^2)  exp(-x^2)  exp(x)/x  exp(exp(x))  x^2*exp(-x^2)  sin(x^2)  cos(exp(x))
# sin(x)/x. Hence the importance of the numerical integration we have seen.
#
# We can also trace the direction field of the solutions of an indefinite integral (see):

#
# See here for the compute of limits of a function
#
# See here for the "study of a function"

# See here for the correct use of x -> x^x and similar function that use ^.
x -> x^x

# It is easy to integrate numerically and graphically functions that are not easy to
# study "by hand". See here for an example.
# An example of integral that you can not study by hand:
F = function(x) sqrt(3*x-x^3); integral(F, 0,3/2)
# 1.720722

# Note 1. If the integral between -Inf and Inf is 0 check whether the integer between 0
# and Inf is finite too: if it does not happen it means that the integral diverges (0
# is output only for a formal simplification).

# Note 2. To integrate a function between a and b if it is not defined in a point c
# between a and b, make the sum of the integral between a and c and that between c and b.
#
# See here  how to calculate areas using integrals

# If F is a function of at least 2 variables, if two of these are u and v,
# derivxy(F,u,v) computes the second partial derivative with respect to u and v,
# derivyx(F,u,v) that with respect to v and u (xy and yx indicate that the derivatives
# are "crossed", even if the variables are not x and y). deriv2(F,u) computes the 2nd
# derivative respect to u.
# If F is a function of x and y, hessian(F,a,b) (see) computes the hessian of F in (a,b).
# (we remember that if in (x,y) all partial derivatives are zero and the hessian is
#  negative/null/positive, then  it is a saddle point / we cannot conlude anything /
#  it is a min/max point if deriv2 respect to x is > / < 0 )

F = function(x,y,z) 5*x^2*y^3+x^3*z
deriv(F,"x"); deriv(F,"z")
# 5*(2*x)*y^3+3*x^2*z    x^3
deriv2(F,"y"); deriv2(F,"z")
# 5*x^2*(3*(2*y))        0
derivxy(F,"x","y"); derivxy(F,"x","z")
# 5*(2*x)*(3*y^2)        3*x^2
#
Q = function(x,y) x^4+6*x^2*y^2+y^4-32*(x^2+y^2)

# deriv(Q,"x") and deriv(Q,"y") are 0 in (0,0).
# To calculate hessian the variables must be x and y
hessian(Q,0,0); hessian(Q,0,4); hessian(Q,2,-2)
# 4096   16384  -8192
deriv2(Q,"x")
# 4 * (3 * x^2) + 6 * 2 * y^2 - 32 * 2
# In (0,0) the hessian in > 0 and deriv2 respect to x is > 0: max point
# And indeed:
MAXF2(Q, -1,1, -1,1)
#  0 0  0    In (0,0) I have a max, and Q(0,0) = 0
#
# The gradient of Q, grad(Q), is x,y -> ( d(Q(x,y)/dx, d(Q(x,y)/dy ):
c( deriv(Q,"x"), deriv(Q,"y") )
# 4*x^3+12*x*y^2-64*x  12*x^2*y+4*y^3-64*y
# For the use of this concept see here.
#
# For Taylor polynomials see here.

#
# For integrals of the functions of two variables we can use the command INTEGRAL. See here.

# For volumes of solids of revolution see here.

# For max and min of functions of two variables see here (commands MINF2, MAXF2).
#
# The mean value (not of a set of values but) of a variable quantity depends on
# which variable one wants to make it dependent on.
# If f is a continuous positive function on [a,b] we assume as mean of f on [a,b] the
# value m such that the area of the drawn rectangle equals the result of integrating f
# on [a,b]:  m = integral(f, a,b) / (b-a)

# For example, the problem:
#     Calculate the average value of the chord of a semi-circle
#          (of radius 1) that is parallel to the diameter
# depends on whether, for example, we identify the position of the chord using
# L, D or A. I can obtain 1, π/2 or 4/π.

FL = function(x) x; I = integral(FL, 0,2) / 2; I
# 1
FD = function(x) 2*sqrt(1-x^2); I = integral(FD, 0,1)/1; I
# 1.570796
pi/2
# 1.570796
FA = function(t) 2*sin(t); I = integral(FA, 0,pi/2) / (pi/2); I
# 1.27324
I*pi
# 4

Other examples of use


text:
## A constrained eigenvalue problem.(English)Zbl 0666.15006

For the problem of finding $$Min(x^ TAx)$$ for a symmetric matrix A subject to $$x^ Tx=1$$ and $$N^ Tx=t$$ theoretical and numerical methods are described. First the linear constraint is removed and then Lagrange multipliers are employed reducing the problem to solve either a secular equation or a quadratic eigenvalue problem.
Reviewer: V.Mehrmann

### MSC:

 15A18 Eigenvalues, singular values, and eigenvectors 65F15 Numerical computation of eigenvalues and eigenvectors of matrices

GQTPAR
Full Text:

### References:

 [1] Draper, N.R., “ridge analysis” of response surfaces, Technometrics, 5, 469-479, (1963) · Zbl 0124.35003 [2] Forsythe, G.E.; Golub, G.H., On the stationary values of a second-degree polynomial on the unit sphere, SIAM J. appl. math., 13, 1050-1068, (1963) · Zbl 0168.03005 [3] Gander, W., Least squares with a quadratic constraint, Numer. math., 36, 291-307, (1981) · Zbl 0437.65031 [4] Golub, G.H., Some modified matrix eigenvalue problems, SIAM rev., 15, 318-334, (1973) · Zbl 0254.65027 [5] Golub, G.H.; Heath, M.; Wahba, G., Generalized cross-validation as a method for choosing a good ridge parameter, Technometrics, 21, 215-223, (1979) · Zbl 0461.62059 [6] Golub, G.H.; Van Loan, C.F., Matrix computations, (1983), John Hopkins U.P Baltimore · Zbl 0559.65011 [7] Moré, J.J., The Levenberg-Marquardt algorithm: implementation and theory, (), 105-116 · Zbl 0372.65022 [8] Moré, J.J.; Sorensen, D.C., Computing a trust region step, SIAM J. sci. statist. comput., 4, 553-572, (1983) · Zbl 0551.65042 [9] Reinsch, Chr.H., Smoothing by spline functions, II, Numer. math., 16, 451-454, (1971) · Zbl 1248.65020 [10] Spjøtvoll, E., A note on a theorem of forsythe and Golub, SIAM J. appl. math., 23, 307-311, (1972) · Zbl 0259.15019 [11] Spjøtvoll, E., Multiple comparison of regression functions, Ann. math. statist., 43, 1076-1088, (1972) · Zbl 0241.62046 [12] von Matt, U., A constrained eigenvalue problem, () · Zbl 0742.65044
This reference list is based on information provided by the publisher or from digital mathematics libraries. Its items are heuristically matched to zbMATH identifiers and may contain data conversion errors. It attempts to reflect the references listed in the original paper as accurately as possible without claiming the completeness or perfect precision of the matching.


text:
96-490 Gerald Teschl
Renormalized Oscillation Theory for Dirac Operators (36K, LaTeX2e) Oct 15, 96
Abstract , Paper (src), View paper (auto. generated ps), Index of related papers

Abstract. Oscillation theory for one-dimensional Dirac operators with separated boundary conditions is investigated. Our main theorem reads: If $\lambda_{0,1}\in \mathbb R$ and if $u,v$ solve the Dirac equation $H u= \lambda_0 u$, $H v= \lambda_1 v$ (in the weak sense) and respectively satisfy the boundary condition on the left/right, then the dimension of the spectral projection $P_{(\lambda_0, \lambda_1)}(H)$ equals the number of zeros of the Wronskian of $u$ and $v$. As an application we establish finiteness of the number of eigenvalues in essential spectral gaps of perturbed periodic Dirac operators.

Files: 96-490.tex


text:
## how to use mort.raw?

Post here to share useful tips and tricks, to ask questions about using your DM42 or to report software-related problems
Walter
Posts: 1610
Joined: Tue May 02, 2017 11:13 am
Location: Close to FRA, Germany

### Re: how to use mort.raw?

HPMike wrote:
Sun May 24, 2020 3:53 pm
DHCrocks wrote:
Sun May 24, 2020 7:09 am
oh, what is DUE variable?
It is a setting to determine if payments are to be made at the end of a payment period (arrears), in which case DUE=0, or if they are made at the beginning of a payment period when due, in which case DUE=1. When payments are made when due, they are lower, since no interest will accrue during that period on the current unpaid balance. In this case, with DUE=0, then PMT=45.8037941842, whereas if DUE=1, then PMT=44.4697030914.
IMO, this DUE is an inartful term here since the payment is due anyway, either at beginning or end of each period. A variable called BEGIN would be more appropriate for what it does.
DM42 SN: 00041 Beta
WP 43S running on this device

HP-35, HP-45, ..., HP-50, WP 34S, WP 31S, DM16L

PierreMengisen
Posts: 118
Joined: Wed Nov 29, 2017 1:38 pm
Location: Neuchâtel CH

### Re: how to use mort.raw?

in French we talk about pre-numerando and post-numerando payment.
These are two modes of treatment. Why not talk about "KOP" (kind of payment)?
Pierre
[TI59 with PC100C; TI-84 Plus CE-T; HP41CV with HP IL loop & 2*82161A DCD & 82162 TP; HP15C; HP28S; DM41; DM41L; DM42; DM41X]

HPMike
Posts: 205
Joined: Fri Jul 21, 2017 11:01 pm
Location: DFW, Texas

### Re: how to use mort.raw?

Walter wrote:
Sun May 24, 2020 5:12 pm
HPMike wrote:
Sun May 24, 2020 3:53 pm
DHCrocks wrote:
Sun May 24, 2020 7:09 am
oh, what is DUE variable?
It is a setting to determine if payments are to be made at the end of a payment period (arrears), in which case DUE=0, or if they are made at the beginning of a payment period when due, in which case DUE=1. When payments are made when due, they are lower, since no interest will accrue during that period on the current unpaid balance. In this case, with DUE=0, then PMT=45.8037941842, whereas if DUE=1, then PMT=44.4697030914.
IMO, this DUE is an inartful term here since the payment is due anyway, either at beginning or end of each period. A variable called BEGIN would be more appropriate for what it does.
BEGIN is not a variable either, since it represents only one case for a binary choice. On my HP 12C, there are two choices, BEGIN for payments made at the beginning of a payment period, and END for payments made at the end of a period. So DUE in fact makes perfect sense for the variable, as it requests the answer to the question of when is the payment due for each payment period. Perhaps the variable name should be WD for “when due.”
DM15L, S/N 00548. DM42, SN: 00159. DM41X, SN: 00973.

Walter
Posts: 1610
Joined: Tue May 02, 2017 11:13 am
Location: Close to FRA, Germany

### Re: how to use mort.raw?

HPMike wrote:
Sun May 24, 2020 7:54 pm
Walter wrote:
Sun May 24, 2020 5:12 pm
HPMike wrote:
Sun May 24, 2020 3:53 pm
It is a setting to determine if payments are to be made at the end of a payment period (arrears), in which case DUE=0, or if they are made at the beginning of a payment period when due, in which case DUE=1. When payments are made when due, they are lower, since no interest will accrue during that period on the current unpaid balance. In this case, with DUE=0, then PMT=45.8037941842, whereas if DUE=1, then PMT=44.4697030914.
IMO, this DUE is an inartful term here since the payment is due anyway, either at beginning or end of each period. A variable called BEGIN would be more appropriate for what it does.
BEGIN is not a variable either, since it represents only one case for a binary choice. On my HP 12C, there are two choices, BEGIN for payments made at the beginning of a payment period, and END for payments made at the end of a period. So DUE in fact makes perfect sense for the variable, as it requests the answer to the question of when is the payment due for each payment period. Perhaps the variable name should be WD for “when due.”
You're right in principle (a variable BEGIN presumes you know the only alternative is end - so if BEGIN = false then end must be true and vice versa). But a variable DUE or even WhenDue by its name asks for a date, or for an answer reading begin or end; none of these answers is true or false. So any intelligent compiler would throw a data type mismatch here.
DM42 SN: 00041 Beta
WP 43S running on this device

HP-35, HP-45, ..., HP-50, WP 34S, WP 31S, DM16L

akaTB
Posts: 521
Joined: Tue May 02, 2017 1:56 pm
Location: Milan, Italy

### Re: how to use mort.raw?

Just to add a little more confusion: in italian DUE means TWO.
Greetings,
Massimo
ajcaton
-+×÷ left is right and right is wrong Casted in gold

Geoff Quickfall
Posts: 326
Joined: Wed May 24, 2017 10:44 pm
Contact:

### Re: how to use mort.raw?

Due = two

Well that is too funny to!

DHCrocks
Posts: 5
Joined: Sat May 23, 2020 1:09 am

### Re: how to use mort.raw?

so sorry but this is my first time using the DM42, I previously had various models of 48, so this it a bit different.

can some one walk me thru the steps? as an example, this is the loan I would like to calculate. $50,000 loan [pv], 60 months [60], 2.5% interest this is what i'm doing [xeq] [mort] 0 [due] 0 [fv] 50,000 [pv] 0.025 60 [n] [pmt][pmt] the result is 1,617.67 which is obviously wrong it should be returning 887.37 thanks for the guidance. DHCrocks Posts: 5 Joined: Sat May 23, 2020 1:09 am ### Re: how to use mort.raw? I just got my DM42 and are not used to using it, I'm coming from using various HP48 for the past 20 years or so. Can someone please walk me thru the steps? I just cant seem to get it to work this is what I try: [xeq] [mort] 0 [due] 0 [fv] 50,000 [pv] 2.5 60 [n] [pmt] [pmt] result is 125,000 which is not correct it should be 887.37. this is for a 50,000 loan, 5 years, 2.5% interest. thank you for the guidance and patients. HPMike Posts: 205 Joined: Fri Jul 21, 2017 11:01 pm Location: DFW, Texas ### Re: how to use mort.raw? DHCrocks wrote: Mon May 25, 2020 8:25 am so sorry but this is my first time using the DM42, I previously had various models of 48, so this it a bit different. can some one walk me thru the steps? as an example, this is the loan I would like to calculate.$50,000 loan [pv], 60 months [60], 2.5% interest

this is what i'm doing
[xeq]
[mort]
0 [due]
0 [fv]
50,000 [pv]
0.025
60 [n]
[pmt][pmt]

the result is 1,617.67 which is obviously wrong it should be returning 887.37

thanks for the guidance.

You forgot to divide the annual interest by twelve to get monthly compounding interest. So, monthly interest = 0.025 / 12 = 0.00208333 will give you the result 887.3680803. Note that interest is entered as decimal value, not a percentage.
DM15L, S/N 00548. DM42, SN: 00159. DM41X, SN: 00973.

DHCrocks
Posts: 5
Joined: Sat May 23, 2020 1:09 am

### Re: how to use mort.raw?

HPMike wrote:
Mon May 25, 2020 2:42 pm

You forgot to divide the annual interest by twelve to get monthly compounding interest. So, monthly interest = 0.025 / 12 = 0.00208333 will give you the result 887.3680803. Note that interest is entered as decimal value, not a percentage.

That’s it! Thank you so much. That’s a strange way of inputting interest. I just hope I can remember it when I actually need to use the program a year from now. lol.


text:
# Formalising Definite Description (Russell)

So, we've been doing a lot of formalising of definite descriptions according to Russell, where we take the Definite Description as a predicate rather than designator.

So, I've formalised sentences like "The major of london ......" with this kind of predicate:

Mx: x is a major of London

But, my question is if I can formalise this even more?

On a recent test I, out of habit, made this key (because I was once told "the more we formalise, the better"):

Mxy: x is a major of y

l: London

But in hindsight I fear this may have been a mistake ...

I translated the sentence "The major of London is happy" as

$$\exists x \Big(Mxl \land \forall y(Myl \rightarrow y=x) \land Hx \Big)$$ but, as said, I fear it that this is instead correct: $$\exists x \Big(Mx \land \forall y(My \rightarrow y=x) \land Hx \Big)$$ where Mx stands for 'x is a major of London'.

Any ideas? I can't find anywhere online that says anything about this. The only examples I can find is of the second style, which is why I am worried.

Thanks :)

A good principle in formalizing is expose just as much logical structure as is relevant to the case in hand.

Suppose you were interested only in the argument "Boris is the mayor of London. If Boris is the mayor of London, then business won't be scared away. So business won't be scared away." All the structure you need to expose to show this argument is valid is revealed by $P, P\to Q \therefore Q$.

Suppose next you are interested in the argument "Boris is blonde, Boris is the current mayor of london, so whoever is the current mayor of london is blonde". Here we do need to expose some more structure, but talk of mayors only appears in the context "is the current mayor of London". Then you might as well treat "is the current mayor of London" as a fused whole, with no locally significant internal structure, and symbolise "Boris is the current mayor of London" as $Mb$.

But suppose that you want to formalize the argument "Boris is the mayor of London, Bill is the mayor of New York; different cities have different mayors, London is not New York, hence Boris is not Bill". Then you will evidently need to expose the relational structure in the claim "Boris is a mayor of London", and formalize that as $Mbl$.

So there is no one "right" formalisation for "Boris is a mayor of London" or "Boris is the mayor of London": how much logical structure to expose in your formalisation will depend on context.

• Thank you for a wonderful answer! Although, we haven't really done any arguments with definite description. For now, all we've done is learn how to write a sentence with definite description. So, we've never really had an important context, other than ".. is happy", "... is bald". But your explanation really made it clear how to do it in different ways, depending on the context. Thank you very much! Jan 12 '14 at 14:51
• Glad to be of help. (If you are going to use this site, do up-tick answers you accept, though -- if only as a helpful pointer to other readers.) Jan 12 '14 at 19:59

When you modified the predicate $M$ to mean "$x$ is a mayor of $y$" (in the second gray box), the notation should have become $Mxy$ instead of $Mx$. Apart from that, your proposed formalization looks fine to me.

• Thanks! That was just a silly mistake. Jan 12 '14 at 14:46


text:
Summarizes a data.frame, somewhat differently than the normal summary function of R. The function is mainly useful as a basic EDA tool on data.frames before they are converted to tasks, but can be used on tasks as well.

Columns can be of type numeric, integer, logical, factor, or character. Characters and logicals will be treated as factors.

summarizeColumns(obj)

## Arguments

obj (data.frame | Task) Input data.

## Value

(data.frame). With columns:

name

Name of column.

type

Data type of column.

na

Number of NAs in column.

disp

Measure of dispersion, for numerics and integers sd is used, for categorical columns the qualitative variation.

mean

Mean value of column, NA for categorical columns.

median

Median value of column, NA for categorical columns.

MAD of column, NA for categorical columns.

min

Minimal value of column, for categorical columns the size of the smallest category.

max

Maximal value of column, for categorical columns the size of the largest category.

nlevs

For categorical columns, the number of factor levels, NA else.

Other eda_and_preprocess: capLargeValues, createDummyFeatures, dropFeatures, mergeSmallFactorLevels, normalizeFeatures, removeConstantFeatures
summarizeColumns(iris)#>           name    type na     mean      disp median     mad  min  max nlevs
#> 5      Species  factor  0       NA 0.6666667     NA      NA 50.0 50.0     3


text:
# Is there a standard for math terminology?

I searched for the definition of "quotient" and I felt there were many answers. Is there a place that defines math terms? Something like a standard?

• I'm afraid the answer is: no, there is no such standard. – Wojowu Dec 20 '15 at 10:06
• There are quite a few things associated with the word quotient: en.wikipedia.org/wiki/Quotient_(disambiguation) – Element118 Dec 20 '15 at 10:07
• You can see ISO 31-11 : "about mathematical signs and symbols, their meanings, verbal equivalents and applications". – Mauro ALLEGRANZA Dec 20 '15 at 10:08

The important thing is context. For example, to the average person quotient means something like $\frac{a}{b}$. To an algebraist there are things called quotient groups.
Mathematical terminology is always changing. A "number" a few thousand years ago might have only meant $1,2,3,...$ where now we have much broader ideas about what constitutes a number.


text:
# Thread: Area of a rectangle

1. ## Area of a rectangle

I am stuck with this question... Can anyone help please:

A rectangle is formed by doubling both the length and the width of a rectangle.

What is the area, in cm2, of this new rectangle?

the length is (a+3) cm
the width is 2 cm

= 2(a+3) x 2squared
= 4 (2a + 6)
= 8a + 24 cmsquared

2. ## Re: Area of a rectangle

Yes, that's correct. If the original rectangle has dimensions L=a+3 cm and W=2, then doubling both gives new dimensions of L= 2(a+3) cm =2a+6 cm, and W = 2(2cm) = 4 cm. The area of the new rectangle is then L x W = (2a+6 cm)(4 cm) = 8a +24 cm^2.

3. ## Re: Area of a rectangle

Originally Posted by Natasha1
I am stuck with this question... Can anyone help please:

A rectangle is formed by doubling both the length and the width of a rectangle.

What is the area, in cm2, of this new rectangle?

the length is (a+3) cm
the width is 2 cm
The length and width of which rectangle? If a+ 3 and 2 are the length and the width of the first rectangle and the new rectangle "is formed by doubling both the length and the width of a rectangle" then the new rectangle has length 2(a+ 3)= 2a+ 6 and width 4.


text:
William Blum's website

# Crowdsourcing Cracklock 64-bit

Websites analytics are sometimes revealing. Among people who come to visit my website there is the LaTeX user category coming to learn about my past work on SumatraPDF. Once in a blue moon somebody may come and peek at my research page if I'm lucky. Some folks get here by accident hoping to find information about my American homonym. But the vast majority visits this site for another reason (...)

Published: 11/26/2015 6:25:21 AM
Tags: Windows, 64-bit, Cracklock

# SyncTeX-related news

Most TeX distributions (TeXLive 2008 and MikTex 2.7) have now been updated to support SyncTeX. The option is activated with the -synctex switch at the command-line as follows.

Published: 10/7/2008 7:44:08 PM
Tags: Windows, SyncTex, SumatraPDF, LaTex

# SyncTeX and SumatraPDF

The author of the pdfsync TeX package, J�r�me Laurens, has developed a new technology called SyncTex that will eventually replace pdfsync. SyncTex is directly integrated in pdftex. This has many benefits: there is no need to load a special package in your .tex document, there is no more incompatibilities with latex packages, and the synchronization is more precise. One disadvantage is that the generated synchronization file can be huge, and for that reason it has an option to zip it.

Published: 6/12/2008 2:00:12 AM
Tags: Windows, SyncTex, SumatraPDF, LaTex

# PDF-Latex synchronization (continued)

This is the dual of the previous post: I have now implemented forward-search in SumatraPDF, i.e. the ability to go from the source .tex file to the corresponding location in the PDF. The communication between the TeX editor and SumatraPDF relies on the DDE protocol. (Most TeX editors such as WinEdt and TeXnicCenter support it.)

Published: 5/15/2008 4:54:47 AM
Tags: Windows, SyncTex, SumatraPDF, LaTex

# Synchronizing PDF files with Latex documents

I am currently in the process of writing up my PhD thesis using Latex. I just realized how spoiled Mac users are when it comes to TeX tool support: nearly all Tex editors support viewing PDF files with source synchronization.

Published: 5/14/2008 2:39:33 AM
Tags: Windows, SyncTex, SumatraPDF, LaTex

# My experience with the XML/Docbook/XSLT toolchain

Recently I have tried to learn how to use the Docbook toolchain to produce documentation. My overall impression is not so good. I have many complaints but I am mainly concerned with its poor performance.

Published: 1/21/2008 5:31:52 AM
Tags: XML, Docbook, XSLT, Windows


text:
+0

I have trouble with piecewise functions.

0
102
1

If$f(x) = \begin{cases} 2x-5 &\quad \text{if } x \ge 3, \\ -x + 5 &\quad \text{if } x < 3, \end{cases}$then for how many values of $x$ is $f(f(x)) = 3$?

Dec 31, 2021

f(f(1)) = f(4) = 3.  The only solution is x = 1, so the number of solutions is $\boxed{1}$.


text:
# including frame comments (as animation on beamer)

I would like to make my frames each time one appear, it would be joined with a comment "frame", how to create such thing on beamer!

Here is a figure about what I hope do! Any proposition or maybe some links help me to do so.

You could simply use columns to place the blocks and notes beside each other - there you can also specify the overlay using <1-> etc.

You should already know how to draw clouds from your previous question, so I leave this part as exercise for you.

\documentclass{beamer}

\begin{document}

\begin{frame}
\begin{columns}
\begin{column}<1->{.6\textwidth}
\begin{block}{title}
content
\end{block}
\end{column}
\begin{column}<1->{.3\textwidth}
Note
\end{column}
\end{columns}
\begin{columns}
\begin{column}<2->{.6\textwidth}
content
\end{column}
\begin{column}<2->{.3\textwidth}
Note
\end{column}
\end{columns}
\begin{columns}
\begin{column}<3->{.6\textwidth}
\begin{exampleblock}{title}
content
\end{exampleblock}
\end{column}
\begin{column}<3->{.3\textwidth}
Note
\end{column}
\end{columns}
\end{frame}

\end{document}


text:
# Introduction

There are countless examples in the literature where deep learning is used to make predictions using protein sequences. In order to use protein sequences as inputs for deep learning models they need to be transformed into numerical descriptors. This can be performed in several different ways depending on the underlying data and the purpose of the model. In this notebook I highlight four different ways to generate such numerical features from protein sequences using; 1) the Amino Acid composition, 2) generation of N-grams, 3) creation of word embeddings and 4) BLOSUM matrices. To illustrate how to use these descriptors I use them to train a deep learning model to predict protein classifications. The models themselves were kept as basic as possible and although reaching good accuracy were not fully optimized.

# Libraries

The below libraries were used to process the data and build the models.

# import the necessary packages
from collections import Counter
from keras.preprocessing.text import Tokenizer
from sklearn.model_selection import train_test_split, KFold
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import MinMaxScaler, LabelBinarizer, OneHotEncoder
from sklearn.feature_extraction.text import CountVectorizer
from sklearn.pipeline import Pipeline
from sklearn import preprocessing
from tensorflow.keras.callbacks import EarlyStopping
from tensorflow.keras.models import Model
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import Dense, Activation, Dropout, Embedding, Flatten, Concatenate, Input, Conv1D, MaxPooling1D, LeakyReLU, LSTM, BatchNormalization
from keras.losses import mean_squared_error
from sklearn.metrics import r2_score, accuracy_score
from keras import backend as K
from keras.utils.np_utils import to_categorical
from sklearn import linear_model
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import cv2
np.random.seed(1337) # for reproducibility

# Data

The dataset was compiled form the Protein Data Bank and posted by SHAHIR on kaggle for people to use to test different modeling approaches. The PDB database is used to visualize and store three-dimensional structural data of proteins. Each structure comes with a comprehensive list of all the experimental details that was used to generate the structure and the characteristics of the protein. It also classifies the protein into a family based on structural and sequence similarity with other proteins.

## Processing

Before using the data several steps were performed to make it more presentable for our deep learning model which is outlined below.

The data consists of two separate files. One file containing all the details about the structure and a second file containing the details on the protein sequence. Both files are loaded and merged together based on Structure ID. The file also contains structural data from DNA and RNA. These will be removed as we are only interested in the protein structures. We are only interested in the protein IDs and the sequences so we drop all other columns.

pdb_data_no_dups = pd.read_csv("Data/pdb_data_no_dups.csv")
#pdb_data_no_dups = pdb_data_no_dups[[""]]
#Merge data on structureId and select only protein structures
data_merged = pdb_data_no_dups.merge(pdb_data_seq, how='inner', on='structureId').query('macromoleculeType_x == "Protein"')
#Select columns we are interested in and remove duplicates
data_merged = data_merged[["structureId", "classification", "sequence"]]
data_merged = data_merged.drop_duplicates(["structureId", "classification", "sequence"])
data_merged = data_merged.reset_index(drop=True)
#Print data
'The data has {} rows and {} columns'.format(data_merged.shape[0], data_merged.shape[1])

## 'The data has 166681 rows and 3 columns'

pd.set_option('display.max_columns', None)
data_merged

##        structureId                            classification  \
## 0             101M                          OXYGEN TRANSPORT
## 1             102L                     HYDROLASE(O-GLYCOSYL)
## 2             102M                          OXYGEN TRANSPORT
## 3             103L                     HYDROLASE(O-GLYCOSYL)
## 4             103M                          OXYGEN TRANSPORT
## ...            ...                                       ...
## 166676        9RSA            HYDROLASE (PHOSPHORIC DIESTER)
## 166677        9RUB                      LYASE(CARBON-CARBON)
## 166678        9WGA                       LECTIN (AGGLUTININ)
## 166679        9XIA  ISOMERASE(INTRAMOLECULAR OXIDOREDUCTASE)
## 166680        9XIM  ISOMERASE(INTRAMOLECULAR OXIDOREDUCTASE)
##
##                                                  sequence
## 1       MNIFEMLRIDEGLRLKIYKDTEGYYTIGIGHLLTKSPSLNAAAKSE...
## 3       MNIFEMLRIDEGLRLKIYKDTEGYYTIGIGHLLTKSPSLNSLDAAK...
## ...                                                   ...
## 166676  KETAAAKFERQHMDSSTSAASSSNYCNQMMKSRNLTKDRCKPVNTF...
## 166677  MDQSSRYVNLALKEEDLIAGGEHVLCAYIMKPKAGYGYVATAAHFA...
## 166678  ERCGEQGSNMECPNNLCCSQYGYCGMGGDYCGKGCQNGACWTSKRC...
## 166679  MNYQPTPEDRFTFGLWTVGWQGRDPFGDATRRALDPVESVQRLAEL...
## 166680  SVQATREDKFSFGLWTVGWQARDAFGDATRTALDPVEAVHKLAEIG...
##
## [166681 rows x 3 columns]

Rows with missing values are removed

data_merged.isna().sum()

## structureId       0
## classification    1
## sequence          3
## dtype: int64

data_merged = data_merged.dropna()

Having a closer look at the data one can see that some of the sequences contain AAs other then the 21 standard ones (B|O|U|X|Z).

#Count all different AAs
sequences = data_merged["sequence"].dropna()
charCounter = CountVectorizer(analyzer="char", ngram_range=(1,1))
charVec = charCounter.fit_transform(sequences.to_numpy())
AA = charCounter.get_feature_names_out()
AAcounts = np.asarray(charVec.sum(axis=0))[0]
print(dict(zip(AA, AAcounts)))

## {'a': 3378009, 'b': 19, 'c': 642515, 'd': 2484519, 'e': 2835673, 'f': 1746745, 'g': 3211970, 'h': 1176659, 'i': 2389749, 'k': 2517051, 'l': 3940616, 'm': 1001217, 'n': 1857711, 'o': 2, 'p': 2048893, 'q': 1659623, 'r': 2136127, 's': 2837198, 't': 2460573, 'u': 50, 'v': 2990099, 'w': 628670, 'x': 177570, 'y': 1546907, 'z': 28}

For the sake of simplicity I decided to remove all those sequences.

#Remove sequences that contain non standard AAs
AA_mask = sequences.str.contains('^[RHKDESTNQCGPAILMFWYV]+') data_merged = data_merged.loc[AA_mask] '{} rows were removed'.format((AA_mask == False).sum())  ## '4377 rows were removed'  If we look at the size distribution of the sequences we can see that quite a large proportion consists of sequences shorter then 20 AAs. seqLengths = data_merged["sequence"].str.len() lengths = Counter(seqLengths.to_list()) plt.clf() # visualize _ = plt.hist(seqLengths.to_list(), bins=120, edgecolor='black') _ = plt.xlim(0, 1200) plt.show()  These are most likely from proteins for which only part of the structure has been solved and it will be hard to predict an accurate classification for those. To get more accurate models I removed all sequences that have a length with less then 20 AAs. #Filter out sequences that are shorter then 20 AAs AA_length_mask = data_merged["sequence"].str.len() >= 20 data_merged = data_merged.loc[AA_length_mask]  There are >4000 different classifications defined. Most of them consisting of far to little sequences to train a model. data_merged.groupby(["classification"])["classification"] \ .count() \ .sort_values(ascending = False)  ## classification ## HYDROLASE 23465 ## TRANSFERASE 16557 ## OXIDOREDUCTASE 14288 ## IMMUNE SYSTEM 7920 ## HYDROLASE/HYDROLASE INHIBITOR 5661 ## ... ## NOVEL SEQUENCE 1 ## NONSTRUCTURAL PROTEIN 1 ## NON-HEME IRON PROTEIN 1 ## NK CELL 1 ## virus/inhibitor 1 ## Name: classification, Length: 4383, dtype: int64  To build the models I focus on the top 10 classifications. #Filter for top10 classifications classes = 10 data_merged['count'] = data_merged.groupby(["classification"])["classification"].transform('count') top_class = data_merged["count"].sort_values(ascending=False).unique()[classes] data_top_classes = data_merged[data_merged["count"] > top_class].sort_values(by=["count"], ascending=False) data_top_classes = data_top_classes.reset_index(drop=True) #Print results data_top_classes.groupby(["classification"])["classification"] \ .count() \ .sort_values(ascending = False)  ## classification ## HYDROLASE 23465 ## TRANSFERASE 16557 ## OXIDOREDUCTASE 14288 ## IMMUNE SYSTEM 7920 ## HYDROLASE/HYDROLASE INHIBITOR 5661 ## LYASE 4476 ## TRANSCRIPTION 4288 ## TRANSPORT PROTEIN 3751 ## SIGNALING PROTEIN 3284 ## VIRAL PROTEIN 2761 ## Name: classification, dtype: int64  ## Transformation To build the models we binarize the classes. See: link for an explanation on choosing the right encoding method. lb = LabelBinarizer() y = lb.fit_transform(data_top_classes["classification"]) y[1:10]  ## array([[1, 0, 0, 0, 0, 0, 0, 0, 0, 0], ## [1, 0, 0, 0, 0, 0, 0, 0, 0, 0], ## [1, 0, 0, 0, 0, 0, 0, 0, 0, 0], ## [1, 0, 0, 0, 0, 0, 0, 0, 0, 0], ## [1, 0, 0, 0, 0, 0, 0, 0, 0, 0], ## [1, 0, 0, 0, 0, 0, 0, 0, 0, 0], ## [1, 0, 0, 0, 0, 0, 0, 0, 0, 0], ## [1, 0, 0, 0, 0, 0, 0, 0, 0, 0], ## [1, 0, 0, 0, 0, 0, 0, 0, 0, 0]])  The sequences are transformed to a numpy array. X = data_top_classes["sequence"].to_numpy()  # Amino Acid composition model As a first most basic model we can use the Amino Acid (AA) composition of each sequence to train our model. To do this we count the number of AAs for each sequence. vectorizer = CountVectorizer(analyzer="char") X_vec = vectorizer.fit_transform(X) X_AAC = X_vec.toarray() feature_names = vectorizer.get_feature_names_out() pd.DataFrame(X_AAC, columns = feature_names)  ## a c d e f g h i k l m n p q r s t v \ ## 0 22 1 23 20 11 16 6 27 29 18 2 16 19 7 6 13 19 14 ## 1 19 4 20 15 8 18 11 17 18 24 8 9 13 4 11 18 17 10 ## 2 23 3 21 23 19 29 7 26 17 23 8 24 22 11 18 31 26 12 ## 3 23 3 21 22 19 29 7 26 17 23 8 24 22 12 18 31 26 12 ## 4 40 0 22 16 9 26 13 16 9 23 5 11 20 13 15 13 28 28 ## ... .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. .. ## 86446 13 4 4 3 4 18 0 5 5 11 2 4 6 8 5 16 8 9 ## 86447 20 5 12 13 18 25 6 10 7 30 4 21 29 17 12 20 30 22 ## 86448 12 2 5 5 5 15 1 2 4 7 3 4 3 8 8 13 9 7 ## 86449 16 5 20 9 20 25 7 10 6 28 3 19 28 20 14 22 20 22 ## 86450 53 3 28 28 17 32 6 30 35 45 7 28 24 19 13 19 24 25 ## ## w y ## 0 6 11 ## 1 1 9 ## 2 11 26 ## 3 11 26 ## 4 4 12 ## ... .. .. ## 86446 2 8 ## 86447 5 9 ## 86448 3 8 ## 86449 4 10 ## 86450 10 17 ## ## [86451 rows x 20 columns]  ## AAC Model For the AA model I use three dense layers with the last layer having 10 units which is equal to the amount of classes that we try to predict. I use early stopping based on the accuracy of the test data to avoid over fitting. To limit the training time I increased the batch size to 128. The model is wrapped into a function such that I can re-use the model for Cross Validation. def Model_AAC(X_train, y_train, X_test, y_test, classes = 10): #Define model model = Sequential() model.add(Dense(16, input_dim=X_train.shape[1], activation="relu")) model.add(Dense(8, activation="relu")) model.add(Dense(classes, activation="softmax")) #Define optimizer opt = Adam(learning_rate=0.01) early_stopping = EarlyStopping(min_delta=0.01, patience=5, monitor='val_accuracy', verbose=0, mode="max", restore_best_weights=True ) #Compile model model.compile(loss='categorical_crossentropy', optimizer=opt, metrics=['accuracy']) model.fit(x=X_train, y=y_train, validation_data=(X_test, y_test), epochs=100, callbacks=[early_stopping], verbose=0, batch_size=128) y_pred = model.predict(X_test) test_acc = accuracy_score(np.argmax(y_test, axis=1), np.argmax(y_pred, axis=1)) return test_acc  ## Valdidation To test the model we check its performance on the test data that we set apart and using 5-fold cross validation. ### Training data To asses the model the data is split into a train set and 20% of the data is held back to evaluate our model. To avoid bias I normalized the data using the MinMaxScaler() (see: link for an explanation on normalization of data). #Generate train and test data X_train, X_test, y_train, y_test = train_test_split(X_AAC, y, test_size=.2, random_state=42, shuffle=True) #Scale data scaler = MinMaxScaler() X_train = scaler.fit_transform(pd.DataFrame(X_train, columns = feature_names)) #Fit and transform data on training data X_test = scaler.transform(pd.DataFrame(X_test, columns = feature_names)) #Transform data using the training parameters for fitting. #Fit data to Model train_acc = Model_AAC(X_train, y_train, X_test, y_test) print("Accuracy-train = {:.3f}".format(train_acc))  ## Accuracy-train = 0.452  ### Cross Validation Estimate accuracy using 5-fold cross validation. #Generate X, y data ACC = [] kf = KFold(n_splits=5, random_state=42, shuffle=True) i = 1 for train_index, test_index in kf.split(X_AAC): X_train, X_test = X_AAC[train_index], X_AAC[test_index] y_train, y_test = y[train_index], y[test_index] scaler = MinMaxScaler() X_train = scaler.fit_transform(pd.DataFrame(X_train)) X_test = scaler.transform(pd.DataFrame(X_test)) pred = Model_AAC(X_train, y_train, X_test, y_test) print("Accuracy-fold-{} = {:.3f}".format(i, pred)) i = i+1 ACC += [pred]  ## Accuracy-fold-1 = 0.435 ## Accuracy-fold-2 = 0.448 ## Accuracy-fold-3 = 0.428 ## Accuracy-fold-4 = 0.443 ## Accuracy-fold-5 = 0.429  print("Accuracy-mean = {:.3f}".format(np.mean(ACC, dtype=object)))  ## Accuracy-mean = 0.437  ## Conclusion Using the AA composition alone we get an accuracy of around 0.4. This is not too bad considering we are trying to predict 10 different classes. This type of encoding can work well for shorter sequences like peptides. However, with this type of encoding any information on the order of the sequences is lost. Especially for longer sequences this is important as many classes are defined by specific domains within the sequence. # N-grams model To preserve some information on the order of AAs in the sequences we can count the occurrences of the most common AA combinations across all the sequences. In the below example we look at the 20 most common 3 letter AA combinations and then count the occurrence of each for all the sequences. vectorizer = CountVectorizer(analyzer="char", ngram_range=(3,3), max_features = 20) X_train_vec = vectorizer.fit_transform(X) feature_names = vectorizer.get_feature_names_out() #Print results pd.DataFrame(X_train_vec.toarray(), columns = feature_names)  ## aaa aag aal ala ale all avl eal ell gla gsg gss hhh laa \ ## 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 ## 1 0 0 0 0 0 0 0 1 0 0 0 0 1 0 ## 2 0 0 0 0 0 0 0 0 0 0 1 0 0 0 ## 3 0 0 0 0 0 0 0 0 0 0 1 0 0 0 ## 4 0 0 1 0 0 0 0 0 0 0 0 0 4 0 ## ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ## 86446 0 0 0 0 0 0 0 0 0 0 1 0 0 0 ## 86447 0 0 0 0 0 0 0 0 0 1 0 0 0 0 ## 86448 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ## 86449 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ## 86450 2 0 2 1 0 0 0 1 0 1 0 0 0 1 ## ## lke lla lle lll sgs ssg ## 0 0 0 0 0 1 0 ## 1 0 0 0 0 0 0 ## 2 0 0 0 0 0 0 ## 3 0 0 0 0 0 0 ## 4 0 1 0 0 0 0 ## ... ... ... ... ... ... ... ## 86446 0 0 0 0 1 0 ## 86447 0 0 0 0 0 1 ## 86448 0 0 0 0 0 0 ## 86449 0 0 0 0 0 0 ## 86450 0 1 0 0 0 0 ## ## [86451 rows x 20 columns]  We use the above principle to count the occurrence of AA combinations with a length of 1 to 3 AAs and use the top 8000 features as input for our model. Note that we split the data in a train and test set first. We then generate all the N-grams based on the training data. This makes sure that we don’t have data leakage when we estimate the accuracy on the test data. The occurrences are then counted for both the test and train data. As the final step we scale the values (see the AA composition model). #Generate train and test data X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.2, random_state=42, shuffle=True) #Generate N-grams and fit on train and test data vectorizer = CountVectorizer(analyzer="char", ngram_range=(1,3), max_features = 8000) X_train_vec = vectorizer.fit_transform(X_train) X_test_vec = vectorizer.transform(X_test) feature_names = vectorizer.get_feature_names_out() #Scale values scaler = MinMaxScaler() X_train_AAN_scaled = scaler.fit_transform(pd.DataFrame(X_train_vec.toarray(), columns = feature_names)) X_test_AAN_scaled = scaler.transform(pd.DataFrame(X_test_vec.toarray(), columns = feature_names)) #Use X_train parameters for fitting. #Print results pd.set_option('display.max_columns', 10) pd.DataFrame(X_train_AAN_scaled, columns = feature_names)  ## a aa aaa aac aad ... yys yyt yyv yyw \ ## 0 0.036585 0.000000 0.000000 0.0 0.0 ... 0.0 0.0 0.000000 0.0 ## 1 0.004878 0.000000 0.000000 0.0 0.0 ... 0.0 0.5 0.000000 0.0 ## 2 0.029268 0.013245 0.000000 0.0 0.2 ... 0.0 0.0 0.000000 0.0 ## 3 0.151220 0.119205 0.040000 0.0 0.0 ... 0.0 0.0 0.000000 0.0 ## 4 0.007317 0.000000 0.000000 0.0 0.0 ... 0.0 0.0 0.000000 0.0 ## ... ... ... ... ... ... ... ... ... ... ... ## 69155 0.080488 0.033113 0.006667 0.0 0.0 ... 0.0 0.5 0.000000 0.0 ## 69156 0.007317 0.000000 0.000000 0.0 0.0 ... 0.0 0.0 0.333333 0.0 ## 69157 0.119512 0.039735 0.000000 0.0 0.2 ... 0.0 0.0 0.000000 0.0 ## 69158 0.007317 0.000000 0.000000 0.0 0.0 ... 0.0 0.0 0.000000 0.0 ## 69159 0.078049 0.013245 0.000000 0.0 0.0 ... 0.0 0.0 0.000000 0.0 ## ## yyy ## 0 0.0 ## 1 0.0 ## 2 0.0 ## 3 0.0 ## 4 0.0 ## ... ... ## 69155 0.0 ## 69156 0.0 ## 69157 0.0 ## 69158 0.0 ## 69159 0.0 ## ## [69160 rows x 8000 columns]  ## Model For comparison we use the same set-up as the Amino Acid composition model. However, more layers and nodes will most likely be able to improve the accuracy since we have far more features. def Model_AAN(X_train, y_train, X_test, y_test, classes = 10): #Define model model = Sequential() model.add(Dense(16, input_dim=X_train.shape[1], activation="relu")) model.add(Dense(8, activation="relu")) model.add(Dense(classes, activation="softmax")) #Define optimizer opt = Adam(learning_rate=0.01) early_stopping = EarlyStopping(min_delta=0.01, patience=5, monitor='val_accuracy', verbose=0, mode="max", restore_best_weights=True ) #Compile model model.compile(loss='categorical_crossentropy', optimizer=opt, metrics=['accuracy']) model.fit(x=X_train, y=y_train, validation_data=(X_test, y_test), epochs=100, callbacks=[early_stopping], verbose=0, batch_size=128) y_pred = model.predict(X_test) test_acc = accuracy_score(np.argmax(y_test, axis=1), np.argmax(y_pred, axis=1)) return test_acc  ## Validation ### Training data The same as with the AA composition model the values are scaled and fit to the model #Fit data to Model train_acc = Model_AAN(X_train_AAN_scaled, y_train, X_test_AAN_scaled, y_test) print("Accuracy-train = {:.3f}".format(train_acc))  ## Accuracy-train = 0.806  ### Cross Validation #Generate X, y data ACC = [] kf = KFold(n_splits=5, random_state=42, shuffle=True) i = 1 for train_index, test_index in kf.split(X): X_train, X_test = X[train_index], X[test_index] y_train, y_test = y[train_index], y[test_index] #Generate N-grams vectorizer = CountVectorizer(analyzer="char", ngram_range=(1,3), max_features = 8000) X_train_vec = vectorizer.fit_transform(X_train) X_test_vec = vectorizer.transform(X_test) feature_names = vectorizer.get_feature_names_out() #Scale values scaler = MinMaxScaler() X_train_scaled = scaler.fit_transform(pd.DataFrame(X_train_vec.toarray(), columns = feature_names)) X_test_scaled = scaler.transform(pd.DataFrame(X_test_vec.toarray(), columns = feature_names)) #Use X_train parameters for fitting. pred = Model_AAN(X_train_scaled, y_train, X_test_scaled, y_test) print("Accuracy-fold-{} = {:.3f}".format(i, pred)) i = i+1 ACC += [pred]  ## Accuracy-fold-1 = 0.807 ## Accuracy-fold-2 = 0.802 ## Accuracy-fold-3 = 0.811 ## Accuracy-fold-4 = 0.794 ## Accuracy-fold-5 = 0.804  print("Accuracy-mean = {:.3f}".format(np.mean(ACC)))  ## Accuracy-mean = 0.804  ## Conclusion Using N-grams we get an accuracy of around 0.8. This is a very big improvement over using only the AA composition. The increase in accuracy is mostly because with this kind of encoding we are capturing the AA motifs that are most common among the sequences. We might be able to increase the accuracy even further by playing around with the length of the words and increasing the complexity of the model by adding more layers and nodes. # Embedding model For this encoding the sequences are first split into AAs and each AA is assigned a number (tokenized). #Tokenize sequences tk = Tokenizer(char_level=True) tk.fit_on_texts(X) seq_tok_all = tk.texts_to_sequences(X) seq_tok_all[1]  ## [18, 15, 10, 6, 1, 4, 12, 20, 1, 9, 7, 13, 16, 2, 16, 10, 1, 17, 7, 6, 7, 8, 3, 8, 4, 3, 4, 4, 7, 12, 5, 6, 2, 6, 12, 10, 10, 7, 5, 1, 9, 11, 5, 3, 11, 13, 1, 8, 16, 10, 1, 13, 8, 17, 17, 17, 16, 7, 17, 8, 3, 3, 13, 1, 6, 1, 9, 7, 11, 16, 3, 2, 9, 4, 10, 3, 5, 2, 18, 7, 9, 7, 11, 10, 12, 3, 10, 7, 18, 2, 1, 9, 7, 3, 7, 9, 19, 18, 14, 2, 3, 17, 6, 4, 17, 4, 18, 7, 8, 12, 3, 17, 8, 9, 3, 17, 10, 5, 1, 16, 14, 12, 3, 5, 11, 2, 10, 14, 8, 3, 7, 8, 18, 14, 5, 1, 5, 20, 3, 9, 1, 14, 6, 3, 8, 12, 9, 15, 18, 1, 2, 5, 1, 15, 9, 10, 8, 5, 1, 12, 7, 7, 8, 5, 10, 16, 20, 3, 17, 6, 16, 8, 1, 5, 13, 5, 9, 14, 2, 1, 5, 1, 6, 12, 13, 13, 6, 4, 1, 15, 5, 16, 2, 2, 17, 4, 2, 6, 1, 11, 5, 9, 9, 1, 12, 8, 10, 12, 8, 8, 4, 9, 18, 6, 9, 2, 20, 13, 12, 14, 1, 11, 5, 5, 13, 8, 7, 10, 11, 11, 2, 1, 11, 10, 12, 6, 2, 2, 7, 6, 2, 6, 2, 1, 3, 10, 10, 11, 9, 2, 9, 7, 7, 14]  The protein sequences have variable lengths while for deep learning models equal size sequences are required. For this reason shorter sequences need to be padded. The longest sequence in the data is around 5000 AAs long therefore by default all other sequences will be padded to the same length. This will result in a very large data set that will take very long to train. For this reason I decided to add a maximum length of 256 AAs. Sequences shorther then 256 AAs will be padded and sequences longer then 256 AAs will be trimmed. max_length = 256 #max_length = len(max(X, key=len)) X_vec = pad_sequences(seq_tok_all, maxlen=max_length, padding='post') #padding='post'  ## Model For this model a Convolutional kernel (Conv1D) is used as this type seemed to work better then a dense layer with this type of sequence representation. Before the data is fed into the model an Embedding layer is used which which transforms the tokenized sequences into a continues numeric vector of length 8. def Model_embed(X_train, y_train, X_test, y_test, max_length=max_length): #Define model model = Sequential() model.add(Embedding(input_dim=len(tk.word_index)+1, output_dim=8, input_length=max_length, trainable=True)) model.add(Conv1D(filters=64, kernel_size=6, padding="same", activation="relu")) model.add(Flatten()) model.add(Dense(10, activation="softmax")) opt = Adam(learning_rate=0.01) early_stopping = EarlyStopping(min_delta=0.01, patience=5, monitor='val_accuracy', verbose=0, mode="max", restore_best_weights=True ) #Compile model model.compile(loss='categorical_crossentropy', optimizer=opt, metrics=['accuracy']) model.fit(x=X_train, y=y_train, validation_data=(X_test, y_test), epochs=100, callbacks=[early_stopping], verbose=0, batch_size=128) y_pred = model.predict(X_test) test_acc = accuracy_score(np.argmax(y_test, axis=1), np.argmax(y_pred, axis=1)) return test_acc  ### Training data X_train, X_test, y_train, y_test = train_test_split(X_vec, y, test_size=.2, random_state=42, shuffle=True) train_acc_embed = Model_embed(X_train, y_train, X_test, y_test, max_length=max_length) print("Accuracy-train = {:.3f}".format(train_acc_embed))  ## Accuracy-train = 0.762  ### Cross Validation #Generate X, y data ACC = [] kf = KFold(n_splits=5, random_state=42, shuffle=True) i = 1 for train_index, test_index in kf.split(X_vec): X_train, X_test = X_vec[train_index], X_vec[test_index] y_train, y_test = y[train_index], y[test_index] pred = Model_embed(X_train, y_train, X_test, y_test, max_length=max_length) print("Accuracy-fold-{} = {:.3f}".format(i, pred)) i = i+1 ACC += [pred]  ## Accuracy-fold-1 = 0.749 ## Accuracy-fold-2 = 0.756 ## Accuracy-fold-3 = 0.757 ## Accuracy-fold-4 = 0.764 ## Accuracy-fold-5 = 0.758  print("Accuracy-mean = {:.3f}".format(np.mean(ACC)))  ## Accuracy-mean = 0.757  ## Conclusion The embedded model reaches an accuracy of around 0.75. This is much better then the AA composition model but still lower then the N-grams model. This is most likely because we loose some of the sequence information due to the padding and trimming that is required to get equal length vectors for our deep learning model. # BLOSUM encoding BLOSUM matrices are widely used in bioinformatics to generate sequence alignments. We can use this encoding also as an input for deep learning. BLOSUM encoding of the sequences was performed using the protr R package using the below code. library(protr) #Function for BLOSUM encoding of the AA sequences. BLOSUM_df <- function(sequences, submat="AABLOSUM62", k=5, lag=7, silent=TRUE){ df <- t(sapply(sequences, function(i) protr::extractBLOSUM(i, submat = submat, k=k, lag=lag, silent=silent))) df <- data.frame(df) rownames(df) <- c() return(df) } #Import data sequences <- pydata_top_classes["sequence"]
X_BLOSUM <-
BLOSUM_df(sequences\$sequence, submat="AABLOSUM62", k=9, lag=18, silent=TRUE) #9 / 18
#write_csv(X_BLOSUM, "BLOSUM_k9_lag19.csv")

BLOSUM sequence matrices generated in R were saved into a .csv file and loaded back into python

BLOSUM_seq_matrixes = pd.read_csv("Data/BLOSUM_k9_lag19.csv")
pd.set_option('display.max_columns', 6)
BLOSUM_seq_matrixes

##        scl1.lag1  scl2.lag1  scl3.lag1  ...  scl8.7.lag18  scl9.7.lag18  \
## 0      -0.002615   0.000283  -0.002937  ...     -0.003019      0.000230
## 1      -0.000948   0.000633   0.003090  ...      0.005362     -0.001396
## 2       0.005699   0.004398   0.003010  ...     -0.002283     -0.003664
## 3       0.005689   0.004419   0.002920  ...     -0.002295     -0.003665
## 4       0.002448   0.008009   0.000728  ...      0.000387      0.002231
## ...          ...        ...        ...  ...           ...           ...
## 86446  -0.004528   0.012028   0.008579  ...     -0.002316     -0.006078
## 86447  -0.002841   0.004932  -0.007188  ...      0.000123      0.001226
## 86448   0.004192   0.007186   0.023998  ...     -0.008741     -0.006935
## 86449  -0.002852   0.005154  -0.003875  ...     -0.000257      0.000208
## 86450  -0.001217   0.003631  -0.001252  ...     -0.000895      0.000950
##
##        scl9.8.lag18
## 0          0.000997
## 1          0.001682
## 2          0.003491
## 3          0.003506
## 4          0.004065
## ...             ...
## 86446      0.000662
## 86447      0.002007
## 86448      0.002186
## 86449      0.002532
## 86450      0.000996
##
## [86451 rows x 1458 columns]

X_BLOSUM = BLOSUM_seq_matrixes.to_numpy()

## Model

def Model_BLOSUM(X_train, y_train, X_test, y_test, classes = 10):
#Define model
model = Sequential()
#Define optimizer

early_stopping = EarlyStopping(min_delta=0.01,
patience=5,
monitor='val_accuracy',
verbose=0,
mode="max",
restore_best_weights=True
)
#Compile model
model.compile(loss='categorical_crossentropy', optimizer=opt, metrics=['accuracy'])

model.fit(x=X_train, y=y_train,
validation_data=(X_test, y_test),
epochs=100,
callbacks=[early_stopping],
verbose=0,
batch_size=128)

y_pred = model.predict(X_test)
test_acc = accuracy_score(np.argmax(y_test, axis=1), np.argmax(y_pred, axis=1))
return test_acc

## Validation

### Training data

X_train, X_test, y_train, y_test = train_test_split(X_BLOSUM, y, test_size=.2, random_state=42, shuffle=True)

train_acc_BLOSUM = Model_BLOSUM(X_train, y_train, X_test, y_test)

print("Accuracy-train = {:.3f}".format(train_acc_BLOSUM))

## Accuracy-train = 0.753

### Cross Validation

ACC = []
kf = KFold(n_splits=5, random_state=42, shuffle=True)
i = 1
for train_index, test_index in kf.split(X_BLOSUM):
X_train, X_test = X_BLOSUM[train_index], X_BLOSUM[test_index]
y_train, y_test = y[train_index], y[test_index]

pred = Model_BLOSUM(X_train, y_train, X_test, y_test)
print("Accuracy-fold-{} = {:.3f}".format(i, pred))
i = i+1

ACC += [pred]

## Accuracy-fold-1 = 0.777
## Accuracy-fold-2 = 0.760
## Accuracy-fold-3 = 0.776
## Accuracy-fold-4 = 0.757
## Accuracy-fold-5 = 0.766

print("Accuracy-mean = {:.3f}".format(np.mean(ACC)))

## Accuracy-mean = 0.767

## Conclusion

The BLOSUM model reaches an accuracy of around 0.75. This is similar to the N-grams model. However, generating the BLOSUM matrixes took much langer then generating the N-grams. Because there are so many features generated a higher accuracy might be reached by increasing layers and nodes.