Sunday, August 21, 2016

Geometric Algebra

The basic building block of information is a variable. Let M be the number of values (objects,elements) (extent) in one variable and N the number of values in another variable. This could be the number of unit lengths on the sides of a rectangle, but it could also be the number of different values in a column of a database table. A joint selection will have MN values to choose from, i.e. the extent of the joint variable is MN. That is the origin of multiplication. It is applied to rectangles in elementary school, but in physics education (i.e. description of the world) we later on encounter products of various kinds, like s = vt and p = mv and so on (s, t, m, v, p being extents of variables). The idea that allows us to do this is the abstract notion of information, and, specifically regarding multiplication, of the number of possibilities to jointly fix values of more variables. If s is the extent (length) of a line swept by a point and t is the extent of another line swept by the pointer on a clock, then v = s ⁄ t is a (reduced, relational) extent, that needs to be multiplied again by t to regain the original extent s. Physicists quite creatively introduce many helper variables of such kind and work with them using the rules of algebra of numbers. Over the history they reached certain obstacles that necessitated the addition of new numbers, of which the most unintuitive one is the addition of the imaginary numbers. But they are used all over physics. So to understand the complex numbers is crucial in order to understand the modern description of nature.

Note

Geometric algebra is a conceptually more intuitive approach to complex numbers and quaternions.

When one does not want to specify the number of combination of jointly fixing more variables, but only fix them to create a context, one uses vectors.

In an n dimensional space two vectors still only define a plane. The possible combinations swept by the vectors is the area of the parallelogram enclosed by the vectors. Let's reduce this to the multiplication - area of rectangle analogy by decomposing the vectors into orthogonal components along () and vertical () to an arbitrary direction in the same plane.

(a + a)(b + b) = ab + ab = ab − ba
  • The components in the same direction don't enclose any area and are 0 therefore
  • The last step shows that the two area terms have opposite direction. If we take first in both cases, i.e. give the area the same direction, then  − 1 encodes the direction.
  • Opposite direction is needed to make the product 0, if both a and b have the same direction, because in this case they would not enclose any area.

Isn't this the best analogy to the product of the sides of a rectangle? It is, but with the difference that a and b here are not perpendicular. When a and b point in the same direction, the product becomes 0. So this product contains both, the rectangle product and the relative direction, more specifically ab − ba = ab(cosθasinθb − cosθbsinθa) = |a||b|sin(θb − θa). This shows the anticommutativity of this product.

For two variables (area and relative direction) one needs two other variables to maintain the information. Therefore this outer product (ab) is accompanied by an inner product (or dot product) ab + ab = ab(cosθacosθb + sinθasinθb) = |a||b|cos(θb − θa).

Note

The outer product corresponds to the cross product usually taught at school. But look for duality for more.

Only both of these two complementary products determine |a||b| and θ. The geometric product is the combination of these products and serves as natural extension of the rectangle product to vectors and as such it is written without symbol. It results in a (scalar + area):

  • Parallel components are multiplied and added to give a scalar.
  • Vertical components are multiplied and subtracted to give a directed area.

So for two vectors we have

ab = ab + ab = (1)/(2)(ab + ba) + (1)/(2)(ab − ba)

The outer product makes a new kind of magnitude (area is different than length) and therefore points to a new direction not part of the other n variables forming the basis of an n-dimensional vector space. In general the vector space is extended by outer products of 2, 3, ..., n vectors, called bivector, trivector or generally grade k vectors or k blades.

Also a scalar is include, because it is linked via the scalar multiplication and it results from the inner product.

There are (nk) different k-vectors. For example in 3D-space there are 3 = (32) plane direction and thus bivectors.

There is only one n-vector, as there is only one scalar. Therefore this is called a pseudoscalar.

The scalar and the k-vectors from grade 1 to grade n make up all we need to describe a fixation of a set of variables and their relative orientations. This is called a multivector. Vector is for the combination of base variables. And multivector also includes the derived variables, i.e. the outer products. Because we have the 0, a base variable is also a vector and a vector is also a multivector.

Algebraic closure, i.e. having the result again in the same set, is important, because we can reuse the same operations or algorithms over again. Normally this is at the cost of losing information. It is an abstraction of information. This information then needs to be stored outside, which normally means less storage space, though. E.g. a number without a unit still needs to be made sense of via the separately stored units. It is tedious to carry along units in all calculations.

