-
Notifications
You must be signed in to change notification settings - Fork 1
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Implement own approximation of terrain #2
Comments
Comparing of OCC B-Spline surface and our B-Spline surface (yellow crosses). There is small difference, but I have to note that approximated points are on OCC surface as well on our surface. When OCC B-Spline surface and SciPy B-Spline surface (yellow corsses) is compared, then the surface is the same: |
This is quite significant difference. So either we do not set data of the OCC B-spline correctly or we have wrong evaluation of the B-spline. Try to compare on simpler B-spline, e.g. the surface constant in one direction. |
Concerning difference maps. What is the color scale? |
There is probably some difference between our and OCC implementation of B-Spline surfaces. I will do investigation of OCC source code tomorrow. Yes, color ramps of difference maps can be misleading, because maximal values are different. I will add information about max values and I will try to use same color ramps. |
Great. On 03/09/2016 11:14 AM, Jiri Hnidek wrote:
Mgr. Jan Brezina, Ph. D. |
I just discovered that it is possible to set number of knot vectors for SciPy too. Thus I compared computation of B-Spline approximation using 100x100 knot vectors
More details:
I'm attaching generated terrains:
I updated QR metod a little (8453 seconds improved to 1453 seconds):
The method: |
Can you set only number of knots or the knot vectors itself? I.e. can SciPy How you do calculation of the Z matrix (do you mean Z vector?)? On Thu, Mar 17, 2016 at 12:48 PM, Jiri Hnidek notifications@github.com
|
When I'm looking at: http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.interpolate.bisplrep.html , then it seems that you can also force SciPy to use non-eqidistant knot vectors ( This is code of computation of Z matrix (useful: http://mathesaurus.sourceforge.net/matlab-numpy.html) print('Computing QR ...')
start_time = time.time()
# Converted from Matlab code: [q r] = qr(B);
q_mat, r_mat = numpy.linalg.qr(b_mat, mode='full')
end_time = time.time()
print('Computed in {0} seconds.'.format(end_time - start_time))
print('Computing Z matrix ...')
start_time = time.time()
# Converted from Matlab code: z = r \ q' * g
z_mat = numpy.linalg.lstsq(r_mat, q_mat.transpose())[0] * g_z_mat
end_time = time.time()
print('Computed in {0} seconds.'.format(end_time - start_time)) This code does not use sparse vector / matrices. |
Must use fact that r_mat is upper triangular to avoid lstsq solve and use just matrix multiplication. To be discussed. |
Consider: Rather use z_vec, g_z_vec ?? |
Sparse version of matrices is very slow for some reason. Sci offers many variants of sparse matrices: I use dictionary based matrices. "dok" matrices are intended for constructing matrices, but constructing is very slow:
Here is simple comparison: Dense matrices:
Sparse matrices:
|
Efficiency of computation with sparse matrices strongly depends on What about the csr format, especially for B? The same performance issues? On 04/07/2016 02:44 PM, Jiri Hnidek wrote:
|
@jirikopal it was computed for 755 points and the number of (u/v)knots was 20. I used "dok" format, because of assembling, but assembling was slow too (due to kron()?). CSC or CSR formats are good candidates for multiplication of matrices. I though that performance improvements of sparse matrix implementation could be good task for you :-). This Matlab code wasn't 100% clear to me: Sm = S > 1e-3;
S =S.*Sm;
r = rank(S);
Sdi = diag(diag(S(1:r,1:r)).^(-1)); %diag(diag()) is some evil matlab trick, right? :-)
Si = sparse(zeros(a,b));
Si(1:r,1:r) = Sdi; % IMHO you can crop some relevant values, when rank != min(a,b) BTW: I noticed big oscillations, when len(uknots) * len(vknots) > len(B). |
Try assembling in COO and then convert to CSR. This is common approach Making SVD etc. for sparse matrices make little sense, these are meant to solve B^t B z = B^t g iteratively. On 04/07/2016 02:44 PM, Jiri Hnidek wrote:
Mgr. Jan Brezina, Ph. D. |
@jbrezmorf I wrote that SVD used dense B matrix in both cases (dense and sparse solution), because scipy nor numpy does not implement SVD for sparse matrices. |
I know, that's why it makes little sense in terms of both memory and performance. Anyway, what is difference between "dense" and "sparse" case for other steps "Si matrix" and "Z matrix"? Could you past a code snippet? |
@jbrezmorf it is in repo in branch JH-svd-approx: https://github.com/GeoMop/bapprox/blob/JH-svd-approx/src/approx/terrain.py#L392 Simplified: print('Creating Si matrix ...')
if sparse is True:
mat_u = scipy.sparse.dok_matrix(mat_u)
mat_v = scipy.sparse.dok_matrix(mat_v)
# Make compatible with Matlab code
mat_v = mat_v.transpose()
# Filter too small values and invert other values
for key, value in enumerate(mat_s):
if value < filter_thresh:
mat_s[key] = 0.0
else:
mat_s[key] = 1.0 / value
size = max(mat_s.shape)
if sparse is True:
mat_si = scipy.sparse.spdiags(mat_s, 0, size, size)
else:
mat_si = numpy.diagflat(mat_s)
print('Computing Z matrix ...')
z_mat = mat_v * (mat_si.transpose() * (mat_u.transpose() * mat_g)) |
Regularization term solve this problem. On 04/11/2016 10:12 AM, Jiri Hnidek wrote:
|
Implement own approximation of terrain data.
The text was updated successfully, but these errors were encountered: