Finish implementation of ellipsoidal tube interpolation

[pgagarinov]

Deadline: end of this semester

Ellipsoidal tubes are difficult to interpolate using the convential methods because a few properties very specific to ellipsoidal tubes, especially tight ellipsoidal tubes have to be preserved. These properties are 1) E_1[t] belongs E_2[t] at all times 2) E_1[t] touches E_2[t] along l(t)

The property 1) is easy to satisfy even by using a linear interpolation. However given that both E_1 and E_2 tubes are known on a certain time grid {t_i} a linear interpolation cannot guarantee that an interpolated tube (which is based on linearly interpolated shape matrix of an ellipsoid) satisfies property 2). This is because EllTube.interp method is based on nearest interpolation which is not suitable for control-synthesis problems. The goal is to implement a proper interpolation that doesn't break both properties. Luckily such interpolation does exist and has been implemented at 70%. Thus we need just to finish the implementation glueing loose ends here and there. The interpolation is based on ODE solution via Runge-Kutta 45 algorithm by Dormand Prince (see articles in https://drive.google.com/open?id=0B5OqL4IVOIC8elN4MHpfbHhTWFU). The idea is that Dormand Prince suggested an interpolation that uses the same coefficients as ODE solver so that this interpolation has the same level of precision as ODE solution. Actually this interpolation is already used by ode45 solver. According to Mathworks Dormand and Prince have sent this interpolation via an email. There is no any theoretical proof that I was able to find though but this is the best that we have. The fact that this interpolation is implemented in an official ode45 solver gives it a certain level of credability.

Anyways, branch "issue_10" contains a partial implementation of the task. I suggest you create a copy of this branch and merge your changes back to issue_10 fromtime to time. Kirill Aksenov will be helping your with this task by implementing a comparison functionality for matrix functions so he will also be merging to issue_10. Once issue_10 is ready you will merge it to master branch. You are responsible for issue_10, not Kirill but your task depend on his work.

Step I: study what has been already done, namely

1) A few classes in gras.ode package. These classes allow to generate a special interpolation object from ode45 solver so that this object can be used for interpolating solution of any ODE.

products+gras+ode+test+mlunit\SuiteBasic.m products+gras+ode+test+mlunit\SuiteOde45Reg.m products+gras+ode\IMatrixSysInterp.m products+gras+ode\MatrixODE45InterpFunc.m products+gras+ode\MatrixODE45ScalarTimeInterpFunc.m products+gras+ode\MatrixSysNearestInterp.m products+gras+ode\MatrixSysODERegInterpSolver.m products+gras+ode\MatrixSysODESolver.m products+gras+ode\MatrixSysReshapeOde45RegInterp.m products+gras+ode\VecOde45RegInterp.m products+gras+ode\ode45reg.m

You need to understand what these classes do by studying tests for them in SuiteBasic and SuiteOde45Reg test classes.

2) Changes to MatrixFunctions in

products+gras+mat products+gras+mat+fcnlib

packages. Matrix function object is an entity that is responsible for a delayed calculation (symbolic). The idea of an ellipsoidal tube interpolation is that instead of calculating ellipsoidal matrices at certain time grid {t_i} we "remember" the formulas themselves in a symbolic form. This way we can then calculate the shape matrix at any given t, not just at time moments from {t_i}. For instance, if you need to multiply matrix A by vector b you write A*b but in our case we do something different, we do opFactory=gras.mat.CompositeMatrixOperations();

aMatrixFunc=... (comes from somewhere, isa(aMatrixFunc,'gras.mat.IMatrixFunction')=true bVecFunc=...(comes from somewhere , isa(bVecFunc,'gras.mat.IMatrixFunction')=true)

multResultFunc=opFactory.rMultiplyByVec(aMatrixFunc, bVecFunc);

then we can calculate a result of Ab on a certain time grid timeVec just by calling multResultFunc.evaluate(timeVec). So all those new matrix function classes added in issue_10 branch to gras.mat and gras.mat.fcnlib package are responsible for different kinds of operations (like , /, indexing etc).

3) Changes to ellipsoidal tube calculation classes.

