Skip to content

Instantly share code, notes, and snippets.

@mathdoodle
Last active October 15, 2017 17:56
Show Gist options
  • Save mathdoodle/7fa42ff59204b222b971a2f28d49c811 to your computer and use it in GitHub Desktop.
Save mathdoodle/7fa42ff59204b222b971a2f28d49c811 to your computer and use it in GitHub Desktop.
Geometric Computer Algebra

Geometric Computer Algebra

Overview

Demonstrates use of the GeoCAS library. GeoCAS provides a multivector type that is parameterized by the field type and supports arbitrary metrics.

By having a parameterized field, the Multivector<T> type can be used for different kinds of numeric and symbolic computation. The most simple parameterization is to use number for the T parameter. This requires using an appropriate adapter for the number field. GeoCAS conveniently provides the NumberFieldAdapter class for this purpose.

Here is an example of how to initialize an algebra for Euclidean 3D space:

const G3 = GeoCAS.algebra([1, 1, 1], new GeoCAS.NumberFieldAdapter())

export const zero = G3.zero;
export const one = G3.one;
export const e1 = G3.unit(0);
export const e2 = G3.unit(1);
export const e3 = G3.unit(2);

If you would like to customize the labels for the basis vectors to be different from the standard 'e1', 'e2', ..., then provide the labels as the last parameter in the algebra generation function:

const G3 = GeoCAS.algebra([1, 1, 1], new GeoCAS.NumberFieldAdapter(), ['i','j','k'])

GeoCAS provides several other field adapters (FieldAdapter<T> type)

import { Geometric3 as G3 } from 'davinci-newton'
// Create shortcuts for some values.
export const e1 = G3.e1
export const e2 = G3.e2
export const e3 = G3.e3
export const meter = G3.meter
export const kilogram = G3.kilogram
export const second = G3.second
export const coulomb = G3.coulomb
export const ampere = G3.ampere
export const kelvin = G3.kelvin
export const mole = G3.mole
export const candela = G3.candela
export const newton = meter * kilogram / (second * second)
export const joule = newton * meter
export const volt = joule / coulomb
/**
* Print the HTML string without a line ending.
*/
export function print(html: string): void {
const element = document.getElementById('info');
if (element) {
element.innerHTML = element.innerHTML + html;
}
}
/**
* Print the HTML string and go to the next line.
*/
export function println(html: string): void {
print(html + '<br>');
}
/**
* Print the value of a variable.
*/
export function printvar(name: string, value: any): void {
println('<b>' + name + '</b> => ' + value);
}
const G = GeoCAS.algebra([+1, -1], new GeoCAS.ComplexFieldAdapter(), ['e_1','e_2'])
export const zero = G.zero;
export const one = G.one;
export const e1 = G.unit(0);
export const e2 = G.unit(1);
const G3 = GeoCAS.algebra([1, 1, 1], new GeoCAS.NumberFieldAdapter(), ['right','up','back'])
export const zero = G3.zero;
export const one = G3.one;
export const e1 = G3.unit(0);
export const e2 = G3.unit(1);
export const e3 = G3.unit(2);
<!doctype html>
<html>
<head>
<style>
/* STYLE-MARKER */
</style>
<script src="https://jspm.io/system.js"></script>
<script type="text/x-mathjax-config">
<!-- I much prefer using the single dollar-sign for inline mathematics -->
MathJax.Hub.Config({tex2jax: {inlineMath: [['$','$'], ['\\(','\\)']]}});
</script>
<script type="text/javascript" async src="https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS_CHTML"></script>
<!-- SCRIPTS-MARKER -->
</head>
<body>
<p id='info'></p>
<script>
// CODE-MARKER
</script>
<script>
System.import('./script.js')
</script>
</body>
</html>
export function abs<T extends {__abs__: () => T}>(arg: T): T {
return arg.__abs__();
}
import orthoFramesToVersor from './orthoFramesToVersor';
export default function() {
describe("G2", function() {
const G2 = GeoCAS.algebra([1, 1], new GeoCAS.ComplexFieldAdapter())
const one = G2.one
const e1 = G2.unit(0)
const e2 = G2.unit(1)
const I = e1 ^ e2
const A = [e1, e2]
const B = [e2, e1]
const V = orthoFramesToVersor(A, B, G2)
it("[e1, e2] => [e2, e1]", function() {
expect(coord(V, one)).toBe(0)
expect(coord(V, e1)).toBeCloseTo(+0.7071, 4)
expect(coord(V, e2)).toBeCloseTo(-0.7071, 4)
expect(coord(V, I)).toBe(0)
})
})
describe("G4", function() {
const G4 = GeoCAS.algebra([1, 1, 1, 1], new GeoCAS.NumberFieldAdapter())
// const one = G4.one
const [e1, e2, e3, e4] = G4.units
// const I = e1 ^ e2
const A = [e1, e2, e3, e4]
const B = [e3, e4, e1, e2]
const V = orthoFramesToVersor(A, B, G4)
it("[e1, e2, e3, e4] => [e3, e4, e1, e2]", function() {
expect(V.extractGrade(0).isZero()).toBeTruthy()
expect(V.extractGrade(1).isZero()).toBeTruthy()
expect(V.extractGrade(2).isZero()).toBeFalsy()
expect(V.extractGrade(2).toString()).toBe("-0.4999999999999999 * e1 ^ e2 - 0.4999999999999999 * e2 ^ e3 + 0.4999999999999999 * e1 ^ e4 - 0.4999999999999999 * e3 ^ e4")
expect(V.scp(~(e1 ^ e2)).scalarCoordinate()).toBeCloseTo(-0.5, 15)
expect(V.scp(~(e2 ^ e3)).scalarCoordinate()).toBeCloseTo(-0.5, 15)
expect(V.scp(~(e1 ^ e4)).scalarCoordinate()).toBeCloseTo(+0.5, 15)
expect(V.scp(~(e3 ^ e4)).scalarCoordinate()).toBeCloseTo(-0.5, 15)
expect(V.extractGrade(3).isZero()).toBeTruthy()
expect(V.extractGrade(4).isZero()).toBeTruthy()
})
})
}
/**
* Helper function to project out the (real) part of the coordinate corresponding to a blade.
* This may be used with the ComplexFieldAdapter.
*/
function coord(X: GeoCAS.Multivector<GeoCAS.Complex>, B: GeoCAS.Multivector<GeoCAS.Complex>): number {
return X.scp(B).scalarCoordinate().x;
}
/**
* Create a type alias to make the algorithm a bit neater.
*/
export type Multivector<T> = GeoCAS.Multivector<T>;
export type Algebra<T> = GeoCAS.Algebra<T>;
const cos = GeoCAS.cosineOfAngleBetweenBlades;
// const squaredNorm = GeoCAS.squaredNorm;
/**
* Determines a versor for two frames related by an orthogonal transform.
*
* A and B are lists of corresponding vectors, representing the frames.
* A is the start frame, B is the final frame.
*
* This is an application of the Cartan-Dieudonne theorem that states:
* In n-D, any orthogonal transformation can be represented as at most n planar reflections.
*/
export default function orthoFramesToVersor<T>(A: Multivector<T>[], B: Multivector<T>[], G: Algebra<T>): Multivector<T> {
return orthoFramesToVersorRecursive(A, B, G.one, G);
}
function orthoFramesToVersorRecursive<T>(A: Multivector<T>[], B: Multivector<T>[], V: Multivector<T>, G: Algebra<T>): Multivector<T> {
if (A.length > 0) {
const j = bestIndex(A, B);
const a = A[j];
const b = B[j];
if (abs(cos(a, b) - G.one) < G.ε) {
return orthoFramesToVersorRecursive(removeAt(A, j), removeAt(B, j), V, G);
}
else {
const e = a - b;
// Don't let irrational number rounding errors from sqrt propagate...
const e2 = e | e
const R = e.direction();
return orthoFramesToVersorRecursive(removeAt(A, j).map(x => -e * x * e / e2), removeAt(B, j), R * V, G);
}
}
else {
return V;
}
}
function abs<T>(x: Multivector<T>): Multivector<T> {
return x.__abs__();
}
/**
* Returns a copy of the array with the element at index removed.
*/
function removeAt<T>(xs: T[], index: number): T[] {
const result: T[] = [];
for (let i = 0; i < xs.length; i++) {
if (i !== index) {
result.push(xs[i]);
}
}
return result;
}
/**
* Determine the best vector for numerical stability.
*/
function bestIndex<T>(A: Multivector<T>[], B: Multivector<T>[]): number {
const N = A.length;
let max: Multivector<T> | undefined = void 0;
let idx = 0;
for (let k = 0; k < N; k++) {
const x = A[k].sub(B[k]);
const squaredNorm = x.scp(x.rev());
if (typeof max === 'undefined' || squaredNorm > max) {
idx = k;
max = squaredNorm;
}
}
return idx;
}
{
"description": "Geometric Computer Algebra",
"dependencies": {
"DomReady": "1.0.0",
"geocas": "1.13.0",
"jasmine": "2.5.2",
"davinci-newton": "0.0.43"
},
"operatorOverloading": true,
"name": "copy-of-units-of measure scratchpad",
"version": "0.1.0",
"keywords": [
"GeoCAS",
"Multivector"
],
"author": "David Geo Holmes"
}
import { zero, one, e1, e2 } from './G11';
import { println } from './extras';
import { abs } from './math';
function show(lhs: string, rhs: any): void {
println(`$$${lhs} = ${rhs}$$`)
}
function direction(arg: GeoCAS.Multivector<any>) {
show(`\\hat{${arg}}`, arg.direction())
}
function magnitude(arg: GeoCAS.Multivector<any>) {
show(`|${arg}|`, abs(arg))
}
/**
* Executed when the DOM is ready...
*/
function main() {
show("zero", zero)
show("one", one)
show("e_1", e1)
direction(e1)
magnitude(e1)
direction(e2)
magnitude(e2)
println(`$$|\\hat{e}_2|\\hat{e}_2 = ${abs(e2) * e2.direction()}$$`)
}
function colorize(arg: any, color: string) {
return "<span style='color:" + color + "'>" + arg + "</span>";
}
// Wait for the DOM to be loaded.
function runMe() {
try {
main();
}
catch (e) {
println(colorize(e.message, 'red'));
}
}
//
DomReady.ready(runMe);
body {
background-color: #000000;
}
#info {
position: absolute;
left: 20px;
top: 20px;
font-size: 26px;
color: #00FF00;
}
<!DOCTYPE html>
<html>
<head>
<!-- STYLES-MARKER -->
<style>
/* STYLE-MARKER */
</style>
<script src='https://jspm.io/system.js'></script>
<!-- SCRIPTS-MARKER -->
</head>
<body>
<script>
// CODE-MARKER
</script>
<script>
System.import('./tests.js')
</script>
</body>
</html>
import orthoFramesToVersor from './orthoFramesToVersor.spec'
window['jasmine'] = jasmineRequire.core(jasmineRequire)
jasmineRequire.html(window['jasmine'])
const env = jasmine.getEnv()
const jasmineInterface = jasmineRequire.interface(window['jasmine'], env)
extend(window, jasmineInterface)
const htmlReporter = new jasmine.HtmlReporter({
env: env,
getContainer: function() { return document.body },
createElement: function() { return document.createElement.apply(document, arguments) },
createTextNode: function() { return document.createTextNode.apply(document, arguments) },
timer: new jasmine.Timer()
})
env.addReporter(htmlReporter)
DomReady.ready(function() {
htmlReporter.initialize()
describe("orthoFramesToVersor", orthoFramesToVersor)
env.execute()
})
/*
* Helper function for extending the properties on objects.
*/
export default function extend<T>(destination: T, source: any): T {
for (let property in source) {
destination[property] = source[property]
}
return destination
}
{
"allowJs": true,
"checkJs": true,
"declaration": true,
"emitDecoratorMetadata": true,
"experimentalDecorators": true,
"jsx": "react",
"module": "system",
"noImplicitAny": true,
"noImplicitReturns": true,
"noImplicitThis": true,
"noUnusedLocals": true,
"noUnusedParameters": true,
"preserveConstEnums": true,
"removeComments": false,
"skipLibCheck": true,
"sourceMap": true,
"strictNullChecks": true,
"suppressImplicitAnyIndexErrors": true,
"target": "es5",
"traceResolution": true
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment