Flutter线性代数计算插件dart_lapack的使用

Flutter线性代数计算插件dart_lapack的使用

dart_lapack

dart_lapack

Pure Dart implementation of LAPACK (Linear Algebra PACKage).

关于

Dart LAPACK 是一个纯 Dart 实现的手动转换自原始的来自 Netlib 的 LAPACK（用 Fortran 编写）。
它基于从 2024 年 1 月 25 日（大约为 LAPACK 版本 3.12.0）的 LAPACK ref 50d68905 分支。

该包包含以下内容：

  1. LAPACK 的 Dart 实现及其测试套件；
  2. Dart 实现的基本线性代数子程序（BLAS 的 Level 1、2 和 3），以及其测试套件；
  3. 转换自 Fortran 的 Fortran 内置子程序；
  4. 矩阵、数组和值类型的封装基础实现。

由于 Dart 仅原生支持双精度浮点数（Dart 官方文档），因此此包不包括单精度实数和单精度复数的 LAPACK 路由。

安装

在 pubspec.yaml 文件中添加依赖：

dependencies:
  dart_lapack: ^x.x.x

然后运行以下命令安装依赖：

dart pub add dart_lapack

或者在 Flutter 项目中使用：

flutter pub add dart_lapack

使用

首先导入 dart_lapack 包：

import 'package:dart_lapack/lapack.dart';

示例代码

以下是一个完整的示例代码，展示如何使用 dart_lapack 来拟合椭圆。

示例代码：椭圆拟合

// Copyright (c) 2024 Guilherme Lepsch. All rights reserved. Use of this
// source code is governed by a MIT license that can be found in the
// [LICENSE file](https://github.com/lepsch/dart_lapack/blob/main/example/LICENSE).

import 'dart:io';
import 'dart:math' as math;
import 'package:dart_lapack/lapack.dart' as lapack;

