Finish implementation of ellipsoidal tube interpolation #39
Comments
|
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? |
|
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. |
|
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. |
|
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. |
|
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? |
|
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. |
|
Next question is about |
These classes belong to 3 categories: solvers (like lreachuncert.ExtIntEllApxBuilder), 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
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 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. |
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? |
|
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. |
|
I have some questions about 3)
What exactly I need to do? I could not find information about changes to
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? |
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
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 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 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 This I believe should do it - the comparison will be performed based on BLOB representation for interp fields. 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. |
|
Thank you very much for the detailed explanation. I have a few questions about this piece of code I added this code to isEqualScalarInternal but after executing this line |
|
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. |
Am I right, that I need to multiply RstTransInterpObj by for each element of list? |
|
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, ... t_n]; This solution has all the good properties of ellipsoidal estimates: both E_{-}[t_i] 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 t_i 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^{ If you use this algorithm E_{-}[t_i] has the same level of precision as E_{-}[t^{ To make this process easier we implemented ode45reg that returns a special interpolation object that can later be used 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 |
|
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.
ans = So I think that this temporary solution will not work. |
|
This will work, just bear with me for a second. Let us consider the following example:
rel.('QArrayInterp') or rel.QArrayInterp returns a cell vector. The is the case in EllTube relation - relRef.(interpFieldName).setComparisonMode(ObjectComparisonMode.BLOB); we use cellfun(@(x)relRef.(interpFieldName),@(x)x.setComparisonMode(ObjectComparisonMode.BLOB));
This is absolutely not true, you can call The idea is that Does this make sense? |
|
Thank you, I will try to do so. |
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
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
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.
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 A*b 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).
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.
Finish changes to problem dynamics classes (see above). Right now only a class for LTI system has been changed.
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).
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.
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
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
An example of such matrix function is the one build via gras.mat.CompositeMatrixOperations.subarrayInit
The text was updated successfully, but these errors were encountered: