OpenCMISS / iron

Source code repository for OpenCMISS-Iron
9 stars 62 forks source link

Nonlinear fitting #148

Open PrasadBabarendaGamage opened 6 years ago

PrasadBabarendaGamage commented 6 years ago

Currently, only linear fitting is implemented:

equationsSetSpecification = [iron.EquationsSetClasses.FITTING,
                             iron.EquationsSetTypes.DATA_FITTING_EQUATION,
                             iron.EquationsSetSubtypes.DATA_POINT_VECTOR_STATIC_FITTING,
                             iron.EquationsSetFittingSmoothingTypes.SOBOLEV_VALUE]
problemSpecification = [iron.ProblemClasses.FITTING,
                        iron.ProblemTypes.DATA_FITTING,
                        iron.ProblemSubtypes.DATA_POINT_VECTOR_STATIC_FITTING]

We would like to implement nonlinear fitting.

Add a new problem subtype:

iron.EquationsSetSubtypes.DATA_POINT_VECTOR_STATIC_NONLINEAR_FITTING

Rename existing linear fitting as:

iron.EquationsSetSubtypes.DATA_POINT_VECTOR_STATIC_LINEAR_FITTING

Proposed implementation steps:

PrasadBabarendaGamage commented 6 years ago

Not too sure of the difference between the fitting problem types below? e.g. static fitting vs data point vector static fitting. Are some of these old - can some be removed?

SELECT CASE(problem%specification(3))
CASE(PROBLEM_STATIC_FITTING_SUBTYPE)
CASE(PROBLEM_STANDARD_DATA_FITTING_SUBTYPE)
CASE(PROBLEM_GENERALISED_DATA_FITTING_SUBTYPE)
CASE(PROBLEM_MAT_PROPERTIES_DATA_FITTING_SUBTYPE)
CASE(PROBLEM_DATA_POINT_VECTOR_STATIC_FITTING_SUBTYPE)
CASE(PROBLEM_DATA_POINT_VECTOR_QUASISTATIC_FITTING_SUBTYPE)
CASE(PROBLEM_VECTOR_DATA_FITTING_SUBTYPE)
CASE(PROBLEM_DIV_FREE_VECTOR_DATA_FITTING_SUBTYPE)
chrispbradley commented 6 years ago

Yes, there are a few old types that haven't been merged/removed. The better one to look at is

equationsSetSpecification = [iron.EquationsSetClasses.FITTING, iron.EquationsSetTypes.DATA_FITTING_EQUATION, iron.EquationsSetSubtypes.DATA_POINT_FITTING, iron.EquationsSetFittingSmoothingTypes.SOBOLEV_VALUE]

and

problemSpecification = [iron.ProblemClasses.FITTING, iron.ProblemTypes.DATA_FITTING, iron.ProblemSubtypes.STATIC_FITTING]

Your equations set will be a fibre fitting as the objective function is quite different. There will need to be the distinction between linear/non-linear etc. You will also need to have something that is specific to using the Diffusion tensor i.e. f^T.D.f

For the problem, you will definatly need to distinguish linear/non-linear static/dynamic. What the particular form of objective function (apart from linear/nonlinear) shouldn't matter.

PrasadBabarendaGamage commented 6 years ago

Thanks @chrispbradley. I'll add in the following:

equationsSetSpecification = [iron.EquationsSetClasses.FITTING, iron.EquationsSetTypes.FIBRE_FITTING_EQUATION, iron.EquationsSetSubtypes.FIBRE_POINT_NONLINEAR_FITTING, iron.EquationsSetFittingSmoothingTypes.SOBOLEV_VALUE]

and

problemSpecification = [iron.ProblemClasses.FITTING, iron.ProblemTypes.FIBRE_FITTING, iron.ProblemSubtypes.STATIC_NONLINEAR_FITTING]

Lets set this up to use the diffusion tensor by default (without specifying additional subtypes) and then discuss how to add in a subtype to select between different types of objective functions - will simplify/speed up implementation.

