few minor changes to run example code in FAIRFEM #5

hariagr · 2017-07-11T03:22:56Z

‘yn’ is not defined when wc is empty. Therefore, xn should be
updated by yRef.
Addpath for add-ons folders during FAIR startup.
solver is defined as a char, not string. Probably, this issue has
recently arises because of structural changes in R2017a.

1. ‘yn’ is not defined when wc is empty. Therefore, xn should be updated by yRef. 2. Addpath for add-ons folders during FAIR startup. 3. solver is defined as a char, not string. Probably, this issue has recently arises because of structural changes in R2017a.

These changes are required to use matrix-free functionality of hyper elastic regularization function in the add-on FEMAPP (Only for 2D). To verify the changes, the file “EFEM_SSD_MBvsMF.m” compares the results from MF and MB.

These are minor changes to run matrixfree functionality of FEMPIRE. But, data fidelity term is still not matrixfree. Prob, this involves little more work. Will come back to this later. Added a file to compare results from MF and MB.

lruthotto · 2017-07-15T06:58:23Z

thanks for this excellent job. Can you please resolve the conflict and clean up the .m~ files?

lruthotto · 2017-07-20T09:20:28Z

add-ons/FAIRFEM/EFEM_SSD_MBvsMF.m

@@ -0,0 +1,38 @@
+% ==================================================================================
+% (c) Lars Ruthotto 2012/07/26, see FAIR.2 and FAIRcopyright.m.


can you please update the description of the examples and use the new header from the other examples? We decided to let Github keep track of the authorship for us.

Please add a brief description about what this example is supposed to do.

Sure. I will do that

lruthotto · 2017-07-20T09:21:04Z

add-ons/FAIRFEM/EFEM_MP_MBvsMF.m

@@ -0,0 +1,38 @@
+% ==================================================================================
+% (c) Lars Ruthotto 2012/07/26, see FAIR.2 and FAIRcopyright.m.


See comment below about the description of the example.

Sure. I will do that.

lruthotto · 2017-07-20T09:24:28Z

add-ons/FAIRFEM/FEMPIREobjFctn.m