In geometric algebra closure is achieved via keeping the different results in a multivector, which stores the number for each grade separately. So for n base variables we have 2n component variables of the multivector (2n numbers) (E.g. 3D it is 8: α, a, b, c, ab, ac, bc, abc). To assert this closure we need to define how two general multivectors multiply and for that we still need to specify how to multiply grade combinations ArBs, not just two vectors. This can be done by reducing vectors to orthonormal components whose geometric product is eiej = δij + (1 − δij)eiej, i.e. we have either a scalar 1 or the wedge product. The result of ArBs depends on which ei each combination of the (ns) and (nr) components have in common and thus annihilate. Those who do not annihilate augment the dimension of the wedge product.

We can also use the following formula only valid for vector times r-vector, but which can be repeatedly applied.

aAr = aAr + aAr = (1)/(2)(aAr − ( − 1)rAra) + (1)/(2)(aAr + ( − 1)rAra)

Here the dot part reduces the grade by one

aAr = a⋅∧ak = ( − 1)k + 1aaki ≠ kai

and the wedge part increases it by one.

Application

The pseudoscalar unit I comprises all base vectors. I2 is -1 for n=2,3 and in general it transforms a r-vector into a (n-r)-vector. This is called a duality transformation. In 3D it can be used to go from a blade to the vector normal to the plane, which is the traditional cross product (v1 × v2 =  − I(v1v2))

For n=2: e1I = e1(e1e2) = e2 and Ie1 = (e1e2)e1 =  − e2. I = e1e2 is the imaginary unit. It makes one dimension into two dimensions: e1(a + bI) = ae1 + be2.

For n=3: e2e3 = i, e1e3 = j, e1e2 = k result in the quaternions with i2 = j2 = k2 = ijk =  − 1.

Reflection on a plane of unit normal vector n:  − nan =  − n(an + an)n =  − n(aI + a) = an − an

Rotation by 2θ is a composition of two reflections: mnanm = RaR. With (nm)2 = (nm − nm)(nm − mn) =  − sin2θ we have

R = nm = nm + nm = cosθ − (mn)/(sinθ)sinθ = e − Bθ

If unit vector a should be rotated to unit vector b, let n be the unit vector between a and b:

a’ = bnanb = RaR R = bn = b(a + b)/(|a + b|) = (1 + ba)/((2(1 + ab)))

These are just elementary applications to demonstrate the usefulness of geometric algebra. See [1] and [2] for a more comprehensive account.

Software

There is a python simpy implementation with a good documentation. Since Simpy version 1.0 it is separately on github. bookGA is a succinct book on geometric algebra.

Here an example:

import sys

from sympy import symbols,sin,cos,sqrt
from galgebra.ga import Ga

(o3d,ex,ey,ez) = Ga.build('e*x|y|z',g=[1,1,1])

(r,th,phi,alpha,beta,gamma) = symbols('r theta phi alpha beta gamma',real=True)
(x_a,y_a,z_a,x_b,y_b,z_b,ab_mag,th_ab) = symbols('x_a y_a z_a x_b y_b z_b ab_mag theta_ab',real=True)

I = ex^ey^ez

a = o3d.mv('a','vector')
b = o3d.mv('b','vector')
c = o3d.mv('c','vector')
ab = a-b

ab_norm = ab/ab_mag

R_ab = cos(th_ab/2) +I*ab_norm*cos(th_ab/2)
R_ab_rev = R_ab.rev()

e__ab_x = R_ab * ex * R_ab_rev
e__ab_y = R_ab * ey * R_ab_rev
e__ab_z = R_ab * ez * R_ab_rev

R_phi = cos(phi/2)-(ex^ey)*sin(phi/2)
R_phi_rev = R_phi.rev()

e_phi = (R_phi * ey * R_phi.rev())

R_th = cos(th/2)+I*e_phi*sin(th/2)
R_th_rev = R_th.rev()

e_r = (R_th*R_phi*ex*R_phi_rev*R_th_rev).trigsimp()
#cos(phi)*cos(theta)*e_x + sin(phi)*cos(theta)*e_y + sin(theta)*e_z
e_th = (R_th*R_phi*ez*R_phi_rev*R_th_rev).trigsimp()
#-sin(theta)*cos(phi)*e_x - sin(phi)*sin(theta)*e_y + cos(theta)*e_z
e_phi = e_phi.trigsimp()
#-sin(phi)*e_x + cos(phi)*e_y
[1]New Foundations of Classical Mechanics, David Hestenes
[2]Geometric Algebra for Physicists, C.J.L.Doran and A.N.Lasenby

Geometric Algebra

The basic building block of information is a variable. Let M be the number of values (objects,elements) (extent) in one variable and N the number of values in another variable. This could be the number of unit lengths on the sides of a rectangle, but it could also be the number of different values in a column of a database table. A joint selection will have MN values to choose from, i.e. the extent of the joint variable is MN. That is the origin of multiplication. It is applied to rectangles in elementary school, but in physics education (i.e. description of the world) we later on encounter products of various kinds, like s = vt and p = mv and so on (s, t, m, v, p being extents of variables). The idea that allows us to do this is the abstract notion of information, and, specifically regarding multiplication, of the number of possibilities to jointly fix values of more variables. If s is the extent (length) of a line swept by a point and t is the extent of another line swept by the pointer on a clock, then v = s ⁄ t is a (reduced, relational) extent, that needs to be multiplied again by t to regain the original extent s. Physicists quite creatively introduce many helper variables of such kind and work with them using the rules of algebra of numbers. Over the history they reached certain obstacles that necessitated the addition of new numbers, of which the most unintuitive one is the addition of the imaginary numbers. But they are used all over physics. So to understand the complex numbers is crucial in order to understand the modern description of nature.

Note

Geometric algebra is a conceptually more intuitive approach to complex numbers and quaternions.

When one does not want to specify the number of combination of jointly fixing more variables, but only fix them to create a context, one uses vectors.

In an n dimensional space two vectors still only define a plane. The possible combinations swept by the vectors is the area of the parallelogram enclosed by the vectors. Let's reduce this to the multiplication - area of rectangle analogy by decomposing the vectors into orthogonal components along () and vertical () to an arbitrary direction in the same plane.

(a + a)(b + b) = ab + ab = ab − ba
  • The components in the same direction don't enclose any area and are 0 therefore
  • The last step shows that the two area terms have opposite direction. If we take first in both cases, i.e. give the area the same direction, then  − 1 encodes the direction.
  • Opposite direction is needed to make the product 0, if both a and b have the same direction, because in this case they would not enclose any area.

Isn't this the best analogy to the product of the sides of a rectangle? It is, but with the difference that a and b here are not perpendicular. When a and b point in the same direction, the product becomes 0. So this product contains both, the rectangle product and the relative direction, more specifically ab − ba = ab(cosθasinθb − cosθbsinθa) = |a||b|sin(θb − θa). This shows the anticommutativity of this product.

For two variables (area and relative direction) one needs two other variables to maintain the information. Therefore this outer product (ab) is accompanied by an inner product (or dot product) ab + ab = ab(cosθacosθb + sinθasinθb) = |a||b|cos(θb − θa).

Note

The outer product corresponds to the cross product usually taught at school. But look for duality for more.

Only both of these two complementary products determine |a||b| and θ. The geometric product is the combination of these products and serves as natural extension of the rectangle product to vectors and as such it is written without symbol. It results in a (scalar + area):

  • Parallel components are multiplied and added to give a scalar.
  • Vertical components are multiplied and subtracted to give a directed area.

So for two vectors we have

ab = ab + ab = (1)/(2)(ab + ba) + (1)/(2)(ab − ba)

The outer product makes a new kind of magnitude (area is different than length) and therefore points to a new direction not part of the other n variables forming the basis of an n-dimensional vector space. In general the vector space is extended by outer products of 2, 3, ..., n vectors, called bivector, trivector or generally grade k vectors or k blades.

Also a scalar is include, because it is linked via the scalar multiplication and it results from the inner product.

There are (nk) different k-vectors. For example in 3D-space there are 3 = (32) plane direction and thus bivectors.

There is only one n-vector, as there is only one scalar. Therefore this is called a pseudoscalar.