chrispbradley commented 6 years ago

Hi Prasad, I would have thought that the greater defining factor of this objective function is that it involves a diffusion tensor. Yes it is nonlinear but that is a consequence of the f^T.D.f rather than that the diffusion tensor is a consequence of the nolinearity. I would think something like

equationsSetSpecification = [iron.EquationsSetClasses.FITTING, iron.EquationsSetTypes.FIBRE_FITTING_EQUATION, iron.EquationsSetSubtypes.DIFFUSION_TENSOR_FIBRE_FITTING, iron.EquationsSetFittingSmoothingTypes.SOBOLEV_VALUE]

???

The problem specification of

problemSpecification = [iron.ProblemClasses.FITTING, iron.ProblemTypes.FIBRE_FITTING, iron.ProblemSubtypes.STATIC_NONLINEAR_FITTING]

seems fine as the nonlinearity is the defining feature of the solver workflow.

PrasadBabarendaGamage commented 6 years ago

Thanks Chris, yup sure, will use that subtype for the equationsSet.

chrispbradley commented 6 years ago

Hi Prasad, did you see https://github.com/OpenCMISS/iron/commit/3c7f4acdfac45fd4adfdf191c6b2e7fe9de46604 ? Just trying to give you enough time to change it before I comment to change it again ;-)

PrasadBabarendaGamage commented 6 years ago

Yup, lets do it as:

equationsSetSpecification = [iron.EquationsSetClasses.FITTING, iron.EquationsSetTypes.DATA_FITTING_EQUATION, iron.EquationsSetSubtypes.DIFFUSION_TENSOR_FIBRE_FITTING, iron.EquationsSetFittingSmoothingTypes.SOBOLEV_VALUE]

Does that mean the GAUSS fitting should also be moved under data fitting?

chrispbradley commented 6 years ago

Hi Prasad, no, Gauss and data point fitting are fundamentally different in that we need a switch to access the "values" we need from either the data point or Gauss point data structures. That is, we could easily define an objective function being the sum of the error over "points" for data or Gauss points. The objective is just f^t.D.f regardless of whether or not it is stored at a Gauss or data point. The example above is for data points but we could easily define

equationsSetSpecification = [iron.EquationsSetClasses.FITTING, iron.EquationsSetTypes.GAUSS_POINT_FITTING_EQUATION, iron.EquationsSetSubtypes.DIFFUSION_TENSOR_FIBRE_FITTING, iron.EquationsSetFittingSmoothingTypes.SOBOLEV_VALUE]

PrasadBabarendaGamage commented 6 years ago

Ah, ok, so when writing up the Equations set we would name it something like DIFFUSION_TENSOR_FIBRE_FITTING and inside there check the EquationsSetTypes to determine the input data (either Gauss point or data points) to prevent duplicating code?

PrasadBabarendaGamage commented 6 years ago

Actually, everything seems to be in Fitting_FiniteElementCalculate (currently it initialises a lot of linear matrices/mappings pointers). Should I add in the nonlinear matrices there too or rename it LinearFitting_FiniteElementCalculate, and add in a NonlinearFitting_FiniteElementCalculate?

chrispbradley commented 6 years ago

Depends what you need to do. If you are doing something that would be the same between data and Gauss but different between fitting objectives then look at the subtype. The type should split into gauss and data point subroutines for the actual equations and so then look at the subsubtype (or whatever it is called).

chrispbradley commented 6 years ago

The Fitting_FiniteElementCalculate will only be called if it is a linear (FEM) equation set. If the equations set is nonlinear then it will call ...ResidualEvaluate ...JacobianEvaluate etc.

PrasadBabarendaGamage commented 6 years ago

Ah of course! Its been a while : )

chrispbradley commented 6 years ago

Which probably do not exist for fitting at the moment and so will need to be added.

PrasadBabarendaGamage commented 6 years ago

Yup, will add them in. Thanks.

chrispbradley commented 6 years ago

Great