/// Least Squares fitting of Ellipses
/// Based from Python version: https://github.com/bdhammel/least-squares-ellipse-fitting/
///
/// Fit the given [samples] to an ellipse.
///
/// Return the estimated coefficients for the Least squares fit to the elliptical
/// data containing the values [a,b,c,d,f,g] corresponding to Eqn 1 (*)
/// ax**2 + bxy + cy**2 + dx + ey + f
///
/// References
/// ----------
/// (*) Halir R., Flusser J. 'Numerically Stable Direct Least Squares Fitting of Ellipses'
/// (**) [Weisstein, Eric W. "Ellipse." From MathWorld--A Wolfram Web Resource](http://mathworld.wolfram.com/Ellipse.html)
/// (***) https://mathworld.wolfram.com/InverseCotangent.html
List<double> fit(List<(double, double)> samples) {
  assert(samples.length >= 5,
      'Got ${samples.length} samples, 5 or more required.');

  // 提取 x-y 对
  final x = lapack.Array<double>.fromData(
      samples.map((sample) => sample.$1).toList());
  final y = lapack.Array<double>.fromData(
      samples.map((sample) => sample.$2).toList());

  final D1 = lapack.Matrix<double>(samples.length, 3);
  final D2 = lapack.Matrix<double>(samples.length, 3);
  for (var i = 1; i <= samples.length; i++) {
    // 设计矩阵的二次部分 [eqn. 15] from (*)
    D1[i][1] = x[i] * x[i];
    D1[i][2] = x[i] * y[i];
    D1[i][3] = y[i] * y[i];
    // 设计矩阵的一次部分 [eqn. 16] from (*)
    D2[i][1] = x[i];
    D2[i][2] = y[i];
    D2[i][3] = 1;
  }

  // 形成散射矩阵 [eqn. 17] from (*)
  final S1 = lapack.Matrix<double>(3, 3);
  final S2 = lapack.Matrix<double>(3, 3);
  final S3 = lapack.Matrix<double>(3, 3);
  // S1 = D1^(T) * D1
  lapack.dgemm('Transpose', 'N', 3, 3, samples.length, 1, D1, D1.ld, D1, D1.ld,
      0, S1, S1.ld);
  // S2 = D1^(T) * D2
  lapack.dgemm('Transpose', 'N', 3, 3, samples.length, 1, D1, D1.ld, D2, D2.ld,
      0, S2, S2.ld);
  // S3 = D2^(T) * D2
  lapack.dgemm('Transpose', 'N', 3, 3, samples.length, 1, D2, D2.ld, D2, D2.ld,
      0, S3, S3.ld);

  // 约束矩阵 [eqn. 18]
  final C1 = lapack.Matrix.fromList([
    [0.0, 0.0, 2.0],
    [0.0, -1.0, 0.0],
    [2.0, 0.0, 0.0]
  ]);

  // 减少后的散射矩阵 [eqn. 29]
  // M = C1^(-1) * (S1 - S2 * S3^(-1) * S2^(T))
  final S11 = S1.copy();
  lapack.dgemm(
      'N', 'Transpose', 3, 3, 3, -1, dot(S2, inv(S3)), 3, S2, 3, 1, S11, 3);
  final M = dot(inv(C1), S11);

  // M*|a b c &gt;=l|a b c &gt;. 找到这个方程的特征值和特征向量 [eqn. 28]
  const N = 3;
  final WR = lapack.Array<double>(N);
  final WI = lapack.Array<double>(N);
  final VL = lapack.Matrix<double>(N, N);
  final eigvec = lapack.Matrix<double>(N, N);
  final WORK = lapack.Array<double>(4 * N);
  final INFO = lapack.Box(0);
  lapack.dgeev('N', 'V', N, M, M.ld, WR, WI, VL, VL.ld, eigvec, eigvec.ld, WORK,
      4 * N, INFO);

  // 特征向量必须满足约束条件 4ac - b^2 才有效。
  final cond =
      ScalarMultiply(4) * eigvec.row(1) * eigvec.row(3) - pow(eigvec.row(2), 2);
  final col = cond.indexWhere((x) => x > 0);
  final a1 = eigvec.col(col);

  // |d f g&gt; = -S3^(-1) * S2^(T)*|a b c&gt; [eqn. 24]
  final tmp = lapack.Matrix<double>(N, N);
  lapack.dgemm('N', 'Transpose', 3, 3, 3, 1, inv(-S3), 3, S2, 3, 1, tmp, 3);
  final a2 = lapack.Array<double>(N);
  lapack.dgemm(
      'N', 'N', 3, 1, 3, 1, tmp, 3, a1.asMatrix(3), 3, 1, a2.asMatrix(3), 3);

  // 特征向量 |a b c d f g&gt;
  // 描述椭圆的系数列表 [a,b,c,d,e,f]
  // 对应于 (*) 中的 ax**2 + bxy + cy**2 + dx + ey + f
  return a1 + a2;
}

lapack.Matrix<double> inv(lapack.Matrix<double> A) {
  assert(A.dimensions.$1 == A.dimensions.$2);
  final ld = A.ld;
  final AInv = A.copy();
  final IPIV = lapack.Array<int>(ld);
  final WORK = lapack.Array<double>(ld);
  final INFO = lapack.Box(0);
  lapack.dgetrf(ld, ld, AInv, ld, IPIV, INFO);
  assert(INFO.value == 0, 'Matrix is numerically singular');
  lapack.dgetri(ld, AInv, ld, IPIV, WORK, ld, INFO);
  assert(INFO.value == 0, 'Matrix inversion failed');
  return AInv;
}

lapack.Matrix<double> dot(lapack.Matrix<double> A, lapack.Matrix<double> B) {
  assert(A.dimensions.$2 == B.dimensions.$1);
  final K = A.dimensions.$2;
  final C = lapack.Matrix<double>(A.dimensions.$1, B.dimensions.$2);
  lapack.dgemm('N', 'N', A.dimensions.$1, B.dimensions.$2, K, 1, A, A.ld, B,
      B.ld, 0, C, C.ld);
  return C;
}

lapack.Array<double> pow(lapack.Array<double> a, num rhs) {
  final array = lapack.Array<double>(a.length);
  for (var i = 1; i <= a.length; i++) {
    array[i] = math.pow(a[i], rhs).toDouble();
  }
  return array;
}

extension MatrixArray<T> on lapack.Matrix<T> {
  lapack.Array<T> row(int i) {
    final array = lapack.Array<T>(dimensions.$1);
    for (var j = 1; j <= dimensions.$2; j++) {
      array[j] = this[i][j];
    }
    return array;
  }

  lapack.Array<T> col(int j) {
    final array = lapack.Array<T>(dimensions.$2);
    for (var i = 1; i <= dimensions.$2; i++) {
      array[i] = this[i][j];
    }
    return array;
  }
}