The scalar and the k-vectors from grade 1 to grade n make up all we need to describe a fixation of a set of variables and their relative orientations. This is called a multivector. Vector is for the combination of base variables. And multivector also includes the derived variables, i.e. the outer products. Because we have the 0, a base variable is also a vector and a vector is also a multivector.

Algebraic closure, i.e. having the result again in the same set, is important, because we can reuse the same operations or algorithms over again. Normally this is at the cost of losing information. It is an abstraction of information. This information then needs to be stored outside, which normally means less storage space, though. E.g. a number without a unit still needs to be made sense of via the separately stored units. It is tedious to carry along units in all calculations.

In geometric algebra closure is achieved via keeping the different results in a multivector, which stores the number for each grade separately. So for n base variables we have 2n component variables of the multivector (2n numbers) (E.g. 3D it is 8: α, a, b, c, ab, ac, bc, abc). To assert this closure we need to define how two general multivectors multiply and for that we still need to specify how to multiply grade combinations ArBs, not just two vectors. This can be done by reducing vectors to orthonormal components whose geometric product is eiej = δij + (1 − δij)eiej, i.e. we have either a scalar 1 or the wedge product. The result of ArBs depends on which ei each combination of the (ns) and (nr) components have in common and thus annihilate. Those who do not annihilate augment the dimension of the wedge product.

We can also use the following formula only valid for vector times r-vector, but which can be repeatedly applied.

aAr = aAr + aAr = (1)/(2)(aAr − ( − 1)rAra) + (1)/(2)(aAr + ( − 1)rAra)

Here the dot part reduces the grade by one

aAr = a⋅∧ak = ( − 1)k + 1aaki ≠ kai

and the wedge part increases it by one.

Application

The pseudoscalar unit I comprises all base vectors. I2 is -1 for n=2,3 and in general it transforms a r-vector into a (n-r)-vector. This is called a duality transformation. In 3D it can be used to go from a blade to the vector normal to the plane, which is the traditional cross product (v1 × v2 =  − I(v1v2))

For n=2: e1I = e1(e1e2) = e2 and Ie1 = (e1e2)e1 =  − e2. I = e1e2 is the imaginary unit. It makes one dimension into two dimensions: e1(a + bI) = ae1 + be2.

For n=3: e2e3 = i, e1e3 = j, e1e2 = k result in the quaternions with i2 = j2 = k2 = ijk =  − 1.

Reflection on a plane of unit normal vector n:  − nan =  − n(an + an)n =  − n(aI + a) = an − an

Rotation by 2θ is a composition of two reflections: mnanm = RaR. With (nm)2 = (nm − nm)(nm − mn) =  − sin2θ we have

R = nm = nm + nm = cosθ − (mn)/(sinθ)sinθ = e − Bθ

If unit vector a should be rotated to unit vector b, let n be the unit vector between a and b:

a’ = bnanb = RaR R = bn = b(a + b)/(|a + b|) = (1 + ba)/((2(1 + ab)))

These are just elementary applications to demonstrate the usefulness of geometric algebra. See [1] and [2] for a more comprehensive account.

Software

There is a python simpy implementation with a good documentation. Since Simpy version 1.0 it is separately on github. bookGA is a succinct book on geometric algebra.

Here an example:

import sys

from sympy import symbols,sin,cos,sqrt
from galgebra.ga import Ga

(o3d,ex,ey,ez) = Ga.build('e*x|y|z',g=[1,1,1])

(r,th,phi,alpha,beta,gamma) = symbols('r theta phi alpha beta gamma',real=True)
(x_a,y_a,z_a,x_b,y_b,z_b,ab_mag,th_ab) = symbols('x_a y_a z_a x_b y_b z_b ab_mag theta_ab',real=True)

I = ex^ey^ez

a = o3d.mv('a','vector')
b = o3d.mv('b','vector')
c = o3d.mv('c','vector')
ab = a-b

ab_norm = ab/ab_mag

R_ab = cos(th_ab/2) +I*ab_norm*cos(th_ab/2)
R_ab_rev = R_ab.rev()

e__ab_x = R_ab * ex * R_ab_rev
e__ab_y = R_ab * ey * R_ab_rev
e__ab_z = R_ab * ez * R_ab_rev

R_phi = cos(phi/2)-(ex^ey)*sin(phi/2)
R_phi_rev = R_phi.rev()

e_phi = (R_phi * ey * R_phi.rev())

R_th = cos(th/2)+I*e_phi*sin(th/2)
R_th_rev = R_th.rev()

e_r = (R_th*R_phi*ex*R_phi_rev*R_th_rev).trigsimp()
#cos(phi)*cos(theta)*e_x + sin(phi)*cos(theta)*e_y + sin(theta)*e_z
e_th = (R_th*R_phi*ez*R_phi_rev*R_th_rev).trigsimp()
#-sin(theta)*cos(phi)*e_x - sin(phi)*sin(theta)*e_y + cos(theta)*e_z
e_phi = e_phi.trigsimp()
#-sin(phi)*e_x + cos(phi)*e_y
[1]New Foundations of Classical Mechanics, David Hestenes
[2]Geometric Algebra for Physicists, C.J.L.Doran and A.N.Lasenby

Thursday, August 18, 2016

NI DAQmx

Data Acquisition with the NI USB-6212

Data Acquisition with the NI USB-6212

There is a C API.

Help for the C API is installed in:

C:\Program Files (x86)\National Instruments\NI-DAQ\Docs\cdaqmx.chm

But the usage is best learned via the examples in:

C:\Users\public\Documents\National Instruments\NI-DAQ\Examples\DAQmx ANSI C

I worked with the PyDAQmx python wrapping. The wrapper's function names are the C ones without the DAQmx prefix. So the same help as for the ANSI C functions can be used. With this python wrapper the tests are a good starting point.

Here are some tests of mine.

Write to analogue out and read back in at the same time

Physically connect ao0 with ai0. The code writes a signal to ao0 and reads it back in via ai0. Dev2 is the device as shown in NI MAX (Measurement & Automation Explorer).

import ctypes
import numpy as np
import time
import matplotlib.pyplot as plt
from PyDAQmx import *

DEV = "Dev2"
f_s = 1000

tAI = Task()
tAI.CreateAIVoltageChan(DEV+"/ai0","Voltage",DAQmx_Val_RSE,-10,10,DAQmx_Val_Volts,None)
tAI.CfgSampClkTiming("",f_s,DAQmx_Val_Rising,DAQmx_Val_FiniteSamps,f_s)

written = int32()
dataO = np.array([sin(24.0*v/f_s)*f_s/v/24.0 for v in range(1,f_s+1)])
tAO = Task()
tAO.CreateAOVoltageChan(DEV+"/ao0","",-10.0,10.0,DAQmx_Val_Volts,None)
tAO.CfgSampClkTiming("",f_s,DAQmx_Val_Rising,DAQmx_Val_ContSamps,f_s)
tAO.WriteAnalogF64(f_s,0,10.0,DAQmx_Val_GroupByChannel,dataO,None,None)
tAO.StartTask()

dataI = np.zeros(f_s)
dataI = np.zeros(f_s)
readI = int32()
tAI.ReadAnalogF64(f_s,10.0,DAQmx_Val_GroupByChannel,dataI,f_s,byref(readI),None)
tAI.StartTask()
time.sleep(1)
tAI.StopTask()
tAI.ClearTask()

tAO.StopTask()
tAO.ClearTask()

plt.plot(dataI)
plt.show()

Reading more channels at the same time

I tried to create two tasks, each with an input channel. When starting them I got the error:

The specified resource is reserved.

Making the two channels part of the same task works.

import ctypes
import numpy as np
import time
import matplotlib.pyplot as plt
from PyDAQmx import *

DEV = "Dev2"
f_s = 1000
t = 1#s
V_MAX = 10
SHUNT = 100
I_MAX = 0.1
timeout = 10.0
Nch = 2

tsk = Task()
tsk.CreateAIVoltageChan(DEV+"/ai0","Voltage",DAQmx_Val_RSE,-V_MAX,V_MAX,DAQmx_Val_Volts,None)
# Taking the previous or the next line away makes the commented lines below work
tsk.CreateAICurrentChan(DEV+"/ai1","Current",DAQmx_Val_Cfg_Default,-I_MAX,I_MAX,DAQmx_Val_Amps,DAQmx_Val_Default,SHUNT,None) #

Ns = t*f_s
tsk.CfgSampClkTiming("",f_s,DAQmx_Val_Rising,DAQmx_Val_FiniteSamps,Ns)