products+gras+ellapx+lreachuncert+probdyn\LReachProblemLTIDynamics.m (this class and many others are responsible for calculating things like A(t), B(t), B(t)P(t)B'(t), C(t)Q(t)C'(t), B(t)p(t), C(t)q(t) as well as a center of ellipsoidal tubes that are solutions to \dot(x)=A(t)x(t)+B(t)p(t), C(t)q(t). In master branch this solution to this ODE is obtained on a certain time grid and then interpolated via a cubic spline. In issue_10 LReachProblemLTIDynamics has been changed to use an interpolated object returned by ODE solver itself. This way cube interpolation is replaced by the interpolation by Dormand-Prince that has the same level of precision as the solution itself. You need to change LReachProblemDynamicsInterp so that it uses the same approach as LReachProblemLTIDynamics. Also please avoid copy-pasting, if you can place the common functionality in LReachProblemDynamicsInterp and LReachProblemLTIDynamics in a new base class - please do so.

products+gras+ellapx+lreachplain\AGoodDirs.m (a base class for calculating a transition matrix X(s,t) and l(t)) products+gras+ellapx+lreachplain\GoodDirsContinuousGen.m (implementations of X(t,s) and l(t) calculations for a generic continuous time systems products+gras+ellapx+lreachplain\GoodDirsContinuousLTI.m (same for time-independent systems) products+gras+ellapx+lreachplain\GoodDirsDiscrete.m (same for discrete-time systems)

The changes to these classes are very similar to the changes to LReachProblemLTIDynamics - a spline-based interpolation is either replaced or complemented with Dormand-Prince interpolation.

An important thing to remember here is that GoodDirs classes do not operate with X(t,s) directly, instead they calculate an normalized X(t,s) denoted as R(s,t) (R(s,t)=X(s,t)/||X(s,t)||). R is suitable for calculating good directions l(t) to the same extent as X because for good directions their norm doesn't matter. But at the same time R has a few advantages - its norm is kept around 1 which gives l(t) with norm also around 1 (more or less). If we use X then for certain A a norm of l_s(t) = X(s,t)'l_s can grow or drop down to zero with an exponential speed.

products+gras+ellapx+lreachuncert\ExtIntEllApxBuilder.m (this class is responsible for calculating both external and internal approximations and building ellipsoidal tubes. In issue_10 branch tube constructors require a few additional inputs associated with an interpolation objects).

For instance, if previously it was sufficient to provide QArray:double[nDims,nDims,nTime](- Q%28t%29), aMat:double[nDims,nTimes](- a%28t%29) along with a few additional inputs to build E(t)=E(Q(t),a(t)) tube, now we also need interpolation objects for both Q(t) and a(t) so that we can calculate E(t) at any arbitrary t \in [t_0,t_1].

Changes to

products+gras+ellapx+smartdb+rels\EllTube.m products+gras+ellapx+smartdb+rels\EllTubeBasic.m

introduced in issue_10 reflect that.

Step II: Finish the implementation.

1) Finish changes to problem dynamics classes (see above). Right now only a class for LTI system has been changed. 2) Finish changes to GoodDirs classes. Changes have been introduced by none of the changes have been thoroughly tested. There are a few tests in products+gras+ellapx+lreachplain+test but you need to extend the coverage by verifying these classes one a few simple linear systems and comparing the results with ethalon results (calculated analytically on a piece of paper).

3) Finish changes to EllTube clases in gras.ellapx.smartdb.rels classes and make sure the tests in gras.ellapx.smartdb.test pass. A potential challenge is a comparison of ellipsoidal tubes. When tube objects do not contain any interpolation objects as properties (like in master branch) comparing tubes is easy and this is successfully done. However in issue_10 branch ellipsoidal tubes contain interpolation object that have dependencies on a lot of interpolation objects. Thus we need to be able to compare interpolation objects with a certain precision. For that purpose we need to be able to compare any object imlpemented in gras.mat and gras.mat.fcnlib package.

I suggest we approach this problem in two steps: a) implementing a stub comparison code that can recognize as equal objects that are 100% equal and not equal otherwise. This can be done via inheriting gras.mat.AMatrixFunctionComparable from modgen.common.obj.HandleObjectCloner object which provides all necessary comparison functionality. All you need to do is overriding isEqualScalarInternal in AMatrixFunctionComparable so that it calls isEqualScalarAsBlobInternal method. This way all matrix function objects will be compared as BLOBs which should be sufficient for fixing all existing tests in gras.ellapx.smartdb.test package.

b) A correct implementation of matrix function comparison will require providing a matrix function class-specific implementation of isEqualScalarInternal method. This will be done by Kirill Aksenov.

Please note that some of the tests that fail in issue_10 right now are because of getRGoodDirCurveInterpObj method missing in AGoodDirs. (R stands for R(s,t) - normalized transition matrix, see above). Instead of getRGoodDirCurveInterpObj you need to use getRGoodDirOneCurveInterpList that you need to implement. Here is an explanation why and how.

fromQArraysInternal in products+gras+ellapx+smartdb+rels\EllTubeBasic.m generates interpolation objects for good curves using subarray method of CompositeMatrixOperations factory:

line 331: STubeData.ltGoodDirInterpObjList{iLDir} = ... compositeMatrixOperationsObj.subarray(... ltGoodDirInterObj,[{1:nRows} {iLDir}]);

This needs to be changed by implementing a method AGoodDirections in analogous to getRGoodDirOneCurveSplineList but returning a list of interpolation objects for each l_i. This method needs to be called

getRGoodDirOneCurveInterpList and EllTubeBasic.fromQArraysInternal should accept a list of these interpolation objects returned by getRGoodDirOneCurveInterpList method instead of a single ltGoodDirInterObj. Once this is done the transformation performed on line 331 won't be needed and STubeData.ltGoodDirInterpObjList can be taken directly from the input.

AGoodDirs.getRGoodDirOneCurveInterpList should be implemented based directly on getRstTransInterpObj method that returns an interpolation objects for R(s,t)' - normalized X(s,t)' matrix. A normalized l(t) could then be taken directly by multiplying R(s,t)' by l_i. This can be done via gras.mat.CompositeMatrixOperations.rMultiplyByVec.

4) Re-cache regression tests in gras.ellapx.uncertcalc.test package. Most of the tests from this package use the calculation results stored in mat files as ethalons Since you have a new structure of EllTube classes comparison won't work so you need to re-cache the results (update ethalons on disk).

Assuming you are fixed all the bugs in ExtIntEllApxBuilder the test re-caching can be done via

gras.ellapx.uncertcalc.test.regr.run_regr_tests('reCache',true) gras.ellapx.uncertcalc.test.regr.run_regr_tests('reCache',true,'nParallelProcesses',8) (for 4-core processor with hyper-threading like i7 4790k for instance)

Note: cache files are located in \products+gras+ellapx+uncertcalc+test+regr+mlunit\TestData\SuiteRegression\testRegression_out

5) Enable the following tests for elltool.reach.ReachContinuos class that are currently disabled because "interp" method of EllTube class doesn't work correctly (it uses a nearest interpolation while it should provide a "true" interpolation).

products+elltool+reach+test+mlunit\ContinuousIsEqualTestCase.m

DISABLED_testIsEqualEnclosedTimeVecs DISABLED_testIsEqualNotEnclosedTimeVecs

products+elltool+reach+test+mlunit\ContinuousReachTestCase.m

DISABLE_testEvolve DISABLE_testCut

Once the test are disable you can make sure that all tests for ReachContinuos pass via running elltool.reach.test.run_cont_tests

6) Remove matrix function that are not used anywhere (some were added just in case and we might not need them any longer). An example of such matrix function is the one build via gras.mat.CompositeMatrixOperations.subarrayInit

[thelastpride]

If I understand correctly, first I need to download a copy of the branch 10 and make a rebase to the latest master. The question is, in favour of which to resolve conflicts?

[pgagarinov]

Correct.

Conflict resolution is something you need to do manually and very carefully line-by-line, there is no simple rule like "always prefer master" or "always prefer the branch" because they both contain valuable changes. If there were such a simple rule the conflict resolution process would be automated but it can't be.

The plan should be a) studying all the changes made in the branch and in master and making sure you understand them. Without understanding what is changed and why you won't be able to resolve conflicts. In many cases conflicts can only be resolved by merging via editing code manually. It is not that hard if you keep the goal in mind - the goal is to have master on top of which all changes made in the branch are applied. What you should not do is to always prefer the branch because the branch can contain mistakes, it contains unfinished implementation and it is out of date.

Instead of resolving conflicts you can just create a new branch and manually introduce the changes made in the branch, it is safer but can be more time-consuming, the choice is yours.