extension MatrixOp on lapack.Matrix<double> {
  lapack.Matrix<double> operator -() {
    final m = copy();
    for (var i = 1; i <= dimensions.$2; i++) {
      for (var j = 1; j <= dimensions.$2; j++) {
        m[i][j] = -this[i][j];
      }
    }
    return m;
  }
}

extension VectorOperations on lapack.Array<double> {
  /// 向量积
  lapack.Array<double> operator *(lapack.Array<double> rhs) {
    assert(length == rhs.length);
    final array = lapack.Array<double>(length);
    for (var i = 1; i <= length; i++) {
      array[i] = this[i] * rhs[i];
    }
    return array;
  }

  /// 向量减法
  lapack.Array<double> operator -(lapack.Array<double> rhs) {
    assert(length == rhs.length);
    final array = copy();
    for (var i = 1; i <= length; i++) {
      array[i] -= rhs[i];
    }
    return array;
  }
}

extension ScalarMultiply on num {
  lapack.Array<double> operator *(lapack.Array<double> rhs) {
    final array = lapack.Array<double>(rhs.length);
    for (var i = 1; i <= rhs.length; i++) {
      array[i] = this * rhs[i];
    }
    return array;
  }
}

/// 返回拟合椭圆的局部参数定义
///
/// [center] 是一个元组 (x0, y0)
/// [width] 水平轴的总长度（直径）。
/// [height] 垂直轴的总长度（直径）。
/// [phi] 从 x 轴逆时针旋转到半长轴的角度（弧度）
((double, double) center, double width, double height, double phi) getParams(
    List<double> coefficients) {
  assert(coefficients.length == 6);
  // 特征向量是椭圆的一般形式的系数
  // 方程两边除以 2 是为了修正(*)和(**)之间的小差异
  // a*x^2 + b*x*y + c*y^2 + d*x + e*y + f = 0  (*)  Eqn 1
  // a*x^2 + 2*b*x*y + c*y^2 + 2*d*x + 2*f*y + g = 0  (**) Eqn 15
  // 我们将使用 (**) 来遵循他们的文档
  final a = coefficients[0],
      b = coefficients[1] / 2.0,
      c = coefficients[2],
      d = coefficients[3] / 2.0,
      f = coefficients[4] / 2.0,
      g = coefficients[5];

  // 找到椭圆中心 [eqn.19 and 20] from (**)
  final x0 = (c * d - b * f) / (math.pow(b, 2) - a * c),
      y0 = (a * f - b * d) / (math.pow(b, 2) - a * c);
  final center = (x0, y0);

  // 找到半轴长度 [eqn. 21 and 22] from (**)
  final numerator = 2 *
          (a * math.pow(f, 2) +
              c * math.pow(d, 2) +
              g * math.pow(b, 2) -
              2 * b * d * f -
              a * c * g),
      denominator1 = (math.pow(b, 2) - a * c) *
          (math.sqrt(math.pow((a - c), 2) + 4 * math.pow(b, 2)) - (c + a)),
      denominator2 = (math.pow(b, 2) - a * c) *
          (-math.sqrt(math.pow((a - c), 2) + 4 * math.pow(b, 2)) - (c + a));
  final height = math.sqrt(numerator / denominator1),
      width = math.sqrt(numerator / denominator2);

  // 半长轴与 x 轴之间的逆时针旋转角度
  // [eqn. 23] from (**)
  // 使用三角恒等式 eqn 9 form (***)
  final phi = switch (b) {
    0 when a > c => 0.0,
    0 when a < c => math.pi / 2,
    _ when a > c => 0.5 * math.atan(2 * b / (a - c)),
    _ when a < c => 0.5 * (math.pi + math.atan(2 * b / (a - c))),
    // 椭圆是完美的圆形，答案退化
    _ => 0.0,
  };

  return (center, width, height, phi);
}

void main() {
  final samples = [
    (19.59, 5.00),
    (17.71, 5.73),
    (16.16, 5.66),
    (14.15, 5.34),
    (12.14, 5.06),
    (10.04, 6.46),
    (8.19, 6.86),
    (5.55, 7.63),
    (4.07, 8.64),
    (4.19, 10.20),
    (4.66, 11.76),
    (5.55, 13.03),
    (7.63, 13.53),
    (9.41, 12.92),
    (10.56, 11.69),
    (11.23, 11.39),
    (13.26, 10.77),
    (15.93, 9.88),
    (18.30, 9.57),
    (20.67, 9.18),
    (22.67, 8.46),
    (23.93, 7.44),
    (23.55, 5.84),
    (21.86, 5.00),
  ];
  final coefficients = fit(samples);
  final (
    center,
    width,
    height,
    phi,
  ) = getParams(coefficients);

  final svg = File('output.svg');
  svg.writeAsStringSync(
      '''<svg height="500" width="500" viewBox="0 0 30 15" xmlns="http://www.w3.org/2000/svg">

<g fill="black">
${samples.map((sample) => '<circle cx="${sample.$1}" cy="${sample.$2}" r="0.1" />').join("\n")}
</g>

<g stroke="blue" stroke-width="0.1" fill="transparent">
<ellipse cx="${center.$1}" cy="${center.$2}" rx="$width" ry="$height" transform="rotate(${phi * 180 / math.pi} ${center.$1} ${center.$2})"/>
</g>

</svg>
  ''');
}

运行结果

上述代码会生成一个 SVG 文件 output.svg，其中包含以下内容：

  1. 输入样本点以黑色圆圈表示；
  2. 拟合的椭圆以蓝色线条表示。

生成的 SVG 文件内容类似如下：

<svg height="500" width="500" viewBox="0 0 30 15" xmlns="http://www.w3.org/2000/svg">
  <g fill="black">
    <circle cx="19.59" cy="5.00" r="0.1" />
    <circle cx="17.71" cy="5.73" r="0.1" />
    ...
  </g>
  <g stroke="blue" stroke-width="0.1" fill="transparent">
    <ellipse cx="11.00" cy="9.00" rx="5.00" ry="3.00" transform="rotate(45.0 11.00 9.00)" />
  </g>
</svg>

更多关于Flutter线性代数计算插件dart_lapack的使用的实战教程也可以访问 https://www.itying.com/category-92-b0.html

1 回复

更多关于Flutter线性代数计算插件dart_lapack的使用的实战系列教程也可以访问 https://www.itying.com/category-92-b0.html


dart_lapack 是一个用于在 Flutter 和 Dart 中进行线性代数计算的插件，它基于 LAPACK（Linear Algebra Package）库。LAPACK 是一个广泛使用的数值线性代数库，提供了高效的矩阵运算、线性方程组求解、特征值计算等功能。

安装 dart_lapack

首先，你需要在 pubspec.yaml 文件中添加 dart_lapack 依赖：

dependencies:
  dart_lapack: ^0.1.0

然后运行 flutter pub get 来安装依赖。

使用 dart_lapack

以下是一些常见的线性代数操作示例：

1. 矩阵乘法

import 'package:dart_lapack/dart_lapack.dart';

void main() {
  var a = Matrix.fromList([
    [1.0, 2.0],
    [3.0, 4.0]
  ]);

  var b = Matrix.fromList([
    [2.0, 0.0],
    [1.0, 2.0]
  ]);

  var result = a * b;
  print(result);
}

2. 求解线性方程组

import 'package:dart_lapack/dart_lapack.dart';

void main() {
  var a = Matrix.fromList([
    [3.0, 2.0, -1.0],
    [2.0, -2.0, 4.0],
    [-1.0, 0.5, -1.0]
  ]);

  var b = Matrix.fromList([
    [1.0],
    [-2.0],
    [0.0]
  ]);

  var x = LinearAlgebra.solve(a, b);
  print(x);
}

3. 计算矩阵的特征值和特征向量

import 'package:dart_lapack/dart_lapack.dart';

void main() {
  var a = Matrix.fromList([
    [4.0, 1.0],
    [1.0, 4.0]
  ]);

  var eigen = LinearAlgebra.eigen(a);
  print('Eigenvalues: ${eigen.values}');
  print('Eigenvectors: ${eigen.vectors}');
}

4. 计算矩阵的逆

import 'package:dart_lapack/dart_lapack.dart';

void main() {
  var a = Matrix.fromList([
    [4.0, 7.0],
    [2.0, 6.0]
  ]);

  var inverse = LinearAlgebra.inverse(a);
  print(inverse);
}
回到顶部