KerMor  0.9
Model order reduction for nonlinear dynamical systems and nonlinear approximation
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ACoreFun.m
Go to the documentation of this file.
1 namespace dscomponents{
2 
3 
4 /* (Autoinserted by mtoc++)
5  * This source code has been filtered by the mtoc++ executable,
6  * which generates code that can be processed by the doxygen documentation tool.
7  *
8  * On the other hand, it can neither be interpreted by MATLAB, nor can it be compiled with a C++ compiler.
9  * Except for the comments, the function bodies of your M-file functions are untouched.
10  * Consequently, the FILTER_SOURCE_FILES doxygen switch (default in our Doxyfile.template) will produce
11  * attached source files that are highly readable by humans.
12  *
13  * Additionally, links in the doxygen generated documentation to the source code of functions and class members refer to
14  * the correct locations in the source code browser.
15  * However, the line numbers most likely do not correspond to the line numbers in the original MATLAB source files.
16  */
17 
18 class ACoreFun
19  :public KerMorObject,
20  public general.AProjectable {
82  public: /* ( setObservable ) */
83 
106  public: /* ( setObservable ) */
107 
127  .sparse<logical> JSparsityPattern = "[]";
151  integer xDim = "[]";
171  integer fDim = "[]";
191  public:
192 
206  public: /* ( Transient ) */
207 
226 
228 
229 
230  public: /* ( Transient ) */
231 
232 
233  ACoreFun(sys) {
234  this = this@KerMorObject;
235  this.registerProps(" CustomProjection ",...
236  " MultiArgumentEvaluations "," JSparsityPattern ",...
237  " TimeDependent "," xDim "," fDim ");
238  this.setSystem(sys);
239  }
240 
241 
242  function setSystem(sys) {
243  this.System= sys;
244  }
245 
246 
247  function target = project(V,W,target) {
248  if nargin < 4
249  if this.CustomProjection
250  error(" The target parameter must be given if custom projection is used. ");
251  end
252  target = this.clone;
253  end
254  target = project@general.AProjectable(this, V, W, target);
255  /* New state space argument size is the number of columns of V. */
256  if isempty(V)
257  target.xDim= this.xDim;
258  else
259  target.xDim= size(V,2);
260  end
261  /* New state space output size is the number of columns of W. */
262  if isempty(W)
263  target.fDim= this.fDim;
264  else
265  target.fDim= size(W,2);
266  end
267  }
296  function fx = evaluate(x,t) {
297  cp = this.CustomProjection;
298  if ~cp && ~isempty(this.V)
299  x = this.V*x;
300  end
301  fx = this.evaluateCoreFun(x, t);
302  if ~cp && ~isempty(this.V)
303  fx = this.W^t*fx;
304  end
305  }
335  function fx = evaluateMulti(colvec<double> x,double t,colvec<double> mu) {
336  singlemu = size(mu,2) == 1;
337  if isempty(mu)
338  mu = [];
339  singlemu = true;
340  end
341  if singlemu
342  this.prepareSimulation(mu);
343  end
344  if isempty(t)
345  t = double.empty(0,size(x,2));
346  elseif length(t) == 1
347  t = ones(1,size(x,2))*t;
348  end
349  m = size(x,2);
350  if m == 0
351  m = length(t);
352  end
353  fx = zeros(this.fDim, m);
354  for idx = 1:size(x,2)
355  if ~singlemu
356  this.prepareSimulation(mu(:, idx));
357  end
358  fx(:,idx) = this.evaluate(x(:,idx), t(:,idx));
359  end
360  }
381  this.mu= mu;
382  }
395  function J = getStateJacobian(x,t) {
396  if ~isempty(this.V)
397  x = this.V*x;
398  end
399  J = this.getStateJacobianImpl(x, t);
400  if ~isempty(this.V)
401  J = this.W^t*(J*this.V);
402  end
403  }
424  function J = getStateJacobianImpl(colvec<double> x,double t) {
425  J = this.getStateJacobianFD(x, t);
426  }
438  function copy = clone(copy) {
439  if nargin == 1 || ~isa(copy," dscomponents.ACoreFun ")
440  error(" Incorrect call to clone. As this class is abstract, a subclass of ACoreFun has to be passed as second argument. ");
441  end
442  copy = clone@general.AProjectable(this, copy);
443  /* Copy local properties */
444  copy.CustomProjection= this.CustomProjection;
445  /* Dont clone the associated system */
446  copy.System= this.System;
447  copy.JSparsityPattern= this.JSparsityPattern;
448  copy.TimeDependent= this.TimeDependent;
449  copy.xDim= this.xDim;
450  copy.fDim= this.fDim;
451  }
452 
453 
454  public: /* ( Sealed ) */ /* ( Transient ) */
455 
456  function [matrix<double>J , dx ] = getStateJacobianFD(x,t,rowvec<integer> partidx) {
457  d = size(x,1);
458  if nargin < 4
459  partidx = 1:this.xDim;
460  end
461  len = length(partidx);
462  X = repmat(x,1,len);
463  dx = ones(d,1)*sqrt(eps(class(x))).*max(abs(x),1);
464 
465  DX = sparse(1:d,1:d,dx,d,d);
466  del = 1:this.xDim;
467  del(partidx) = [];
468  DX(:,del) = [];
469 
470  /* Evaluate makes use of built-in multi-argument evaluation */
471  J = (this.evaluateMulti(X+DX,t,this.mu) - repmat(this.evaluate(x,t),1,len))*diag(1./dx);
472  /* Create sparse matrix if pattern is set */
473  if ~isempty(this.JSparsityPattern) && nargin < 4
474  [i, j] = find(this.JSparsityPattern);
475  J = sparse(i,j,J(logical(this.JSparsityPattern)),this.fDim,this.xDim);
476  end
477  }
504  public: /* ( Abstract ) */ /* ( Transient ) */
505 
506  virtual function fx = evaluateCoreFun(colvec<double> x,double t) = 0;
526  /* % Getter & Setter */
527  public: /* ( Transient ) */
528 
529 
530 #if 0 //mtoc++: 'set.CustomProjection'
531 function CustomProjection(value) {
532  if ~islogical(value)
533  error(" Property must be logical/boolean. either true or false ");
534  end
535  this.CustomProjection= value;
536  }
537 
538 #endif
539 
540 
541 
542 #if 0 //mtoc++: 'set.JSparsityPattern'
543 function JSparsityPattern(value) {
544  if ~isempty(value) && ~issparse(value)
545  error(" JSparsityPattern must be a sparse matrix. ");
546  end
547  this.JSparsityPattern= value;
548  }
549 
550 #endif
551 
552 
553 
554 #if 0 //mtoc++: 'set.xDim'
555 function xDim(value) {
556  if ~isempty(value) && ~isposintscalar(value)
557  error(" xDim must be a positive integer. ");
558  end
559  this.xDim= value;
560  }
561 
562 #endif
563 
564 
565 
566 #if 0 //mtoc++: 'set.fDim'
567 function fDim(value) {
568  if ~isempty(value) && ~isposintscalar(value)
569  error(" fDim must be a positive integer. ");
570  end
571  this.fDim= value;
572  }
573 
574 #endif
575 
576 
577  function res = test_MultiArgEval(mudim) {
578  if nargin < 2
579  mudim = 100;
580  end
581  x = rand(this.xDim,200);
582  mui = rand(mudim,200);
583  fxm = this.evaluateMulti(x,1:200,mui);
584  fxs = zeros(size(x));
585  for idx = 1:size(x,2)
586  this.prepareSimulation(mui(:, idx));
587  fxs(:,idx) = this.evaluate(x(:,idx), idx);
588  end
589  err = sum((fxm-fxs).^2,1);
590  res = all(err < eps);
591  if ~res
592  plot(err);
593  end
594  }
605  function logicalres = test_Jacobian(matrix<double> xa,rowvec<double> ta,matrix<double> mua) {
606 
607  if nargin < 2
608  xa = rand(this.xDim,20);
609  ta = 1:20;
610  mua = rand(20,20);
611  else
612  if nargin < 4 || isempty(mua)
613  mua = double.empty(0,size(xa,2));
614  if nargin < 3 || isempty(ta)
615  ta = double.empty(0,size(xa,2));
616  end
617  end
618  end
619  d = size(xa,1);
620  perstep = floor((256*1024^2)/(8*d)); /* 256MB chunks */
621 
622  steps = ceil(d/perstep);
623  gpi = steps == 1 && size(xa,2) > 1;
624  if gpi
625  pi = ProcessIndicator(" Comparing %d %dx%d jacobians with finite differences ",...
626  size(xa,2),false,size(xa,2),this.fDim,this.xDim);
627  end
628  if size(mua,2) == 1
629  mua = repmat(mua,1,size(xa,2));
630  end
631  res = true;
632  JSP = this.JSparsityPattern;
633  haspattern = ~isempty(JSP);
634  zero = ~JSP;
635  for i = 1:size(xa,2)
636  x = xa(:,i); t = ta(:,i);
637  this.prepareSimulation(mua(:,i))
638  Jc = this.getStateJacobian(x, t);
639 
640  /* % Numerical jacobian */
641  if ~gpi && steps > 1
642  pi = ProcessIndicator(" Comparing %dx%d jacobian with finite differences over %d blocks of size %d ",...
643  steps,false,this.fDim,this.xDim,steps,min(d,perstep));
644  end
645  for k = 1:steps
646  pos = (k-1)*perstep+1:min(d,k*perstep);
647  J = this.getStateJacobianFD(x, t, pos);
648  if steps > 1
649  pi.step;
650  end
651 
652  pm = PlotManager;
653  pm.LeaveOpen= true;
654  if haspattern && any(J(zero(:,pos)))
655  warning(" Jacobian Pattern mismatch between given and finite-difference version! ");
656  ax = pm.nextPlot(," Nonzero positions ");/* #ok */
657 
658  errpos = zero(:,pos) & logical(J);
659  spy(errpos);
660  [rows,columns] = find(errpos)/* #ok */
661 
662  end
663 
664  iszero = J == 0;
665  diff = abs(J-Jc(:,pos));
666  reldiff = abs((J-Jc(:,pos))./J);
667  abserr = max(diff(:));
668  if KerMor.App.Verbose > 1
669  di = diff(~iszero);
670  r = reldiff(~iszero);
671  [r, idx] = sort(r," descend ");
672  di = di(idx);
673  pos = find(r < 1e-2,1);
674  if ~isempty(pos) && pos > 1
675  ax = pm.nextPlot(,...
676  sprintf(" Relative errors larger than 1e-2\nCorresponding absolute difference "));/* #ok */
677 
678  ax = plotyy(1:pos,r(1:pos),1:pos,di(1:pos)," plot ");
679  axis(ax," tight ");
680  end
681 
682  pm2 = PlotManager(false,1,2);
683  pm2.LeaveOpen= true;
684  ax = pm2.nextPlot(," FD-Jacobian ");/* #ok */
685 
686  spy(J);
687  ax = pm2.nextPlot(," Analytic Jacobian ");/* #ok */
688 
689  spy(Jc);
690  ax = pm.nextPlot(,...
691  sprintf(" Points greater than ten percent\nof the max abs error %g ",abserr));/* #ok */
692 
693  spy(diff>abserr/10);
694  pm.done;
695  pm2.done;
696  end
697 
698 /* reltol = 1e-5;
699  *
700  * % Continue with pointwise relative error
701  * pointwisereldiff = abs(diff./J);
702  * % Remove those where the computed finite difference
703  * % jacobian is exactly zero
704  * pointwisereldiff(iszero) = 0;
705  * % Remove difference values, for which the finite
706  * % difference jacobian has values smaller than
707  * % 5e-7 and the corresponding values of the
708  * % given jacobian are even 4 orders of magnitude
709  * % smaller
710  * hlp = (abs(Jc) < 1e-4*abs(J)) & (abs(J) < 1e-6);
711  * pointwisereldiff(hlp) = 0;
712  *
713  * [pointwisereldiff, idx] = sort(pointwisereldiff(:),'descend');
714  * num = length(find(pointwisereldiff > reltol));
715  * maxpwrelerr = pointwisereldiff(1);
716  *
717  * dispnum = min(20,num);
718  * % dispnum = num;
719  * [posi,posj] = ind2sub(size(J), idx(1:dispnum));
720  * indic = sprintf('%d, ',posi);
721  * indjc = sprintf('%d, ',posj);
722  * vals = sprintf('%g, ',pointwisereldiff(1:dispnum));
723  * Jvals = sprintf('%g, ',J(idx(1:dispnum)));
724  * Jcvals = sprintf('%g, ',full(Jc(idx(1:dispnum))));
725  * maxJval = max(abs(J(idx(1:num))));
726  * % Check if/where the pointwise relative error is greater
727  * % than the tolerance
728  * if maxpwrelerr >= reltol;
729  * if abserr > 1e-6
730  * res = false;
731  * fprintf('Failed. Max absolute error %g >= %g\n',abserr,1e-6);
732  * elseif maxJval > reltol
733  * res = false;
734  * fprintf('Failed. Some points with relative error >= %g have jacobian absolute values >= %g \n',reltol,maxJval);
735  * end
736  * fprintf(2,'Max absolute error %g, max relative %g, %d are >= %g pointwise tolerance (all <= %g <= %g).\n',abserr,maxpwrelerr,num,reltol,maxJval, reltol);
737  * fprintf('Max errors: %s\nJFD: %s\nJ: %s\nRows: %s\nCols: %s\n',vals,Jvals,Jcvals,indic,indjc);
738  * end */
739  end
740  if ~gpi && steps > 1
741  pi.stop;
742  end
743  if gpi
744  pi.step;
745  end
746  end
747  if gpi
748  pi.stop;
749  end
750  }
771  protected: /* ( Static ) */ /* ( Transient ) */
772 
773  static function obj = loadobj(obj,from) {
774  if nargin == 2
775  obj.CustomProjection= from.CustomProjection;
776  obj.JSparsityPattern= from.JSparsityPattern;
777  obj.TimeDependent= from.TimeDependent;
778  obj.xDim= from.xDim;
779  obj.fDim= from.fDim;
780  obj = loadobj@general.AProjectable(obj, from);
781  obj = loadobj@KerMorObject(obj, from);
782  else
783  obj = loadobj@general.AProjectable(obj);
784  obj = loadobj@KerMorObject(obj);
785  end
786  }
787 
805 };
806 }
807 
808 
809 
integer fDim
The current output dimension of the function.
Definition: ACoreFun.m:171
A matlab matrix.
CustomProjection
Set this property if the projection process is customized by overriding the default project method...
Definition: ACoreFun.m:108
function J = getStateJacobian(x, t)
Default implementation of jacobian matrix evaluation via finite differences.
Definition: ACoreFun.m:395
Interface for all components that can be projected.
Definition: AProjectable.m:18
function prepareSimulation(colvec< double > mu)
A method that allows parameter-dependent computations to be performed before a simulation using this ...
Definition: ACoreFun.m:380
logical TimeDependent
Flag that indicates if the ACoreFun is (truly) time-dependent.
Definition: ACoreFun.m:84
Base class for any KerMor class.
Definition: KerMorObject.m:17
sort
ort the handle objects in any array in ascending or descending order.
function registerProps(varargin)
Call this method at any class that defines DPCM observed properties.
Definition: DPCMObject.m:125
function res = test_MultiArgEval(mudim)
Convenience function that tests if a custom MultiArgumentEvaluation works as if called with single ar...
Definition: ACoreFun.m:577
An integer value.
PlotManager: Small class that allows the same plots generated by some script to be either organized a...
Definition: PlotManager.m:17
A boolean value.
Basic interface for all dynamical system's core functions Inherits the AProjectable interface...
Definition: ACoreFun.m:18
logical LeaveOpen
Flag indicating if the plots should be left open once the PlotManager is deleted (as variable) ...
Definition: PlotManager.m:213
KerMorObject()
Constructs a new KerMor object.
Definition: KerMorObject.m:55
integer xDim
The current state space dimension of the function's argument .
Definition: ACoreFun.m:151
virtual function fx = evaluateCoreFun(colvec< double > x,double t)
Actual method used to evaluate the dynamical sytems' core function.
models.BaseFirstOrderSystem System
The system associated with the current ACoreFun.
Definition: ACoreFun.m:193
static function obj = loadobj(obj, from)
Definition: ACoreFun.m:773
V
The matrix of the biorthogonal pair .
Definition: AProjectable.m:61
colvec< double > mu
The current model parameter mu for evaluations. Will not be persisted as only valid for runtime durin...
Definition: ACoreFun.m:208
#define X(i, j)
function copy = clone(copy)
The interface method with returns a copy of the current class instance.
Definition: ACoreFun.m:438
function fx = evaluate(x, t)
Evaluates the f-approximation. Depending on a possible projection and the CustomProjection-property t...
Definition: ACoreFun.m:296
Speed test * all(1:3)
sparse< logical > JSparsityPattern
Sparsity pattern for the jacobian matrix.
Definition: ACoreFun.m:127
function fx = evaluateMulti(colvec< double > x,double t,colvec< double > mu)
Evaluates this function on multiple locations and maybe multiple times and parameters.
Definition: ACoreFun.m:335
function res = isposintscalar(value)
isposintscalar: Backwards-compatibility function for matlab versions greater than 2012a ...
function [ matrix< double > J , dx ] = getStateJacobianFD(x, t,rowvec< integer > partidx)
Implementation of jacobian matrix evaluation via finite differences.
Definition: ACoreFun.m:456
Global configuration class for all KerMor run-time settings.
Definition: KerMor.m:17
function logical res = test_Jacobian(matrix< double > xa,rowvec< double > ta,matrix< double > mua)
Tests the custom provided jacobian matrix against the default finite difference computed one...
Definition: ACoreFun.m:605
static function KerMor theinstance = App()
The singleton KerMor instance.
Definition: KerMor.m:910
W
The matrix of the biorthogonal pair .
Definition: AProjectable.m:72
function J = getStateJacobianImpl(colvec< double > x,double t)
Default implementation of state jacobians. uses finite differences.
Definition: ACoreFun.m:424
Base class for all KerMor first-order dynamical systems.
function target = project(V, W, target)
Sets the protected matrices that can be utilized on core function evaluations after projection...
Definition: ACoreFun.m:247
ProcessIndicator: A simple class that indicates process either via waitbar or text output...
function setSystem(sys)
Definition: ACoreFun.m:242