# #A third output channel on a separate task does not work (1)
# written = int32()
# dataAO = np.array([sin(24.0*v/f_s)*f_s/v/24.0 for v in range(1,f_s+1)])
# tAO = Task()
# tAO.CreateAOVoltageChan(DEV+"/ao0","",-10.0,10.0,DAQmx_Val_Volts,None)
# tAO.CfgSampClkTiming("",f_s,DAQmx_Val_Rising,DAQmx_Val_ContSamps,f_s)
# tAO.WriteAnalogF64(f_s,0,10.0,DAQmx_Val_GroupByChannel,dataAO,None,None)
# tAO.StartTask()

read = int32()
data = np.zeros(Nch*Ns)
tsk.ReadAnalogF64(Ns,timeout,DAQmx_Val_GroupByChannel,data,Nch*Ns,byref(read),None)
tsk.StartTask()
time.sleep(1)

# #A third output channel on a separate task does not work (2)
#tAO.StopTask()
#tAO.ClearTask()

tsk.StopTask()

data.shape = (Ns,Nch)
plt.plot(data)
plt.show()

Put StartTask() after WriteAnalogF64(). ReadAnalogF64() will work properly also without calling StartTask(), but by calling it one keeps control of starting and stopping.

DAQmx_Val_GroupByChannel makes blocks by channel in the buffer. data.shape = (Ns,Nch) takes account of that.

I used the above samples (for t > 1) to verify that a frequency is there in a rather noise signal.

def has_frequency(data, f, f_s):
"""data is sampled over a multiple of the period, e.g. t = int(10/f+1)

The code uses heuristic parameters for a special context.
"""
    data = data.transpose()[0]
    N = data.size
    ssf = sps.medfilt(data,29) #filter first
    amplitudes = np.abs(np.fft.rfft(ssf))*2/N #then do an FFT
    f_idx = int(f*N/f_s)
    allthere = []
    #rectangular signal
    #check that in several modes the same frequency is present
    for m in range(1,4):
        Nwindow = 9
        rng = np.arange(-Nwindow,Nwindow,dtype=np.int)
        amps = amplitudes[m*f_idx+rng]
        imax_in_rng = m*f_idx + rng[np.argmax(amps)]
        there = abs(m*f_idx - imax_in_rng) <= 1 #allow for f tolerance
        allthere.append(there)
    return all(allthere)

Writing and reading from digital lines

StartTask() is not needed, as it is called automatically.

class digital_in_out:
    def __init__(self):
        self.tDO = Task()
        self.DO_names = "FootEnable,HandAbort,HandStart,DI,RC1,RC2,RC3,RC4".split(',')
        self.tDO.CreateDOChan(DEV+"/port0/line0:7",self.DO_names,DAQmx_Val_ChanForAllLines)
        self.do = np.array([0]*8,np.uint8)
        self.tDI = Task()
        self.DI_names = "DO1,DO2,Stim,Ready,Leadoff,Current,Power".split(',')
        self.tDI.CreateDIChan(DEV+"/port2/line0:6",,DAQmx_Val_ChanForAllLines)
        self.di = np.array([0]*7,np.uint8)
    def __end__(self):
        self.tDO.StopTask()
        self.tDO.ClearTask()
        self.tDI.StopTask()
        self.tDI.ClearTask()
    def set(self, what):
        for w in what:
            self.do[self.DO_names.index(w)] = 1
        self.tDO.WriteDigitalLines(1,1,10.0,DAQmx_Val_GroupByChannel,self.do,None,None);
    def clear(self, what):
        for w in what:
            self.do[self.DO_names.index(w)] = 0
        self.tDO.WriteDigitalLines(1,1,10.0,DAQmx_Val_GroupByChannel,self.do,None,None);
    def press(self, what):
        self.set(what)
        time.sleep(0.2)
        self.clear(what)
    def read(self, what):
        readDI = int32()
        sampleDI = int32()
        self.tDI.ReadDigitalLines(1,10.0,DAQmx_Val_GroupByChannel,self.di,len(self.DI_names),byref(sampleDI),byref(readDI),None)
        return int(self.di[self.DI_names.index(what)])

dio = digital_in_out()
dio.set("RC1 RC2".split())
dio.press(["HandStart"])
dio.read("Stim")
#...