search

xjiang4 / ellipsoids

Automatically exported from code.google.com/p/ellipsoids
Other
0 stars 1 forks source link

Implement an interpolation technique based on ODE solution #126

Open GoogleCodeExporter opened 8 years ago

GoogleCodeExporter commented 8 years ago 
Currently gras.ode.ode45reg uses ntrp45 function to interpolation ODE solution 
for an arbitrary set of potions based on adaptive step-size solution. We need 
to implement an interpolation class that uses the same logic as in ntrp45 to 
build ODE solution on an arbitrary time grid. This class should implement 
gras.mat.IMatrixFunction interface and its instance should be returned by 
gras.ode.ode45reg function. The returned instance should be stored within 
EllTube objects for use in interp method (currently this method uses nearest 
point interpolation). Once this is done testCut and testEvolve methods in  
elltool.reach.test.mlunit.ContinuousReachTestCase will start to pass again. 
Right now they are disabled.

Original issue reported on code.google.com by heartofm...@gmail.com on 18 Aug 2013 at 2:49

GoogleCodeExporter commented 8 years ago 
[deleted comment]
GoogleCodeExporter commented 8 years ago 
1) Study gras.ode.ode45reg function. 
The function is a modified version of ode45 Runge-Kutta solver which solves a 
regularized ode of the following form:

\dot(y)(t)=F(y,t)+m(t) where y(t) and m(t) are N-dimensional vectors. This ODE 
is defined in domain defined by fRegOde function. An example
of such domain is all(eig(Y(t))>0) where Y(t) is a matrix composed from y(t) 
vector in the following way:

[y_1(t),...,y_n(t); y_n_1(t),...,y_2n(t);....;.... y_n*n(t)], N=n*n;

basically this corresponds to matrix equation

\dot(Y(t))=F(Y(t),t)+M(t) transformed into \dot(y)(t)=F(y,t)+m(t).

The idea is to solve such ODE choosing M(t) dynamically such that Y(t) stays 
within all(eig(Y(t))>epsilon) domain.

Line 186 

[isStrictViol,yInterimNewRegVec]=fOdeReg(tStep,yInterimNewVec);

in ode45reg calls fOdeReg function for each step. For the example above 
isStrictViol would be true 

all(eig(Y(t))>0)==false and false otherwise. yInterimNewRegVec = yInterimNewVec 
+h * M(t) where M(t) is chosen such that

all(eig(yInterimNewRegVec)>eps)=true.

Thus basically fOdeReg performs a dynamic regularization. However, if y(t) is 
within the system definition domain 
then fOdeReg does nothing and yInterimNewRegVec= yInterimNewVec.

To better understand how this function works you need to study the tests in 
gras.ode.test.mlunit.SuiteBasic (see below),
especially testNNDedTriuSpline test.

a) Study the tests located in gras.ode.test.mlunit.SuiteBasic, particularly 
testReg test.

To run all the tests use

res=gras.ode.test.run_tests

Also you can use 

res=mlunitext.runtestcase('gras.ode.test.mlunit.SuiteBasic','testParams',{'gras.
ode.ode45reg','ode45','marker','ode45'})

to run all tests for ode45reg

and

res=mlunitext.runtestcase('gras.ode.test.mlunit.SuiteBasic','testReg','testParam
s',{'gras.ode.ode45reg','ode45','marker','ode45'})

You can use 

res.getErrorFailMessage()

to see the failures

to run a specific test (testReg in this case).

b) Study the following part of the code in ode45reg function.

    nsteps = nsteps + 1;
    switch outputAt
        case 'SolverSteps'        % computed points, no refinement
            nout_new = 1;
            tout_new = tnew;
            yout_new = ynew;
        case 'RefinedSteps'       % computed points, with refinement
            tref = t + (tnew-t)*S;
            nout_new = refine;
            tout_new = [tref, tnew];
            yout_new = [ntrp45(tref,t,y,h,f,fOdeReg), ynew];
        case 'RequestedPoints'    % output only at tspan points
            nout_new =  0;
            tout_new = [];
            yout_new = [];
            while next <= ntspan
                if tnew < tspan(next)
                    break;
                end
                nout_new = nout_new + 1;
                tout_new = [tout_new, tspan(next)];
                if tspan(next) == tnew
                    yout_new = [yout_new, ynew];
                else
                    yout_new = [yout_new, ntrp45(tspan(next),t,y,h,f,fOdeReg)];
                end
                next = next + 1;
            end
    end

This part performs an interpolation of ODE solution calculated on adaptive 
grid. The interpolation is done dynamically
in the process of solving ODE. For instance, you ask ode45reg to calculate a 
solution on [1,2,3,4,5] time grid
while the solver find a solution on [1,0.1,0.3332,0.56,..,1.99,2.233,..., 
0.456] time grid. The question is - how do we find y(t) at
t=1,2 etc. The answer is in ntrp45 which uses f(..) values corresponding to 
different Runge-Kutta steps at last time
preceeding the requested time (for t=2 the last preceeding time would be 1.99 
for the example above) and uses Runge-Kutta formulas
to find y(t) at t=2. Moreover y(2) would have the same order of precision.

2) Once you do 1) you need to modify ode45reg in a such way that basically 
stores all data necessary for ntrp45 to 
calculate ODE solution at any arbitrary grid.