-    dJ = reshape(Mesh.mfPi(reshape(dD',[],Mesh.dim),'C'),1,[]) + dS;
+    %P  = @(x) reshape(Mesh.mfPi(reshape(x,[],dim),'C'),[],1);
+    %dTcmod = P((sdiag(det)*dT)')' + sdiag(Tc)*dDet;
+    P = kron(speye(dim),Mesh.PC);


this could be a costly operation (which will be done in each call!). In principle, P is only used once and we have matrixFree code for it.

I asssume dT below is sparse matrix of size n by dim*n, where n=prod(m). See if, supplying 'matrixFree',1 in the interpolation would render this a dense n by dim array. Then you can apply P or P' to each column in a matrix free fashion.

This should enable you to avoid the sdiag as well. Can you give this a try?

Certainly. I will work on it. Thanks for suggestion.

I have modified the code. Please see the new commit.

lruthotto · 2017-07-20T09:26:05Z

Hey Hari, I resolved the conflicts and added three problems. The one concerning the documentation are critical. See if you also want to address the one regarding the performance.

Matrixfree implementations of - derivative of determinants, - diagonal of hessian matrix, and - Hessian matrix

hariagr · 2017-08-13T12:06:42Z

Hi Lars,

I was busy with my work. Sorry for delay.

Now, I have modified the matrixfree implementation. Please have a look. Let me know if something else needs to be done.

Thanks,
Hari

lruthotto · 2017-08-31T15:08:53Z

Hi Hari,

I have one comment before we can merge this (already did two minor changes to this). My goal is to finish this soon, but it's important to get things right, I think.

I noticed that your code for computing the diagonal of the Hessian is only for 2D problems. It would be great if we could extend this to 3D, where matrix-free methods are most attractive. I believe that our mesh classes provide the necessary functionality already. See one change I did to average quantities from the nodes to the barycenters of the elements. Can you help me modify the remainder of this code accordingly and write a small test file for 3D? I have push access to your repo so we can work on this together.

Thanks a lot,
Lars

* 'master' of https://github.com/hariagr/FAIR.m: another minor edit two minor changes

Implement matrix free code to evaluate determinant. Added a script “verifyMF_getdeterminant” to verify results

Modified HyperElasticFEM function. Work to do: MatrixFree code for diagonal elements of Hessian matrix.

lruthotto · 2017-09-09T15:47:48Z

Hi Hari!

Thanks for your hard work and patience. It's getting close I think. One final thing that I see is to make the function getDiag in FEMPIREobjFctn work for 3D problems. If you take a look at EFEM_FEMPIRE_Mice3D the code runs until the Hessian solve. I think we need a switch statement here and implement the 2D and 3D separately.

Best,
Lars

MatrixFree code for diagonal of Hessian matrix is implemented. We may have an issue in the implementation. - if we set only alphaVolume parameter to check the accuracy of Matrixfree code with MB code, there is a slight error in the calculation. You can run hyperelastic.m to see the error norm. It appears that the differences are high only for few nodes. These nodes varies with different yc/yn values. Generally, d2Svolume has very high values. Small numerical errors in least decimal points are probably showing up. These are few quick observations.

Most of the pieces required for MatrixFree code are completed now. Run EFEM_FEMPIRE_Mice3D.m. I have not done rigorous testing so far. Specially for 3D code, few of the calculations are being done multiple times. May be we can restructure the code little bit to reduce computational time.

lruthotto · 2017-09-24T18:57:20Z

Thanks, Hari. This works for me! Great job !

lruthotto · 2017-09-24T19:05:33Z

Do you want me to merge this right now, or should I wait for performance improvements? I've also noticed that the matrix free code does not run for non mass-preserving problems since there is no tailored solver.

On a different note: Did you compare your code for computing the diagonal of the distance term with getDiagVolume in hyperElasticFEM.m ? Therein, I think we do not recompute things that are not needed. Might help in improving the efficiency for the objective function as well. One idea could be to pull this code out of the regularizer and call it from both functions. What do you think?

hariagr · 2017-09-25T08:53:06Z

Hi Lars,

One of the major change could be to store gradient of basis functions, i.e., variable "dphi" from the function getBasisGradient, in the MESH structure. These calculations does not depends on the current estimate "yc", therefore we can compute gradient once and store them in the MESH structure.
We can pull out cofactor3D function from both regularizer and distance/objective function, and call it wherever required.

I can mainly see these changes right now to improve efficiency.

I didn't check the code for non mass-preserving problems. I will try to look at it asap. There is no hurry from my side for merging the code. Probably, these changes should not take too much time. We can wait and merge once we are done with our changes.

Hari

lruthotto · 2017-09-25T13:02:16Z

Hi Hari,

storing the dphi in the mesh structure is a good idea I think. There is a way in MATLAB to have dependent fields (i.e., that are computed once and then stored). If we can get this to work, this is the way to go I think. I'm not sure if this actually works in our setup or we have to return the mesh object in every function call.
Do you mean extracting the getDiagVolume part or really the cofactor3D? In any case, collecting code that does the same thing will be beneficial in the long run. So we should do it.

Moved getBasisGradient and Cofactor3D functions in a separate script file. Gradient of basis functions are now computed only once in the entire function call, by moving these calculations to MESH structure main file. NOTE: There is a slight difference between estimates coming out from MatrixFree and Non-MatrixFree code, in the 3D case.

hariagr · 2017-09-27T12:07:42Z

Hi Lars,

I have changed the code. dphi is now stored in the mesh structure. cofactor3D code is separated out.

In the last comment I wrote that the Matrixfree code (MF) and non-matrifree (MB) code results are not matching. But, that is not right. They are matching!! By mistake, I used different regularization parameter for MB and MF code during verification.

Hari

Minor changes in FEMobjFctn.m for MatrixFree functionality and added a PCG solver function.

hariagr · 2017-10-17T21:44:20Z

Hi Lars,

In the MLIRFEM script, image pixel values on a coarse grid are computed at two steps, i.e.,

In the first step, the pixel values are defined on the cell-centered points of a finer grid ( i.e. a rectangular grid) and then restricted to a coarser grid by averaging adjacent cells (see getMultilevel.m).
Now, computed pixel values are estimated on barycentre points of triangular elements of FEM mesh.

Instead of this, we can probably define pixels values directly on the FEM triangular mesh and compute approximation on a coarse FEM mesh. Probably, this may results in a more accurate approximation.

I am not sure how much beneficial it would be in view of registration accuracy. Probably, you can comment on it.

Thanks,
Hari

herrinj · 2018-07-16T19:47:04Z

Checked out Hari's branch today. I have some updates to hyperElasticFEM.m to add. Specifically, I will push code that has matrix-based and matrix-free implementations with support for 2D and 3D problems. This includes function calls to extract the diagonal of the Hessian in the matrix case for both 2D and 3D. Will upload as soon as I have everything compatible with the branch.

lruthotto · 2018-07-25T16:08:56Z

Perfect. Thank you both for your help!

…

On Mon, Jul 16, 2018 at 3:47 PM, James Herring ***@***.***> wrote: Checked out Hari's branch today. I have some updates to hyperElasticFEM.m to add. Specifically, I will push code that has matrix-based and matrix-free implementations with support for 2D and 3D problems. This includes function calls to extract the diagonal of the Hessian in the matrix case for both 2D and 3D. Will upload as soon as I have everything compatible with the branch. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#5 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AEDYB75JIQ01wz3_CdeV7odVmnapEBjbks5uHO24gaJpZM4OTrbw> .

hariagr added 3 commits July 11, 2017 05:16

Matrix-free FEMAPP - SSD

e889f4d

These changes are required to use matrix-free functionality of hyper elastic regularization function in the add-on FEMAPP (Only for 2D). To verify the changes, the file “EFEM_SSD_MBvsMF.m” compares the results from MF and MB.

Matrixfree code - MP

2827a84

These are minor changes to run matrixfree functionality of FEMPIRE. But, data fidelity term is still not matrixfree. Prob, this involves little more work. Will come back to this later. Added a file to compare results from MF and MB.

hariagr and others added 4 commits July 17, 2017 15:53

Delete MLIRFEM.m~

4fa6cee

Delete GaussNewton.m~

dab3d1f

*.m~ ignore

6157cf3

Merge branch 'master' into master

ae2199d

lruthotto reviewed Jul 20, 2017

View reviewed changes

Improvements in the MatrixFree code

2d156f5

Matrixfree implementations of - derivative of determinants, - diagonal of hessian matrix, and - Hessian matrix

hariagr and others added 3 commits August 14, 2017 14:17

deleted files

4df61b1

two minor changes

4ce857a

another minor edit

5d236e9

Hari Om Aggraval added 3 commits September 3, 2017 09:22

Merge branch 'master' of https://github.com/hariagr/FAIR.m

1da7ca4

* 'master' of https://github.com/hariagr/FAIR.m: another minor edit two minor changes

3D - Matrixfree code ... in progress

6952a8b

Implement matrix free code to evaluate determinant. Added a script “verifyMF_getdeterminant” to verify results

MatrixFree code for 3D... in progress

f3e9d3f

Modified HyperElasticFEM function. Work to do: MatrixFree code for diagonal elements of Hessian matrix.

Hari Om Aggraval added 2 commits September 15, 2017 17:19

MatrixFree Functionality for non-mass preserving case

925f778

Minor changes in FEMobjFctn.m for MatrixFree functionality and added a PCG solver function.

herrinj mentioned this pull request Jul 17, 2018

hyperElasticFEM matrix-free 3D #8

Open

		@@ -0,0 +1,38 @@
		% ==================================================================================
		% (c) Lars Ruthotto 2012/07/26, see FAIR.2 and FAIRcopyright.m.

few minor changes to run example code in FAIRFEM #5

Are you sure you want to change the base?

few minor changes to run example code in FAIRFEM #5

Uh oh!

Conversation

hariagr commented Jul 11, 2017

Uh oh!

lruthotto commented Jul 15, 2017

Uh oh!

lruthotto Jul 20, 2017

Choose a reason for hiding this comment

Uh oh!

hariagr Jul 20, 2017

Choose a reason for hiding this comment

Uh oh!

lruthotto Jul 20, 2017

Choose a reason for hiding this comment

Uh oh!

hariagr Jul 20, 2017

Choose a reason for hiding this comment

Uh oh!

lruthotto Jul 20, 2017

Choose a reason for hiding this comment

Uh oh!

hariagr Jul 20, 2017

Choose a reason for hiding this comment

Uh oh!

hariagr Aug 13, 2017

Choose a reason for hiding this comment

Uh oh!

lruthotto commented Jul 20, 2017

Uh oh!

hariagr commented Aug 13, 2017

Uh oh!

lruthotto commented Aug 31, 2017

Uh oh!

lruthotto commented Sep 9, 2017

Uh oh!

lruthotto commented Sep 24, 2017

Uh oh!

lruthotto commented Sep 24, 2017

Uh oh!

hariagr commented Sep 25, 2017

Uh oh!

lruthotto commented Sep 25, 2017

Uh oh!

hariagr commented Sep 27, 2017

Uh oh!

hariagr commented Oct 17, 2017

Uh oh!

herrinj commented Jul 16, 2018

Uh oh!

lruthotto commented Jul 25, 2018 via email

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants