|
| 1 | +#!/usr/bin/env python3 |
| 2 | +# -*- coding: utf-8 -*- |
| 3 | +""" |
| 4 | +Micro-benchmarks for spatialmath: base functions and classes for SO(3), |
| 5 | +SE(3), quaternions and twists, plus a NumPy baseline. Run from a checkout:: |
| 6 | +
|
| 7 | + python benchmarks/benchmark_smtb.py |
| 8 | +
|
| 9 | +The output starts with a summary of the machine and package versions, so a |
| 10 | +pasted table is self-describing. |
| 11 | +
|
| 12 | +Created on Fri Apr 10 14:22:36 2020 |
| 13 | +
|
| 14 | +@author: Peter Corke |
| 15 | +""" |
| 16 | + |
| 17 | +import os |
| 18 | +import platform |
| 19 | +import subprocess |
| 20 | +import timeit as _timeit |
| 21 | + |
| 22 | +import numpy as np |
| 23 | +from ansitable import ANSITable, Column |
| 24 | + |
| 25 | +import spatialmath |
| 26 | + |
| 27 | +N = 10_000 |
| 28 | +REPEATS = 5 |
| 29 | + |
| 30 | +table = None |
| 31 | + |
| 32 | + |
| 33 | +def cpu_info() -> str: |
| 34 | + """Best-effort, portable one-line CPU description |
| 35 | +
|
| 36 | + :return: CPU name and core count, plus clock speed if available |
| 37 | + :rtype: str |
| 38 | +
|
| 39 | + No new hard dependency: psutil is used for clock speed only if already |
| 40 | + installed. Some platforms (e.g. Apple Silicon) don't expose a single |
| 41 | + meaningful clock speed, so a missing or nonsensical reading is omitted. |
| 42 | + """ |
| 43 | + system = platform.system() |
| 44 | + name = None |
| 45 | + |
| 46 | + if system == "Darwin": |
| 47 | + try: |
| 48 | + name = subprocess.check_output( |
| 49 | + ["sysctl", "-n", "machdep.cpu.brand_string"], text=True |
| 50 | + ).strip() |
| 51 | + except Exception: |
| 52 | + pass |
| 53 | + elif system == "Linux": |
| 54 | + try: |
| 55 | + with open("/proc/cpuinfo") as f: |
| 56 | + for line in f: |
| 57 | + if line.lower().startswith("model name"): |
| 58 | + name = line.split(":", 1)[1].strip() |
| 59 | + break |
| 60 | + except Exception: |
| 61 | + pass |
| 62 | + elif system == "Windows": |
| 63 | + name = platform.processor() or None |
| 64 | + |
| 65 | + if not name: |
| 66 | + name = platform.processor() or platform.machine() or "unknown CPU" |
| 67 | + |
| 68 | + info = f"{name} ({os.cpu_count() or '?'} cores)" |
| 69 | + |
| 70 | + try: |
| 71 | + import psutil |
| 72 | + |
| 73 | + freq = psutil.cpu_freq() |
| 74 | + # real clock speeds are hundreds to thousands of MHz; some platforms |
| 75 | + # report bogus single-digit values instead of raising |
| 76 | + if freq and freq.max and freq.max > 100: |
| 77 | + info += f", {freq.max:.0f} MHz" |
| 78 | + except Exception: |
| 79 | + pass |
| 80 | + |
| 81 | + return info |
| 82 | + |
| 83 | + |
| 84 | +def print_machine_summary() -> None: |
| 85 | + print(f"CPU: {cpu_info()}") |
| 86 | + print(f"OS: {platform.platform()}") |
| 87 | + print(f"Python: {platform.python_version()}") |
| 88 | + print(f"numpy: {np.__version__}") |
| 89 | + print(f"spatialmath: {spatialmath.__version__}") |
| 90 | + print(f"Timing: min of {REPEATS} repeats x {N} calls") |
| 91 | + |
| 92 | + |
| 93 | +def new_table(): |
| 94 | + return ANSITable( |
| 95 | + Column("Operation", headalign="^"), |
| 96 | + Column("Time (μs)", headalign="^", fmt="{:.2f}"), |
| 97 | + border="thick", |
| 98 | + ) |
| 99 | + |
| 100 | + |
| 101 | +def timeit(stmt, label, setup): |
| 102 | + # min-of-repeats suppresses system jitter (GC pauses, OS scheduling) |
| 103 | + # much better than a single long run of the same total work |
| 104 | + t = min(_timeit.repeat(stmt=stmt, setup=setup, number=N, repeat=REPEATS)) |
| 105 | + table.row(label, t / N * 1e6) |
| 106 | + |
| 107 | + |
| 108 | +def section(title): |
| 109 | + # print the previous table (if any), then start a fresh one under a |
| 110 | + # plain header naming the category of test that follows |
| 111 | + global table |
| 112 | + if table is not None: |
| 113 | + table.print() |
| 114 | + print(f"\n{title}\n") |
| 115 | + table = new_table() |
| 116 | + |
| 117 | + |
| 118 | +print_machine_summary() |
| 119 | + |
| 120 | +# ------------------------------------------------------------------------- # |
| 121 | +transforms_setup = ''' |
| 122 | +from spatialmath import SE3 |
| 123 | +from spatialmath import base |
| 124 | +
|
| 125 | +import numpy as np |
| 126 | +from collections import namedtuple |
| 127 | +Rt = namedtuple('Rt', 'R t') |
| 128 | +X1 = SE3.Rand() |
| 129 | +X2 = SE3.Rand() |
| 130 | +T1 = X1.A |
| 131 | +T2 = X2.A |
| 132 | +R1 = base.t2r(T1) |
| 133 | +R2 = base.t2r(T2) |
| 134 | +t1 = base.transl(T1) |
| 135 | +t2 = base.transl(T2) |
| 136 | +Rt1 = Rt(R1, t1) |
| 137 | +Rt2 = Rt(R2, t2) |
| 138 | +v = np.r_[1,2,3] |
| 139 | +v2 = np.r_[1,2,3, 1] |
| 140 | +''' |
| 141 | + |
| 142 | +section("SO(3) / SE(3) base functions") |
| 143 | + |
| 144 | +timeit('base.getvector(0.2)', "base.getvector(x)", transforms_setup) |
| 145 | +timeit('base.rotx(0.2, unit="rad")', "base.rotx(x)", transforms_setup) |
| 146 | +timeit('base.trotx(0.2, unit="rad")', "base.trotx(x)", transforms_setup) |
| 147 | +timeit('base.t2r(T1)', "base.t2r(T1)", transforms_setup) |
| 148 | +timeit('base.r2t(R1)', "base.r2t(R1)", transforms_setup) |
| 149 | +timeit('T1 @ T2', "T1 @ T2 (4x4)", transforms_setup) |
| 150 | +timeit( |
| 151 | + 'T1[:3,:3] @ T2[:3,:3] + T1[:3,:3] @ T2[:3,3]', |
| 152 | + "T1 @ T2 decomposed (slice)", |
| 153 | + transforms_setup, |
| 154 | +) |
| 155 | +timeit( |
| 156 | + '(Rt1.R @ Rt2.R, Rt1.R @ Rt2.t)', |
| 157 | + "T1 @ T2 decomposed (namedtuple)", |
| 158 | + transforms_setup, |
| 159 | +) |
| 160 | +timeit('base.trinv(T1)', "base.trinv(T1)", transforms_setup) |
| 161 | +timeit( |
| 162 | + '(Rt1.R.T, -Rt1.R.T @ Rt1.t)', "T1 inverse decomposed (R,t)", transforms_setup |
| 163 | +) |
| 164 | +timeit('np.linalg.inv(T1)', "np.linalg.inv(T1)", transforms_setup) |
| 165 | +timeit('T1 @ v2', "T1 @ v2 (4,4)*(4,)", transforms_setup) |
| 166 | + |
| 167 | +# ------------------------------------------------------------------------- # |
| 168 | +section("SE3 class") |
| 169 | + |
| 170 | +timeit('SE3()', "SE3()", transforms_setup) |
| 171 | +timeit('SE3.Rx(0.2)', "SE3.Rx(x)", transforms_setup) |
| 172 | +timeit('SE3(T1)', "SE3(T1)", transforms_setup) |
| 173 | +timeit('SE3(T1, check=False)', "SE3(T1, check=False)", transforms_setup) |
| 174 | +timeit('SE3([T1], check=False)', "SE3([T1], check=False)", transforms_setup) |
| 175 | +timeit('T1[:3,:3]', "T1[:3,:3] (raw slice)", transforms_setup) |
| 176 | +timeit('X1.A', "X1.A (property)", transforms_setup) |
| 177 | +timeit('X1 * X2', "X1 * X2", transforms_setup) |
| 178 | +timeit('X1 @ X2', "X1 @ X2 (normalized)", transforms_setup) |
| 179 | +timeit('X1.inv()', "X1.inv()", transforms_setup) |
| 180 | +timeit('X1 * v', "X1 * v", transforms_setup) |
| 181 | +timeit('a = X1.log()', "X1.log()", transforms_setup) |
| 182 | +timeit('SE3.RPY([0.1, 0.2, 0.3])', "SE3.RPY(rpy)", transforms_setup) |
| 183 | +timeit('X1.rpy()', "X1.rpy()", transforms_setup) |
| 184 | +timeit('SE3.Eul([0.1, 0.2, 0.3])', "SE3.Eul(eul)", transforms_setup) |
| 185 | +timeit('X1.eul()', "X1.eul()", transforms_setup) |
| 186 | +timeit('X1.UnitQuaternion()', "X1.UnitQuaternion()", transforms_setup) |
| 187 | + |
| 188 | +# ------------------------------------------------------------------------- # |
| 189 | +quat_setup = ''' |
| 190 | +from spatialmath import base |
| 191 | +from spatialmath import UnitQuaternion, SO3 |
| 192 | +import numpy as np |
| 193 | +q1 = base.qrand() |
| 194 | +q2 = base.qrand() |
| 195 | +v = np.r_[1,2,3] |
| 196 | +R1 = SO3.Rand().R |
| 197 | +Q1 = UnitQuaternion.Rx(0.2) |
| 198 | +Q2 = UnitQuaternion.Ry(0.3) |
| 199 | +''' |
| 200 | + |
| 201 | +section("Quaternion base functions") |
| 202 | + |
| 203 | +timeit('a = base.qqmul(q1,q2)', "base.qqmul(q1, q2)", quat_setup) |
| 204 | +timeit('a = base.qvmul(q1,v)', "base.qvmul(q1, v)", quat_setup) |
| 205 | + |
| 206 | +# ------------------------------------------------------------------------- # |
| 207 | +section("UnitQuaternion class") |
| 208 | + |
| 209 | +timeit('a = UnitQuaternion()', "UnitQuaternion()", quat_setup) |
| 210 | +timeit('a = UnitQuaternion.Rx(0.2)', "UnitQuaternion.Rx(x)", quat_setup) |
| 211 | +timeit('a = UnitQuaternion(R1)', "UnitQuaternion(R1)", quat_setup) |
| 212 | +timeit('a = Q1.R', "Q1.R (property)", quat_setup) |
| 213 | +timeit('a = Q1 * Q2', "Q1 * Q2", quat_setup) |
| 214 | +timeit('a = Q1 * v', "Q1 * v", quat_setup) |
| 215 | +timeit('a = Q1.SE3()', "Q1.SE3()", quat_setup) |
| 216 | + |
| 217 | +# ------------------------------------------------------------------------- # |
| 218 | +twist_setup = ''' |
| 219 | +from spatialmath import SE3, Twist3 |
| 220 | +from spatialmath import base |
| 221 | +import numpy as np |
| 222 | +from math import cos |
| 223 | +S1 = Twist3(SE3.Rand()) |
| 224 | +S2 = Twist3(SE3.Rand()) |
| 225 | +X1 = SE3.Rand() |
| 226 | +T1 = X1.A |
| 227 | +A1 = X1.Ad() |
| 228 | +se3 = S1.skewa() |
| 229 | +s = np.r_[1,2,3,4,5,6] |
| 230 | +v = np.r_[1,2,3] |
| 231 | +''' |
| 232 | + |
| 233 | +section("Twist / exponential-map base functions") |
| 234 | + |
| 235 | +timeit('a = base.skew(v)', "base.skew(v)", twist_setup) |
| 236 | +timeit('a = base.skewa(s)', "base.skewa(s)", twist_setup) |
| 237 | +timeit('a = base.vexa(se3)', "base.vexa(se3)", twist_setup) |
| 238 | +timeit('a = base.trlog(T1)', "base.trlog(T1)", twist_setup) |
| 239 | +timeit('a = base.trlog(T1, twist=True)', "base.trlog(T1, twist=True)", twist_setup) |
| 240 | +timeit('a = base.trexp(se3)', "base.trexp(se3)", twist_setup) |
| 241 | +timeit('a = base.rodrigues(v)', "base.rodrigues(v)", twist_setup) |
| 242 | +timeit('a = A1 @ s', "A1 @ s (6,6)*(6,)", twist_setup) |
| 243 | +timeit('a = cos(0.3)', "math.cos(x)", twist_setup) |
| 244 | +timeit('a = np.cos(0.3)', "np.cos(x)", twist_setup) |
| 245 | + |
| 246 | +# ------------------------------------------------------------------------- # |
| 247 | +section("Twist3 class") |
| 248 | + |
| 249 | +# NB: Twist3 * Twist3 and Twist3.Ad() both round-trip through the SE3 |
| 250 | +# exponential/logarithm maps (trexp/trlog above), which is why they cost |
| 251 | +# several times more than a single trexp or trlog call - see |
| 252 | +# claude-notes/twist3-timing-investigation.md for the profiling detail. |
| 253 | + |
| 254 | +timeit('a = Twist3()', "Twist3()", twist_setup) |
| 255 | +timeit('a = Twist3(X1)', "Twist3(X1) (via log)", twist_setup) |
| 256 | +timeit('a = S1.inv()', "S1.inv()", twist_setup) |
| 257 | +timeit('a = S1 * S2', "S1 * S2 (product of exponentials)", twist_setup) |
| 258 | +timeit('a = S1.Ad()', "S1.Ad() (via SE3 + tr2adjoint)", twist_setup) |
| 259 | +timeit('a = S1.exp(1)', "S1.exp() -> SE3", twist_setup) |
| 260 | + |
| 261 | +# ------------------------------------------------------------------------- # |
| 262 | +misc_setup = """ |
| 263 | +from spatialmath import base |
| 264 | +import numpy as np |
| 265 | +s = np.r_[1.0,2,3,4,5,6] |
| 266 | +s3 = np.r_[1.0,2,3] |
| 267 | +a = np.r_[1.0, 2.0, 3.0] |
| 268 | +b = np.r_[-5.0, 4.0, 3.0] |
| 269 | +
|
| 270 | +A = np.random.randn(6,6) |
| 271 | +As = (A + A.T) / 2 |
| 272 | +bb = np.random.randn(6) |
| 273 | +""" |
| 274 | + |
| 275 | +section("NumPy linear-algebra baseline") |
| 276 | + |
| 277 | +timeit("c = np.linalg.inv(As)", "np.linalg.inv(As)", misc_setup) |
| 278 | +timeit("c = np.linalg.pinv(As)", "np.linalg.pinv(As)", misc_setup) |
| 279 | +timeit("c = np.linalg.solve(As, bb)", "np.linalg.solve(As, b)", misc_setup) |
| 280 | +timeit("c = np.cross(a,b)", "np.cross(a, b)", misc_setup) |
| 281 | +timeit("c = base.cross(a,b)", "base.cross(a, b)", misc_setup) |
| 282 | +timeit("a = np.inner(s,s).sum()", "np.inner(s, s).sum()", misc_setup) |
| 283 | +timeit("a = np.linalg.norm(s) ** 2", "np.linalg.norm(s) ** 2", misc_setup) |
| 284 | +timeit("a = base.normsq(s)", "base.normsq(s)", misc_setup) |
| 285 | +timeit("a = (s ** 2).sum()", "(s ** 2).sum()", misc_setup) |
| 286 | +timeit("a = np.sum(s ** 2)", "np.sum(s ** 2)", misc_setup) |
| 287 | +timeit("a = np.linalg.norm(s)", "np.linalg.norm(s) [R6]", misc_setup) |
| 288 | +timeit("a = base.norm(s)", "base.norm(s) [R6]", misc_setup) |
| 289 | +timeit("a = np.linalg.norm(s3)", "np.linalg.norm(s3) [R3]", misc_setup) |
| 290 | +timeit("a = base.norm(s3)", "base.norm(s3) [R3]", misc_setup) |
| 291 | + |
| 292 | +table.print() |
0 commit comments