Instead of using [t,yVec,dyRegMat]=ode45reg call a user would use

[t,yVec,dyRegMat,interpObj]=ode45reg call where interpObj corresponds to an 
interpolation object.

A user could then use interpObj to calculate ode solution at any requested 
point vector timeVec
 using interpObj.evaluate(timeVec) call. 

interpObj should be of gras.ode.VecOde45RegInterp class. The class should have 
evaluate method that would be used like this:

[yMat,yRegMat]=interpObj.evaluate(timeVec) where yMat would contain the y(t) at 
requested points in columns and 
yRegMat - m(t).

Internally interpObj would keep all the derivatives necessary for calculating 
ODE solution at any requested time grid. 

You would have to write a test that verifies that your implementation is 
correct. I suggest the following test:

a) call 
[t,yMat,dyRegMat,interpObj]=ode45reg(fOdeDeriv,fOdeReg,timeVec,y0Vec,...) once 
for some fOdeDeriv, fOdeReg,timeVec and y0Vec.
b) fixInterpObj=interpObj
c) vary timeVec by considering different time vectors such that timeVec(1) and 
timeVec(end) remain unchanged. For each 
timeVec call [t,yMat,dyRegMat,~]=ode45reg(fOdeDeriv,fOdeReg,timeVec,y0Vec,...) 
for the same fOdeDeriv, fOdeReg and y0Vec.
and compare yMat and yRegMat with 
[yInterpMat,yRegInterpMat]=fixInterpObj.evaluate(timeVec)

3) Study gras.ode.MatrixSysODESolver class and the tests for it located in 
gras.interp.test.mlunit.SuiteBasic (testMatrixSysODESolver test)

You need to slightly modify this class so that it it returns varargout 
arguments that exceed nEqs * nFuncs directly from the solver function
stored in fSolveFunc field of the class. New functionality needs to be covered 
by a simple test.

4) Implement gras.ode.MatrixSysODEInterpSolver class inherited from 
gras.ode.MatrixSysODEInterpSolver. This class 
should use the last output of solve method of base class (which will be 
interpObj) and create an instance of new
gras.ode.MatrixSysReshapeOde45RegInterp. This class in its constructor should 
just accept an object of 
gras.ode.VecOde45RegInterp class (returned
by MatrixSysODESolver) along with reshape sizes, known to 
MatrixSysODEInterpSolver, store this object along with size vectors
 somewhere in its private fields and in its evaluate method just do reshape to transform y(t) and m(t) vectors
into Y_i(t) and M_i(t) matrices for both outputs of 
gras.ode.VecOde45RegInterp.evaluate method (here i stands 
for a number of matrix equation within a system of matrix equations. When we 
build both internal and external approximations
in parallel we find a solution to a system of two equations 
\dot{X}_{-}(t)=f_{-}(X_{-}(t),t) and \dot{X}_{+}(t)=f_{+}(X_{-}(t),t)
in parallel, where X_{-}(t) stands for a configuration matrix of internal 
ellipsoid while X_{+}(t) - for a configuration matrix of
external ellipsoid.)

Thus, for example for a system of 2 matrix equations a call of 
gras.ode.MatrixSysODEInterpSolver.solve method would return 

[timeVec,intApxArr,extApxArr,intRegArr,extRegArr,interpObj]=matrixSysOdeInterpSo
lvObj.evaluate(timeVec);

where interpObj is an object of gras.ode.MatrixSysReshapeOde45RegInterp class 
(which you need to implement)

Calling evaluate method for interpObj object 

[interpIntApxArr,interpExtApxArr,interpIntRegArr,interpExtRegArr]=interpObj.eval
uate(timeVec)

you would get interpIntApxArr,interpExtApxArr,interpIntRegArr,interpExtRegArr 
that are equal to 

intApxArr,extApxArr,intRegArr,extRegArr

However the main idea of having interpObj is that you can call evaluate with 
much dense timeVec without actually solving ODE.

5) Study the code of the following methods in 
gras.ellapx.lreachuncert.ExtIntEllApxBuilder class

calcRegEllApxMatrix
calcEllApxMatrixDeriv

These methods correspond to fOdeDeriv,fOdeReg inputs of 
gras.ode.MatrixSysReshapeOde45RegInterp class and
define ODE system for X_{-}(t) and X_{+}(t) (internal and external ellipsoidal 
approximations).  