[thelastpride]

Where can I learn what exactly do classes LReachProblemDynamicsInterp and LReachProblemLTIDynamics and what is the difference between them? I could not find any help on the subject in ET.

[pgagarinov]

The source code is your main source because you are a developer, not an end-user in this case. These classes as we can guess from their names represent a dynamics of linear control systems for different kinds of A(t). LReachProblemDynamicsInterp corresponds to a generic A(t) which is interpolated on a certain time grid while LReachProblemLTIDynamics corresponds to A(t)=A and uses a matrix exponent to calculate transition matrices, l(t), B(t)*P(t)B(t)' and some other stuff used in ellipsoidal approximations. Just go through the code of these classes, look at the interfaces and everything will be clear.

[thelastpride]

We Kirill spent a lot of time to rebase. As a result, we found that the branch issue_10 itself does not pass a very large number of tests (for example gras.test.run_tests gived us << FAILED >> || TESTS: 167, FAILURES: 9, ERRORS: 29). Does is ean that we first need to fix the branch, and then make a rebase, or we can try to first make a rebase, and then try to fix errors?

[pgagarinov]

I think it is safer to finish re-basing first and then take it from there. It is only natural to have that many failures just because implementation is not finished. However after rebase you need to make sure that you did not revert any useful changes made in master or introduced some misprints that broke the code. Finally, when you rebase it is always easier to squash all commits made in the branch into a single commit (two commits into one) and only then rebase. Such an approach can lower a number of conflict resolutions.

[thelastpride]

Next question is about +gras+ellapx+lreachplain +gras+ellapx+lreachuncert Am I right that all huncert classes inherited from plain classes? Could you explain a little bit what they are doing? I'm sorry, but I have not found any documentation on this subject, and code review did not help me.

[pgagarinov]

Am I right that all huncert classes inherited from plain classes? Yes

Could you explain a little bit what they are doing?

These classes belong to 3 categories: solvers (like lreachuncert.ExtIntEllApxBuilder), problem definition classes (like LReachContProblemDef) that are just containers for system parameters and problem dynamics classes that are responsible for calculating B_P_B', X(t,s), C_Q_C' and other auxilary stuff used for ellipsoidal approximations.

The topology of the toolbox is such that ReachContinuous and ReachDiscrete classes are just high-level wrapper around classes from gras.ellapx. The latter do most of calculation work, especially solver classes. Actually, we have a separate command-line tool for calculating ellipsoidal tubes by using solvers from gras.ellapx directly:

gras.ellapx.uncertcalc.listconf 'advanced' 'advancednodisp' 'default' 'ellDemo3test' 'test' 'test2d' 'test2dbad' 'test2dinternal' 'test2dnodisp' 'test3d' 'uosc8' 'uosc8full' 'utest2d' 'x2dtest'

gras.ellapx.uncertcalc.editconf('advanced ') [SRunProp, SRunAuxProp]=gras.ellapx.uncertcalc.run('advanced');

The last call returns a structure with ellipsoidal tubes (relations) where tuples represent individual ellipsoidal tubes.

I think the best way for you to gain a deeper understanding of the classes in gras.ellapx is to follow a source code of gras.ellapx.uncertcalc.run. Just run it as specified above and follow a code execution line by line in a debugger, it is actually very simple. Then compare the way "run" calls the solver classes and ReachContinuous calls them. You will see a lot of similarities.

The reason for having low-level classes is functionality separate, the high level ReachContinous sometimes easier to use because it provides ease-to-use methods and allows for plotting reach sets (as a union or intersection of ellipsoidal tubes) while gras.ellapx.uncertcalc.run does something much more simple - just calculates a set of ellipsoidal tubes with specific parameters but does this task very well.

Also please have a look at the tests and examples in gras.ellapx.smartdb.test - they will give you a good understanding how relations (tables) of ellipsoidal tubes work.

[thelastpride]

I`m sorry, there is one more commit. I accidentally called him "Issue #10" and message about it now here https://github.com/SystemAnalysisDpt-CMC-MSU/ellipsoids/issues/10

[thelastpride]

2) Finish changes to GoodDirs classes. Changes have been introduced by none of the changes have been thoroughly tested. There are a few tests in products+gras+ellapx+lreachplain+test but you need to extend the coverage by verifying these classes one a few simple linear systems and comparing the results with ethalon results (calculated analytically on a piece of paper).

In +gras+ellapx+lreachplain+test there are 4 examples for simple systems. Am I right that GoodDirsTestCase now testing only process of creation classes without checking them?

[pgagarinov]

Yes, you are right, this is one of the reasons this task is unfinished. As mentioned in the task description (see above) you need to enhance the test coverage for GoodDirs classes.

[thelastpride]

I have some questions about 3)

Finish changes to EllTube clases in gras.ellapx.smartdb.rels

What exactly I need to do? I could not find information about changes to products+gras+ellapx+smartdb+rels\EllTube.m products+gras+ellapx+smartdb+rels\EllTubeBasic.m in issue_10.

A potential challenge is a comparison of ellipsoidal tubes.

What should I do, if Kirill had already completed his method? Where should I paste the code? In tests or classes in? Could you give me an example of a test of the toolbox, which fails due to a problem of comparison matrices?

[pgagarinov]

What exactly I need to do? I could not find information about changes to

Well, it is overly optimistic to expect a documentation for changed made to the source code. But I understand that you might get a bit frustrated when you look at what has been changed without knowing anything why. So let me explain the idea. Initially EllTube-like classes were designed to represet an ellipsoidal tube on a fixed time grid. These classes are like tables or relations of ellipsoidal tube properties (tube are stored in tuples and different properties of the tubes - in fields or columns). Internally inside EllTube-like object we have SData structure with a separate field for each property. There are also SIsValueNull and SIsNull structure for storing NULL information for certain properties but this is not relevant in the context of this task. You can understand the way these classes work by studying tests in gras.ellapx.smartdb.rels.test.SuiteEllTube test case in master. Just run various tests and see how stuff works.

Now about the changes. The fact that the tube information is only stored on a fixed grid is not good for for control synthesis (as I explained on a seminar). This is because for control synthesis we need to know E_{-}[t](internal ellipsoidal approximation) at any t, not just t_i from timeVec (which is a field on EllTube and other classes). To solve this problem we decided to take the following approach: introduce additional fields called "QArrayInterpObjList" (should be called QArrayInterpObj i.e without "List", please rename it) that represents an interpolation object for Q(t). There should an obvious consistency i.e. QArrayInterpObj.evaluate(t_i) should be equal to QArray(:,:,timeVec(i)). We also introduce similar objects for both a(t) (consistent with aMat) and l(t) (consistent with ltGoodDirMat). These properties complement QArray, aMat and ltGoodDirMat fields and are used (at least shoudl be used) in "interp" method which is a key method which doesn't work as needed even right now in master. In master this method is implemented via "nearest" interpolation which means that if we call interp with certain t a closest t_i from timeVec is found and ellipsoidal properties for t_i are returned(i.e. Q(t_i),q()).

What should I do, if Kirill had already completed his method?

isEqual in Ellipsoidal tubes is implemented in such way that "isEqual" for interp objects (that are fields of EllTube class) are called automatically. By default "isEqual" for interp objects compares handles (i.e. points that are different in different objects). This is the reason all the tests in SuiteEllTube fail right now in your branch:

(object arrays #1 and #2):(element #1):(SData):(1).QArrayInterpObjList-->{1}values are different (1).ltGoodDirInterpObjList-->{1}values are different

Once Kirill overrides isEqual comparison will be perfomred via Kirill's method automatically, no changes are needed on your part. But meanwhile you need to have a default implementation that performs a comparison with an absolute precision (as described in a)). This can be done by comparing interp objects as BLOBs. You can either ask Kirill to provide a temporary implementation that does the following:

            leftBlobVec=getByteStreamFromArray(leftObj);
            rightBlobVec=getByteStreamFromArray(rightObj);
            isOk=isequal(leftBlobVec,rightBlobVec);

This will set isOk to true only if leftObj and rightObj are absolutely equal and to false otherwise.

But I suggest you make a temporary hack on your end so that you can work in parallel in the same branch.

Go to +gras+ellapx+smartdb+rels\ATypifiedAdjustedRel.m and in isEqualScalarInternal method add the following lines at the beginning of the method:

    import modgen.common.obj.ObjectComparisonMode;
    %

    relRef=self;
    interpObjectFieldNameList={'QArrayInterpObj',<list all interp object field names here>};
    nInterpFields=numel(interpObjectFieldNameList);
    isInterpFieldPresentVec=relRef.isField(interpObjectFieldNameList);
    for iInterpField=1:nInterpFields
        if isInterpFieldPresentVec(iInterpField)
            interpFieldName=interpObjectFieldNameList{iInterpField};
            relRef.(interpFieldName).setComparisonMode(ObjectComparisonMode.BLOB);
            end
        end
    end
    %
    %% IMPORTANT - now do the same for otherRel input of isEqualScalarInternal as well. i.each
    relRef=otherRel;
    %% execute the same code as above.

This I believe should do it - the comparison will be performed based on BLOB representation for interp fields. Once Kirill finishes his part you just remove this code.

So using my explanation please study the changes to EllTube-like classes carefully line by line (via Git log) and then figure out where they are either incorrect or not finished. To do that you need to study internals of EllTube-like classes by looking at their source code and tests.

[thelastpride]

Thank you very much for the detailed explanation. I have a few questions about this piece of code import modgen.common.obj.ObjectComparisonMode;

   %

    relRef=self;
    interpObjectFieldNameList={'QArrayInterpObj',<list all interp object field names here>};
    nInterpFields=numel(interpObjectFieldNameList);
    isInterpFieldPresentVec=relRef.isField(interpObjectFieldNameList);
    for iInterpField=1:nInterpFields
        if isInterpFieldPresentVec(iInterpField)
            interpFieldName=interpObjectFieldNameList{iInterpField};
            relRef.(interpFieldName).setComparisonMode(ObjectComparisonMode.BLOB);
            end
        end
    end
    %
    %% IMPORTANT - now do the same for otherRel input of isEqualScalarInternal as well. i.each
    relRef=otherRel;
    %% execute the same code as above.

I added this code to isEqualScalarInternal but after executing this line relRef.(interpFieldName).setComparisonMode(ObjectComparisonMode.Blob); an error occurs. Struct contents reference from a non-struct array object. I do not understand what it means and where it comes from.

[pgagarinov]

I suspect that this is because you need to write self.SData.(interpFieldName) instead of self.(interpFieldName). If it fails again - use debugger to stop execution on these lines and see which line fails and why. I just gave you a hint and my code snippet can sometime contain misprints. In such cases you are expected to understand the idea and make the corrections all by yourself, not to ask me - "hey, you code is wrong, fix it please for me." I hope you understand. Thanks.

[thelastpride]

A normalized l(t) could then be taken directly by multiplying R(s,t)' by l_i. This can be done via gras.mat.CompositeMatrixOperations.rMultiplyByVec.

Am I right, that I need to multiply RstTransInterpObj by

ltGoodDirOneCurveSplineList

for each element of list?

[pgagarinov]

A direction of your though is right but the fact that you asked this question means that you didn't understand the global idea behind this task. Please make sure you understand the following.

The whole idea of "interp" objects is to eliminate an imprecision of interpolation (that is already used in master branch). And there are things that we just cannot interpolate. For instance for two ellipsoidal tubes E{-}[t] and E{-}[t] calculated by ode45 (via some higher-level wrapper classes that support solving matrix-valued ODEs) we calculate (in master) these tubes on a fixed time grid timeVec=[t_1,... ,t_i, ... tn]; This solution has all the good properties of ellipsoidal estimates: both E{-}[ti] and E{+}[t_i] touches a reachability domain for i, both are included in each other for each i and both via union(intersection) given an exact representation of the reachability domain (in theory, in practice depends on a number of l_k).

BUT, once you start to deviate from ti BY !!!CUBIC SPLINE INTERPOLATION!!! all these properties are lost. This is because an interpolated solution at t* is not the same as the exact solution at t* (as if you included t* into timeVec at the beginning). This is very inconvenient because you cannot perform a precise control synthesis as E{-}[t] can be no longer inside E{+}[t]. All in all, everything breaks at t* just because cubic splines do not know anything about ODE. Thus we need a different approach which is to use the method implemented by ode45 itself. This method is to use ODE-aware interpolation. ode45 uses this interpolation (out of the box, implemented by Mathworks) to calculate E{-}[t*] via the following steps:

a) calculate E_{-}[t^{~}_k] where t^{~}k are knots for a adaptive time grid b) use special ODE-specific coefficients and special Runge-Kutta-based interpolation schema to calculate E{-}[t_i]

If you use this algorithm E_{-}[ti] has the same level of precision as E{-}[t^{~}_k] and properties of ellipsoidal approximation DO NOT BREAK at t_i even if t_i doesn't match any of t^{~}_k.

To make this process easier we implemented ode45reg that returns a special interpolation object that can later be used as an ordinary cubic spline-based interpolation object to calculate interpolated values via "evaluate" method. Around ode45reg we implemented a few higher-level wrappers that make things much easier for matrix-valued systems.

Now, back to your question. The fact that you sugget to combine ODE-aware interpolation object with cubic spline-based interpolation object breaks the whole idea described above. The rule of thumb is - apples to apples, splines to splines, ODE-aware interpolants to ODE-aware interpolants. This you should use ltGoodDirOneCurveInterpList instead of ltGoodDirOneCurveSplineList.

There is only one exception in this rule - we can still use splines to approximate values of A(t), C(t) etc. that define our systems because this doesn't break anything, it just lowers a precision which we can incresea indefinitely just by increasing a number of cube spline knots. What everything that happens after A(t) should use EITHER splines (for backward-compatibility and high-performance mode when we want to know our ellipsoidal approximation only at t_i ) OR ODE-aware interpolants. And you never mix them with each other.

[thelastpride]

After working with the debugger, I realized that the parameter setComparisonMode(ObjectComparisonMode.Blob) can be set only for the entire object as a whole, but not for an individual of his field.

testProjectAdv(self) calls isEqual for ATypifiedAdjustedRel. I tried to set the comparison mode as varargin{indObj}.setComparisonMode(ObjectComparisonMode.Blob); in ATypifiedAdjustedRel.isEqual After that HandleObjectCloner sets flags as case ObjectComparisonMode.Blob isAsBlobDefault=true; isAsHandleDefault=false; But the comparison is still not working correctly. (object arrays #1 and #2):blob sizes are different in testProjectAdv(self) rel1Proj=rel.projectStatic(projMat); rel2Proj=rel.projectStatic({projMat}); [isPos,reportStr]=rel1Proj.isEqual(rel2Proj); But `isequal(getByteStreamFromArray(rel1Proj),getByteStreamFromArray(rel2Proj))

ans =

 1`

So I think that this temporary solution will not work. Maybe I'm wrong, and is necessary set the mode setComparisonMode(ObjectComparisonMode.Blob) for QArrayInterpObjList, aMatInterpObj and ltGoodDirInterpObjList separately?

[pgagarinov]

This will work, just bear with me for a second. Let us consider the following example:

SData=struct('QArrayInterp',{{1;2}},'aMatInter',{{3;4}}); rel=smartdb.relations.DynamicRelation(SData);

rel.('QArrayInterp') or rel.QArrayInterp returns a cell vector. The is the case in EllTube relation - QArrayInterp is not a field of EllTube object in a normal sense, it is a dynamic property that points to a specific field inside of self.SData. Anyways, what is inside this field is MatrixFunction object that Kirill needs to inherit from HandleObjectCloner. Thus self.(QArrayInterp) is a cell array of objects inherited from HandleObjectCloner. My example needs to be modified a bit so that instead of

relRef.(interpFieldName).setComparisonMode(ObjectComparisonMode.BLOB);

we use

cellfun(@(x)relRef.(interpFieldName),@(x)x.setComparisonMode(ObjectComparisonMode.BLOB));

setComparisonMode(ObjectComparisonMode.Blob) can be set only for the entire object as a whole, >but not for an individual of his field.

This is absolutely not true, you can call setComparisonMode for any object inherited from HandleObjectCloner even when this object is a field or is an element of a field of another object. If self.QArrayInterp is a cell array of MatrixFunction classes inherited from HandleObjectCloner then you can call setComparisonMode.

The idea is that isEqual method of EllTube classes is smart enough to compare its fields recursively so if you change the way the fields are compared (by calling setComparisonMode) the behavior of EllTube.isEqual will change automatically so that "interp" fields are compared as BLOBs while other fields will be compared as before.

Does this make sense?

[thelastpride]

Thank you, I will try to do so.