- html - 出于某种原因,IE8 对我的 Sass 文件中继承的 html5 CSS 不友好?
- JMeter 在响应断言中使用 span 标签的问题
- html - 在 :hover and :active? 上具有不同效果的 CSS 动画
- html - 相对于居中的 html 内容固定的 CSS 重复背景?
我在 cython 中编写了以下 Matrix 类
用于矩阵求逆和其他一些线性代数运算。我尝试使用LU 分解 来计算矩阵的逆。代码速度很好。我试图实现 this code在 cython
中。我已经检查了我的每一行代码并与给定的代码进行了几次比较,但我仍然返回错误的答案。
矩阵.pyx
from libcpp.vector cimport vector
import cython
cimport cython
import numpy as np
cimport numpy as np
import ctypes
from libc.math cimport log, exp, pow, fabs
from libc.stdint cimport *
from libcpp.string cimport string
from libc.stdio cimport *
from libcpp cimport bool
cdef extern from "<iterator>" namespace "std" nogil:
cdef cppclass iterator[Category, T, Distance, Pointer, Reference]:
pass
cdef cppclass output_iterator_tag:
pass
cdef cppclass input_iterator_tag:
pass
cdef cppclass forward_iterator_tag(input_iterator_tag):
pass
cdef extern from "<algorithm>" namespace "std" nogil:
void fill [ForwardIterator, T](ForwardIterator, ForwardIterator, T& )
cdef class Matrix:
def __cinit__(self, size_t rows=0, size_t columns=0, bint Identity=False, bint ones=False):
self._rows=rows
self._columns=columns
self.matrix=new vector[double]()
self.matrix.resize(rows*columns)
if Identity:
self._IdentityMatrix()
if ones:
self._fillWithOnes()
def __dealloc__(self):
del self.matrix
property rows:
def __get__(self):
return self._rows
def __set__(self, size_t x):
self._rows = x
property columns:
def __get__(self):
return self._columns
def __set__(self, size_t y):
self._columns = y
cpdef double getVal(self, size_t r, size_t c):
return self.matrix[0][r*self._columns+c]
cpdef void setVal(self, size_t r, size_t c, double v):
self.matrix[0][r*self._columns+c] = v
@cython.boundscheck(False)
@cython.wraparound(False)
cdef void _fillWithOnes(self):
fill(self.matrix.begin(),self.matrix.end(),1.)
cdef void _IdentityMatrix(self):
cdef size_t i
if (self._rows!=self._columns):
raise Exception('In order to generate identity matrix, the number of rows and columns must be equal')
else:
for i from 0 <= i <self._columns:
self.setVal(i,i,1.0)
@cython.boundscheck(False)
@cython.wraparound(False)
cpdef Matrix Inv(self):
cdef Matrix A_inverse = Matrix(self._rows,self._columns)
cdef MatrixList LU = ludcmp(self)
cdef Matrix A = LU.get(0)
cdef Matrix indx = LU.get(1)
cdef Matrix d = LU.get(2)
cdef double det = d.getVal(0,0)
cdef int i, j
cdef np.ndarray[np.float64_t, ndim=2] L = np.zeros((self._rows,self._columns),dtype=np.float64)
cdef np.ndarray[np.float64_t, ndim=2] U = np.zeros((self._rows,self._columns),dtype=np.float64)
cdef Matrix col = Matrix(self._rows,1)
for i from 0 <= i < self._rows:
for j from 0 <= j < self._columns:
if (j>i):
U[i,j]=A.getVal(i,j)
L[i,j]=0
elif (j<i):
U[i,j]=0
L[i,j]=A.getVal(i,j)
else:
U[i,j]=A.getVal(i,j)
L[i,j]=1
print "product of a lower triangular matrix L and an upper triangular matrix U:", np.dot(L, U)
for i from 0 <= i < self._rows:
det*= A.getVal(i,i)
for j from 0 <= j < self._columns:
if (i==j):
col.setVal(j,0,1)
col=lubksb(A, indx, col)
for j from 0 <= j < self._columns:
A_inverse.setVal(j,i,col.getVal(j,0))
print "determinant of matrix %.4f"%(det)
return A_inverse
cdef class MatrixList:
def __cinit__(self):
self.inner = []
cdef void append(self, Matrix a):
self.inner.append(a)
cdef Matrix get(self, int i):
return <Matrix> self.inner[i]
def __len__(self):
return len(self.inner)
@cython.boundscheck(False)
@cython.wraparound(False)
cdef Matrix lubksb(Matrix a, Matrix indx, Matrix b):
cdef int n = a.rows
cdef int i, ip, j
cdef int ii = 0
cdef double su
for i from 0 <= i < n:
ip = <int>indx.getVal(i,0)
su = b.getVal(ip,0)
b.setVal(ip,0, b.getVal(i,0))
if (ii):
for j from ii <= j < (i-1):
su -= a.getVal(i,j) * b.getVal(j,0)
elif (su):
ii = i
b.setVal(i, 0, su)
for i from n > i >= 0:
su = b.getVal(i,0)
for j from (i+1) <= j < n:
su -= a.getVal(i,j) * b.getVal(j,0)
b.setVal(i, 0, su/a.getVal(i,i))
return b
@cython.boundscheck(False)
@cython.wraparound(False)
cdef MatrixList ludcmp(Matrix a):
#Given a matrix a_{nxn}, this routine replaces it by the LU decomposition of a row-wise permutation of itself.
cdef MatrixList LU = MatrixList()
cdef int n = a.rows
cdef int i, j, k, imax
cdef double big, dum, su, temp
cdef Matrix vv = Matrix(n,1)
cdef Matrix indx = Matrix(n,1) #an output vector that records the row permutation effected by the partial pivoting
cdef Matrix d = Matrix(1,1, ones= True) #an output as +1 or -1 depending on whether the number of row interchanges was even or odd, respectively
cdef double TINY = 1.1e-16
for i from 0 <= i < n:
big = 0.0
for j from 0 <= j < n:
temp=fabs(a.getVal(i,j))
if (temp > big):
big=temp
if (big ==0.0):
raise Exception("ERROR! ludcmp: Singular matrix\n")
vv.setVal(i,0,1.0/big)
for j from 0 <= j < n:
for i from 0 <= i < j:
su = a.getVal(i,j)
for k from 0 <= k < i:
su -= a.getVal(i,k)*a.getVal(k,j)
a.setVal(i,j,su)
big=0.0
for i from j<= i< n:
su = a.getVal(i,j)
for k from 0 <= k < j:
su -= a.getVal(i,k)*a.getVal(k,j)
a.setVal(i, j, su)
dum=vv.getVal(i,0)*fabs(su )
if (dum >= big):
big=dum
imax=i
if (j != imax):
for k from 0 <= k < n:
dum = a.getVal(imax,k)
a.setVal(imax, k, a.getVal(j,k))
a.setVal(j,k, dum)
d.setVal(0, 0, -d.getVal(0,0))
vv.setVal(imax, 0, vv.getVal(j, 0))
indx.setVal(j, 0, imax)
if (a.getVal(j,j) == 0.0):
a.setVal(j,j, TINY)
if (j != (n-1)):
dum=1.0/a.getVal(j,j)
for i from (j+1)<= i <n:
a.setVal(i,j, a.getVal(i,j)*dum)
LU.append(<Matrix>a)
LU.append(<Matrix>indx)
LU.append(<Matrix>d)
return LU
矩阵.pxd
from libcpp.vector cimport vector
cdef class MatrixList:
cdef list inner
cdef void append(self, Matrix a)
cdef Matrix get(self, int i)
cdef class Matrix:
cdef vector[double] *matrix
cdef size_t _rows
cdef size_t _columns
cdef bint Identity
cdef bint ones
cpdef double getVal(self, size_t r, size_t c)
cpdef void setVal(self, size_t r, size_t c, double v)
cpdef Matrix transpose(self)
cdef void _IdentityMatrix(self)
cdef void _fillWithOnes(self)
cpdef Matrix Inv(self)
cdef Matrix lubksb(Matrix a, Matrix indx, Matrix b)
cdef MatrixList ludcmp(Matrix a)
如果您能帮助我们找到错误,我们将不胜感激。
示例:
import numpy
from matrix import Matrix
from numpy.linalg import inv
import timeit
import numpy as np
r=numpy.random.random((100, 100))
d=Matrix(r.shape[0],r.shape[1])
for i in range(d.rows):
for j in range(d.columns):
d.setVal(i,j,r[i,j])
start = timeit.default_timer()
x=d.Inv()
stop = timeit.default_timer()
print "LU decomposition:", stop - start
最佳答案
我发现我在 lubksb
函数中犯了一些非常小的错误,通过修复这些错误我得到了正确的答案。这是固定代码:
@cython.boundscheck(False)
@cython.wraparound(False)
cdef Matrix lubksb(Matrix a, Matrix indx, Matrix b):
cdef int n = a.rows
cdef int i, ip, j
cdef int ii = 0
cdef double su
for i from 0 <= i < n:
ip = <int>indx.getVal(i,0)
su = b.getVal(ip,0)
b.setVal(ip,0, b.getVal(i,0))
if (ii>=0):
for j from ii <= j <= (i-1):
su -= a.getVal(i,j) * b.getVal(j,0)
elif (su):
ii = i
b.setVal(i, 0, su)
for i from n >= i >= 0:
su = b.getVal(i,0)
for j from (i+1) <= j < n:
su -= a.getVal(i,j) * b.getVal(j,0)
if (a.getVal(i,i)==0.0):
a.setVal(i,i, 1.1e-16)
b.setVal(i, 0, su/a.getVal(i,i))
return b
关于python - 使用 LU 分解求逆矩阵,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/48491962/
要在 R 中绘制正态分布曲线,我们可以使用: (x = seq(-4,4, length=100)) y = dnorm(x) plot(x, y) 如 dnorm将 y 计算为 x 的函数,R 是否
@XmlTransient 阻止将 JavaBeans 属性映射到 XML 表示。是否存在与此相反的情况,这意味着即使 WebService 未使用的方法也会被映射?如果这不可能,是否存在解决方法?
我有以下键数组: var keys = [{userId: "333"}, {userId: "334"}] 这个对象数组: var users = [ {id: "333", firstName:
我正在寻找将字符串转换为类型的通用方法。 例如: class SomeThing { public void Add(T value) { //... } pub
我看到了this question , 并弹出这个想法。 有没有一种在 PHP 中执行此操作的有效方法? 编辑 有演示最好吗? 最佳答案 你可以使用 pear 包 Math_Matrix为此。 关于矩
如何在 python 中求逆矩阵?我自己实现了它,但它是纯 python,我怀疑那里有更快的模块可以做到这一点。 最佳答案 你应该看看 numpy如果您进行矩阵操作。这是一个主要用C语言编写的模块,比
是否有比使用 IF ELSE 构造更简单的方法来反转 bool 值? 通常我会使用! bool 值前面。但这在 Navision 中不起作用 最佳答案 您可以使用 NOT 关键字代替 !。 关于nav
假设我有一个对象响应。现在我想检查一个 bool 变量,success,在 Response 下并做一个早期返回是 response 不成功。 if(response == null || !resp
任何人都可以提供/引用多维行主要顺序的“索引->偏移”*转换的倒数。此外,(伪)代码将不胜感激。 http://en.wikipedia.org/wiki/Row-major_order 举个例子,简
我有一个看起来像这样的系统: z1 = 5*x1 + x2*cos(x3) z2 = x1*sin(x3) + 3*x2 z3 = 3*x1 - 2*x2 这是微分方程组的变换(只是为了提供一些背景信
我正在使用org.apache.commons.math3.transform类FastFourierTransformer,我现在尝试在真实数据集上应用FFT,并应用逆FFT来获取原始数据集。我的问
关闭。这个问题不符合Stack Overflow guidelines .它目前不接受答案。 这个问题似乎与 help center 中定义的范围内的编程无关。 . 关闭 9 年前。 Improve
背景 我需要使用已知的累积分布函数 (CDF) 从相当复杂的概率密度函数 (PDF) 中随机采样,并且我正在尝试使用 inverse transform sampling 。这应该很容易做到,因为我有
是否有任何 System.identityHashCode (object) 的逆函数能够从 System.identityHashCode (object) 的结果中提供对象的值? 最佳答案 Sys
有没有办法在mysql中获取group by语句的逆?我的用例是删除所有重复项。 假设我的表格如下所示: ID | columnA | ... 1 | A 2 | A 3 | A 4
我有一个查询,它给我一个公司列表(tblprov)及其相应的类别(tblrubro) 两个表通过查找表 (tblprovxrubro) 相关 SELECT p.id, p.name, r.idCat,
我有一个 jpg 图像,在矩形中有一个圆形物体,我想使圆形物体的环境透明... (本例去除红色区域) 借助这个iOS make part of an UIImage transparent和“UIBe
我想知道是否可以在不需要临时数组的情况下通过 Cholesky 分解获得矩阵的逆。截至目前,我可以在不使用临时数组的情况下进行 cholesky 分解,但从那里我还没有想出一种方法来获得原始矩阵的逆矩
是否可以在 Angular 中使用逆$watch? 我的问题 我使用 Angular-translate,并且我想对每个缺少的翻译使用 $http.put 。但我收到此错误: "10 $digets(
我正在执行 radix-2 dif 逆 fft。我正在使用共轭和缩放的属性来返回结果。我共轭我的输入 vector ,执行常规 radix-2 fft(不是 ifft),共轭结果,然后按 1.0/N
我是一名优秀的程序员,十分优秀!