You need to write a test for MatrixSysReshapeOde45RegInterp based on this 
system and make sure that for the 
interpolated X_{-}(t) and X_{+}(t) the property X_{-}(t)<= X_{+}(t) and 
\rho(l(t)|X_{-}(t))=\rho(l(t)|X_{+}(t))
holds.

Original comment by heartofm...@gmail.com on 16 Sep 2013 at 5:25

GoogleCodeExporter commented 8 years ago 
Step 5) updated - instead of writing the test for 
MatrixSysReshapeOde45RegInterp directly you need to update the following 
classes 

gras.ellapx.lreachuncert.ExtIntEllApxBuilder (line 208)
gras.ellapx.lreachplain.ATightEllApxBuilder (line 75)

so that the interpObj object returned by solved is saved as an additional field 
in EllTube relation. The relations are created via 

gras.ellapx.smartdb.rels.EllTube.fromQArrays
and
gras.ellapx.smartdb.rels.EllTube.fromQMArrays 

methods. These methods need to be refactored so that they accept a list of 
interpObj objects (for each tube) as an additional argument. 

these objects should be stored into a relation as "tubeInterpObj" field that 
needs to be defined in gras.ellapx.smartdb.F class. Once you define it as 
TUBE_INTERP_OBJ you should add the field declaration int 
gras.ellapx.smartdb.rels.EllTubeBasic as

FCODE_TUBE_INTERP_OBJ

Once all this is done you need to make getInterpDataInternal method in 
EllTubeBasic class use interpObj objects from the new field as opposed to using 
the "nearest" interpolation as right now.  As the final accord enable testCut 
and testEvolve tests in  elltool.reach.test.mlunit.ContinuousReachTestCase and 
make sure that they pass.

Original comment by heartofm...@gmail.com on 20 Sep 2013 at 11:00

GoogleCodeExporter commented 8 years ago

Original comment by vadimdan...@gmail.com on 27 Sep 2013 at 6:07

GoogleCodeExporter commented 8 years ago 
1) SuiteBasic.testInterp would be called twice from gras.ode.test.run_tests, 
one time for ode45reg and one time for ode113reg. The second call would fail so 
we need to create a separate test case "SuiteOde45Reg" next to SuiteBasic, move 
the test there and modify gras.ode.test.run_tests to call the tests from the 
new test case class.

2) The test code produce code checking warnings (ttt and nTimePoints variables 
seem to be unused). Please make sure that you fix all such problems and see a 
green square in the right upper corner of editor window.

3) absTol=relTol=1e-14 is too low. What I mean is using 1e-14 for comparing the 
solutions, not for using such a low value for ode45reg. Please use 1e-8 for 
finding the solution and 1e-14 for comparing interpObj and ode45reg

4) Please run "run_tests" as the final step to make sure that all tests from 
gras.ode package pass

Original comment by heartofm...@gmail.com on 30 Sep 2013 at 12:04

GoogleCodeExporter commented 8 years ago

Original comment by vadimdan...@gmail.com on 1 Oct 2013 at 4:49

GoogleCodeExporter commented 8 years ago 
Please make SuiteOde45Reg test case you've created inherit SuiteBasic so that 
we have just two test cases - SuiteBasic for ode113reg and SuiteOde45Reg for 
ode45reg. All tests from SuiteBasic would be inherited by SuiteOde45Reg. 
run_tests would then call just 2 test cases, not 3.

Original comment by heartofm...@gmail.com on 1 Oct 2013 at 5:12

GoogleCodeExporter commented 8 years ago 
Also, please

1) Move abstol into compare function by removing abstol as an input argument. 
Also, please rename it into COMPARE_TOL constant.
2) Rename absTol into ABS_TOL as this is a constant.
3) Remove copy-pasting inside check function by created a nested function 
inside compare.
4) Are you sure that by doing 

                tVaryVec = sin(pi*tBeginVec/2);
                tVaryVec(2) = tBeginVec(2);

you get monotonous values in tVaryVec?

Original comment by heartofm...@gmail.com on 2 Oct 2013 at 7:04

GoogleCodeExporter commented 8 years ago

Original comment by vadimdan...@gmail.com on 4 Oct 2013 at 7:49

GoogleCodeExporter commented 8 years ago

Original comment by heartofm...@gmail.com on 5 Oct 2013 at 4:09

GoogleCodeExporter commented 8 years ago

Original comment by vadimdan...@gmail.com on 7 Oct 2013 at 5:24

GoogleCodeExporter commented 8 years ago

Original comment by heartofm...@gmail.com on 8 Oct 2013 at 1:12

GoogleCodeExporter commented 8 years ago

Original comment by vadimdan...@gmail.com on 8 Oct 2013 at 1:46

GoogleCodeExporter commented 8 years ago 
A few minor problems and questions.

1) The following size desciption can't be true, it rather char[1,] than 
char[1,1] which would correspond to 1 symbol.

39  +   % dataType: char[1,1] - data type, result of the
    40  +   % function superiorfloat()

2) nIter - can we change this name to "nKnots" to highlight the fact that 
VecOde45RegInterp doesn't know anything
about where SData came from?

3)  47  +   % iNext: double[1,1] - next iteration index - for it just doesn't 
make sense to pass this into
the constructor and iNext changes within the class itself. This is rather an 
internal property of the class than
an input parameter. Would you agree?

4) Please change nEquations into nSysDims to stress that we have a single 
vectorial equation and nSysDims is its dimensionality
5) You still use "i" as a counter which is a bad name as "i" is Matlab built-in 
constant for imaginary unit.

6) "tCVec: cell[1,nIter] - nodes of adaptive time grid" - your type 
specification doesn't tell what is inside of those
cells. You need to write cell[1,nIter] of <type>[<size>] where "of" can be 
specified more than once if you have nested cells.
7) What is the difference between "new adaptive gird" and just "adaptive" grid" 
and what is their role in interpolation? It should be
clear from their description.

8) "fCVec: cell[1,nIter] - values ??of the derivatives" what do you have 
question marks here?

Original comment by heartofm...@gmail.com on 8 Oct 2013 at 2:59

GoogleCodeExporter commented 8 years ago

Original comment by vadimdan...@gmail.com on 9 Oct 2013 at 4:41

GoogleCodeExporter commented 8 years ago 
Good job, please strictly follow the reintegration algorithm from 
http://code.google.com/p/ellipsoids/wiki/Issues_workflow_and_branching_policy 
to reintegrate.

Original comment by heartofm...@gmail.com on 10 Oct 2013 at 4:17

GoogleCodeExporter commented 8 years ago 
I think you can reintegrate now

Original comment by heartofm...@gmail.com on 16 Oct 2013 at 6:04

GoogleCodeExporter commented 8 years ago

Original comment by vadimdan...@gmail.com on 17 Oct 2013 at 1:48

GoogleCodeExporter commented 8 years ago

Original comment by heartofm...@gmail.com on 21 Oct 2013 at 7:42

GoogleCodeExporter commented 8 years ago

Original comment by vadimdan...@gmail.com on 24 Oct 2013 at 3:47

GoogleCodeExporter commented 8 years ago

Original comment by heartofm...@gmail.com on 26 Oct 2013 at 6:14

GoogleCodeExporter commented 8 years ago 
About the interface.

MatrixSysReshapeOde45RegInterp class you already implemented is designed to 
perform interpolation using ODE solution. But sometimes (in the tests for 
EllTube methods) you might not have any ODE to provide you with 
MatrixSysReshapeOde45RegInterp object. Thus you need to provide an ability to 
construct MatrixSysReshapeOde45RegInterp -like object with the same interface 
which requires just matrix data at certain points and doesn't require ODE.

I suggest that you implement gras.ode.MatrixSysNearestInterp class that 
implements gras.ode.IMatrixSysInterp interface. The interface should contain 
just one abstract method - evaluate. MatrixSysReshapeOde45RegInterp should be 
made to inherit the same interface. Constructor of MatrixSysNearestInterp class 
should accept a list of arrays of matrix values at timeVec grid for each 
equation/function. For 2 functions and 2 equations the list would contain 4 
elements with nDims x nDims x nTimePoints arrays in each. Internally evaluate 
method would perform the nearest interpolation as in interpInternal method.

Moreover, in addition to QArray and MArray you would have to interpolate aMat 
(centers) and ltGoodDirMat (good directions). For that you could still use 
MatrixSysODERegInterpSolver which would replace MatrixODESolver classes in 
gras.ellapx.lreachplain.GoodDirsContinuousGen (for good directions) and in 
Dynamics objects in +gras\+ellapx\+lreachplain\+probdyn and 
+gras\+ellapx\+lreachuncert\+probdyn packages (Anastasia knows about these 
objects, if you can ask her as well)

Original comment by heartofm...@gmail.com on 1 Nov 2013 at 3:13

GoogleCodeExporter commented 8 years ago 
Please delete your old branch

Original comment by heartofm...@gmail.com on 12 Nov 2013 at 10:18

GoogleCodeExporter commented 8 years ago

Original comment by vadimdan...@gmail.com on 24 Nov 2013 at 1:33

GoogleCodeExporter commented 8 years ago 
Replied in email

Original comment by heartofm...@gmail.com on 27 Nov 2013 at 1:57

GoogleCodeExporter commented 8 years ago

Original comment by vadimdan...@gmail.com on 3 Dec 2013 at 5:32

GoogleCodeExporter commented 8 years ago 
Details were sent via email

Original comment by heartofm...@gmail.com on 25 Dec 2013 at 4:59