Menu
▾
▴
octave-cvsupdate — CVS update notifications
This list is closed, nobody may subscribe to it.
| 2001 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
|
Jul
(10) |
Aug
(5) |
Sep
(3) |
Oct
(41) |
Nov
(41) |
Dec
(33) |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 2002 |
Jan
(75) |
Feb
(10) |
Mar
(170) |
Apr
(174) |
May
(66) |
Jun
(11) |
Jul
(10) |
Aug
(44) |
Sep
(73) |
Oct
(28) |
Nov
(139) |
Dec
(52) |
| 2003 |
Jan
(35) |
Feb
(93) |
Mar
(62) |
Apr
(10) |
May
(55) |
Jun
(70) |
Jul
(37) |
Aug
(16) |
Sep
(56) |
Oct
(31) |
Nov
(57) |
Dec
(83) |
| 2004 |
Jan
(85) |
Feb
(67) |
Mar
(27) |
Apr
(37) |
May
(75) |
Jun
(85) |
Jul
(160) |
Aug
(68) |
Sep
(104) |
Oct
(25) |
Nov
(39) |
Dec
(23) |
| 2005 |
Jan
(10) |
Feb
(45) |
Mar
(43) |
Apr
(19) |
May
(108) |
Jun
(31) |
Jul
(41) |
Aug
(23) |
Sep
(65) |
Oct
(58) |
Nov
(44) |
Dec
(54) |
| 2006 |
Jan
(96) |
Feb
(27) |
Mar
(69) |
Apr
(59) |
May
(67) |
Jun
(35) |
Jul
(13) |
Aug
(461) |
Sep
(160) |
Oct
(399) |
Nov
(32) |
Dec
(72) |
| 2007 |
Jan
(316) |
Feb
(305) |
Mar
(318) |
Apr
(54) |
May
(194) |
Jun
(173) |
Jul
(282) |
Aug
(91) |
Sep
(227) |
Oct
(365) |
Nov
(168) |
Dec
(18) |
| 2008 |
Jan
(71) |
Feb
(111) |
Mar
(155) |
Apr
(173) |
May
(70) |
Jun
(67) |
Jul
(55) |
Aug
(83) |
Sep
(32) |
Oct
(68) |
Nov
(80) |
Dec
(29) |
| 2009 |
Jan
(46) |
Feb
(18) |
Mar
(95) |
Apr
(76) |
May
(140) |
Jun
(98) |
Jul
(84) |
Aug
(123) |
Sep
(94) |
Oct
(131) |
Nov
(142) |
Dec
(125) |
| 2010 |
Jan
(128) |
Feb
(158) |
Mar
(172) |
Apr
(134) |
May
(94) |
Jun
(84) |
Jul
(32) |
Aug
(127) |
Sep
(167) |
Oct
(109) |
Nov
(69) |
Dec
(78) |
| 2011 |
Jan
(39) |
Feb
(58) |
Mar
(52) |
Apr
(47) |
May
(56) |
Jun
(76) |
Jul
(55) |
Aug
(54) |
Sep
(165) |
Oct
(255) |
Nov
(328) |
Dec
(263) |
| 2012 |
Jan
(82) |
Feb
(147) |
Mar
(400) |
Apr
(216) |
May
(209) |
Jun
(160) |
Jul
(86) |
Aug
(141) |
Sep
(156) |
Oct
(6) |
Nov
|
Dec
|
| 2015 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
|
Jul
(1) |
Aug
|
Sep
(1) |
Oct
|
Nov
(1) |
Dec
(2) |
| 2016 |
Jan
|
Feb
(2) |
Mar
(2) |
Apr
(1) |
May
(1) |
Jun
(2) |
Jul
(1) |
Aug
(1) |
Sep
|
Oct
|
Nov
(1) |
Dec
|
| 2019 |
Jan
|
Feb
|
Mar
(1) |
Apr
|
May
|
Jun
|
Jul
|
Aug
|
Sep
|
Oct
|
Nov
|
Dec
|
| 2021 |
Jan
|
Feb
|
Mar
|
Apr
(3) |
May
(4) |
Jun
(8) |
Jul
(2) |
Aug
(5) |
Sep
(9) |
Oct
|
Nov
|
Dec
|
|
From: <par...@us...> - 2012-05-07 12:40:00
|
Revision: 10373
http://octave.svn.sourceforge.net/octave/?rev=10373&view=rev
Author: paramaniac
Date: 2012-05-07 12:39:54 +0000 (Mon, 07 May 2012)
Log Message:
-----------
control-devel: fix arx initial conditions
Modified Paths:
--------------
trunk/octave-forge/extra/control-devel/inst/arx.m
Modified: trunk/octave-forge/extra/control-devel/inst/arx.m
===================================================================
--- trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-07 04:39:42 UTC (rev 10372)
+++ trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-07 12:39:54 UTC (rev 10373)
@@ -22,13 +22,16 @@
Y = dat.y{1};
U = dat.u{1};
Ts = dat.tsam{1};
+
+ n = max (na, nb);
## avoid warning: toeplitz: column wins anti-diagonal conflict
PhiY = toeplitz (Y(1:end-1, :), [Y(1, :); zeros(na-1, 1)]);
PhiU = toeplitz (U(1:end-1, :), [U(1, :); zeros(nb-1, 1)]);
Phi = [-PhiY, PhiU];
+ Phi = Phi(n:end, :)
- ## Theta = Phi \ Y(2:end, :); # naive formula
+ ## Theta = Phi \ Y(n+1:end, :); # naive formula
## solve linear least squares problem by pseudoinverse
## the pseudoinverse is computed by singular value decomposition
@@ -42,7 +45,7 @@
V = V(:, 1:r);
S = S(1:r);
U = U(:, 1:r);
- Theta = V * (S .\ (U' * Y(2:end, :))); # U' is the conjugate transpose
+ Theta = V * (S .\ (U' * Y(n+1:end, :))); # U' is the conjugate transpose
A = [1; Theta(1:na)]; # a0 = 1, a1 = Theta(1), an = Theta(n)
B = [0; Theta(na+1:end)]; # b0 = 0 (leading zero required by filt)
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-05-07 04:39:48
|
Revision: 10372
http://octave.svn.sourceforge.net/octave/?rev=10372&view=rev
Author: paramaniac
Date: 2012-05-07 04:39:42 +0000 (Mon, 07 May 2012)
Log Message:
-----------
control: handle special case num=0 or den=1 when displaying z^-1
Modified Paths:
--------------
trunk/octave-forge/main/control/inst/@tf/display.m
trunk/octave-forge/main/control/inst/filt.m
Modified: trunk/octave-forge/main/control/inst/@tf/display.m
===================================================================
--- trunk/octave-forge/main/control/inst/@tf/display.m 2012-05-07 03:27:54 UTC (rev 10371)
+++ trunk/octave-forge/main/control/inst/@tf/display.m 2012-05-07 04:39:42 UTC (rev 10372)
@@ -66,9 +66,11 @@
tfp = isa (num, "tfpoly");
- if (tfp && num == 0)
+ if (num == tfpoly (0))
str = [" ", name, ": 0"];
- elseif (tfp && den == 1)
+ elseif ((tfp && den == 1) || (! tfp && isequal (den, 1)))
+ ## elseif (den == tfpoly (1)) doesn't work because it
+ ## would mistakingly accept non-tfpoly denominators like [0, 1]
str = [" ", name, ": "];
numstr = tfpoly2str (num, tfvar);
str = [str, numstr];
Modified: trunk/octave-forge/main/control/inst/filt.m
===================================================================
--- trunk/octave-forge/main/control/inst/filt.m 2012-05-07 03:27:54 UTC (rev 10371)
+++ trunk/octave-forge/main/control/inst/filt.m 2012-05-07 04:39:42 UTC (rev 10372)
@@ -78,11 +78,13 @@
switch (nargin)
case 0 # filt ()
sys = tf ();
+ ## sys.inv = true;
return;
case 1 # filt (sys), filt (matrix)
if (isa (num, "lti") || is_real_matrix (num))
sys = tf (num);
+ ## sys.inv = true; # would be a problem for continuous-time LTI models
return;
else
print_usage ();
@@ -115,7 +117,9 @@
den = cellfun (@postpad, den, lmax, "uniformoutput", false);
## use standard tf constructor
- ## sys is stored and displayed in standard z form, not z^-1
+ ## sys is stored in standard z form, not z^-1
+ ## so we can mix it with regular transfer function models
+ ## property "inv", true displays sys in z^-1 form
sys = tf (num, den, tsam, "inv", true, varargin{:});
endswitch
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-05-07 03:28:01
|
Revision: 10371
http://octave.svn.sourceforge.net/octave/?rev=10371&view=rev
Author: paramaniac
Date: 2012-05-07 03:27:54 +0000 (Mon, 07 May 2012)
Log Message:
-----------
control: handle LTI property 'inv' in get & set routines
Modified Paths:
--------------
trunk/octave-forge/main/control/inst/@tf/__get__.m
trunk/octave-forge/main/control/inst/@tf/__set__.m
Modified: trunk/octave-forge/main/control/inst/@tf/__get__.m
===================================================================
--- trunk/octave-forge/main/control/inst/@tf/__get__.m 2012-05-06 20:20:20 UTC (rev 10370)
+++ trunk/octave-forge/main/control/inst/@tf/__get__.m 2012-05-07 03:27:54 UTC (rev 10371)
@@ -1,4 +1,4 @@
-## Copyright (C) 2009, 2010 Lukas F. Reichlin
+## Copyright (C) 2009, 2010, 2012 Lukas F. Reichlin
##
## This file is part of LTI Syncope.
##
@@ -20,7 +20,7 @@
## Author: Lukas Reichlin <luk...@gm...>
## Created: October 2009
-## Version: 0.2
+## Version: 0.3
function val = __get__ (sys, prop)
@@ -34,6 +34,9 @@
case {"tfvar", "variable"}
val = sys.tfvar;
+ case "inv"
+ val = sys.inv;
+
otherwise
error ("tf: get: invalid property name");
endswitch
Modified: trunk/octave-forge/main/control/inst/@tf/__set__.m
===================================================================
--- trunk/octave-forge/main/control/inst/@tf/__set__.m 2012-05-06 20:20:20 UTC (rev 10370)
+++ trunk/octave-forge/main/control/inst/@tf/__set__.m 2012-05-07 03:27:54 UTC (rev 10371)
@@ -43,10 +43,10 @@
endif
case "inv"
- if (islogical (val))
+ if (isscalar (val))
sys.inv = logical (val);
else
- error ("tf: set: property 'inv' must be a logical");
+ error ("tf: set: property 'inv' must be a scalar logical");
endif
otherwise
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-05-06 20:20:26
|
Revision: 10370
http://octave.svn.sourceforge.net/octave/?rev=10370&view=rev
Author: paramaniac
Date: 2012-05-06 20:20:20 +0000 (Sun, 06 May 2012)
Log Message:
-----------
control: remove completed task from projects file (transfer function variable z^-1)
Modified Paths:
--------------
trunk/octave-forge/main/control/devel/PROJECTS
Modified: trunk/octave-forge/main/control/devel/PROJECTS
===================================================================
--- trunk/octave-forge/main/control/devel/PROJECTS 2012-05-06 20:17:43 UTC (rev 10369)
+++ trunk/octave-forge/main/control/devel/PROJECTS 2012-05-06 20:20:20 UTC (rev 10370)
@@ -91,8 +91,6 @@
Transfer Function Models:
------------------------
- * Support "z^-1" as tfvar (@tf and @tfpoly)
-
* Find a transfer function to state-space conversion routine which
supports individual denominators for each input-output channel
(@tf/__sys2ss__.m). Support for non-proper descriptor state-space
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-05-06 20:17:49
|
Revision: 10369
http://octave.svn.sourceforge.net/octave/?rev=10369&view=rev
Author: paramaniac
Date: 2012-05-06 20:17:43 +0000 (Sun, 06 May 2012)
Log Message:
-----------
control: fix display of transfer function variable z^-1
Modified Paths:
--------------
trunk/octave-forge/main/control/inst/filt.m
trunk/octave-forge/main/control/inst/tfpoly2str.m
Modified: trunk/octave-forge/main/control/inst/filt.m
===================================================================
--- trunk/octave-forge/main/control/inst/filt.m 2012-05-06 19:36:54 UTC (rev 10368)
+++ trunk/octave-forge/main/control/inst/filt.m 2012-05-06 20:17:43 UTC (rev 10369)
@@ -57,9 +57,9 @@
##
## Transfer function 'H' from input 'u1' to output ...
##
-## 3 z
-## y1: -------------
-## z^2 + 4 z + 2
+## 3 z^-1
+## y1: -------------------
+## 1 + 4 z^-1 + 2 z^-2
##
## Sampling time: unspecified
## Discrete-time model.
Modified: trunk/octave-forge/main/control/inst/tfpoly2str.m
===================================================================
--- trunk/octave-forge/main/control/inst/tfpoly2str.m 2012-05-06 19:36:54 UTC (rev 10368)
+++ trunk/octave-forge/main/control/inst/tfpoly2str.m 2012-05-06 20:17:43 UTC (rev 10369)
@@ -27,13 +27,21 @@
function str = tfpoly2str (p, tfvar = "x")
+ ## TODO: simplify this ugly code
+
str = "";
lp = numel (p);
- if (lp > 0)
- ## first element (lowest order 0)
- a = p(1);
+ if (lp > 0) # first element (lowest order)
+ idx = find (p); # first non-zero element
+ if (isempty (idx))
+ str = "0";
+ return;
+ else
+ idx = idx(1);
+ endif
+ a = p(idx);
if (a < 0)
cs = "-";
@@ -41,11 +49,18 @@
cs = "";
endif
- str = [cs, num2str(abs (a), 4)];
+ if (idx == 1)
+ str = [cs, num2str(abs (a), 4)];
+ else
+ if (abs (a) == 1)
+ str = [cs, __variable__(tfvar, idx-1)];
+ else
+ str = [cs, __coefficient__(a), " ", __variable__(tfvar, idx-1)];
+ endif
+ endif
- if (lp > 1)
- ## remaining elements of higher order
- for k = 2 : lp
+ if (lp > idx) # remaining elements of higher order
+ for k = idx+1 : lp
a = p(k);
if (a != 0)
@@ -84,6 +99,6 @@
function str = __variable__ (tfvar, n)
- str = [tfvar, "^-", num2str(n)];
+ str = [tfvar, "^-", num2str(n)];
endfunction
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-05-06 19:37:01
|
Revision: 10368
http://octave.svn.sourceforge.net/octave/?rev=10368&view=rev
Author: paramaniac
Date: 2012-05-06 19:36:54 +0000 (Sun, 06 May 2012)
Log Message:
-----------
control: support transfer function variable z^-1
Modified Paths:
--------------
trunk/octave-forge/main/control/inst/@tf/__set__.m
trunk/octave-forge/main/control/inst/@tf/__sys_group__.m
trunk/octave-forge/main/control/inst/@tf/display.m
trunk/octave-forge/main/control/inst/@tf/tf.m
trunk/octave-forge/main/control/inst/filt.m
Added Paths:
-----------
trunk/octave-forge/main/control/inst/tfpoly2str.m
Modified: trunk/octave-forge/main/control/inst/@tf/__set__.m
===================================================================
--- trunk/octave-forge/main/control/inst/@tf/__set__.m 2012-05-06 15:53:35 UTC (rev 10367)
+++ trunk/octave-forge/main/control/inst/@tf/__set__.m 2012-05-06 19:36:54 UTC (rev 10368)
@@ -1,4 +1,4 @@
-## Copyright (C) 2009 Lukas F. Reichlin
+## Copyright (C) 2009, 2012 Lukas F. Reichlin
##
## This file is part of LTI Syncope.
##
@@ -20,7 +20,7 @@
## Author: Lukas Reichlin <luk...@gm...>
## Created: October 2009
-## Version: 0.1
+## Version: 0.2
function sys = __set__ (sys, prop, val)
@@ -42,6 +42,13 @@
error ("tf: set: invalid transfer function variable");
endif
+ case "inv"
+ if (islogical (val))
+ sys.inv = logical (val);
+ else
+ error ("tf: set: property 'inv' must be a logical");
+ endif
+
otherwise
error ("tf: set: invalid property name");
Modified: trunk/octave-forge/main/control/inst/@tf/__sys_group__.m
===================================================================
--- trunk/octave-forge/main/control/inst/@tf/__sys_group__.m 2012-05-06 15:53:35 UTC (rev 10367)
+++ trunk/octave-forge/main/control/inst/@tf/__sys_group__.m 2012-05-06 19:36:54 UTC (rev 10368)
@@ -1,4 +1,4 @@
-## Copyright (C) 2009 Lukas F. Reichlin
+## Copyright (C) 2009, 2012 Lukas F. Reichlin
##
## This file is part of LTI Syncope.
##
@@ -22,7 +22,7 @@
## Author: Lukas Reichlin <luk...@gm...>
## Created: September 2009
-## Version: 0.1
+## Version: 0.2
function retsys = __sys_group__ (sys1, sys2)
@@ -61,4 +61,8 @@
retsys.tfvar = sys1.tfvar;
endif
-endfunction
\ No newline at end of file
+ if (sys1.inv && sys2.inv)
+ retsys.inv = true;
+ endif
+
+endfunction
Modified: trunk/octave-forge/main/control/inst/@tf/display.m
===================================================================
--- trunk/octave-forge/main/control/inst/@tf/display.m 2012-05-06 15:53:35 UTC (rev 10367)
+++ trunk/octave-forge/main/control/inst/@tf/display.m 2012-05-06 19:36:54 UTC (rev 10368)
@@ -1,4 +1,4 @@
-## Copyright (C) 2009, 2010, 2011 Lukas F. Reichlin
+## Copyright (C) 2009, 2010, 2011, 2012 Lukas F. Reichlin
##
## This file is part of LTI Syncope.
##
@@ -20,7 +20,7 @@
## Author: Lukas Reichlin <luk...@gm...>
## Created: September 2009
-## Version: 0.3
+## Version: 0.4
function display (sys)
@@ -31,12 +31,19 @@
[outname, p] = __labels__ (outname, "y");
disp ("");
+
+ if (sys.inv)
+ [num, den] = filtdata (sys);
+ else
+ num = sys.num;
+ den = sys.den;
+ endif
for nu = 1 : m
disp (["Transfer function '", sysname, "' from input '", inname{nu}, "' to output ..."]);
disp ("");
for ny = 1 : p
- __disp_frac__ (sys.num{ny, nu}, sys.den{ny, nu}, sys.tfvar, outname{ny});
+ __disp_frac__ (num{ny, nu}, den{ny, nu}, sys.tfvar, outname{ny});
endfor
endfor
@@ -56,10 +63,12 @@
function __disp_frac__ (num, den, tfvar, name)
MAX_LEN = 12; # max length of output name
+
+ tfp = isa (num, "tfpoly");
- if (num == 0)
+ if (tfp && num == 0)
str = [" ", name, ": 0"];
- elseif (den == 1)
+ elseif (tfp && den == 1)
str = [" ", name, ": "];
numstr = tfpoly2str (num, tfvar);
str = [str, numstr];
Modified: trunk/octave-forge/main/control/inst/@tf/tf.m
===================================================================
--- trunk/octave-forge/main/control/inst/@tf/tf.m 2012-05-06 15:53:35 UTC (rev 10367)
+++ trunk/octave-forge/main/control/inst/@tf/tf.m 2012-05-06 19:36:54 UTC (rev 10368)
@@ -1,4 +1,4 @@
-## Copyright (C) 2009, 2010, 2011 Lukas F. Reichlin
+## Copyright (C) 2009, 2010, 2011, 2012 Lukas F. Reichlin
##
## This file is part of LTI Syncope.
##
@@ -119,7 +119,7 @@
## Author: Lukas Reichlin <luk...@gm...>
## Created: September 2009
-## Version: 0.2.1
+## Version: 0.3
function sys = tf (num = {}, den = {}, varargin)
@@ -196,7 +196,8 @@
tfdata = struct ("num", {num},
"den", {den},
- "tfvar", tfvar); # struct for tf-specific data
+ "tfvar", tfvar,
+ "inv", false); # struct for tf-specific data
ltisys = lti (p, m, tsam); # parent class for general lti data
Modified: trunk/octave-forge/main/control/inst/filt.m
===================================================================
--- trunk/octave-forge/main/control/inst/filt.m 2012-05-06 15:53:35 UTC (rev 10367)
+++ trunk/octave-forge/main/control/inst/filt.m 2012-05-06 19:36:54 UTC (rev 10368)
@@ -116,7 +116,7 @@
## use standard tf constructor
## sys is stored and displayed in standard z form, not z^-1
- sys = tf (num, den, tsam, varargin{:});
+ sys = tf (num, den, tsam, "inv", true, varargin{:});
endswitch
endfunction
Added: trunk/octave-forge/main/control/inst/tfpoly2str.m
===================================================================
--- trunk/octave-forge/main/control/inst/tfpoly2str.m (rev 0)
+++ trunk/octave-forge/main/control/inst/tfpoly2str.m 2012-05-06 19:36:54 UTC (rev 10368)
@@ -0,0 +1,89 @@
+## Copyright (C) 2012 Lukas F. Reichlin
+##
+## This file is part of LTI Syncope.
+##
+## LTI Syncope is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## LTI Syncope is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{str} =} tfpoly2str (@var{p})
+## @deftypefnx {Function File} {@var{str} =} tfpoly2str (@var{p}, @var{tfvar})
+## Return the string of a polynomial with string @var{tfvar} as variable.
+## @end deftypefn
+
+## Author: Lukas Reichlin <luk...@gm...>
+## Created: May 2012
+## Version: 0.1
+
+function str = tfpoly2str (p, tfvar = "x")
+
+ str = "";
+
+ lp = numel (p);
+
+ if (lp > 0)
+ ## first element (lowest order 0)
+ a = p(1);
+
+ if (a < 0)
+ cs = "-";
+ else
+ cs = "";
+ endif
+
+ str = [cs, num2str(abs (a), 4)];
+
+ if (lp > 1)
+ ## remaining elements of higher order
+ for k = 2 : lp
+ a = p(k);
+
+ if (a != 0)
+ if (a < 0)
+ cs = " - ";
+ else
+ cs = " + ";
+ endif
+
+ if (abs (a) == 1)
+ str = [str, cs, __variable__(tfvar, k-1)];
+ else
+ str = [str, cs, __coefficient__(a), " ", __variable__(tfvar, k-1)];
+ endif
+ endif
+ endfor
+
+ endif
+ endif
+
+endfunction
+
+
+function str = __coefficient__ (a)
+
+ b = abs (a);
+
+ if (b == 1)
+ str = "";
+ else
+ str = num2str (b, 4);
+ endif
+
+endfunction
+
+
+function str = __variable__ (tfvar, n)
+
+ str = [tfvar, "^-", num2str(n)];
+
+endfunction
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-05-06 15:53:41
|
Revision: 10367
http://octave.svn.sourceforge.net/octave/?rev=10367&view=rev
Author: paramaniac
Date: 2012-05-06 15:53:35 +0000 (Sun, 06 May 2012)
Log Message:
-----------
control-devel: use filt instead of tf
Modified Paths:
--------------
trunk/octave-forge/extra/control-devel/inst/arx.m
Modified: trunk/octave-forge/extra/control-devel/inst/arx.m
===================================================================
--- trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-06 14:32:55 UTC (rev 10366)
+++ trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-06 15:53:35 UTC (rev 10367)
@@ -28,13 +28,13 @@
PhiU = toeplitz (U(1:end-1, :), [U(1, :); zeros(nb-1, 1)]);
Phi = [-PhiY, PhiU];
- Theta = Phi \ Y(2:end, :) # naive formula
+ ## Theta = Phi \ Y(2:end, :); # naive formula
## solve linear least squares problem by pseudoinverse
## the pseudoinverse is computed by singular value decomposition
## M = U S V* ---> M+ = V S+ U*
## Th = Ph \ Y = Ph+ Y
- ## Th = V S+ U' Y, S+ = 1 ./ diag (S)
+ ## Th = V S+ U* Y, S+ = 1 ./ diag (S)
[U, S, V] = svd (Phi, 0); # 0 for "economy size" decomposition, U overwrites input U
S = diag (S); # extract main diagonal
@@ -42,11 +42,11 @@
V = V(:, 1:r);
S = S(1:r);
U = U(:, 1:r);
- Theta = V * (S .\ (U' * Y(2:end, :))) # U' is the conjugate transpose
+ Theta = V * (S .\ (U' * Y(2:end, :))); # U' is the conjugate transpose
- A = [1; Theta(1:na)]; # a0 = 1, a1 = Theta(1), an = Theta(n)
- B = Theta(na+1:end); # b0 = 0 (leading zeros are removed by tf, no need to add one)
+ A = [1; Theta(1:na)]; # a0 = 1, a1 = Theta(1), an = Theta(n)
+ B = [0; Theta(na+1:end)]; # b0 = 0 (leading zero required by filt)
- sys = tf ({B, 1}, {A, A}, Ts);
+ sys = filt ({B, 1}, {A, A}, Ts);
endfunction
\ No newline at end of file
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <as...@us...> - 2012-05-06 14:33:01
|
Revision: 10366
http://octave.svn.sourceforge.net/octave/?rev=10366&view=rev
Author: asnelt
Date: 2012-05-06 14:32:55 +0000 (Sun, 06 May 2012)
Log Message:
-----------
Relaxed condition for invalid rows in mnpdf and mnrnd so that larger numbers of
categories are supported.
Modified Paths:
--------------
trunk/octave-forge/main/statistics/NEWS
trunk/octave-forge/main/statistics/inst/mnpdf.m
trunk/octave-forge/main/statistics/inst/mnrnd.m
Modified: trunk/octave-forge/main/statistics/NEWS
===================================================================
--- trunk/octave-forge/main/statistics/NEWS 2012-05-05 13:03:58 UTC (rev 10365)
+++ trunk/octave-forge/main/statistics/NEWS 2012-05-06 14:32:55 UTC (rev 10366)
@@ -1,3 +1,9 @@
+Summary of important user-visible changes for statistics 1.1.3:
+-------------------------------------------------------------------
+
+ ** The functions mnpdf and mnrnd are now also usable for greater numbers
+ of categories for which the rows do not exactly sum to 1.
+
Summary of important user-visible changes for statistics 1.1.2:
-------------------------------------------------------------------
Modified: trunk/octave-forge/main/statistics/inst/mnpdf.m
===================================================================
--- trunk/octave-forge/main/statistics/inst/mnpdf.m 2012-05-05 13:03:58 UTC (rev 10365)
+++ trunk/octave-forge/main/statistics/inst/mnpdf.m 2012-05-06 14:32:55 UTC (rev 10366)
@@ -115,7 +115,9 @@
t = x .* log (p);
t(x == 0) = 0;
y = exp (gammaln (n+1) - sum (gammaln (x+1), 2) + sum (t, 2));
- y(sum (p, 2) != 1) = NaN;
+ # Set invalid rows to NaN
+ k = (abs (sum (p, 2) - 1) > 1e-6);
+ y(k) = NaN;
endfunction
Modified: trunk/octave-forge/main/statistics/inst/mnrnd.m
===================================================================
--- trunk/octave-forge/main/statistics/inst/mnrnd.m 2012-05-05 13:03:58 UTC (rev 10365)
+++ trunk/octave-forge/main/statistics/inst/mnrnd.m 2012-05-06 14:32:55 UTC (rev 10366)
@@ -135,23 +135,24 @@
# Upper bounds of categories
ub = cumsum (p, 2);
+ # Make sure that the greatest upper bound is 1
+ gub = ub(:, end);
+ ub(:, end) = 1;
# Lower bounds of categories
lb = [zeros(sz(1), 1) ub(:, 1:(end-1))];
# Draw multinomial samples
x = zeros (sz);
for i = 1:sz(1)
- # Check if the probabilities sum to 1
- if (ub(i, end) == 1)
- # Draw uniform random numbers
- r = repmat (rand (n(i), 1), 1, sz(2));
- # Compare the random numbers of r to the cumulated probabilities of p and
- # count the number of samples for each category
- x(i, :) = sum (r <= repmat (ub(i, :), n(i), 1) & r > repmat (lb(i, :), n(i), 1), 1);
- else
- x(i, :) = NaN;
- endif
+ # Draw uniform random numbers
+ r = repmat (rand (n(i), 1), 1, sz(2));
+ # Compare the random numbers of r to the cumulated probabilities of p and
+ # count the number of samples for each category
+ x(i, :) = sum (r <= repmat (ub(i, :), n(i), 1) & r > repmat (lb(i, :), n(i), 1), 1);
endfor
+ # Set invalid rows to NaN
+ k = (abs (gub - 1) > 1e-6);
+ x(k, :) = NaN;
endfunction
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-05-05 13:04:06
|
Revision: 10365
http://octave.svn.sourceforge.net/octave/?rev=10365&view=rev
Author: paramaniac
Date: 2012-05-05 13:03:58 +0000 (Sat, 05 May 2012)
Log Message:
-----------
control-devel: add some debug code
Modified Paths:
--------------
trunk/octave-forge/extra/control-devel/devel/GlassFurnace.m
trunk/octave-forge/extra/control-devel/devel/PowerPlant.m
trunk/octave-forge/extra/control-devel/src/devel_slicot_functions.cc
Added Paths:
-----------
trunk/octave-forge/extra/control-devel/devel/CDplayer_a.m
trunk/octave-forge/extra/control-devel/devel/Destillation_a.m
trunk/octave-forge/extra/control-devel/devel/GlassFurnace_a.m
trunk/octave-forge/extra/control-devel/devel/PowerPlant_a.m
trunk/octave-forge/extra/control-devel/devel/ident_a.m
trunk/octave-forge/extra/control-devel/src/slident_a.cc
trunk/octave-forge/extra/control-devel/src/slident_b.cc
trunk/octave-forge/extra/control-devel/src/slident_c.cc
Added: trunk/octave-forge/extra/control-devel/devel/CDplayer_a.m
===================================================================
--- trunk/octave-forge/extra/control-devel/devel/CDplayer_a.m (rev 0)
+++ trunk/octave-forge/extra/control-devel/devel/CDplayer_a.m 2012-05-05 13:03:58 UTC (rev 10365)
@@ -0,0 +1,67 @@
+%{
+Contributed by:
+ Favoreel
+ KULeuven
+ Departement Electrotechniek ESAT/SISTA
+Kardinaal Mercierlaan 94
+B-3001 Leuven
+Belgium
+ wou...@es...
+Description:
+ Data from the mechanical construction of a CD player arm.
+ The inputs are the forces of the mechanical actuators
+ while the outputs are related to the tracking accuracy of the arm.
+ The data was measured in closed loop, and then through a two-step
+ procedure converted to open loop equivalent data
+ The inputs are highly colored.
+Sampling:
+Number:
+ 2048
+Inputs:
+ u: forces of the mechanical actuators
+Outputs:
+ y: tracking accuracy of the arm
+References:
+ We are grateful to R. de Callafon of the
+ Mechanical Engineering Systems and Control group of Delft, who
+ provided us with these data.
+
+ - Van Den Hof P., Schrama R.J.P., An Indirect Method for Transfer
+ Function Estimation From Closed Loop Data. Automatica, Vol. 29,
+ no. 6, pp. 1523-1527, 1993.
+
+Properties:
+Columns:
+ Column 1: input u1
+ Column 2: input u2
+ Column 1: output y1
+ Column 2: output y2
+Category:
+ mechanical systems
+
+%}
+
+close all, clc
+
+load CD_player_arm-1.dat
+U=CD_player_arm_1(:,1:2);
+Y=CD_player_arm_1(:,3:4);
+
+
+dat = iddata (Y, U)
+
+[sys, x0] = ident_a (dat, 15, 8) % s=15, n=8
+
+
+[y, t] = lsim (sys, U, [], x0);
+
+err = norm (Y - y, 1) / norm (Y, 1)
+
+figure (1)
+p = columns (Y);
+for k = 1 : p
+ subplot (2, 1, k)
+ plot (t, Y(:,k), 'b', t, y(:,k), 'r')
+endfor
+
+
Added: trunk/octave-forge/extra/control-devel/devel/Destillation_a.m
===================================================================
--- trunk/octave-forge/extra/control-devel/devel/Destillation_a.m (rev 0)
+++ trunk/octave-forge/extra/control-devel/devel/Destillation_a.m 2012-05-05 13:03:58 UTC (rev 10365)
@@ -0,0 +1,83 @@
+%{
+This file describes the data in the destill.dat file.
+1. Contributed by:
+ Peter Van Overschee
+ K.U.Leuven - ESAT - SISTA
+ K. Mercierlaan 94
+ 3001 Heverlee
+ Pet...@es...
+2. Process/Description:
+ Data of a simulation (not real !) related to the identification
+ of an ethane-ethylene destillationcolumn. The series consists of 4
+ series:
+ U_dest, Y_dest: without noise (original series)
+ U_dest_n10, Y_dest_n10: 10 percent additive white noise
+ U_dest_n20, Y_dest_n20: 20 percent additive white noise
+ U_dest_n30, Y_dest_n30: 30 percent additive white noise
+3. Sampling time
+ 15 min.
+4. Number of samples:
+ 90 samples
+5. Inputs:
+ a. ratio between the reboiler duty and the feed flow
+ b. ratio between the reflux rate and the feed flow
+ c. ratio between the distillate and the feed flow
+ d. input ethane composition
+ e. top pressure
+6. Outputs:
+ a. top ethane composition
+ b. bottom ethylene composition
+ c. top-bottom differential pressure.
+7. References:
+ R.P. Guidorzi, M.P. Losito, T. Muratori, The range error test in the
+ structural identification of linear multivariable systems,
+ IEEE transactions on automatic control, Vol AC-27, pp 1044-1054, oct.
+ 1982.
+8. Known properties/peculiarities
+
+9. Some MATLAB-code to retrieve the data
+ !gunzip destill.dat.Z
+ load destill.dat
+ U=destill(:,1:20);
+ Y=destill(:,21:32);
+ U_dest=U(:,1:5);
+ U_dest_n10=U(:,6:10);
+ U_dest_n20=U(:,11:15);
+ U_dest_n30=U(:,16:20);
+ Y_dest=Y(:,1:3);
+ Y_dest_n10=Y(:,4:6);
+ Y_dest_n20=Y(:,7:9);
+ Y_dest_n30=Y(:,10:12);
+%}
+
+ close all, clc
+
+load destill.dat
+U=destill(:,1:20);
+Y=destill(:,21:32);
+U_dest=U(:,1:5);
+U_dest_n10=U(:,6:10);
+U_dest_n20=U(:,11:15);
+U_dest_n30=U(:,16:20);
+Y_dest=Y(:,1:3);
+Y_dest_n10=Y(:,4:6);
+Y_dest_n20=Y(:,7:9);
+Y_dest_n30=Y(:,10:12);
+
+
+dat = iddata (Y_dest, U_dest)
+
+[sys, x0] = ident_a (dat, 5, 4) % s=5, n=4
+
+
+[y, t] = lsim (sys, U_dest, [], x0);
+%[y, t] = lsim (sys, U_dest);
+
+err = norm (Y_dest - y, 1) / norm (Y_dest, 1)
+
+figure (1)
+%plot (t, Y_dest, 'b')
+plot (t, Y_dest, 'b', t, y, 'r')
+legend ('y measured', 'y simulated', 'location', 'southeast')
+
+
Modified: trunk/octave-forge/extra/control-devel/devel/GlassFurnace.m
===================================================================
--- trunk/octave-forge/extra/control-devel/devel/GlassFurnace.m 2012-05-05 06:48:41 UTC (rev 10364)
+++ trunk/octave-forge/extra/control-devel/devel/GlassFurnace.m 2012-05-05 13:03:58 UTC (rev 10365)
@@ -37,8 +37,9 @@
%}
-clear all, close all, clc
+close all, clc
+
load glassfurnace.dat
T=glassfurnace(:,1);
U=glassfurnace(:,2:4);
Added: trunk/octave-forge/extra/control-devel/devel/GlassFurnace_a.m
===================================================================
--- trunk/octave-forge/extra/control-devel/devel/GlassFurnace_a.m (rev 0)
+++ trunk/octave-forge/extra/control-devel/devel/GlassFurnace_a.m 2012-05-05 13:03:58 UTC (rev 10365)
@@ -0,0 +1,66 @@
+%{
+This file describes the data in the glassfurnace.dat file.
+1. Contributed by:
+ Peter Van Overschee
+ K.U.Leuven - ESAT - SISTA
+ K. Mercierlaan 94
+ 3001 Heverlee
+ Pet...@es...
+2. Process/Description:
+ Data of a glassfurnace (Philips)
+3. Sampling time
+
+4. Number of samples:
+ 1247 samples
+5. Inputs:
+ a. heating input
+ b. cooling input
+ c. heating input
+6. Outputs:
+ a. 6 outputs from temperature sensors in a cross section of the
+ furnace
+7. References:
+ a. Van Overschee P., De Moor B., N4SID : Subspace Algorithms for
+ the Identification of Combined Deterministic-Stochastic Systems,
+ Automatica, Special Issue on Statistical Signal Processing and Control,
+ Vol. 30, No. 1, 1994, pp. 75-93
+ b. Van Overschee P., "Subspace identification : Theory,
+ Implementation, Application" , Ph.D. Thesis, K.U.Leuven, February 1995.
+8. Known properties/peculiarities
+
+9. Some MATLAB-code to retrieve the data
+ !gunzip glassfurnace.dat.Z
+ load glassfurnace.dat
+ T=glassfurnace(:,1);
+ U=glassfurnace(:,2:4);
+ Y=glassfurnace(:,5:10);
+
+%}
+
+ close all, clc
+
+load glassfurnace.dat
+T=glassfurnace(:,1);
+U=glassfurnace(:,2:4);
+Y=glassfurnace(:,5:10);
+
+
+dat = iddata (Y, U)
+
+[sys, x0] = ident_a (dat, 10, 5) % s=10, n=5
+
+
+[y, t] = lsim (sys, U, [], x0);
+
+err = norm (Y - y, 1) / norm (Y, 1)
+
+figure (1)
+p = columns (Y);
+for k = 1 : p
+ subplot (3, 2, k)
+ plot (t, Y(:,k), 'b', t, y(:,k), 'r')
+endfor
+%title ('DaISy: Glass Furnace')
+%legend ('y measured', 'y simulated', 'location', 'southeast')
+
+
Modified: trunk/octave-forge/extra/control-devel/devel/PowerPlant.m
===================================================================
--- trunk/octave-forge/extra/control-devel/devel/PowerPlant.m 2012-05-05 06:48:41 UTC (rev 10364)
+++ trunk/octave-forge/extra/control-devel/devel/PowerPlant.m 2012-05-05 13:03:58 UTC (rev 10365)
@@ -40,7 +40,7 @@
%}
-clear all, close all, clc
+close all, clc
load powerplant.dat
U=powerplant(:,1:5);
Added: trunk/octave-forge/extra/control-devel/devel/PowerPlant_a.m
===================================================================
--- trunk/octave-forge/extra/control-devel/devel/PowerPlant_a.m (rev 0)
+++ trunk/octave-forge/extra/control-devel/devel/PowerPlant_a.m 2012-05-05 13:03:58 UTC (rev 10365)
@@ -0,0 +1,81 @@
+%{
+This file describes the data in the powerplant.dat file.
+1. Contributed by:
+ Peter Van Overschee
+ K.U.Leuven - ESAT - SISTA
+ K. Mercierlaan 94
+ 3001 Heverlee
+ Pet...@es...
+2. Process/Description:
+ data of a power plant (Pont-sur-Sambre (France)) of 120 MW
+3. Sampling time
+ 1228.8 sec
+4. Number of samples:
+ 200 samples
+5. Inputs:
+ 1. gas flow
+ 2. turbine valves opening
+ 3. super heater spray flow
+ 4. gas dampers
+ 5. air flow
+6. Outputs:
+ 1. steam pressure
+ 2. main stem temperature
+ 3. reheat steam temperature
+7. References:
+ a. R.P. Guidorzi, P. Rossi, Identification of a power plant from normal
+ operating records. Automatic control theory and applications (Canada,
+ Vol 2, pp 63-67, sept 1974.
+ b. Moonen M., De Moor B., Vandenberghe L., Vandewalle J., On- and
+ off-line identification of linear state-space models, International
+ Journal of Control, Vol. 49, Jan. 1989, pp.219-232
+8. Known properties/peculiarities
+
+9. Some MATLAB-code to retrieve the data
+ !gunzip powerplant.dat.Z
+ load powerplant.dat
+ U=powerplant(:,1:5);
+ Y=powerplant(:,6:8);
+ Yr=powerplant(:,9:11);
+
+%}
+
+ close all, clc
+
+load powerplant.dat
+U=powerplant(:,1:5);
+Y=powerplant(:,6:8);
+Yr=powerplant(:,9:11);
+
+inname = {'gas flow',
+ 'turbine valves opening',
+ 'super heater spray flow',
+ 'gas dampers',
+ 'air flow'};
+
+outname = {'steam pressure',
+ 'main steam temperature',
+ 'reheat steam temperature'};
+
+tsam = 1228.8;
+
+dat = iddata (Y, U, tsam, 'outname', outname, 'inname', inname)
+
+[sys, x0] = ident_a (dat, 10, 8) % s=10, n=8
+
+
+[y, t] = lsim (sys, U, [], x0);
+
+err = norm (Y - y, 1) / norm (Y, 1)
+
+figure (1)
+p = columns (Y);
+for k = 1 : p
+ subplot (3, 1, k)
+ plot (t, Y(:,k), 'b', t, y(:,k), 'r')
+endfor
+%title ('DaISy: Power Plant')
+%legend ('y measured', 'y simulated', 'location', 'southeast')
+
+st = isstable (sys)
+
Added: trunk/octave-forge/extra/control-devel/devel/ident_a.m
===================================================================
--- trunk/octave-forge/extra/control-devel/devel/ident_a.m (rev 0)
+++ trunk/octave-forge/extra/control-devel/devel/ident_a.m 2012-05-05 13:03:58 UTC (rev 10365)
@@ -0,0 +1,61 @@
+function [sys, x0] = ident_a (dat, s = [], nn = [])
+
+n=nn;
+
+ %nobr = 15;
+ meth = 2; % 2 % geht: meth/alg 1/1,
+ alg = 0; % 0 % geht nicht: meth/alg 0/1
+ jobd = 1;
+ batch = 3;
+ conct = 1;
+ ctrl = 0; %1;
+ rcond = 0.0;
+ tol = -1.0; % 0;
+
+ [ns, l, m, e] = size (dat);
+
+ if (isempty (s) && isempty (n))
+ nsmp = ns(1);
+ nobr = fix ((nsmp+1)/(2*(m+l+1)));
+ ctrl = 0; # confirm system order estimate
+ n = 0;
+ % nsmp >= 2*(m+l+1)*nobr - 1
+ % nobr <= (nsmp+1)/(2*(m+l+1))
+ elseif (isempty (s))
+ s = min (2*n, n+10);
+ nsmp = ns(1);
+ nobr = fix ((nsmp+1)/(2*(m+l+1)));
+ nobr = min (nobr, s);
+ ctrl = 1; # no confirmation
+ elseif (isempty (n))
+ nobr = s;
+ ctrl = 0; # confirm system order estimate
+ n = 0;
+ else # s & n non-empty
+ nsmp = ns(1);
+ nobr = fix ((nsmp+1)/(2*(m+l+1)));
+ if (s > nobr)
+ error ("ident: s > nobr");
+ endif
+ nobr = s;
+ ctrl = 1;
+ ## TODO: specify n for IB01BD
+ endif
+
+ %nsmp = ns(1)
+ %nobr = fix ((nsmp+1)/(2*(m+l+1)))
+ % nsmp >= 2*(m+l+1)*nobr - 1
+ % nobr <= (nsmp+1)/(2*(m+l+1))
+%nobr = 10
+ [r, sv, n] = slident_a (dat.y{1}, dat.u{1}, nobr, n, meth, alg, jobd, batch, conct, ctrl, rcond, tol);
+n
+n = nn;
+ [a, b, c, d, q, ry, s, k] = slident_b (dat.y{1}, dat.u{1}, nobr, n, meth, alg, jobd, batch, conct, ctrl, rcond, tol, \
+ r, sv, n);
+
+ x0 = slident_c (dat.y{1}, dat.u{1}, nobr, n, meth, alg, jobd, batch, conct, ctrl, rcond, tol, \
+ a, b, c, d);
+
+ sys = ss (a, b, c, d, -1);
+
+endfunction
Modified: trunk/octave-forge/extra/control-devel/src/devel_slicot_functions.cc
===================================================================
--- trunk/octave-forge/extra/control-devel/src/devel_slicot_functions.cc 2012-05-05 06:48:41 UTC (rev 10364)
+++ trunk/octave-forge/extra/control-devel/src/devel_slicot_functions.cc 2012-05-05 13:03:58 UTC (rev 10365)
@@ -1,3 +1,6 @@
// #include "slib01ad.cc" // preprocess the input-output data
#include "slident.cc" // system identification
#include "slitest.cc" // debug system identification
+#include "slident_a.cc"
+#include "slident_b.cc"
+#include "slident_c.cc"
\ No newline at end of file
Added: trunk/octave-forge/extra/control-devel/src/slident_a.cc
===================================================================
--- trunk/octave-forge/extra/control-devel/src/slident_a.cc (rev 0)
+++ trunk/octave-forge/extra/control-devel/src/slident_a.cc 2012-05-05 13:03:58 UTC (rev 10365)
@@ -0,0 +1,422 @@
+/*
+
+Copyright (C) 2012 Lukas F. Reichlin
+
+This file is part of LTI Syncope.
+
+LTI Syncope is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+LTI Syncope is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>.
+
+SLICOT system identification
+Uses SLICOT IB01AD, IB01BD and IB01CD by courtesy of NICONET e.V.
+<http://www.slicot.org>
+
+Author: Lukas Reichlin <luk...@gm...>
+Created: March 2012
+Version: 0.1
+
+*/
+
+#include <octave/oct.h>
+#include <f77-fcn.h>
+#include "common.h"
+
+extern "C"
+{
+ int F77_FUNC (ib01ad, IB01AD)
+ (char& METH, char& ALG, char& JOBD,
+ char& BATCH, char& CONCT, char& CTRL,
+ int& NOBR, int& M, int& L,
+ int& NSMP,
+ double* U, int& LDU,
+ double* Y, int& LDY,
+ int& N,
+ double* R, int& LDR,
+ double* SV,
+ double& RCOND, double& TOL,
+ int* IWORK,
+ double* DWORK, int& LDWORK,
+ int& IWARN, int& INFO);
+
+ int F77_FUNC (ib01bd, IB01BD)
+ (char& METH, char& JOB, char& JOBCK,
+ int& NOBR, int& N, int& M, int& L,
+ int& NSMPL,
+ double* R, int& LDR,
+ double* A, int& LDA,
+ double* C, int& LDC,
+ double* B, int& LDB,
+ double* D, int& LDD,
+ double* Q, int& LDQ,
+ double* RY, int& LDRY,
+ double* S, int& LDS,
+ double* K, int& LDK,
+ double& TOL,
+ int* IWORK,
+ double* DWORK, int& LDWORK,
+ bool* BWORK,
+ int& IWARN, int& INFO);
+
+ int F77_FUNC (ib01cd, IB01CD)
+ (char& JOBX0, char& COMUSE, char& JOB,
+ int& N, int& M, int& L,
+ int& NSMP,
+ double* A, int& LDA,
+ double* B, int& LDB,
+ double* C, int& LDC,
+ double* D, int& LDD,
+ double* U, int& LDU,
+ double* Y, int& LDY,
+ double* X0,
+ double* V, int& LDV,
+ double& TOL,
+ int* IWORK,
+ double* DWORK, int& LDWORK,
+ int& IWARN, int& INFO);
+}
+
+// PKG_ADD: autoload ("slident_a", "devel_slicot_functions.oct");
+DEFUN_DLD (slident_a, args, nargout,
+ "-*- texinfo -*-\n\
+Slicot IB01AD Release 5.0\n\
+No argument checking.\n\
+For internal use only.")
+{
+ int nargin = args.length ();
+ octave_value_list retval;
+
+ if (nargin != 12)
+ {
+ print_usage ();
+ }
+ else
+ {
+////////////////////////////////////////////////////////////////////////////////////
+// SLICOT IB01AD - preprocess the input-output data //
+////////////////////////////////////////////////////////////////////////////////////
+
+ // arguments in
+ char meth;
+ char alg;
+ char jobd;
+ char batch;
+ char conct;
+ char ctrl;
+ char metha;
+ char jobda; // ??? unused
+
+ Matrix y = args(0).matrix_value ();
+ Matrix u = args(1).matrix_value ();
+ int nobr = args(2).int_value ();
+ int nuser = args(3).int_value ();
+
+ const int imeth = args(4).int_value ();
+ const int ialg = args(5).int_value ();
+ const int ijobd = args(6).int_value ();
+ const int ibatch = args(7).int_value ();
+ const int iconct = args(8).int_value ();
+ const int ictrl = args(9).int_value ();
+
+ double rcond = args(10).double_value ();
+ double tol = args(11).double_value ();
+ double tolb = args(10).double_value (); // tolb = rcond
+
+
+ switch (imeth)
+ {
+ case 0:
+ meth = 'M';
+ metha = 'M';
+ break;
+ case 1:
+ meth = 'N';
+ metha = 'N';
+ break;
+ case 2:
+ meth = 'C';
+ metha = 'N'; // no typo here
+ break;
+ default:
+ error ("slib01ad: argument 'meth' invalid");
+ }
+
+ switch (ialg)
+ {
+ case 0:
+ alg = 'C';
+ break;
+ case 1:
+ alg = 'F';
+ break;
+ case 2:
+ alg = 'Q';
+ break;
+ default:
+ error ("slib01ad: argument 'alg' invalid");
+ }
+
+ if (meth == 'C')
+ jobd = 'N';
+ else if (ijobd == 0)
+ jobd = 'M';
+ else
+ jobd = 'N';
+
+ switch (ibatch)
+ {
+ case 0:
+ batch = 'F';
+ break;
+ case 1:
+ batch = 'I';
+ break;
+ case 2:
+ batch = 'L';
+ break;
+ case 3:
+ batch = 'O';
+ break;
+ default:
+ error ("slib01ad: argument 'batch' invalid");
+ }
+
+ if (iconct == 0)
+ conct = 'C';
+ else
+ conct = 'N';
+
+ if (ictrl == 0)
+ ctrl = 'C';
+ else
+ ctrl = 'N';
+
+
+ int m = u.columns (); // m: number of inputs
+ int l = y.columns (); // l: number of outputs
+ int nsmp = y.rows (); // nsmp: number of samples
+ // y.rows == u.rows is checked by iddata class
+ // TODO: check minimal nsmp size
+
+ if (batch == 'O')
+ {
+ if (nsmp < 2*(m+l+1)*nobr - 1)
+ error ("slident: require NSMP >= 2*(M+L+1)*NOBR - 1");
+ }
+ else
+ {
+ if (nsmp < 2*nobr)
+ error ("slident: require NSMP >= 2*NOBR");
+ }
+
+ int ldu;
+
+ if (m == 0)
+ ldu = 1;
+ else // m > 0
+ ldu = nsmp;
+
+ int ldy = nsmp;
+
+ // arguments out
+ int n;
+ int ldr;
+
+ if (metha == 'M' && jobd == 'M')
+ ldr = max (2*(m+l)*nobr, 3*m*nobr);
+ else if (metha == 'N' || (metha == 'M' && jobd == 'N'))
+ ldr = 2*(m+l)*nobr;
+ else
+ error ("slib01ad: could not handle 'ldr' case");
+
+ Matrix r (ldr, 2*(m+l)*nobr);
+ ColumnVector sv (l*nobr);
+
+ // workspace
+ int liwork;
+
+ if (metha == 'N') // if METH = 'N'
+ liwork = (m+l)*nobr;
+ else if (alg == 'F') // if METH = 'M' and ALG = 'F'
+ liwork = m+l;
+ else // if METH = 'M' and ALG = 'C' or 'Q'
+ liwork = 0;
+
+ // TODO: Handle 'k' for DWORK
+
+ int ldwork;
+ int ns = nsmp - 2*nobr + 1;
+
+ if (alg == 'C')
+ {
+ if (batch == 'F' || batch == 'I')
+ {
+ if (conct == 'C')
+ ldwork = (4*nobr-2)*(m+l);
+ else // (conct == 'N')
+ ldwork = 1;
+ }
+ else if (metha == 'M') // && (batch == 'L' || batch == 'O')
+ {
+ if (conct == 'C' && batch == 'L')
+ ldwork = max ((4*nobr-2)*(m+l), 5*l*nobr);
+ else if (jobd == 'M')
+ ldwork = max ((2*m-1)*nobr, (m+l)*nobr, 5*l*nobr);
+ else // (jobd == 'N')
+ ldwork = 5*l*nobr;
+ }
+ else // meth == 'N' && (batch == 'L' || batch == 'O')
+ {
+ ldwork = 5*(m+l)*nobr + 1;
+ }
+ }
+ else if (alg == 'F')
+ {
+ if (batch != 'O' && conct == 'C')
+ ldwork = (m+l)*2*nobr*(m+l+3);
+ else if (batch == 'F' || batch == 'I') // && conct == 'N'
+ ldwork = (m+l)*2*nobr*(m+l+1);
+ else // (batch == 'L' || '0' && conct == 'N')
+ ldwork = (m+l)*4*nobr*(m+l+1)+(m+l)*2*nobr;
+ }
+ else // (alg == 'Q')
+ {
+ // int ns = nsmp - 2*nobr + 1;
+
+ if (ldr >= ns && batch == 'F')
+ {
+ ldwork = 4*(m+l)*nobr;
+ }
+ else if (ldr >= ns && batch == 'O')
+ {
+ if (metha == 'M')
+ ldwork = max (4*(m+l)*nobr, 5*l*nobr);
+ else // (meth == 'N')
+ ldwork = 5*(m+l)*nobr + 1;
+ }
+ else if (conct == 'C' && (batch == 'I' || batch == 'L'))
+ {
+ ldwork = 4*(nobr+1)*(m+l)*nobr;
+ }
+ else // if ALG = 'Q', (BATCH = 'F' or 'O', and LDR < NS), or (BATCH = 'I' or 'L' and CONCT = 'N')
+ {
+ ldwork = 6*(m+l)*nobr;
+ }
+ }
+
+/*
+IB01AD.f Lines 438-445
+C FURTHER COMMENTS
+C
+C For ALG = 'Q', BATCH = 'O' and LDR < NS, or BATCH <> 'O', the
+C calculations could be rather inefficient if only minimal workspace
+C (see argument LDWORK) is provided. It is advisable to provide as
+C much workspace as possible. Almost optimal efficiency can be
+C obtained for LDWORK = (NS+2)*(2*(M+L)*NOBR), assuming that the
+C cache size is large enough to accommodate R, U, Y, and DWORK.
+*/
+
+// warning ("==================== ldwork before: %d =====================", ldwork);
+// ldwork = (ns+2)*(2*(m+l)*nobr);
+ldwork = max (ldwork, (ns+2)*(2*(m+l)*nobr));
+// ldwork *= 3;
+// warning ("==================== ldwork after: %d =====================", ldwork);
+
+
+/*
+IB01AD.f Lines 291-195:
+c the workspace used for alg = 'q' is
+c ldrwrk*2*(m+l)*nobr + 4*(m+l)*nobr,
+c where ldrwrk = ldwork/(2*(m+l)*nobr) - 2; recommended
+c value ldrwrk = ns, assuming a large enough cache size.
+c for good performance, ldwork should be larger.
+
+somehow ldrwrk and ldwork must have been mixed up here
+
+*/
+
+
+ OCTAVE_LOCAL_BUFFER (int, iwork, liwork);
+ OCTAVE_LOCAL_BUFFER (double, dwork, ldwork);
+
+ // error indicators
+ int iwarn = 0;
+ int info = 0;
+
+
+ // SLICOT routine IB01AD
+ F77_XFCN (ib01ad, IB01AD,
+ (metha, alg, jobd,
+ batch, conct, ctrl,
+ nobr, m, l,
+ nsmp,
+ u.fortran_vec (), ldu,
+ y.fortran_vec (), ldy,
+ n,
+ r.fortran_vec (), ldr,
+ sv.fortran_vec (),
+ rcond, tol,
+ iwork,
+ dwork, ldwork,
+ iwarn, info));
+
+
+ if (f77_exception_encountered)
+ error ("ident: exception in SLICOT subroutine IB01AD");
+
+ static const char* err_msg[] = {
+ "0: OK",
+ "1: a fast algorithm was requested (ALG = 'C', or 'F') "
+ "in sequential data processing, but it failed; the "
+ "routine can be repeatedly called again using the "
+ "standard QR algorithm",
+ "2: the singular value decomposition (SVD) algorithm did "
+ "not converge"};
+
+ static const char* warn_msg[] = {
+ "0: OK",
+ "1: the number of 100 cycles in sequential data "
+ "processing has been exhausted without signaling "
+ "that the last block of data was get; the cycle "
+ "counter was reinitialized",
+ "2: a fast algorithm was requested (ALG = 'C' or 'F'), "
+ "but it failed, and the QR algorithm was then used "
+ "(non-sequential data processing)",
+ "3: all singular values were exactly zero, hence N = 0 "
+ "(both input and output were identically zero)",
+ "4: the least squares problems with coefficient matrix "
+ "U_f, used for computing the weighted oblique "
+ "projection (for METH = 'N'), have a rank-deficient "
+ "coefficient matrix",
+ "5: the least squares problem with coefficient matrix "
+ "r_1 [6], used for computing the weighted oblique "
+ "projection (for METH = 'N'), has a rank-deficient "
+ "coefficient matrix"};
+
+
+ error_msg ("ident", info, 2, err_msg);
+ warning_msg ("ident", iwarn, 5, warn_msg);
+
+
+ // resize
+ int rs = 2*(m+l)*nobr;
+ r.resize (rs, rs);
+
+
+ // return values
+ retval(0) = r;
+ retval(1) = sv;
+ retval(2) = octave_value (n);
+ }
+
+ return retval;
+}
Added: trunk/octave-forge/extra/control-devel/src/slident_b.cc
===================================================================
--- trunk/octave-forge/extra/control-devel/src/slident_b.cc (rev 0)
+++ trunk/octave-forge/extra/control-devel/src/slident_b.cc 2012-05-05 13:03:58 UTC (rev 10365)
@@ -0,0 +1,669 @@
+/*
+
+Copyright (C) 2012 Lukas F. Reichlin
+
+This file is part of LTI Syncope.
+
+LTI Syncope is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+LTI Syncope is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>.
+
+SLICOT system identification
+Uses SLICOT IB01AD, IB01BD and IB01CD by courtesy of NICONET e.V.
+<http://www.slicot.org>
+
+Author: Lukas Reichlin <luk...@gm...>
+Created: March 2012
+Version: 0.1
+
+*/
+
+#include <octave/oct.h>
+#include <f77-fcn.h>
+#include "common.h"
+
+extern "C"
+{
+ int F77_FUNC (ib01ad, IB01AD)
+ (char& METH, char& ALG, char& JOBD,
+ char& BATCH, char& CONCT, char& CTRL,
+ int& NOBR, int& M, int& L,
+ int& NSMP,
+ double* U, int& LDU,
+ double* Y, int& LDY,
+ int& N,
+ double* R, int& LDR,
+ double* SV,
+ double& RCOND, double& TOL,
+ int* IWORK,
+ double* DWORK, int& LDWORK,
+ int& IWARN, int& INFO);
+
+ int F77_FUNC (ib01bd, IB01BD)
+ (char& METH, char& JOB, char& JOBCK,
+ int& NOBR, int& N, int& M, int& L,
+ int& NSMPL,
+ double* R, int& LDR,
+ double* A, int& LDA,
+ double* C, int& LDC,
+ double* B, int& LDB,
+ double* D, int& LDD,
+ double* Q, int& LDQ,
+ double* RY, int& LDRY,
+ double* S, int& LDS,
+ double* K, int& LDK,
+ double& TOL,
+ int* IWORK,
+ double* DWORK, int& LDWORK,
+ bool* BWORK,
+ int& IWARN, int& INFO);
+
+ int F77_FUNC (ib01cd, IB01CD)
+ (char& JOBX0, char& COMUSE, char& JOB,
+ int& N, int& M, int& L,
+ int& NSMP,
+ double* A, int& LDA,
+ double* B, int& LDB,
+ double* C, int& LDC,
+ double* D, int& LDD,
+ double* U, int& LDU,
+ double* Y, int& LDY,
+ double* X0,
+ double* V, int& LDV,
+ double& TOL,
+ int* IWORK,
+ double* DWORK, int& LDWORK,
+ int& IWARN, int& INFO);
+}
+
+// PKG_ADD: autoload ("slident_b", "devel_slicot_functions.oct");
+DEFUN_DLD (slident_b, args, nargout,
+ "-*- texinfo -*-\n\
+Slicot IB01AD Release 5.0\n\
+No argument checking.\n\
+For internal use only.")
+{
+ int nargin = args.length ();
+ octave_value_list retval;
+
+ if (nargin != 15)
+ {
+ print_usage ();
+ }
+ else
+ {
+////////////////////////////////////////////////////////////////////////////////////
+// SLICOT IB01AD - preprocess the input-output data //
+////////////////////////////////////////////////////////////////////////////////////
+
+ // arguments in
+ char meth;
+ char alg;
+ char jobd;
+ char batch;
+ char conct;
+ char ctrl;
+ char metha;
+ char jobda; // ??? unused
+
+ Matrix y = args(0).matrix_value ();
+ Matrix u = args(1).matrix_value ();
+ int nobr = args(2).int_value ();
+ int nuser = args(3).int_value ();
+
+ const int imeth = args(4).int_value ();
+ const int ialg = args(5).int_value ();
+ const int ijobd = args(6).int_value ();
+ const int ibatch = args(7).int_value ();
+ const int iconct = args(8).int_value ();
+ const int ictrl = args(9).int_value ();
+
+ double rcond = args(10).double_value ();
+ double tol = args(11).double_value ();
+ double tolb = args(10).double_value (); // tolb = rcond
+
+ Matrix r = args(12).matrix_value ();
+ Matrix sv = args(13).matrix_value ();
+ int n = args(14).int_value ();
+
+
+ switch (imeth)
+ {
+ case 0:
+ meth = 'M';
+ metha = 'M';
+ break;
+ case 1:
+ meth = 'N';
+ metha = 'N';
+ break;
+ case 2:
+ meth = 'C';
+ metha = 'N'; // no typo here
+ break;
+ default:
+ error ("slib01ad: argument 'meth' invalid");
+ }
+
+ switch (ialg)
+ {
+ case 0:
+ alg = 'C';
+ break;
+ case 1:
+ alg = 'F';
+ break;
+ case 2:
+ alg = 'Q';
+ break;
+ default:
+ error ("slib01ad: argument 'alg' invalid");
+ }
+
+ if (meth == 'C')
+ jobd = 'N';
+ else if (ijobd == 0)
+ jobd = 'M';
+ else
+ jobd = 'N';
+
+ switch (ibatch)
+ {
+ case 0:
+ batch = 'F';
+ break;
+ case 1:
+ batch = 'I';
+ break;
+ case 2:
+ batch = 'L';
+ break;
+ case 3:
+ batch = 'O';
+ break;
+ default:
+ error ("slib01ad: argument 'batch' invalid");
+ }
+
+ if (iconct == 0)
+ conct = 'C';
+ else
+ conct = 'N';
+
+ if (ictrl == 0)
+ ctrl = 'C';
+ else
+ ctrl = 'N';
+
+
+ int m = u.columns (); // m: number of inputs
+ int l = y.columns (); // l: number of outputs
+ int nsmp = y.rows (); // nsmp: number of samples
+ // y.rows == u.rows is checked by iddata class
+ // TODO: check minimal nsmp size
+
+ if (batch == 'O')
+ {
+ if (nsmp < 2*(m+l+1)*nobr - 1)
+ error ("slident: require NSMP >= 2*(M+L+1)*NOBR - 1");
+ }
+ else
+ {
+ if (nsmp < 2*nobr)
+ error ("slident: require NSMP >= 2*NOBR");
+ }
+
+ int ldu;
+
+ if (m == 0)
+ ldu = 1;
+ else // m > 0
+ ldu = nsmp;
+
+ int ldy = nsmp;
+
+ // arguments out
+ //int n;
+ int ldr;
+
+ if (metha == 'M' && jobd == 'M')
+ ldr = max (2*(m+l)*nobr, 3*m*nobr);
+ else if (metha == 'N' || (metha == 'M' && jobd == 'N'))
+ ldr = 2*(m+l)*nobr;
+ else
+ error ("slib01ad: could not handle 'ldr' case");
+
+ //Matrix r (ldr, 2*(m+l)*nobr);
+ //ColumnVector sv (l*nobr);
+
+ // workspace
+ int liwork;
+
+ if (metha == 'N') // if METH = 'N'
+ liwork = (m+l)*nobr;
+ else if (alg == 'F') // if METH = 'M' and ALG = 'F'
+ liwork = m+l;
+ else // if METH = 'M' and ALG = 'C' or 'Q'
+ liwork = 0;
+
+ // TODO: Handle 'k' for DWORK
+
+ int ldwork;
+ int ns = nsmp - 2*nobr + 1;
+
+ if (alg == 'C')
+ {
+ if (batch == 'F' || batch == 'I')
+ {
+ if (conct == 'C')
+ ldwork = (4*nobr-2)*(m+l);
+ else // (conct == 'N')
+ ldwork = 1;
+ }
+ else if (metha == 'M') // && (batch == 'L' || batch == 'O')
+ {
+ if (conct == 'C' && batch == 'L')
+ ldwork = max ((4*nobr-2)*(m+l), 5*l*nobr);
+ else if (jobd == 'M')
+ ldwork = max ((2*m-1)*nobr, (m+l)*nobr, 5*l*nobr);
+ else // (jobd == 'N')
+ ldwork = 5*l*nobr;
+ }
+ else // meth == 'N' && (batch == 'L' || batch == 'O')
+ {
+ ldwork = 5*(m+l)*nobr + 1;
+ }
+ }
+ else if (alg == 'F')
+ {
+ if (batch != 'O' && conct == 'C')
+ ldwork = (m+l)*2*nobr*(m+l+3);
+ else if (batch == 'F' || batch == 'I') // && conct == 'N'
+ ldwork = (m+l)*2*nobr*(m+l+1);
+ else // (batch == 'L' || '0' && conct == 'N')
+ ldwork = (m+l)*4*nobr*(m+l+1)+(m+l)*2*nobr;
+ }
+ else // (alg == 'Q')
+ {
+ // int ns = nsmp - 2*nobr + 1;
+
+ if (ldr >= ns && batch == 'F')
+ {
+ ldwork = 4*(m+l)*nobr;
+ }
+ else if (ldr >= ns && batch == 'O')
+ {
+ if (metha == 'M')
+ ldwork = max (4*(m+l)*nobr, 5*l*nobr);
+ else // (meth == 'N')
+ ldwork = 5*(m+l)*nobr + 1;
+ }
+ else if (conct == 'C' && (batch == 'I' || batch == 'L'))
+ {
+ ldwork = 4*(nobr+1)*(m+l)*nobr;
+ }
+ else // if ALG = 'Q', (BATCH = 'F' or 'O', and LDR < NS), or (BATCH = 'I' or 'L' and CONCT = 'N')
+ {
+ ldwork = 6*(m+l)*nobr;
+ }
+ }
+
+/*
+IB01AD.f Lines 438-445
+C FURTHER COMMENTS
+C
+C For ALG = 'Q', BATCH = 'O' and LDR < NS, or BATCH <> 'O', the
+C calculations could be rather inefficient if only minimal workspace
+C (see argument LDWORK) is provided. It is advisable to provide as
+C much workspace as possible. Almost optimal efficiency can be
+C obtained for LDWORK = (NS+2)*(2*(M+L)*NOBR), assuming that the
+C cache size is large enough to accommodate R, U, Y, and DWORK.
+*/
+
+// warning ("==================== ldwork before: %d =====================", ldwork);
+// ldwork = (ns+2)*(2*(m+l)*nobr);
+ldwork = max (ldwork, (ns+2)*(2*(m+l)*nobr));
+// ldwork *= 3;
+// warning ("==================== ldwork after: %d =====================", ldwork);
+
+
+/*
+IB01AD.f Lines 291-195:
+c the workspace used for alg = 'q' is
+c ldrwrk*2*(m+l)*nobr + 4*(m+l)*nobr,
+c where ldrwrk = ldwork/(2*(m+l)*nobr) - 2; recommended
+c value ldrwrk = ns, assuming a large enough cache size.
+c for good performance, ldwork should be larger.
+
+somehow ldrwrk and ldwork must have been mixed up here
+
+*/
+
+#if 0
+ OCTAVE_LOCAL_BUFFER (int, iwork, liwork);
+ OCTAVE_LOCAL_BUFFER (double, dwork, ldwork);
+
+ // error indicators
+ int iwarn = 0;
+ int info = 0;
+
+
+ // SLICOT routine IB01AD
+ F77_XFCN (ib01ad, IB01AD,
+ (metha, alg, jobd,
+ batch, conct, ctrl,
+ nobr, m, l,
+ nsmp,
+ u.fortran_vec (), ldu,
+ y.fortran_vec (), ldy,
+ n,
+ r.fortran_vec (), ldr,
+ sv.fortran_vec (),
+ rcond, tol,
+ iwork,
+ dwork, ldwork,
+ iwarn, info));
+
+
+ if (f77_exception_encountered)
+ error ("ident: exception in SLICOT subroutine IB01AD");
+
+ static const char* err_msg[] = {
+ "0: OK",
+ "1: a fast algorithm was requested (ALG = 'C', or 'F') "
+ "in sequential data processing, but it failed; the "
+ "routine can be repeatedly called again using the "
+ "standard QR algorithm",
+ "2: the singular value decomposition (SVD) algorithm did "
+ "not converge"};
+
+ static const char* warn_msg[] = {
+ "0: OK",
+ "1: the number of 100 cycles in sequential data "
+ "processing has been exhausted without signaling "
+ "that the last block of data was get; the cycle "
+ "counter was reinitialized",
+ "2: a fast algorithm was requested (ALG = 'C' or 'F'), "
+ "but it failed, and the QR algorithm was then used "
+ "(non-sequential data processing)",
+ "3: all singular values were exactly zero, hence N = 0 "
+ "(both input and output were identically zero)",
+ "4: the least squares problems with coefficient matrix "
+ "U_f, used for computing the weighted oblique "
+ "projection (for METH = 'N'), have a rank-deficient "
+ "coefficient matrix",
+ "5: the least squares problem with coefficient matrix "
+ "r_1 [6], used for computing the weighted oblique "
+ "projection (for METH = 'N'), has a rank-deficient "
+ "coefficient matrix"};
+
+
+ error_msg ("ident", info, 2, err_msg);
+ warning_msg ("ident", iwarn, 5, warn_msg);
+
+
+ // resize
+ int rs = 2*(m+l)*nobr;
+ r.resize (rs, rs);
+
+ if (nuser > 0)
+ {
+ if (nuser < nobr)
+ {
+ n = nuser;
+ // warning ("ident: nuser (%d) < nobr (%d), n = nuser", nuser, nobr);
+ }
+ else
+ error ("ident: 'nuser' invalid");
+ }
+#endif
+////////////////////////////////////////////////////////////////////////////////////
+// SLICOT IB01BD - estimating system matrices, Kalman gain, and covariances //
+////////////////////////////////////////////////////////////////////////////////////
+
+ // arguments in
+ char job = 'A';
+ char jobck = 'K';
+
+ // TODO: if meth == 'C', which meth should be taken for IB01AD.f, 'M' or 'N'?
+
+ int nsmpl = nsmp;
+
+ if (nsmpl < 2*(m+l)*nobr)
+ error ("slident: nsmpl (%d) < 2*(m+l)*nobr (%d)", nsmpl, nobr);
+
+ // arguments out
+ int lda = max (1, n);
+ int ldc = max (1, l);
+ int ldb = max (1, n);
+ int ldd = max (1, l);
+ int ldq = n; // if JOBCK = 'C' or 'K'
+ int ldry = l; // if JOBCK = 'C' or 'K'
+ int lds = n; // if JOBCK = 'C' or 'K'
+ int ldk = n; // if JOBCK = 'K'
+
+ Matrix a (lda, n);
+ Matrix c (ldc, n);
+ Matrix b (ldb, m);
+ Matrix d (ldd, m);
+
+ Matrix q (ldq, n);
+ Matrix ry (ldry, l);
+ Matrix s (lds, l);
+ Matrix k (ldk, l);
+
+ // workspace
+ int liwork_b;
+ int liw1;
+ int liw2;
+
+ liw1 = max (n, m*nobr+n, l*nobr, m*(n+l));
+ liw2 = n*n; // if JOBCK = 'K'
+ liwork_b = max (liw1, liw2);
+
+ int ldwork_b;
+ int ldw1;
+ int ldw2;
+ int ldw3;
+/*
+ if (meth == 'M')
+ {
+ int ldw1a = max (2*(l*nobr-l)*n+2*n, (l*nobr-l)*n+n*n+7*n);
+ int ldw1b = max (2*(l*nobr-l)*n+n*n+7*n,
+ (l*nobr-l)*n+n+6*m*nobr,
+ (l*nobr-l)*n+n+max (l+m*nobr, l*nobr + max (3*l*nobr+1, m)));
+ ldw1 = max (ldw1a, ldw1b);
+
+ int aw;
+
+ if (m == 0 || job == 'C')
+ aw = n + n*n;
+ else
+ aw = 0;
+
+ ldw2 = l*nobr*n + max ((l*nobr-l)*n+aw+2*n+max(5*n,(2*m+l)*nobr+l), 4*(m*nobr+n)+1, m*nobr+2*n+l );
+ }
+ else if (meth == 'N')
+ {
+ ldw1 = l*nobr*n + max ((l*nobr-l)*n+2*n+(2*m+l)*nobr+l,
+ 2*(l*nobr-l)*n+n*n+8*n,
+ n+4*(m*nobr+n)+1,
+ m*nobr+3*n+l);
+
+ if (m == 0 || job == 'C')
+ ldw2 = 0;
+ else
+ ldw2 = l*nobr*n+m*nobr*(n+l)*(m*(n+l)+1)+ max ((n+l)*(n+l), 4*m*(n+l)+1);
+
+ }
+ else // (meth == 'C')
+ {
+ int ldw1a = max (2*(l*nobr-l)*n+2*n, (l*nobr-l)*n+n*n+7*n);
+ int ldw1b = l*nobr*n + max ((l*nobr-l)*n+2*n+(2*m+l)*nobr+l,
+ 2*(l*nobr-l)*n+n*n+8*n,
+ n+4*(m*nobr+n)+1,
+ m*nobr+3*n+l);
+
+ ldw1 = max (ldw1a, ldw1b);
+
+ ldw2 = l*nobr*n+m*nobr*(n+l)*(m*(n+l)+1)+ max ((n+l)*(n+l), 4*m*(n+l)+1);
+
+ }
+*/
+
+ int ldw1ax = max (2*(l*nobr-l)*n+2*n, (l*nobr-l)*n+n*n+7*n);
+ int ldw1bx = max (2*(l*nobr-l)*n+n*n+7*n,
+ (l*nobr-l)*n+n+6*m*nobr,
+ (l*nobr-l)*n+n+max (l+m*nobr, l*nobr + max (3*l*nobr+1, m)));
+ int ldw1x = max (ldw1ax, ldw1bx);
+
+ int aw;
+
+ if (m == 0 || job == 'C')
+ aw = n + n*n;
+ else
+ aw = 0;
+
+ int ldw2x = l*nobr*n + max ((l*nobr-l)*n+aw+2*n+max(5*n,(2*m+l)*nobr+l), 4*(m*nobr+n)+1, m*nobr+2*n+l );
+
+
+
+ int ldw1y = l*nobr*n + max ((l*nobr-l)*n+2*n+(2*m+l)*nobr+l,
+ 2*(l*nobr-l)*n+n*n+8*n,
+ n+4*(m*nobr+n)+1,
+ m*nobr+3*n+l);
+ int ldw2y;
+ if (m == 0 || job == 'C')
+ int ldw2y = 0;
+ else
+ int ldw2y = l*nobr*n+m*nobr*(n+l)*(m*(n+l)+1)+ max ((n+l)*(n+l), 4*m*(n+l)+1);
+
+
+ int ldw1az = max (2*(l*nobr-l)*n+2*n, (l*nobr-l)*n+n*n+7*n);
+ int ldw1bz = l*nobr*n + max ((l*nobr-l)*n+2*n+(2*m+l)*nobr+l,
+ 2*(l*nobr-l)*n+n*n+8*n,
+ n+4*(m*nobr+n)+1,
+ m*nobr+3*n+l);
+
+ int ldw1z = max (ldw1az, ldw1bz);
+
+ int ldw2z = l*nobr*n+m*nobr*(n+l)*(m*(n+l)+1)+ max ((n+l)*(n+l), 4*m*(n+l)+1);
+
+
+ ldw1 = max (ldw1x, ldw1y, ldw1z);
+ ldw2 = max (ldw2x, ldw2y, ldw2z);
+
+
+
+ ldw3 = max(4*n*n + 2*n*l + l*l + max (3*l, n*l), 14*n*n + 12*n + 5);
+ ldwork_b = max (ldw1, ldw2, ldw3);
+
+ //
+ ldwork_b *= 3;
+
+ OCTAVE_LOCAL_BUFFER (int, iwork_b, liwork_b);
+ OCTAVE_LOCAL_BUFFER (double, dwork_b, ldwork_b);
+ OCTAVE_LOCAL_BUFFER (bool, bwork, 2*n);
+
+
+ // error indicators
+ int iwarn_b = 0;
+ int info_b = 0;
+
+
+ // SLICOT routine IB01BD
+ F77_XFCN (ib01bd, IB01BD,
+ (meth, job, jobck,
+ nobr, n, m, l,
+ nsmpl,
+ r.fortran_vec (), ldr,
+ a.fortran_vec (), lda,
+ c.fortran_vec (), ldc,
+ b.fortran_vec (), ldb,
+ d.fortran_vec (), ldd,
+ q.fortran_vec (), ldq,
+ ry.fortran_vec (), ldry,
+ s.fortran_vec (), lds,
+ k.fortran_vec (), ldk,
+ tolb,
+ iwork_b,
+ dwork_b, ldwork_b,
+ bwork,
+ iwarn_b, info_b));
+
+
+ if (f77_exception_encountered)
+ error ("ident: exception in SLICOT subroutine IB01BD");
+
+ static const char* err_msg_b[] = {
+ "0: OK",
+ "1: error message not specified",
+ "2: the singular value decomposition (SVD) algorithm did "
+ "not converge",
+ "3: a singular upper triangular matrix was found",
+ "4: matrix A is (numerically) singular in discrete-"
+ "time case",
+ "5: the Hamiltonian or symplectic matrix H cannot be "
+ "reduced to real Schur form",
+ "6: the real Schur form of the Hamiltonian or "
+ "symplectic matrix H cannot be appropriately ordered",
+ "7: the Hamiltonian or symplectic matrix H has less "
+ "than N stable eigenvalues",
+ "8: the N-th order system of linear algebraic "
+ "equations, from which the solution matrix X would "
+ "be obtained, is singular to working precision",
+ "9: the QR algorithm failed to complete the reduction "
+ "of the matrix Ac to Schur canonical form, T",
+ "10: the QR algorithm did not converge"};
+
+ static const char* warn_msg_b[] = {
+ "0: OK",
+ "1: warning message not specified",
+ "2: warning message not specified",
+ "3: warning message not specified",
+ "4: a least squares problem to be solved has a "
+ "rank-deficient coefficient matrix",
+ "5: the computed covariance matrices are too small. "
+ "The problem seems to be a deterministic one; the "
+ "gain matrix is set to zero"};
+
+
+ error_msg ("ident", info_b, 10, err_msg_b);
+ warning_msg ("ident", iwarn_b, 5, warn_msg_b);
+
+ // resize
+ a.resize (n, n);
+ c.resize (l, n);
+ b.resize (n, m);
+ d.resize (l, m);
+
+ q.resize (n, n);
+ ry.resize (l, l);
+ s.resize (n, l);
+ k.resize (n, l);
+
+
+
+ // return values
+ retval(0) = a;
+ retval(1) = b;
+ retval(2) = c;
+ retval(3) = d;
+
+ retval(4) = q;
+ retval(5) = ry;
+ retval(6) = s;
+ retval(7) = k;
+ }
+
+ return retval;
+}
Added: trunk/octave-forge/extra/control-devel/src/slident_c.cc
===================================================================
--- trunk/octave-forge/extra/control-devel/src/slident_c.cc (rev 0)
+++ trunk/octave-forge/extra/control-devel/src/slident_c.cc 2012-05-05 13:03:58 UTC (rev 10365)
@@ -0,0 +1,473 @@
+/*
+
+Copyright (C) 2012 Lukas F. Reichlin
+
+This file is part of LTI Syncope.
+
+LTI Syncope is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+LTI Syncope is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>.
+
+SLICOT system identification
+Uses SLICOT IB01AD, IB01BD and IB01CD by courtesy of NICONET e.V.
+<http://www.slicot.org>
+
+Author: Lukas Reichlin <luk...@gm...>
+Created: March 2012
+Version: 0.1
+
+*/
+
+#include <octave/oct.h>
+#include <f77-fcn.h>
+#include "common.h"
+
+extern "C"
+{
+ int F77_FUNC (ib01ad, IB01AD)
+ (char& METH, char& ALG, char& JOBD,
+ char& BATCH, char& CONCT, char& CTRL,
+ int& NOBR, int& M, int& L,
+ int& NSMP,
+ double* U, int& LDU,
+ double* Y, int& LDY,
+ int& N,
+ double* R, int& LDR,
+ double* SV,
+ double& RCOND, double& TOL,
+ int* IWORK,
+ double* DWORK, int& LDWORK,
+ int& IWARN, int& INFO);
+
+ int F77_FUNC (ib01bd, IB01BD)
+ (char& METH, char& JOB, char& JOBCK,
+ int& NOBR, int& N, int& M, int& L,
+ int& NSMPL,
+ double* R, int& LDR,
+ double* A, int& LDA,
+ double* C, int& LDC,
+ double* B, int& LDB,
+ double* D, int& LDD,
+ double* Q, int& LDQ,
+ double* RY, int& LDRY,
+ double* S, int& LDS,
+ double* K, int& LDK,
+ double& TOL,
+ int* IWORK,
+ double* DWORK, int& LDWORK,
+ bool* BWORK,
+ int& IWARN, int& INFO);
+
+ int F77_FUNC (ib01cd, IB01CD)
+ (char& JOBX0, char& COMUSE, char& JOB,
+ int& N, int& M, int& L,
+ int& NSMP,
+ double* A, int& LDA,
+ double* B, int& LDB,
+ double* C, int& LDC,
+ double* D, int& LDD,
+ double* U, int& LDU,
+ double* Y, int& LDY,
+ double* X0,
+ double* V, int& LDV,
+ double& TOL,
+ int* IWORK,
+ double* DWORK, int& LDWORK,
+ int& IWARN, int& INFO);
+}
+
+// PKG_ADD: autoload ("slident_c", "devel_slicot_functions.oct");
+DEFUN_DLD (slident_c, args, nargout,
+ "-*- texinfo -*-\n\
+Slicot IB01AD Release 5.0\n\
+No argument checking.\n\
+For internal use only.")
+{
+ int nargin = args.length ();
+ octave_value_list retval;
+
+ if (nargin != 16)
+ {
+ print_usage ();
+ }
+ else
+ {
+////////////////////////////////////////////////////////////////////////////////////
+// SLICOT IB01AD - preprocess the input-output data //
+////////////////////////////////////////////////////////////////////////////////////
+
+ // arguments in
+ char meth;
+ char alg;
+ char jobd;
+ char batch;
+ char conct;
+ char ctrl;
+ char metha;
+ char jobda; // ??? unused
+
+ Matrix y = args(0).matrix_value ();
+ Matrix u = args(1).matrix_value ();
+ int nobr = args(2).int_value ();
+ int nuser = args(3).int_value ();
+
+ const int imeth = args(4).int_value ();
+ const int ialg = args(5).int_value ();
+ const int ijobd = args(6).int_value ();
+ const int ibatch = args(7).int_value ();
+ const int iconct = args(8).int_value ();
+ const int ictrl = args(9).int_value ();
+
+ double rcond = args(10).double_value ();
+ double tol = args(11).double_value ();
+ double tolb = args(10).double_value (); // tolb = rcond
+
+ Matrix a = args(12).matrix_value ();
+ Matrix b = args(13).matrix_value ();
+ Matrix c = args(14).matrix_value ();
+ Matrix d = args(15).matrix_value ();
+
+
+ switch (imeth)
+ {
+ case 0:
+ meth = 'M';
+ metha = 'M';
+ break;
+ case 1:
+ meth = 'N';
+ metha = 'N';
+ break;
+ case 2:
+ meth = 'C';
+ metha = 'N'; // no typo here
+ break;
+ default:
+ error ("slib01ad: argument 'meth' invalid");
+ }
+
+ switch (ialg)
+ {
+ case 0:
+ alg = 'C';
+ break;
+ case 1:
+ alg = 'F';
+ break;
+ case 2:
+ alg = 'Q';
+ break;
+ default:
+ error ("slib01ad: argument 'alg' invalid");
+ }
+
+ if (meth == 'C')
+ jobd = 'N';
+ else if (ijobd == 0)
+ jobd = 'M';
+ else
+ jobd = 'N';
+
+ switch (ibatch)
+ {
+ case 0:
+ batch = 'F';
+ break;
+ case 1:
+ batch = 'I';
+ break;
+ case 2:
+ batch = 'L';
+ break;
+ case 3:
+ batch = 'O';
+ break;
+ default:
+ error ("slib01ad: argument 'batch' invalid");
+ }
+
+ if (iconct == 0)
+ conct = 'C';
+ else
+ conct = 'N';
+
+ if (ictrl == 0)
+ ctrl = 'C';
+ else
+ ctrl = 'N';
+
+
+ int m = u.columns (); // m: number of inputs
+ int l = y.columns (); // l: number of outputs
+ int nsmp = y.rows (); // nsmp: number of samples
+ // y.rows == u.rows is checked by iddata class
+ // TODO: check minimal nsmp size
+
+ if (batch == 'O')
+ {
+ if (nsmp < 2*(m+l+1)*nobr - 1)
+ error ("slident: require NSMP >= 2*(M+L+1)*NOBR - 1");
+ }
+ else
+ {
+ if (nsmp < 2*nobr)
+ error ("slident: require NSMP >= 2*NOBR");
+ }
+
+ int ldu;
+
+ if (m == 0)
+ ldu = 1;
+ else // m > 0
+ ldu = nsmp;
+
+ int ldy = nsmp;
+
+ // arguments out
+ int n;
+ int ldr;
+
+ if (metha == 'M' && jobd == 'M')
+ ldr = max (2*(m+l)*nobr, 3*m*nobr);
+ else if (metha == 'N' || (metha == 'M' && jobd == 'N'))
+ ldr = 2*(m+l)*nobr;
+ else
+ error ("slib01ad: could not handle 'ldr' case");
+
+ Matrix r (ldr, 2*(m+l)*nobr);
+ ColumnVector sv (l*nobr);
+
+ // workspace
+ int liwork;
+
+ if (metha == 'N') // if METH = 'N'
+ liwork = (m+l)*nobr;
+ else if (alg == 'F') // if METH = 'M' and ALG = 'F'
+ liwork = m+l;
+ else // if METH = 'M' and ALG = 'C' or 'Q'
+ liwork = 0;
+
+ // TODO: Handle 'k' for DWORK
+
+ int ldwork;
+ int ns = nsmp - 2*nobr + 1;
+
+ if (alg == 'C')
+ {
+ if (batch == 'F' || batch == 'I')
+ {
+ if (conct == 'C')
+ ldwork = (4*nobr-2)*(m+l);
+ else // (conct == 'N')
+ ldwork = 1;
+ }
+ else if (metha == 'M') // && (batch == 'L' || batch == 'O')
+ {
+ if (conct == 'C' && batch == 'L')
+ ldwork = max ((4*nobr-2)*(m+l), 5*l*nobr);
+ else if (jobd == 'M')
+ ldwork = max ((2*m-1)*nobr, (m+l)*nobr, 5*l*nobr);
+ else // (jobd == 'N')
+ ldwork = 5*l*nobr;
+ }
+ else // meth == 'N' && (batch == 'L' || batch == 'O')
+ {
+ ldwork = 5*(m+l)*nobr + 1;
+ }
+ }
+ else if (alg == 'F')
+ {
+ if (batch != 'O' && conct == 'C')
+ ldwork = (m+l)*2*nobr*(m+l+3);
+ else if (batch == 'F' || batch == 'I') // && conct == 'N'
+ ldwork = (m+l)*2*nobr*(m+l+1);
+ else // (batch == 'L' || '0' && conct == 'N')
+ ldwork = (m+l)*4*nobr*(m+l+1)+(m+l)*2*nobr;
+ }
+ else // (alg == 'Q')
+ {
+ // int ns = nsmp - 2*nobr + 1;
+
+ if (ldr >= ns && batch == 'F')
+ {
+ ldwork = 4*(m+l)*nobr;
+ }
+ else if (ldr >= ns && batch == 'O')
+ {
+ if (metha == 'M')
+ ldwork = max (4*(m+l)*nobr, 5*l*nobr);
+ else // (meth == 'N')
+ ldwork = 5*(m+l)*nobr + 1;
+ }
+ else if (conct == 'C' && (batch == 'I' || batch == 'L'))
+ {
+ ldwork = 4*(nobr+1)*(m+l)*nobr;
+ }
+ else // if ALG = 'Q', (BATCH = 'F' or 'O', and LDR < NS), or (BATCH = 'I' or 'L' and CONCT = 'N')
+ {
+ ldwork = 6*(m+l)*nobr;
+ }
+ }
+
+/*
+IB01AD.f Lines 438-445
+C FURTHER COMMENTS
+C
+C For ALG = 'Q', BATCH = 'O' and LDR < NS, or BATCH <> 'O', the
+C calculations could be rather inefficient if only minimal workspace
+C (see argument LDWORK) is provided. It is advisable to provide as
+C much workspace as possible. Almost optimal efficiency can be
+C obtained for LDWORK = (NS+2)*(2*(M+L)*NOBR), assuming tha...
[truncated message content] |
|
From: <par...@us...> - 2012-05-05 06:48:48
|
Revision: 10364
http://octave.svn.sourceforge.net/octave/?rev=10364&view=rev
Author: paramaniac
Date: 2012-05-05 06:48:41 +0000 (Sat, 05 May 2012)
Log Message:
-----------
control: explain difference between ctranspose and transpose more clearly
Modified Paths:
--------------
trunk/octave-forge/main/control/inst/@lti/ctranspose.m
trunk/octave-forge/main/control/inst/@lti/transpose.m
Modified: trunk/octave-forge/main/control/inst/@lti/ctranspose.m
===================================================================
--- trunk/octave-forge/main/control/inst/@lti/ctranspose.m 2012-05-05 00:05:09 UTC (rev 10363)
+++ trunk/octave-forge/main/control/inst/@lti/ctranspose.m 2012-05-05 06:48:41 UTC (rev 10364)
@@ -16,7 +16,16 @@
## along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
-## Conjugate transpose of LTI objects. Used by Octave for "sys'".
+## Conjugate transpose or pertransposition of LTI objects.
+## Used by Octave for "sys'".
+## For a transfer-function matrix G, G' denotes the conjugate
+## of G given by G.'(-s) for a continuous-time system or G.'(1/z)
+## for a discrete-time system.
+## The frequency response of the pertransposition of G is the
+## Hermitian (conjugate) transpose of G(jw), i.e.
+## freqresp (G', w) = freqresp (G, w)'.
+## @strong{WARNING:} Do @strong{NOT} use this for dual problems,
+## use the transpose "sys.'" (note the dot) instead.
## Author: Lukas Reichlin <luk...@gm...>
## Created: May 2012
Modified: trunk/octave-forge/main/control/inst/@lti/transpose.m
===================================================================
--- trunk/octave-forge/main/control/inst/@lti/transpose.m 2012-05-05 00:05:09 UTC (rev 10363)
+++ trunk/octave-forge/main/control/inst/@lti/transpose.m 2012-05-05 06:48:41 UTC (rev 10364)
@@ -17,6 +17,8 @@
## -*- texinfo -*-
## Transpose of LTI objects. Used by Octave for "sys.'".
+## Useful for dual problems, i.e. controllability and observability
+## or designing estimator gains with @command{lqr} and @command {place}.
## Author: Lukas Reichlin <luk...@gm...>
## Created: February 2010
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <mtm...@us...> - 2012-05-05 00:05:15
|
Revision: 10363
http://octave.svn.sourceforge.net/octave/?rev=10363&view=rev
Author: mtmiller
Date: 2012-05-05 00:05:09 +0000 (Sat, 05 May 2012)
Log Message:
-----------
comm: marcumq function has been moved to the signal package
Modified Paths:
--------------
trunk/octave-forge/main/comm/INDEX
trunk/octave-forge/main/comm/NEWS
Removed Paths:
-------------
trunk/octave-forge/main/comm/inst/marcumq.m
Modified: trunk/octave-forge/main/comm/INDEX
===================================================================
--- trunk/octave-forge/main/comm/INDEX 2012-05-05 00:04:13 UTC (rev 10362)
+++ trunk/octave-forge/main/comm/INDEX 2012-05-05 00:05:09 UTC (rev 10363)
@@ -177,4 +177,3 @@
vec2mat
qfunc
qfuncinv
- marcumq
Modified: trunk/octave-forge/main/comm/NEWS
===================================================================
--- trunk/octave-forge/main/comm/NEWS 2012-05-05 00:04:13 UTC (rev 10362)
+++ trunk/octave-forge/main/comm/NEWS 2012-05-05 00:05:09 UTC (rev 10363)
@@ -1,6 +1,8 @@
Summary of important user-visible changes for communications 1.1.1:
-------------------------------------------------------------------
+ ** The function `marcumq' has been moved to the signal package
+
** The following functions are new:
convenc
Deleted: trunk/octave-forge/main/comm/inst/marcumq.m
===================================================================
--- trunk/octave-forge/main/comm/inst/marcumq.m 2012-05-05 00:04:13 UTC (rev 10362)
+++ trunk/octave-forge/main/comm/inst/marcumq.m 2012-05-05 00:05:09 UTC (rev 10363)
@@ -1,49 +0,0 @@
-## Copyright (C) 2007 Sylvain Pelissier <syl...@gm...>
-##
-## This program is free software; you can redistribute it and/or modify it under
-## the terms of the GNU General Public License as published by the Free Software
-## Foundation; either version 3 of the License, or (at your option) any later
-## version.
-##
-## This program is distributed in the hope that it will be useful, but WITHOUT
-## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-## details.
-##
-## You should have received a copy of the GNU General Public License along with
-## this program; if not, see <http://www.gnu.org/licenses/>.
-
-## -*- texinfo -*-
-## @deftypefn {Function File} {[@var{Q}] =} marcumq (@var{a}, @var{b}, @var{m})
-## Compute the Marcum Q function.
-## @end deftypefn
-
-function Q = marcumq(a,b,m)
-
- if nargin < 2
- usage(" marcumq(a,b)\n\tmarcumq(a,b,m)");
- end
-
- if nargin < 3
- m = 1;
- end
-
- l = max([size(a,1) size(b,1) size(m,1)]);
- c = max([size(a,2) size(b,2) size(m,2)]);
- a = padarray(a,[l-size(a,1) c-size(a,2)],'replicate','post')
- b = padarray(b,[l-size(b,1) c-size(b,2)],'replicate','post');
- m = padarray(m,[l-size(m,1) c-size(m,2)],'replicate','post');
-
- if (min(a(:)) < 0 || min(b(:)) < 0)
- error("a and b must be positive");
- end
-
- for n = 1:l
- for o = 1:c
- k = 1-m(n,o):100;
- v = (a(n,o)./b(n,o)).^k.*besseli(k,a(n,o).*b(n,o));
- v = sum(v);
- Q(n,o) = v.*exp(-(a(n,o).^2+b(n,o).^2)/2);
- end
- end
-endfunction
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <mtm...@us...> - 2012-05-05 00:04:19
|
Revision: 10362
http://octave.svn.sourceforge.net/octave/?rev=10362&view=rev
Author: mtmiller
Date: 2012-05-05 00:04:13 +0000 (Sat, 05 May 2012)
Log Message:
-----------
signal: add marcumq to function list
Modified Paths:
--------------
trunk/octave-forge/main/signal/INDEX
Modified: trunk/octave-forge/main/signal/INDEX
===================================================================
--- trunk/octave-forge/main/signal/INDEX 2012-05-04 23:57:42 UTC (rev 10361)
+++ trunk/octave-forge/main/signal/INDEX 2012-05-05 00:04:13 UTC (rev 10362)
@@ -155,6 +155,7 @@
Utility
buffer
fracshift
+ marcumq
wkeep
wrev
zerocrossing
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <mtm...@us...> - 2012-05-04 23:57:49
|
Revision: 10361
http://octave.svn.sourceforge.net/octave/?rev=10361&view=rev
Author: mtmiller
Date: 2012-05-04 23:57:42 +0000 (Fri, 04 May 2012)
Log Message:
-----------
marcumq: update license to GPLv3+, revert to padarray until tablify is available, and clean up help
Modified Paths:
--------------
trunk/octave-forge/main/signal/inst/marcumq.m
Modified: trunk/octave-forge/main/signal/inst/marcumq.m
===================================================================
--- trunk/octave-forge/main/signal/inst/marcumq.m 2012-05-04 23:46:50 UTC (rev 10360)
+++ trunk/octave-forge/main/signal/inst/marcumq.m 2012-05-04 23:57:42 UTC (rev 10361)
@@ -1,70 +1,68 @@
## Copyright (C) 2012 Robert T. Short <rt...@ie...>
##
-## This is free software; you can redistribute it and/or modify
-## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2 of the License, or
-## (at your option) any later version.
+## This program is free software; you can redistribute it and/or modify it under
+## the terms of the GNU General Public License as published by the Free Software
+## Foundation; either version 3 of the License, or (at your option) any later
+## version.
##
-## This software is distributed in the hope that it will be useful,
-## but WITHOUT ANY WARRANTY; without even the implied warranty of
-## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-## GNU General Public License for more details.
+## This program is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+## details.
##
-## You should have received a copy of the GNU General Public License
-## along with this program; If not, see <http://www.gnu.org/licenses/>.
+## You should have received a copy of the GNU General Public License along with
+## this program; if not, see <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
-## @deftypefn {Function File} {[@var{Q}] =} marcumq (@var{a}, @var{b}, @var{m}, @var{tol}=eps)
+## @deftypefn {Function File} {@var{Q} = } marcumq (@var{a}, @var{b})
+## @deftypefnx {Function File} {@var{Q} = } marcumq (@var{a}, @var{b}, @var{m})
+## @deftypefnx {Function File} {@var{Q} = } marcumq (@var{a}, @var{b}, @var{m}, @var{tol})
##
-## @noindent
## Compute the generalized Marcum Q function of order @var{M} with
-## noncentrality parameter @var{a} and argument @var{b}. An optional
-## relative tolerance @var{tol} may be included.
+## noncentrality parameter @var{a} and argument @var{b}. If the order
+## @var{M} is omitted it defaults to 1. An optional relative tolerance
+## @var{tol} may be included, the default is @code{eps}.
##
-## @noindent
## If the input arguments are commensurate vectors, this function
-## will produce a table of values (see tablify).
+## will produce a table of values.
##
-## @noindent
## This function computes Marcum's Q function using the infinite
## Bessel series, truncated when the relative error is less than
## the specified tolerance. The accuracy is limited by that of
## the Bessel functions, so reducing the tolerance is probably
## not useful.
##
-## @noindent
## Reference: Marcum, "Tables of Q Functions", Rand Corporation.
##
-## @noindent
## Reference: R.T. Short, "Computation of Noncentral Chi-squared
## and Rice Random Variables", www.phaselockedsystems.com/publications
##
-## @seealso{imarcumq,tablify}
-##
## @end deftypefn
-## Author: Robert T. Short (rt...@ie...)
-## Description: Generalized Marcum Q function.
-function [Q] = marcumq(a,b,M=1,tol=eps)
+function Q = marcumq(a,b,M=1,tol=eps)
if ( (nargin<2) || (nargin>4) )
- usage('[Q] = marcumq(a,b,M=1,tol=eps)');
+ print_usage();
end
if ( any(a<0) || any(b<0) )
- error('Parameters to marcumq must be positive');
+ error("Parameters to marcumq must be positive");
end
if ( any(M<0) || any(floor(M)~=M) )
error("M must be a positive integer");
end
- [a,b,M] = tablify(a,b,M);
+ nr = max([size(a,1) size(b,1) size(M,1)]);
+ nc = max([size(a,2) size(b,2) size(M,2)]);
+ a = padarray(a, [nr - size(a,1) nc - size(a,2)], "replicate", "post");
+ b = padarray(b, [nr - size(b,1) nc - size(b,2)], "replicate", "post");
+ M = padarray(M, [nr - size(M,1) nc - size(M,2)], "replicate", "post");
Q = arrayfun(@mq, a,b,M,tol);
end
% Subfunction to compute the actual Marcum Q function.
-function [Q] = mq(a,b,M,tol)
+function Q = mq(a,b,M,tol)
% Special cases.
if (b==0)
Q = 1;
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <mtm...@us...> - 2012-05-04 23:46:59
|
Revision: 10360
http://octave.svn.sourceforge.net/octave/?rev=10360&view=rev
Author: mtmiller
Date: 2012-05-04 23:46:50 +0000 (Fri, 04 May 2012)
Log Message:
-----------
marcumq: improved rewrite added to signal package, submitted by Robert T. Short <rt...@ie...>. Closes artifact #3522119
Modified Paths:
--------------
trunk/octave-forge/main/signal/NEWS
Added Paths:
-----------
trunk/octave-forge/main/signal/inst/marcumq.m
Modified: trunk/octave-forge/main/signal/NEWS
===================================================================
--- trunk/octave-forge/main/signal/NEWS 2012-05-04 09:09:56 UTC (rev 10359)
+++ trunk/octave-forge/main/signal/NEWS 2012-05-04 23:46:50 UTC (rev 10360)
@@ -4,6 +4,10 @@
signal-x.y.z Release Date: yyyy-mm-dd Release Manager:
===============================================================================
+ ** The function `marcumq' was imported from the communications package and has
+ been completely rewritten to improve performance and fix computational
+ errors.
+
** Package is no longer automatically loaded
** The functions `__ellip_ws' and `__ellip_ws_min' have been removed (they
Added: trunk/octave-forge/main/signal/inst/marcumq.m
===================================================================
--- trunk/octave-forge/main/signal/inst/marcumq.m (rev 0)
+++ trunk/octave-forge/main/signal/inst/marcumq.m 2012-05-04 23:46:50 UTC (rev 10360)
@@ -0,0 +1,391 @@
+## Copyright (C) 2012 Robert T. Short <rt...@ie...>
+##
+## This is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2 of the License, or
+## (at your option) any later version.
+##
+## This software is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{Q}] =} marcumq (@var{a}, @var{b}, @var{m}, @var{tol}=eps)
+##
+## @noindent
+## Compute the generalized Marcum Q function of order @var{M} with
+## noncentrality parameter @var{a} and argument @var{b}. An optional
+## relative tolerance @var{tol} may be included.
+##
+## @noindent
+## If the input arguments are commensurate vectors, this function
+## will produce a table of values (see tablify).
+##
+## @noindent
+## This function computes Marcum's Q function using the infinite
+## Bessel series, truncated when the relative error is less than
+## the specified tolerance. The accuracy is limited by that of
+## the Bessel functions, so reducing the tolerance is probably
+## not useful.
+##
+## @noindent
+## Reference: Marcum, "Tables of Q Functions", Rand Corporation.
+##
+## @noindent
+## Reference: R.T. Short, "Computation of Noncentral Chi-squared
+## and Rice Random Variables", www.phaselockedsystems.com/publications
+##
+## @seealso{imarcumq,tablify}
+##
+## @end deftypefn
+
+## Author: Robert T. Short (rt...@ie...)
+## Description: Generalized Marcum Q function.
+function [Q] = marcumq(a,b,M=1,tol=eps)
+
+ if ( (nargin<2) || (nargin>4) )
+ usage('[Q] = marcumq(a,b,M=1,tol=eps)');
+ end
+ if ( any(a<0) || any(b<0) )
+ error('Parameters to marcumq must be positive');
+ end
+ if ( any(M<0) || any(floor(M)~=M) )
+ error("M must be a positive integer");
+ end
+
+ [a,b,M] = tablify(a,b,M);
+
+ Q = arrayfun(@mq, a,b,M,tol);
+
+end
+
+% Subfunction to compute the actual Marcum Q function.
+function [Q] = mq(a,b,M,tol)
+ % Special cases.
+ if (b==0)
+ Q = 1;
+ N = 0;
+ return;
+ end
+ if (a==0)
+ k = 0:(M-1);
+ Q = exp(-b^2/2)*sum(b.^(2*k)./(2.^k .* factorial(k)));
+ N = 0;
+ return;
+ end
+
+ % The basic iteration. If a<b compute Q_M, otherwise
+ % compute 1-Q_M.
+ k = M;
+ z = a*b;
+ t = 1;
+ k = 0;
+ if (a<b)
+ s = +1;
+ c = 0;
+ x = a/b;
+ d = x;
+ S = besseli(0,z,1);
+ for k=1:M-1
+ t = (d+1/d)*besseli(k,z,1);
+ S = S + t;
+ d = d*x;
+ end
+ N=k++;
+ else
+ s = -1;
+ c = 1;
+ x = b/a;
+ k = M;
+ d = x^M;
+ S = 0;
+ N = 0;
+ end
+
+ do
+ t = d*besseli(abs(k),z,1);
+ S = S + t;
+ d = d*x;
+ N = k++;
+ until (abs(t/S)<tol)
+ Q = c + s*exp( -(a-b)^2/2 )*S;
+end
+
+% Tests for number and validity of arguments.
+%
+%!error
+%! fail(marcumq(1))
+%!error
+%! fail(marcumq(-1,1,1,1,1))
+%!error
+%! fail(marcumq(-1,1))
+%!error
+%! fail(marcumq(1,-1))
+%!error
+%! fail(marcumq(1,1,-1))
+%!error
+%! fail(marcumq(1,1,1.1))
+
+% Notes on tests and accuracy.
+% -----------------------------------
+% The numbers used as the reference (Q) in the tables below are
+% from J.I. Marcum, "Table of Q Functions", Rand Technical Report
+% RM-339, 1950/1/1.
+%
+% There is one discrepency in the tables. Marcum has
+% Q(14.00,17.10) = 0.001078
+% while we compute
+% Q(14.00,17.10) = 0.0010785053 = 0.001079
+% This is obviously a non-problem.
+%
+% As further tests, I created several different versions of the
+% Q function computation, including a Bessel series expansion and
+% numerical integration. All of them agree to with 10^(-16).
+
+%!test
+%! a = [0.00;0.05;1.00;2.00;3.00;4.00;5.00;6.00;7.00;8.00;9.00;10.00;
+%! 11.00;12.00;13.00;14.00;15.00;16.00;17.00;18.00;19.00;20.00;
+%! 21.00;22.00;23.00;24.00];
+%! b = [0.000000, 0.100000, 1.100000, 2.100000, 3.100000, 4.100000];
+%! Q = [1.000000, 0.995012, 0.546074, 0.110251, 0.008189, 0.000224;
+%! 1.000000, 0.995019, 0.546487, 0.110554, 0.008238, 0.000226;
+%! 1.000000, 0.996971, 0.685377, 0.233113, 0.034727, 0.002092;
+%! 1.000000, 0.999322, 0.898073, 0.561704, 0.185328, 0.027068;
+%! 1.000000, 0.999944, 0.985457, 0.865241, 0.526735, 0.169515;
+%! 1.000000, 0.999998, 0.999136, 0.980933, 0.851679, 0.509876;
+%! 1.000000, 1.000000, 0.999979, 0.998864, 0.978683, 0.844038;
+%! 1.000000, 1.000000, 1.000000, 0.999973, 0.998715, 0.977300;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 0.999969, 0.998618;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999966;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000];
+%! q = marcumq(a,b);
+%! assert(q,Q,1e-6);
+
+%!test
+%! a = [0.00;0.05;1.00;2.00;3.00;4.00;5.00;6.00;7.00;8.00;9.00;10.00;
+%! 11.00;12.00;13.00;14.00;15.00;16.00;17.00;18.00;19.00;20.00;
+%! 21.00;22.00;23.00;24.00];
+%! b = [5.100000, 6.100000, 7.100000, 8.100000, 9.100000, 10.10000];
+%! Q = [0.000002, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
+%! 0.000002, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
+%! 0.000049, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
+%! 0.001606, 0.000037, 0.000000, 0.000000, 0.000000, 0.000000;
+%! 0.024285, 0.001420, 0.000032, 0.000000, 0.000000, 0.000000;
+%! 0.161412, 0.022812, 0.001319, 0.000030, 0.000000, 0.000000;
+%! 0.499869, 0.156458, 0.021893, 0.001256, 0.000028, 0.000000;
+%! 0.839108, 0.493229, 0.153110, 0.021264, 0.001212, 0.000027;
+%! 0.976358, 0.835657, 0.488497, 0.150693, 0.020806, 0.001180;
+%! 0.998549, 0.975673, 0.833104, 0.484953, 0.148867, 0.020458;
+%! 0.999965, 0.998498, 0.975152, 0.831138, 0.482198, 0.147437;
+%! 1.000000, 0.999963, 0.998458, 0.974742, 0.829576, 0.479995;
+%! 1.000000, 1.000000, 0.999962, 0.998426, 0.974411, 0.828307;
+%! 1.000000, 1.000000, 1.000000, 0.999961, 0.998400, 0.974138;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 0.999960, 0.998378;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999960;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000];
+%! q = marcumq(a,b);
+%! assert(q,Q,1e-6);
+
+
+%!test
+%! a = [0.00;0.05;1.00;2.00;3.00;4.00;5.00;6.00;7.00;8.00;9.00;10.00;
+%! 11.00;12.00;13.00;14.00;15.00;16.00;17.00;18.00;19.00;20.00;
+%! 21.00;22.00;23.00;24.00];
+%! b = [11.10000, 12.10000, 13.10000, 14.10000, 15.10000, 16.10000];
+%! Q = [0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
+%! 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
+%! 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
+%! 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
+%! 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
+%! 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
+%! 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
+%! 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
+%! 0.000026, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
+%! 0.001155, 0.000026, 0.000000, 0.000000, 0.000000, 0.000000;
+%! 0.020183, 0.001136, 0.000025, 0.000000, 0.000000, 0.000000;
+%! 0.146287, 0.019961, 0.001120, 0.000025, 0.000000, 0.000000;
+%! 0.478193, 0.145342, 0.019778, 0.001107, 0.000024, 0.000000;
+%! 0.827253, 0.476692, 0.144551, 0.019625, 0.001096, 0.000024;
+%! 0.973909, 0.826366, 0.475422, 0.143881, 0.019494, 0.001087;
+%! 0.998359, 0.973714, 0.825607, 0.474333, 0.143304, 0.019381;
+%! 0.999959, 0.998343, 0.973546, 0.824952, 0.473389, 0.142803;
+%! 1.000000, 0.999959, 0.998330, 0.973400, 0.824380, 0.472564;
+%! 1.000000, 1.000000, 0.999958, 0.998318, 0.973271, 0.823876;
+%! 1.000000, 1.000000, 1.000000, 0.999958, 0.998307, 0.973158;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 0.999957, 0.998297;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999957;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000];
+%! q = marcumq(a,b);
+%! assert(q,Q,1e-6);
+
+
+%!test
+%! a = [0.00;0.05;1.00;2.00;3.00;4.00;5.00;6.00;7.00;8.00;9.00;10.00;
+%! 11.00;12.00;13.00;14.00;15.00;16.00;17.00;18.00;19.00;20.00;
+%! 21.00;22.00;23.00;24.00];
+%! b = [17.10000, 18.10000, 19.10000];
+%! Q = [0.000000, 0.000000, 0.000000;
+%! 0.000000, 0.000000, 0.000000;
+%! 0.000000, 0.000000, 0.000000;
+%! 0.000000, 0.000000, 0.000000;
+%! 0.000000, 0.000000, 0.000000;
+%! 0.000000, 0.000000, 0.000000;
+%! 0.000000, 0.000000, 0.000000;
+%! 0.000000, 0.000000, 0.000000;
+%! 0.000000, 0.000000, 0.000000;
+%! 0.000000, 0.000000, 0.000000;
+%! 0.000000, 0.000000, 0.000000;
+%! 0.000000, 0.000000, 0.000000;
+%! 0.000000, 0.000000, 0.000000;
+%! 0.000000, 0.000000, 0.000000;
+%! 0.000024, 0.000000, 0.000000;
+%! 0.001078, 0.000024, 0.000000;
+%! 0.019283, 0.001071, 0.000023;
+%! 0.142364, 0.019197, 0.001065;
+%! 0.471835, 0.141976, 0.019121;
+%! 0.823429, 0.471188, 0.141630;
+%! 0.973056, 0.823030, 0.470608;
+%! 0.998289, 0.972965, 0.822671;
+%! 0.999957, 0.998281, 0.972883;
+%! 1.000000, 0.999957, 0.998274;
+%! 1.000000, 1.000000, 0.999956;
+%! 1.000000, 1.000000, 1.000000];
+%! q = marcumq(a,b);
+%! assert(q,Q,1e-6);
+
+% The tests for M>1 were generating from Marcum's tables by
+% using the formula
+% Q_M(a,b) = Q(a,b) + exp(-(a-b)^2/2)*sum_{k=1}^{M-1}(b/a)^k*exp(-ab)*I_k(ab)
+%
+%!test
+%! M = 2;
+%! a = [0.00;0.05;1.00;2.00;3.00;4.00;5.00;6.00;7.00;8.00;9.00;10.00;11.00;12.00;13.00;14.00;15.00;16.00;17.00;18.00;19.00;20.00;21.00;22.00;23.00;24.000000];
+%!
+%! b = [ 0.00, 0.10, 2.10, 7.10, 12.10, 17.10];
+%! Q = [1.000000, 0.999987, 0.353353, 0.000000, 0.000000, 0.000000;
+%! 1.000000, 0.999988, 0.353687, 0.000000, 0.000000, 0.000000;
+%! 1.000000, 0.999992, 0.478229, 0.000000, 0.000000, 0.000000;
+%! 1.000000, 0.999999, 0.745094, 0.000001, 0.000000, 0.000000;
+%! 1.000000, 1.000000, 0.934771, 0.000077, 0.000000, 0.000000;
+%! 1.000000, 1.000000, 0.992266, 0.002393, 0.000000, 0.000000;
+%! 1.000000, 1.000000, 0.999607, 0.032264, 0.000000, 0.000000;
+%! 1.000000, 1.000000, 0.999992, 0.192257, 0.000000, 0.000000;
+%! 1.000000, 1.000000, 1.000000, 0.545174, 0.000000, 0.000000;
+%! 1.000000, 1.000000, 1.000000, 0.864230, 0.000040, 0.000000;
+%! 1.000000, 1.000000, 1.000000, 0.981589, 0.001555, 0.000000;
+%! 1.000000, 1.000000, 1.000000, 0.998957, 0.024784, 0.000000;
+%! 1.000000, 1.000000, 1.000000, 0.999976, 0.166055, 0.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 0.509823, 0.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 0.846066, 0.000032;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 0.978062, 0.001335;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 0.998699, 0.022409;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 0.999970, 0.156421;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.495223;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.837820;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.976328;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.998564;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999966;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000];
+%! q = marcumq(a,b,M);
+%! assert(q,Q,1e-6);
+
+%!test
+%! M = 5;
+%! a = [0.00;0.05;1.00;2.00;3.00;4.00;5.00;6.00;7.00;8.00;9.00;10.00;11.00;12.00;13.00;14.00;15.00;16.00;17.00;18.00;19.00;20.00;21.00;22.00;23.00;24.000000];
+%!
+%! b = [ 0.00, 0.10, 2.10, 7.10, 12.10, 17.10];
+%! Q = [1.000000, 1.000000, 0.926962, 0.000000, 0.000000, 0.000000;
+%! 1.000000, 1.000000, 0.927021, 0.000000, 0.000000, 0.000000;
+%! 1.000000, 1.000000, 0.947475, 0.000001, 0.000000, 0.000000;
+%! 1.000000, 1.000000, 0.980857, 0.000033, 0.000000, 0.000000;
+%! 1.000000, 1.000000, 0.996633, 0.000800, 0.000000, 0.000000;
+%! 1.000000, 1.000000, 0.999729, 0.011720, 0.000000, 0.000000;
+%! 1.000000, 1.000000, 0.999990, 0.088999, 0.000000, 0.000000;
+%! 1.000000, 1.000000, 1.000000, 0.341096, 0.000000, 0.000000;
+%! 1.000000, 1.000000, 1.000000, 0.705475, 0.000002, 0.000000;
+%! 1.000000, 1.000000, 1.000000, 0.933009, 0.000134, 0.000000;
+%! 1.000000, 1.000000, 1.000000, 0.993118, 0.003793, 0.000000;
+%! 1.000000, 1.000000, 1.000000, 0.999702, 0.045408, 0.000000;
+%! 1.000000, 1.000000, 1.000000, 0.999995, 0.238953, 0.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 0.607903, 0.000001;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 0.896007, 0.000073;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 0.987642, 0.002480;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 0.999389, 0.034450;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 0.999988, 0.203879;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.565165;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.876284;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.984209;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999165;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999983;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000];
+%! q = marcumq(a,b,M);
+%! assert(q,Q,1e-6);
+
+%!test
+%! M = 10;
+%! a = [0.00;0.05;1.00;2.00;3.00;4.00;5.00;6.00;7.00;8.00;9.00;10.00;11.00;12.00;13.00;14.00;15.00;16.00;17.00;18.00;19.00;20.00;21.00;22.00;23.00;24.000000];
+%!
+%! b = [ 0.00, 0.10, 2.10, 7.10, 12.10, 17.10];
+%! Q = [1.000000, 1.000000, 0.999898, 0.000193, 0.000000, 0.000000;
+%! 1.000000, 1.000000, 0.999897, 0.000194, 0.000000, 0.000000;
+%! 1.000000, 1.000000, 0.999931, 0.000416, 0.000000, 0.000000;
+%! 1.000000, 1.000000, 0.999980, 0.002377, 0.000000, 0.000000;
+%! 1.000000, 1.000000, 0.999997, 0.016409, 0.000000, 0.000000;
+%! 1.000000, 1.000000, 0.999999, 0.088005, 0.000000, 0.000000;
+%! 1.000000, 1.000000, 1.000000, 0.302521, 0.000000, 0.000000;
+%! 1.000000, 1.000000, 1.000000, 0.638401, 0.000000, 0.000000;
+%! 1.000000, 1.000000, 1.000000, 0.894322, 0.000022, 0.000000;
+%! 1.000000, 1.000000, 1.000000, 0.984732, 0.000840, 0.000000;
+%! 1.000000, 1.000000, 1.000000, 0.998997, 0.014160, 0.000000;
+%! 1.000000, 1.000000, 1.000000, 0.999972, 0.107999, 0.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 0.391181, 0.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 0.754631, 0.000004;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 0.951354, 0.000266;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 0.995732, 0.006444;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 0.999843, 0.065902;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 0.999998, 0.299616;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.676336;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.925312;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.992390;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999679;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999995;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
+%! 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000];
+%! q = marcumq(a,b,M);
+%! assert(q,Q,1e-6);
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <jpi...@us...> - 2012-05-04 09:10:03
|
Revision: 10359
http://octave.svn.sourceforge.net/octave/?rev=10359&view=rev
Author: jpicarbajal
Date: 2012-05-04 09:09:56 +0000 (Fri, 04 May 2012)
Log Message:
-----------
geometry: polygons2d
Modified Paths:
--------------
trunk/octave-forge/main/geometry/INDEX
trunk/octave-forge/main/geometry/NEWS
Modified: trunk/octave-forge/main/geometry/INDEX
===================================================================
--- trunk/octave-forge/main/geometry/INDEX 2012-05-04 09:06:45 UTC (rev 10358)
+++ trunk/octave-forge/main/geometry/INDEX 2012-05-04 09:09:56 UTC (rev 10359)
@@ -117,13 +117,22 @@
cbezier2poly
2D Polygons
curvature
+ distancePolygons
+ distancePointPolygon
+ distancePointPolyline
drawPolygon
+ expandPolygon
+ medialAxisConvex
oc_polybool
parametrize
+ polygonLoops
+ polygonSelfIntersections
+ polylineSelfIntersections
reversePolyline
reversePolygon
simplifypolygon
simplifypolyline
+ splitPolygons
supportFunction
2D Piecewise polynomial shapes
polygon2shape
Modified: trunk/octave-forge/main/geometry/NEWS
===================================================================
--- trunk/octave-forge/main/geometry/NEWS 2012-05-04 09:06:45 UTC (rev 10358)
+++ trunk/octave-forge/main/geometry/NEWS 2012-05-04 09:09:56 UTC (rev 10359)
@@ -26,6 +26,10 @@
- supportFunction.m: Compute support function of a polygon.
+ - distancePointPolygon.m, distancePointPolyline.m, distancePolygons.m ,
+ expandPolygon.m, medialAxisConvex.m, polygonLoops.m, polygonSelfIntersections.m
+ polylineSelfIntersections.m, splitPolygons.m
+
* Changed functions:
- distancePointEdge.m: Now the function computes the distance between all points
and all edges. A third optional argument provides
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <jpi...@us...> - 2012-05-04 09:06:57
|
Revision: 10358
http://octave.svn.sourceforge.net/octave/?rev=10358&view=rev
Author: jpicarbajal
Date: 2012-05-04 09:06:45 +0000 (Fri, 04 May 2012)
Log Message:
-----------
geometry: polygons2d
Added Paths:
-----------
trunk/octave-forge/main/geometry/inst/polygons2d/distancePointPolygon.m
trunk/octave-forge/main/geometry/inst/polygons2d/distancePointPolyline.m
trunk/octave-forge/main/geometry/inst/polygons2d/distancePolygons.m
trunk/octave-forge/main/geometry/inst/polygons2d/expandPolygon.m
trunk/octave-forge/main/geometry/inst/polygons2d/medialAxisConvex.m
trunk/octave-forge/main/geometry/inst/polygons2d/polygonLoops.m
trunk/octave-forge/main/geometry/inst/polygons2d/polygonSelfIntersections.m
trunk/octave-forge/main/geometry/inst/polygons2d/polylineSelfIntersections.m
trunk/octave-forge/main/geometry/inst/polygons2d/splitPolygons.m
Added: trunk/octave-forge/main/geometry/inst/polygons2d/distancePointPolygon.m
===================================================================
--- trunk/octave-forge/main/geometry/inst/polygons2d/distancePointPolygon.m (rev 0)
+++ trunk/octave-forge/main/geometry/inst/polygons2d/distancePointPolygon.m 2012-05-04 09:06:45 UTC (rev 10358)
@@ -0,0 +1,51 @@
+%% Copyright (C) 2003-2011 David Legland <dav...@gr...>
+%% Copyright (C) 2012 Adapted to Octave by Juan Pablo Carbajal <car...@if...>
+%% All rights reserved.
+%%
+%% Redistribution and use in source and binary forms, with or without
+%% modification, are permitted provided that the following conditions are met:
+%%
+%% 1 Redistributions of source code must retain the above copyright notice,
+%% this list of conditions and the following disclaimer.
+%% 2 Redistributions in binary form must reproduce the above copyright
+%% notice, this list of conditions and the following disclaimer in the
+%% documentation and/or other materials provided with the distribution.
+%%
+%% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS''
+%% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+%% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+%% ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR
+%% ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+%% DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+%% SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+%% CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+%% OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+%% OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+%%
+%% The views and conclusions contained in the software and documentation are
+%% those of the authors and should not be interpreted as representing official
+%% policies, either expressed or implied, of the copyright holders.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{dist} = } distancePointPolygon (@var{point},@var{poly})
+## Compute shortest distance between a point and a polygon
+##
+## @seealso{polygons2d, points2d, distancePointPolyline, distancePointEdge, projPointOnPolyline}
+## @end deftypefn
+
+function varargout = distancePointPolygon(point, poly)
+
+ % eventually copy first point at the end to ensure closed polygon
+ if sum(poly(end, :)==poly(1,:))~=2
+ poly = [poly; poly(1,:)];
+ end
+
+ % call to distancePointPolyline
+ minDist = distancePointPolyline(point, poly);
+
+ % process output arguments
+ if nargout<=1
+ varargout{1} = minDist;
+ end
+
+endfunction
Property changes on: trunk/octave-forge/main/geometry/inst/polygons2d/distancePointPolygon.m
___________________________________________________________________
Added: svn:executable
+ *
Added: trunk/octave-forge/main/geometry/inst/polygons2d/distancePointPolyline.m
===================================================================
--- trunk/octave-forge/main/geometry/inst/polygons2d/distancePointPolyline.m (rev 0)
+++ trunk/octave-forge/main/geometry/inst/polygons2d/distancePointPolyline.m 2012-05-04 09:06:45 UTC (rev 10358)
@@ -0,0 +1,71 @@
+## Copyright (C) 2003-2011 David Legland <dav...@gr...>
+## Copyright (C) 2012 Adapted to Octave by Juan Pablo Carbajal <car...@if...>
+## All rights reserved.
+##
+## Redistribution and use in source and binary forms, with or without
+## modification, are permitted provided that the following conditions are met:
+##
+## 1 Redistributions of source code must retain the above copyright notice,
+## this list of conditions and the following disclaimer.
+## 2 Redistributions in binary form must reproduce the above copyright
+## notice, this list of conditions and the following disclaimer in the
+## documentation and/or other materials provided with the distribution.
+##
+## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS''
+## AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+## IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+## ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR
+## ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+## DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+## SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+## CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+##
+## The views and conclusions contained in the software and documentation are
+## those of the authors and should not be interpreted as representing official
+## policies, either expressed or implied, of the copyright holders.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{dist} = } distancePointPolyline (@var{point},@var{poly})
+## Compute shortest distance between a point and a polyline
+## Example:
+## @example
+## pt1 = [30 20];
+## pt2 = [30 5];
+## poly = [10 10;50 10;50 50;10 50];
+## distancePointPolyline([pt1;pt2], poly)
+## ans =
+## 10
+## 5
+## @end example
+##
+## @seealso{polygons2d, points2d,distancePointEdge, projPointOnPolyline}
+## @end deftypefn
+
+function varargout = distancePointPolyline(point, poly, varargin)
+
+ # number of points
+ Np = size(point, 1);
+
+ # allocate memory for result
+ minDist = inf * ones(Np, 1);
+
+ # process each point
+ for p = 1:Np
+ # construct the set of edges
+ edges = [poly(1:end-1, :) poly(2:end, :)];
+
+ # compute distance between current each point and all edges
+ dist = distancePointEdge(point(p, :), edges);
+
+ # update distance if necessary
+ minDist(p) = min(dist);
+ end
+
+ # process output arguments
+ if nargout<=1
+ varargout{1} = minDist;
+ end
+
+endfunction
Property changes on: trunk/octave-forge/main/geometry/inst/polygons2d/distancePointPolyline.m
___________________________________________________________________
Added: svn:executable
+ *
Added: trunk/octave-forge/main/geometry/inst/polygons2d/distancePolygons.m
===================================================================
--- trunk/octave-forge/main/geometry/inst/polygons2d/distancePolygons.m (rev 0)
+++ trunk/octave-forge/main/geometry/inst/polygons2d/distancePolygons.m 2012-05-04 09:06:45 UTC (rev 10358)
@@ -0,0 +1,43 @@
+## Copyright (C) 2003-2011 David Legland <dav...@gr...>
+## Copyright (C) 2012 Adapted to Octave by Juan Pablo Carbajal <car...@if...>
+## All rights reserved.
+##
+## Redistribution and use in source and binary forms, with or without
+## modification, are permitted provided that the following conditions are met:
+##
+## 1 Redistributions of source code must retain the above copyright notice,
+## this list of conditions and the following disclaimer.
+## 2 Redistributions in binary form must reproduce the above copyright
+## notice, this list of conditions and the following disclaimer in the
+## documentation and/or other materials provided with the distribution.
+##
+## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS''
+## AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+## IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+## ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR
+## ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+## DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+## SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+## CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+##
+## The views and conclusions contained in the software and documentation are
+## those of the authors and should not be interpreted as representing official
+## policies, either expressed or implied, of the copyright holders.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{dist} = } distancePolygons (@var{poly1},@var{poly2})
+## Compute the shortest distance between 2 polygons
+##
+## @end deftypefn
+function dist = distancePolygons(poly1, poly2)
+
+ # compute distance of each vertex of a polygon to the other polygon
+ dist1 = min(distancePointPolygon(poly1, poly2));
+ dist2 = min(distancePointPolygon(poly2, poly1));
+
+ # keep the minimum of the two distances
+ dist = min(dist1, dist2);
+
+endfunction
Property changes on: trunk/octave-forge/main/geometry/inst/polygons2d/distancePolygons.m
___________________________________________________________________
Added: svn:executable
+ *
Added: trunk/octave-forge/main/geometry/inst/polygons2d/expandPolygon.m
===================================================================
--- trunk/octave-forge/main/geometry/inst/polygons2d/expandPolygon.m (rev 0)
+++ trunk/octave-forge/main/geometry/inst/polygons2d/expandPolygon.m 2012-05-04 09:06:45 UTC (rev 10358)
@@ -0,0 +1,85 @@
+## Copyright (C) 2003-2011 David Legland <dav...@gr...>
+## Copyright (C) 2012 Adapted to Octave by Juan Pablo Carbajal <car...@if...>
+## All rights reserved.
+##
+## Redistribution and use in source and binary forms, with or without
+## modification, are permitted provided that the following conditions are met:
+##
+## 1 Redistributions of source code must retain the above copyright notice,
+## this list of conditions and the following disclaimer.
+## 2 Redistributions in binary form must reproduce the above copyright
+## notice, this list of conditions and the following disclaimer in the
+## documentation and/or other materials provided with the distribution.
+##
+## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS''
+## AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+## IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+## ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR
+## ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+## DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+## SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+## CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+##
+## The views and conclusions contained in the software and documentation are
+## those of the authors and should not be interpreted as representing official
+## policies, either expressed or implied, of the copyright holders.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{loops} = } expandPolygon (@var{poly}, @var{dist})
+## Expand a polygon by a given (signed) distance
+##
+## Associates to each edge of the polygon @var{poly} the parallel line located
+## at distance @var{dist} from the current edge, and compute intersections with
+## neighbor parallel lines. The resulting polygon is simplified to remove
+## inner "loops", and can be disconnected.
+## The result is a cell array, each cell containing a simple linear ring.
+##
+## This is a kind of dilation, but behaviour on corners is different.
+## This function keeps angles of polygons, but there is no direct relation
+## between length of 2 polygons.
+##
+## It is also possible to specify negative distance, and get all points
+## inside the polygon. If the polygon is convex, the result equals
+## morphological erosion of polygon by a ball with radius equal to the
+## given distance.
+##
+## @seealso{polygons2d}
+## @end deftypefn
+
+function loops = expandPolygon(poly, dist)
+
+ # eventually copy first point at the end to ensure closed polygon
+ if sum(poly(end, :)==poly(1,:))~=2
+ poly = [poly; poly(1,:)];
+ end
+
+ # number of vertices of the polygon
+ N = size(poly, 1)-1;
+
+ # find lines parallel to polygon edges located at distance DIST
+ lines = zeros(N, 4);
+ for i=1:N
+ side = createLine(poly(i,:), poly(i+1,:));
+ lines(i, 1:4) = parallelLine(side, dist);
+ end
+
+ # compute intersection points of consecutive lines
+ lines = [lines;lines(1,:)];
+ poly2 = zeros(N, 2);
+ for i=1:N
+ poly2(i,1:2) = intersectLines(lines(i,:), lines(i+1,:));
+ end
+
+ # split result polygon into set of loops (simple polygons)
+ loops = polygonLoops(poly2);
+
+ # keep only loops whose distance to original polygon is correct
+ distLoop = zeros(length(loops), 1);
+ for i=1:length(loops)
+ distLoop(i) = distancePolygons(loops{i}, poly);
+ end
+ loops = loops(abs(distLoop-abs(dist))<abs(dist)/1000);
+
+endfunction
Property changes on: trunk/octave-forge/main/geometry/inst/polygons2d/expandPolygon.m
___________________________________________________________________
Added: svn:executable
+ *
Added: trunk/octave-forge/main/geometry/inst/polygons2d/medialAxisConvex.m
===================================================================
--- trunk/octave-forge/main/geometry/inst/polygons2d/medialAxisConvex.m (rev 0)
+++ trunk/octave-forge/main/geometry/inst/polygons2d/medialAxisConvex.m 2012-05-04 09:06:45 UTC (rev 10358)
@@ -0,0 +1,153 @@
+## Copyright (C) 2003-2011 David Legland <dav...@gr...>
+## Copyright (C) 2012 Adapted to Octave by Juan Pablo Carbajal <car...@if...>
+## All rights reserved.
+##
+## Redistribution and use in source and binary forms, with or without
+## modification, are permitted provided that the following conditions are met:
+##
+## 1 Redistributions of source code must retain the above copyright notice,
+## this list of conditions and the following disclaimer.
+## 2 Redistributions in binary form must reproduce the above copyright
+## notice, this list of conditions and the following disclaimer in the
+## documentation and/or other materials provided with the distribution.
+##
+## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS''
+## AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+## IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+## ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR
+## ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+## DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+## SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+## CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+##
+## The views and conclusions contained in the software and documentation are
+## those of the authors and should not be interpreted as representing official
+## policies, either expressed or implied, of the copyright holders.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{n} @var{e}] = } medialAxisConvex (@var{polygon})
+## Compute medial axis of a convex polygon
+##
+## @var{polygon} is given as a set of points [x1 y1;x2 y2 ...], returns
+## the medial axis of the polygon as a graph.
+## @var{n} is a set of nodes. The first elements of @var{n} are the vertices of the
+## original polygon.
+## @var{e} is a set of edges, containing indices of source and target nodes.
+## Edges are sorted according to order of creation. Index of first vertex
+## is lower than index of last vertex, i.e. edges always point to newly
+## created nodes.
+##
+## Notes:
+## - Is not fully implemented, need more development (usually crashes for
+## polygons with more than 6-7 points...)
+## - Works only for convex polygons.
+## - Complexity is not optimal: this algorithm is O(n*log n), but linear
+## algorithms exist.
+##
+## @seealso{polygons2d, bisector}
+## @end deftypefn
+
+function [nodes, edges] = medialAxisConvex(points)
+
+ # eventually remove the last point if it is the same as the first one
+ if points(1,:) == points(end, :)
+ nodes = points(1:end-1, :);
+ else
+ nodes = points;
+ end
+
+ # special case of triangles:
+ # compute directly the gravity center, and simplify computation.
+ if size(nodes, 1)==3
+ nodes = [nodes; mean(nodes, 1)];
+ edges = [1 4;2 4;3 4];
+ return
+ end
+
+ # number of nodes, and also of initial rays
+ N = size(nodes, 1);
+
+ # create ray of each vertex
+ rays = zeros(N, 4);
+ rays(1, 1:4) = bisector(nodes([2 1 N], :));
+ rays(N, 1:4) = bisector(nodes([1 N N-1], :));
+ for i=2:N-1
+ rays(i, 1:4) = bisector(nodes([i+1, i, i-1], :));
+ end
+
+ # add indices of edges producing rays (indices of first vertex, second
+ # vertex is obtained by adding one modulo N).
+ rayEdges = [[N (1:N-1)]' (1:N)'];
+
+ pint = intersectLines(rays, rays([2:N 1], :));
+ #ti = linePosition(pint, rays);
+ #ti = min(linePosition(pint, rays), linePosition(pint, rays([2:N 1], :)));
+ ti = distancePointLine(pint, ...
+ createLine(points([N (1:N-1)]', :), points((1:N)', :)));
+
+ # create list of events.
+ # terms are : R1 R2 X Y t0
+ # R1 and R2 are indices of involved rays
+ # X and Y is coordinate of intersection point
+ # t0 is position of point on rays
+ events = sortrows([ (1:N)' [2:N 1]' pint ti], 5);
+
+ # initialize edges
+ edges = zeros(0, 2);
+
+
+ # -------------------
+ # process each event until there is no more
+
+ # start after index of last vertex, and process N-3 intermediate rays
+ for i=N+1:2*N-3
+ # add new node at the rays intersection
+ nodes(i,:) = events(1, 3:4);
+
+ # add new couple of edges
+ edges = [edges; events(1,1) i; events(1,2) i];
+
+ # find the two edges creating the new emanating ray
+ n1 = rayEdges(events(1, 1), 1);
+ n2 = rayEdges(events(1, 2), 2);
+
+ # create the new ray
+ line1 = createLine(nodes(n1, :), nodes(mod(n1,N)+1, :));
+ line2 = createLine(nodes(mod(n2,N)+1, :), nodes(n2, :));
+ ray0 = bisector(line1, line2);
+
+ # set its origin to emanating point
+ ray0(1:2) = nodes(i, :);
+
+ # add the new ray to the list
+ rays = [rays; ray0];
+ rayEdges(size(rayEdges, 1)+1, 1:2) = [n1 n2];
+
+ # find the two neighbour rays
+ ind = sum(ismember(events(:,1:2), events(1, 1:2)), 2)==0;
+ ir = unique(events(ind, 1:2));
+ ir = ir(~ismember(ir, events(1,1:2)));
+
+ # create new intersections
+ pint = intersectLines(ray0, rays(ir, :));
+ #ti = min(linePosition(pint, ray0), linePosition(pint, rays(ir, :))) + events(1,5);
+ ti = distancePointLine(pint, line1);
+
+ # remove all events involving old intersected rays
+ ind = sum(ismember(events(:,1:2), events(1, 1:2)), 2)==0;
+ events = events(ind, :);
+
+ # add the newly formed events
+ events = [events; ir(1) i pint(1,:) ti(1); ir(2) i pint(2,:) ti(2)];
+
+ # and sort them according to 'position' parameter
+ events = sortrows(events, 5);
+ end
+
+ # centroid computation for last 3 rays
+ nodes = [nodes; mean(events(:, 3:4))];
+ edges = [edges; [unique(events(:,1:2)) ones(3, 1)*(2*N-2)]];
+
+endfunction
Property changes on: trunk/octave-forge/main/geometry/inst/polygons2d/medialAxisConvex.m
___________________________________________________________________
Added: svn:executable
+ *
Added: trunk/octave-forge/main/geometry/inst/polygons2d/polygonLoops.m
===================================================================
--- trunk/octave-forge/main/geometry/inst/polygons2d/polygonLoops.m (rev 0)
+++ trunk/octave-forge/main/geometry/inst/polygons2d/polygonLoops.m 2012-05-04 09:06:45 UTC (rev 10358)
@@ -0,0 +1,143 @@
+## Copyright (C) 2003-2011 David Legland <dav...@gr...>
+## Copyright (C) 2012 Adapted to Octave by Juan Pablo Carbajal <car...@if...>
+## All rights reserved.
+##
+## Redistribution and use in source and binary forms, with or without
+## modification, are permitted provided that the following conditions are met:
+##
+## 1 Redistributions of source code must retain the above copyright notice,
+## this list of conditions and the following disclaimer.
+## 2 Redistributions in binary form must reproduce the above copyright
+## notice, this list of conditions and the following disclaimer in the
+## documentation and/or other materials provided with the distribution.
+##
+## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS''
+## AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+## IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+## ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR
+## ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+## DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+## SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+## CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+##
+## The views and conclusions contained in the software and documentation are
+## those of the authors and should not be interpreted as representing official
+## policies, either expressed or implied, of the copyright holders.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{loops} = } polygonLoops (@var{poly})
+## Divide a possibly self-intersecting polygon into a set of simple loops
+##
+## @var{poly} is a polygone defined by a series of vertices,
+## @var{loops} is a cell array of polygons, containing the same vertices of the
+## original polygon, but no loop self-intersect, and no couple of loops
+## intersect each other.
+##
+## Example:
+## @example
+## poly = [0 0;0 10;20 10;20 20;10 20;10 0];
+## loops = polygonLoops(poly);
+## figure(1); hold on;
+## drawPolygon(loops);
+## polygonArea(loops)
+## @end example
+##
+## @seealso{polygons2d, polygonSelfIntersections}
+## @end deftypefn
+
+function loops = polygonLoops(poly)
+
+ ## Initialisations
+
+ # compute intersections
+ [inters pos1 pos2] = polygonSelfIntersections(poly);
+
+ # case of a polygon without self-intersection
+ if isempty(inters)
+ loops = {poly};
+ return;
+ end
+
+ # array for storing loops
+ loops = cell(0, 1);
+
+ positions = sortrows([pos1 pos2;pos2 pos1]);
+
+
+ ## First loop
+
+ pos0 = 0;
+ loop = polygonSubcurve(poly, pos0, positions(1, 1));
+ pos = positions(1, 2);
+ positions(1, :) = [];
+
+ while true
+ # index of next intersection point
+ ind = find(positions(:,1)>pos, 1, 'first');
+
+ # if not found, break
+ if isempty(ind)
+ break;
+ end
+
+ # add portion of curve
+ loop = [loop;polygonSubcurve(poly, pos, positions(ind, 1))]; ##ok<AGROW>
+
+ # look for next intersection point
+ pos = positions(ind, 2);
+ positions(ind, :) = [];
+ end
+
+ # add the last portion of curve
+ loop = [loop;polygonSubcurve(poly, pos, pos0)];
+
+ # remove redundant vertices
+ loop(sum(loop(1:end-1,:) == loop(2:end,:) ,2)==2, :) = [];
+ if sum(diff(loop([1 end], :))==0)==2
+ loop(end, :) = [];
+ end
+
+ # add current loop to the list of loops
+ loops{1} = loop;
+
+
+ ## Other loops
+
+ Nl = 1;
+ while ~isempty(positions)
+
+ loop = [];
+ pos0 = positions(1, 2);
+ pos = positions(1, 2);
+
+ while true
+ # index of next intersection point
+ ind = find(positions(:,1)>pos, 1, 'first');
+
+ # add portion of curve
+ loop = [loop;polygonSubcurve(poly, pos, positions(ind, 1))]; ##ok<AGROW>
+
+ # look for next intersection point
+ pos = positions(ind, 2);
+ positions(ind, :) = [];
+
+ # if not found, break
+ if pos==pos0
+ break;
+ end
+ end
+
+ # remove redundant vertices
+ loop(sum(loop(1:end-1,:) == loop(2:end,:) ,2)==2, :) = []; ##ok<AGROW>
+ if sum(diff(loop([1 end], :))==0)==2
+ loop(end, :) = []; ##ok<AGROW>
+ end
+
+ # add current loop to the list of loops
+ Nl = Nl + 1;
+ loops{Nl} = loop;
+ end
+
+endfunction
Property changes on: trunk/octave-forge/main/geometry/inst/polygons2d/polygonLoops.m
___________________________________________________________________
Added: svn:executable
+ *
Added: trunk/octave-forge/main/geometry/inst/polygons2d/polygonSelfIntersections.m
===================================================================
--- trunk/octave-forge/main/geometry/inst/polygons2d/polygonSelfIntersections.m (rev 0)
+++ trunk/octave-forge/main/geometry/inst/polygons2d/polygonSelfIntersections.m 2012-05-04 09:06:45 UTC (rev 10358)
@@ -0,0 +1,87 @@
+## Copyright (C) 2003-2011 David Legland <dav...@gr...>
+## Copyright (C) 2012 Adapted to Octave by Juan Pablo Carbajal <car...@if...>
+## All rights reserved.
+##
+## Redistribution and use in source and binary forms, with or without
+## modification, are permitted provided that the following conditions are met:
+##
+## 1 Redistributions of source code must retain the above copyright notice,
+## this list of conditions and the following disclaimer.
+## 2 Redistributions in binary form must reproduce the above copyright
+## notice, this list of conditions and the following disclaimer in the
+## documentation and/or other materials provided with the distribution.
+##
+## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS''
+## AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+## IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+## ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR
+## ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+## DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+## SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+## CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+##
+## The views and conclusions contained in the software and documentation are
+## those of the authors and should not be interpreted as representing official
+## policies, either expressed or implied, of the copyright holders.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{pts} = } polygonSelfIntersections (@var{poly})
+## @deftypefnx {Function File} {[@var{pts} @var{pos1} @var{pos2}] = } polygonSelfIntersections (@var{poly})
+## Find-self intersection points of a polygon
+##
+## Return the position of self intersection points
+##
+## Also return the 2 positions of each intersection point (the position
+## when meeting point for first time, then position when meeting point
+## for the second time).
+##
+## Example
+## @example
+## # use a '8'-shaped polygon
+## poly = [10 0;0 0;0 10;20 10;20 20;10 20];
+## polygonSelfIntersections(poly)
+## ans =
+## 10 10
+## @end example
+##
+## @seealso{polygons2d, polylineSelfIntersections}
+## @end deftypefn
+
+function varargout = polygonSelfIntersections(poly, varargin)
+
+ tol = 1e-14;
+
+ # ensure the last point equals the first one
+ if sum(abs(poly(end, :)-poly(1,:)) < tol) ~= 2
+ poly = [poly; poly(1,:)];
+ end
+
+ # compute intersections by calling algo for polylines
+ [points pos1 pos2] = polylineSelfIntersections(poly);
+
+ # It may append that first vertex of polygon is detected as intersection,
+ # the following tries to detect this
+ n = size(poly, 1) - 1;
+ inds = (pos1 == 0 & pos2 == n) | (pos1 == n & pos2 == 0);
+ points(inds, :) = [];
+ pos1(inds) = [];
+ pos2(inds) = [];
+
+ # remove multiple intersections
+ [points I J] = unique(points, 'rows', 'first'); ##ok<NASGU>
+ pos1 = pos1(I);
+ pos2 = pos2(I);
+
+
+ ## Post-processing
+
+ # process output arguments
+ if nargout <= 1
+ varargout = {points};
+ elseif nargout == 3
+ varargout = {points, pos1, pos2};
+ end
+
+endfunction
Property changes on: trunk/octave-forge/main/geometry/inst/polygons2d/polygonSelfIntersections.m
___________________________________________________________________
Added: svn:executable
+ *
Added: trunk/octave-forge/main/geometry/inst/polygons2d/polylineSelfIntersections.m
===================================================================
--- trunk/octave-forge/main/geometry/inst/polygons2d/polylineSelfIntersections.m (rev 0)
+++ trunk/octave-forge/main/geometry/inst/polygons2d/polylineSelfIntersections.m 2012-05-04 09:06:45 UTC (rev 10358)
@@ -0,0 +1,153 @@
+## Copyright (C) 2003-2011 David Legland <dav...@gr...>
+## Copyright (C) 2012 Adapted to Octave by Juan Pablo Carbajal <car...@if...>
+## All rights reserved.
+##
+## Redistribution and use in source and binary forms, with or without
+## modification, are permitted provided that the following conditions are met:
+##
+## 1 Redistributions of source code must retain the above copyright notice,
+## this list of conditions and the following disclaimer.
+## 2 Redistributions in binary form must reproduce the above copyright
+## notice, this list of conditions and the following disclaimer in the
+## documentation and/or other materials provided with the distribution.
+##
+## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS''
+## AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+## IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+## ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR
+## ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+## DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+## SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+## CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+##
+## The views and conclusions contained in the software and documentation are
+## those of the authors and should not be interpreted as representing official
+## policies, either expressed or implied, of the copyright holders.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{pts} = } polylineSelfIntersections (@var{poly})
+## Find self-intersections points of a polyline
+##
+## Return the position of self intersection points
+##
+## Also return the 2 positions of each intersection point (the position
+## when meeting point for first time, then position when meeting point
+## for the second time).
+##
+## Example
+## @example
+## # use a gamma-shaped polyline
+## poly = [0 0;0 10;20 10;20 20;10 20;10 0];
+## polylineSelfIntersections(poly)
+## ans =
+## 10 10
+##
+## # use a 'S'-shaped polyline
+## poly = [10 0;0 0;0 10;20 10;20 20;10 20];
+## polylineSelfIntersections(poly)
+## ans =
+## 10 10
+## @end example
+##
+## @seealso{polygons2d, intersectPolylines, polygonSelfIntersections}
+## @end deftypefn
+
+function varargout = polylineSelfIntersections(poly, varargin)
+
+ ## Initialisations
+
+ # determine whether the polyline is closed
+ closed = false;
+ if ~isempty(varargin)
+ closed = varargin{1};
+ if ischar(closed)
+ if strcmp(closed, 'closed')
+ closed = true;
+ elseif strcmp(closed, 'open')
+ closed = false;
+ end
+ end
+ end
+
+ # if polyline is closed, ensure the last point equals the first one
+ if closed
+ if sum(abs(poly(end, :)-poly(1,:))<1e-14)~=2
+ poly = [poly; poly(1,:)];
+ end
+ end
+
+ # arrays for storing results
+ points = zeros(0, 2);
+ pos1 = zeros(0, 1);
+ pos2 = zeros(0, 1);
+
+ # number of vertices
+ Nv = size(poly, 1);
+
+
+ ## Main processing
+
+ # index of current intersection
+ ip = 0;
+
+ # iterate over each couple of edge ( (N-1)*(N-2)/2 iterations)
+ for i=1:Nv-2
+ # create first edge
+ edge1 = [poly(i, :) poly(i+1, :)];
+ for j=i+2:Nv-1
+ # create second edge
+ edge2 = [poly(j, :) poly(j+1, :)];
+
+ # check conditions on bounding boxes
+ if min(edge1([1 3]))>max(edge2([1 3]))
+ continue;
+ end
+ if max(edge1([1 3]))<min(edge2([1 3]))
+ continue;
+ end
+ if min(edge1([2 4]))>max(edge2([2 4]))
+ continue;
+ end
+ if max(edge1([2 4]))<min(edge2([2 4]))
+ continue;
+ end
+
+ # compute intersection point
+ inter = intersectEdges(edge1, edge2);
+
+ if sum(isfinite(inter))==2
+ # add point to the list
+ ip = ip + 1;
+ points(ip, :) = inter;
+
+ # also compute positions on the polyline
+ pos1(ip, 1) = i+edgePosition(inter, edge1)-1;
+ pos2(ip, 1) = j+edgePosition(inter, edge2)-1;
+ end
+ end
+ end
+
+ # if polyline is closed, the first vertex was found as an intersection, so
+ # we need to remove it
+ if closed
+ dist = distancePoints(points, poly(1,:));
+ [minDist ind] = min(dist); ##ok<ASGLU>
+ points(ind,:) = [];
+ pos1(ind) = [];
+ pos2(ind) = [];
+ end
+
+ ## Post-processing
+
+ # process output arguments
+ if nargout<=1
+ varargout{1} = points;
+ elseif nargout==3
+ varargout{1} = points;
+ varargout{2} = pos1;
+ varargout{3} = pos2;
+ end
+
+endfunction
Property changes on: trunk/octave-forge/main/geometry/inst/polygons2d/polylineSelfIntersections.m
___________________________________________________________________
Added: svn:executable
+ *
Added: trunk/octave-forge/main/geometry/inst/polygons2d/splitPolygons.m
===================================================================
--- trunk/octave-forge/main/geometry/inst/polygons2d/splitPolygons.m (rev 0)
+++ trunk/octave-forge/main/geometry/inst/polygons2d/splitPolygons.m 2012-05-04 09:06:45 UTC (rev 10358)
@@ -0,0 +1,64 @@
+## Copyright (C) 2003-2011 David Legland <dav...@gr...>
+## Copyright (C) 2012 Adapted to Octave by Juan Pablo Carbajal <car...@if...>
+## All rights reserved.
+##
+## Redistribution and use in source and binary forms, with or without
+## modification, are permitted provided that the following conditions are met:
+##
+## 1 Redistributions of source code must retain the above copyright notice,
+## this list of conditions and the following disclaimer.
+## 2 Redistributions in binary form must reproduce the above copyright
+## notice, this list of conditions and the following disclaimer in the
+## documentation and/or other materials provided with the distribution.
+##
+## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS''
+## AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+## IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+## ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR
+## ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+## DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+## SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+## CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+##
+## The views and conclusions contained in the software and documentation are
+## those of the authors and should not be interpreted as representing official
+## policies, either expressed or implied, of the copyright holders.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{polygons} = } splitPolygons (@var{polygon})
+## Convert a NaN separated polygon list to a cell array of polygons.
+##
+## @var{polygon} is a N-by-2 array of points, with possibly couples of NaN values.
+## The functions separates each component separated by NaN values, and
+## returns a cell array of polygons.
+##
+## @seealso{polygons2d}
+## @end deftypefn
+function polygons = splitPolygons(polygon)
+
+ if iscell(polygon)
+ # case of a cell array
+ polygons = polygon;
+
+ elseif sum(isnan(polygon(:)))==0
+ # single polygon -> no break
+ polygons = {polygon};
+
+ else
+ # find indices of NaN couples
+ inds = find(sum(isnan(polygon), 2)>0);
+
+ # number of polygons
+ N = length(inds)+1;
+ polygons = cell(N, 1);
+
+ # iterate over NaN-separated regions to create new polygon
+ inds = [0;inds;size(polygon, 1)+1];
+ for i=1:N
+ polygons{i} = polygon((inds(i)+1):(inds(i+1)-1), :);
+ end
+ end
+
+endfunction
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <cd...@us...> - 2012-05-03 21:27:03
|
Revision: 10355
http://octave.svn.sourceforge.net/octave/?rev=10355&view=rev
Author: cdf
Date: 2012-05-03 17:03:10 +0000 (Thu, 03 May 2012)
Log Message:
-----------
added block-jacobi iterative solver
Added Paths:
-----------
trunk/octave-forge/extra/secs3d/inst/Utilities/block_jacobi.m
Added: trunk/octave-forge/extra/secs3d/inst/Utilities/block_jacobi.m
===================================================================
--- trunk/octave-forge/extra/secs3d/inst/Utilities/block_jacobi.m (rev 0)
+++ trunk/octave-forge/extra/secs3d/inst/Utilities/block_jacobi.m 2012-05-03 17:03:10 UTC (rev 10355)
@@ -0,0 +1,58 @@
+## Copyright (C) 2012 Carlo de Falco
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING. If not, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+##
+## @deftypefn {Function File} @
+## {[@var{x}, @var{res}, @var{nit}]} = @
+## block_jacobi (@var{A}, @var{b}, @var{x0}, @var{nblocks}, @var{maxit}, @var{tol})
+## Solve the linear system @code{A * x = b} using the block-Jacobi iterative algorithm.
+## The algorithm is implemented as a preconditioned Richardson iteration with a block-diagonal
+## preconditioner.
+## @var{nblocks} is the number of diagonal blocks to be extracted from the matrix @var{A},
+## @var{maxit} the maximum number of iterations.
+## If requested, the output parameter @var{res} contains the value of the residual @code{b - A*x}
+## at the last iteration and @var{nit} the number of iterations that have been performed.
+## @end deftypefn
+
+function [x, res, it] = block_jacobi (A, b, x0, nblocks, maxit, tol)
+
+ m = rows (A);
+ if nblocks > m
+ error ("block_jacobi: cannot partition a %d size matrix in %d blocks", m, nblocks)
+ endif
+
+ lblocks = ceil (m / nblocks);
+ ib = 1; first(ib) = 1; last(ib) = lblocks;
+ while (last(ib) < m)
+ ++ib;
+ first(ib) = last(ib-1) + 1;
+ last(ib) = min ([m first(ib)+lblocks-1]);
+ endwhile
+ nblocks = ib;
+
+ x = x0;
+ for it = 1:maxit
+ pres = res = b - A*x;
+ if (norm (res, inf) < tol); break; endif
+ for ib = 1:nblocks
+ range = first(ib):last(ib);
+ pres(range) = diag (diag (A(range, range))) \ res(range);
+ endfor
+ x += pres;
+ endfor
+
+endfunction
\ No newline at end of file
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-05-03 20:16:25
|
Revision: 10357
http://octave.svn.sourceforge.net/octave/?rev=10357&view=rev
Author: paramaniac
Date: 2012-05-03 19:54:40 +0000 (Thu, 03 May 2012)
Log Message:
-----------
control: add ctranspose for LTI models
Modified Paths:
--------------
trunk/octave-forge/main/control/INDEX
Added Paths:
-----------
trunk/octave-forge/main/control/inst/@frd/__ctranspose__.m
trunk/octave-forge/main/control/inst/@lti/ctranspose.m
trunk/octave-forge/main/control/inst/@ss/__ctranspose__.m
trunk/octave-forge/main/control/inst/@tf/__ctranspose__.m
trunk/octave-forge/main/control/inst/@tfpoly/conj_ct.m
trunk/octave-forge/main/control/inst/@tfpoly/conj_dt.m
Modified: trunk/octave-forge/main/control/INDEX
===================================================================
--- trunk/octave-forge/main/control/INDEX 2012-05-03 17:31:31 UTC (rev 10356)
+++ trunk/octave-forge/main/control/INDEX 2012-05-03 19:54:40 UTC (rev 10357)
@@ -107,6 +107,7 @@
fwcfconred
spaconred
Overloaded Operators
+ @lti/ctranspose
@lti/horzcat
@lti/inv
@lti/minus
Added: trunk/octave-forge/main/control/inst/@frd/__ctranspose__.m
===================================================================
--- trunk/octave-forge/main/control/inst/@frd/__ctranspose__.m (rev 0)
+++ trunk/octave-forge/main/control/inst/@frd/__ctranspose__.m 2012-05-03 19:54:40 UTC (rev 10357)
@@ -0,0 +1,35 @@
+## Copyright (C) 2012 Lukas F. Reichlin
+##
+## This file is part of LTI Syncope.
+##
+## LTI Syncope is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## LTI Syncope is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## Conjugate transpose of FRD objects.
+
+## Author: Lukas Reichlin <luk...@gm...>
+## Created: May 2012
+## Version: 0.1
+
+function sys = __ctranspose__ (sys)
+
+ [p, m, l] = size (sys.H);
+
+ H = mat2cell (sys.H, p, m, ones (1, l))(:);
+
+ H = cellfun (@ctranspose, H, "uniformoutput", false);
+
+ sys.H = cat (3, H{:});
+
+endfunction
\ No newline at end of file
Added: trunk/octave-forge/main/control/inst/@lti/ctranspose.m
===================================================================
--- trunk/octave-forge/main/control/inst/@lti/ctranspose.m (rev 0)
+++ trunk/octave-forge/main/control/inst/@lti/ctranspose.m 2012-05-03 19:54:40 UTC (rev 10357)
@@ -0,0 +1,39 @@
+## Copyright (C) 2012 Lukas F. Reichlin
+##
+## This file is part of LTI Syncope.
+##
+## LTI Syncope is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## LTI Syncope is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## Conjugate transpose of LTI objects. Used by Octave for "sys'".
+
+## Author: Lukas Reichlin <luk...@gm...>
+## Created: May 2012
+## Version: 0.1
+
+function sys = ctranspose (sys)
+
+ if (nargin != 1) # prevent sys = ctranspose (sys1, sys2, sys3, ...)
+ error ("lti: ctranspose: this is an unary operator");
+ endif
+
+ [p, m] = size (sys);
+ ct = isct (sys);
+
+ sys = __ctranspose__ (sys, ct);
+
+ sys.inname = repmat ({""}, p, 1);
+ sys.outname = repmat ({""}, m, 1);
+
+endfunction
\ No newline at end of file
Added: trunk/octave-forge/main/control/inst/@ss/__ctranspose__.m
===================================================================
--- trunk/octave-forge/main/control/inst/@ss/__ctranspose__.m (rev 0)
+++ trunk/octave-forge/main/control/inst/@ss/__ctranspose__.m 2012-05-03 19:54:40 UTC (rev 10357)
@@ -0,0 +1,54 @@
+## Copyright (C) 2012 Lukas F. Reichlin
+##
+## This file is part of LTI Syncope.
+##
+## LTI Syncope is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## LTI Syncope is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## Conjugate transpose of SS models.
+
+## Author: Lukas Reichlin <luk...@gm...>
+## Created: May 2012
+## Version: 0.1
+
+function sys = __ctranspose__ (sys, ct)
+
+ a = sys.a;
+ b = sys.b;
+ c = sys.c;
+ d = sys.d;
+ e = sys.e;
+
+ if (ct) # continuous-time
+ sys.a = -a.';
+ sys.b = -c.';
+ sys.c = b.';
+ sys.d = d.';
+ sys.e = e.';
+ sys.stname = repmat ({""}, rows (a), 1);
+ else # discrete-time
+ [n, m] = size (b);
+ p = rows (c);
+ if (isempty (e))
+ e = eye (n);
+ endif
+ sys.a = blkdiag (e.', eye (p));
+ sys.b = [zeros(n, p); -eye(p)];
+ sys.c = [b.', zeros(m, p)];
+ sys.d = d.';
+ sys.e = [a.', c.'; zeros(p, n+p)];
+ sys.stname = repmat ({""}, n+p, 1);
+ endif
+
+endfunction
Added: trunk/octave-forge/main/control/inst/@tf/__ctranspose__.m
===================================================================
--- trunk/octave-forge/main/control/inst/@tf/__ctranspose__.m (rev 0)
+++ trunk/octave-forge/main/control/inst/@tf/__ctranspose__.m 2012-05-03 19:54:40 UTC (rev 10357)
@@ -0,0 +1,42 @@
+## Copyright (C) 2012 Lukas F. Reichlin
+##
+## This file is part of LTI Syncope.
+##
+## LTI Syncope is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## LTI Syncope is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## Conjugate transpose of TF models.
+
+## Author: Lukas Reichlin <luk...@gm...>
+## Created: May 2012
+## Version: 0.1
+
+function sys = __ctranspose__ (sys, ct)
+
+ num = sys.num;
+ den = sys.den;
+
+ if (ct) # continuous-time
+ num = cellfun (@conj_ct, num, "uniformoutput", false);
+ den = cellfun (@conj_ct, den, "uniformoutput", false);
+ else # discrete-time
+ num = cellfun (@conj_dt, num, "uniformoutput", false);
+ den = cellfun (@conj_dt, den, "uniformoutput", false);
+ ## TODO: shall I make "den" a monic polynomial?
+ endif
+
+ sys.num = num.';
+ sys.den = den.';
+
+endfunction
\ No newline at end of file
Added: trunk/octave-forge/main/control/inst/@tfpoly/conj_ct.m
===================================================================
--- trunk/octave-forge/main/control/inst/@tfpoly/conj_ct.m (rev 0)
+++ trunk/octave-forge/main/control/inst/@tfpoly/conj_ct.m 2012-05-03 19:54:40 UTC (rev 10357)
@@ -0,0 +1,30 @@
+## Copyright (C) 2012 Lukas F. Reichlin
+##
+## This file is part of LTI Syncope.
+##
+## LTI Syncope is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## LTI Syncope is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## Conjugate of continuous-time polynomial. Replace s by -s.
+## For internal use only.
+
+## Author: Lukas Reichlin <luk...@gm...>
+## Created: May 2012
+## Version: 0.1
+
+function p = conj_ct (p)
+
+ p.poly(2:2:end) = -p.poly(2:2:end);
+
+endfunction
Added: trunk/octave-forge/main/control/inst/@tfpoly/conj_dt.m
===================================================================
--- trunk/octave-forge/main/control/inst/@tfpoly/conj_dt.m (rev 0)
+++ trunk/octave-forge/main/control/inst/@tfpoly/conj_dt.m 2012-05-03 19:54:40 UTC (rev 10357)
@@ -0,0 +1,30 @@
+## Copyright (C) 2012 Lukas F. Reichlin
+##
+## This file is part of LTI Syncope.
+##
+## LTI Syncope is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## LTI Syncope is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## Conjugate of discrete-time polynomial. Replace z by 1/z.
+## For internal use only.
+
+## Author: Lukas Reichlin <luk...@gm...>
+## Created: May 2012
+## Version: 0.1
+
+function p = conj_dt (p)
+
+ p.poly = fliplr (p.poly);
+
+endfunction
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <cd...@us...> - 2012-05-03 20:14:37
|
Revision: 10356
http://octave.svn.sourceforge.net/octave/?rev=10356&view=rev
Author: cdf
Date: 2012-05-03 17:31:31 +0000 (Thu, 03 May 2012)
Log Message:
-----------
add test
Modified Paths:
--------------
trunk/octave-forge/extra/secs3d/inst/Utilities/block_jacobi.m
Modified: trunk/octave-forge/extra/secs3d/inst/Utilities/block_jacobi.m
===================================================================
--- trunk/octave-forge/extra/secs3d/inst/Utilities/block_jacobi.m 2012-05-03 17:03:10 UTC (rev 10355)
+++ trunk/octave-forge/extra/secs3d/inst/Utilities/block_jacobi.m 2012-05-03 17:31:31 UTC (rev 10356)
@@ -50,9 +50,21 @@
if (norm (res, inf) < tol); break; endif
for ib = 1:nblocks
range = first(ib):last(ib);
- pres(range) = diag (diag (A(range, range))) \ res(range);
+ pres(range) = A(range, range) \ res(range);
endfor
x += pres;
endfor
-
-endfunction
\ No newline at end of file
+
+endfunction
+
+%!test
+%! msh = bim3c_mesh_properties (msh3m_structured_mesh (0:.1:1, 0:.1:1, 0:.1:1, 1, 1:6));
+%! x = msh.p(1, :).';
+%! A = bim3a_advection_diffusion (msh, 1, x);
+%! b = bim3a_rhs (msh, 1, 1);
+%! dnodes = bim3c_unknowns_on_faces (msh, [1 2]);
+%! inodes = setdiff (1:numel(x), dnodes);
+%! w = u = zeros (size (x));
+%! u(inodes) = A(inodes, inodes) \ b(inodes);
+%! [w(inodes), ~, nit] = block_jacobi (A(inodes, inodes), b(inodes), b(inodes), 5, 400, 1e-9);
+%! assert (u, w, 1e-7);
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <mtm...@us...> - 2012-05-03 04:11:33
|
Revision: 10354
http://octave.svn.sourceforge.net/octave/?rev=10354&view=rev
Author: mtmiller
Date: 2012-05-03 04:11:25 +0000 (Thu, 03 May 2012)
Log Message:
-----------
comm: fix hdf5 library detection
Modified Paths:
--------------
trunk/octave-forge/main/comm/src/Makefile
Modified: trunk/octave-forge/main/comm/src/Makefile
===================================================================
--- trunk/octave-forge/main/comm/src/Makefile 2012-05-02 05:51:55 UTC (rev 10353)
+++ trunk/octave-forge/main/comm/src/Makefile 2012-05-03 04:11:25 UTC (rev 10354)
@@ -1,6 +1,6 @@
sinclude Makeconf
-HDF5_LIBS := $(shell grep "\#define OCTAVE_CONF_HDF5_LIBS" $(shell $(MKOCTFILE) -p OCTINCLUDEDIR)/oct-conf.h | sed "s/^.*LIBS //" )
+HDF5_LIBS := $(shell grep "\#define OCTAVE_CONF_HDF5_LIBS" $(shell $(MKOCTFILE) -p OCTINCLUDEDIR)/oct-conf.h | sed 's/^.*LIBS //;s/"//g' )
GALOISTARGET = gf.oct
GALOISSOURCES = galois.cc galois-def.cc galoisfield.cc gf.cc \
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-05-02 05:52:01
|
Revision: 10353
http://octave.svn.sourceforge.net/octave/?rev=10353&view=rev
Author: paramaniac
Date: 2012-05-02 05:51:55 +0000 (Wed, 02 May 2012)
Log Message:
-----------
quaternion: doc and test for vector = diag (matrix)
Modified Paths:
--------------
trunk/octave-forge/main/quaternion/inst/@quaternion/diag.m
Modified: trunk/octave-forge/main/quaternion/inst/@quaternion/diag.m
===================================================================
--- trunk/octave-forge/main/quaternion/inst/@quaternion/diag.m 2012-05-01 20:43:31 UTC (rev 10352)
+++ trunk/octave-forge/main/quaternion/inst/@quaternion/diag.m 2012-05-02 05:51:55 UTC (rev 10353)
@@ -1,4 +1,4 @@
-## Copyright (C) 2010 Lukas F. Reichlin
+## Copyright (C) 2010, 2012 Lukas F. Reichlin
##
## This program is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
@@ -22,11 +22,13 @@
## If it is negative, it is placed on the -K-th sub-diagonal.
## The default value of K is 0, and the vector is placed
## on the main diagonal.
+## Given a matrix argument, instead of a vector, @command{diag}
+## extracts the @var{K}-th diagonal of the matrix.
## @end deftypefn
## Author: Lukas Reichlin <luk...@gm...>
## Created: May 2010
-## Version: 0.1
+## Version: 0.2
function a = diag (a, b = 0)
@@ -42,7 +44,7 @@
endfunction
-%!shared R, S
+%!shared R, S, T, U
%! Q = quaternion (2, 3, 4, 5);
%! R = diag ([Q, Q, Q]);
%! W = diag ([2, 2, 2]);
@@ -50,4 +52,7 @@
%! Y = diag ([4, 4, 4]);
%! Z = diag ([5, 5, 5]);
%! S = quaternion (W, X, Y, Z);
+%! T = diag (R);
+%! U = [Q; Q; Q];
%!assert (R == S);
+%!assert (T == U);
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-05-01 20:43:37
|
Revision: 10352
http://octave.svn.sourceforge.net/octave/?rev=10352&view=rev
Author: paramaniac
Date: 2012-05-01 20:43:31 +0000 (Tue, 01 May 2012)
Log Message:
-----------
control: don't compute inverse explicitly (where possible) for state-space inversion
Modified Paths:
--------------
trunk/octave-forge/main/control/inst/@ss/__sys_inverse__.m
Modified: trunk/octave-forge/main/control/inst/@ss/__sys_inverse__.m
===================================================================
--- trunk/octave-forge/main/control/inst/@ss/__sys_inverse__.m 2012-05-01 20:31:07 UTC (rev 10351)
+++ trunk/octave-forge/main/control/inst/@ss/__sys_inverse__.m 2012-05-01 20:43:31 UTC (rev 10352)
@@ -1,4 +1,4 @@
-## Copyright (C) 2009, 2010 Lukas F. Reichlin
+## Copyright (C) 2009, 2010, 2012 Lukas F. Reichlin
##
## This file is part of LTI Syncope.
##
@@ -20,7 +20,7 @@
## Author: Lukas Reichlin <luk...@gm...>
## Created: October 2009
-## Version: 0.2
+## Version: 0.3
function sys = __sys_inverse__ (sys)
@@ -49,12 +49,12 @@
else # proper ss
- di = inv (d);
+ bid = b / d;
- sys.a = a - b * di * c;
- sys.b = -b * di;
- sys.c = di * c;
- sys.d = di;
+ sys.a = a - bid * c;
+ sys.b = -bid;
+ sys.c = d \ c;
+ sys.d = inv (d);
endif
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-05-01 20:31:13
|
Revision: 10351
http://octave.svn.sourceforge.net/octave/?rev=10351&view=rev
Author: paramaniac
Date: 2012-05-01 20:31:07 +0000 (Tue, 01 May 2012)
Log Message:
-----------
control-devel: compute pseudoinverse by SVD instead of \ operator
Modified Paths:
--------------
trunk/octave-forge/extra/control-devel/inst/arx.m
Added Paths:
-----------
trunk/octave-forge/extra/control-devel/devel/test_arx.m
Removed Paths:
-------------
trunk/octave-forge/extra/control-devel/inst/test_arx.m
Copied: trunk/octave-forge/extra/control-devel/devel/test_arx.m (from rev 10337, trunk/octave-forge/extra/control-devel/inst/test_arx.m)
===================================================================
--- trunk/octave-forge/extra/control-devel/devel/test_arx.m (rev 0)
+++ trunk/octave-forge/extra/control-devel/devel/test_arx.m 2012-05-01 20:31:07 UTC (rev 10351)
@@ -0,0 +1,10 @@
+u = [ 0; 0.5; 1; 1; 1; 1; 1 ];
+y = [ 0; 0; 0.25; 0.62; 0.81; 0.90; 0.95 ];
+
+dat = iddata (y, u)
+
+sys = arx (dat, 1, 1)
+
+
+ysim = lsim (sys(1,1), u);
+
Modified: trunk/octave-forge/extra/control-devel/inst/arx.m
===================================================================
--- trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-01 20:26:49 UTC (rev 10350)
+++ trunk/octave-forge/extra/control-devel/inst/arx.m 2012-05-01 20:31:07 UTC (rev 10351)
@@ -28,8 +28,22 @@
PhiU = toeplitz (U(1:end-1, :), [U(1, :); zeros(nb-1, 1)]);
Phi = [-PhiY, PhiU];
- Theta = Phi \ Y(2:end, :);
+ Theta = Phi \ Y(2:end, :) # naive formula
+ ## solve linear least squares problem by pseudoinverse
+ ## the pseudoinverse is computed by singular value decomposition
+ ## M = U S V* ---> M+ = V S+ U*
+ ## Th = Ph \ Y = Ph+ Y
+ ## Th = V S+ U' Y, S+ = 1 ./ diag (S)
+
+ [U, S, V] = svd (Phi, 0); # 0 for "economy size" decomposition, U overwrites input U
+ S = diag (S); # extract main diagonal
+ r = sum (S > eps*S(1));
+ V = V(:, 1:r);
+ S = S(1:r);
+ U = U(:, 1:r);
+ Theta = V * (S .\ (U' * Y(2:end, :))) # U' is the conjugate transpose
+
A = [1; Theta(1:na)]; # a0 = 1, a1 = Theta(1), an = Theta(n)
B = Theta(na+1:end); # b0 = 0 (leading zeros are removed by tf, no need to add one)
Deleted: trunk/octave-forge/extra/control-devel/inst/test_arx.m
===================================================================
--- trunk/octave-forge/extra/control-devel/inst/test_arx.m 2012-05-01 20:26:49 UTC (rev 10350)
+++ trunk/octave-forge/extra/control-devel/inst/test_arx.m 2012-05-01 20:31:07 UTC (rev 10351)
@@ -1,10 +0,0 @@
-u = [ 0; 0.5; 1; 1; 1; 1; 1 ];
-y = [ 0; 0; 0.25; 0.62; 0.81; 0.90; 0.95 ];
-
-dat = iddata (y, u)
-
-sys = arx (dat, 1, 1)
-
-
-ysim = lsim (sys(1,1), u);
-
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-05-01 20:26:55
|
Revision: 10350
http://octave.svn.sourceforge.net/octave/?rev=10350&view=rev
Author: paramaniac
Date: 2012-05-01 20:26:49 +0000 (Tue, 01 May 2012)
Log Message:
-----------
control-devel: use -1 (unspecified) as default sampling time for iddata sets
Modified Paths:
--------------
trunk/octave-forge/extra/control-devel/inst/@iddata/iddata.m
trunk/octave-forge/extra/control-devel/inst/__adjust_iddata_tsam__.m
Modified: trunk/octave-forge/extra/control-devel/inst/@iddata/iddata.m
===================================================================
--- trunk/octave-forge/extra/control-devel/inst/@iddata/iddata.m 2012-05-01 17:16:10 UTC (rev 10349)
+++ trunk/octave-forge/extra/control-devel/inst/@iddata/iddata.m 2012-05-01 20:26:49 UTC (rev 10350)
@@ -41,7 +41,7 @@
## where @var{e} denotes the number of experiments
## and n(i) the individual number of samples for each experiment.
## @item tsam
-## Sampling time.
+## Sampling time. If not specified, default value -1 (unspecified) is taken.
## @item @dots{}
## Optional pairs of properties and values.
## @end table
Modified: trunk/octave-forge/extra/control-devel/inst/__adjust_iddata_tsam__.m
===================================================================
--- trunk/octave-forge/extra/control-devel/inst/__adjust_iddata_tsam__.m 2012-05-01 17:16:10 UTC (rev 10349)
+++ trunk/octave-forge/extra/control-devel/inst/__adjust_iddata_tsam__.m 2012-05-01 20:26:49 UTC (rev 10350)
@@ -18,7 +18,7 @@
## -*- texinfo -*-
## Check whether tsam is a e-by-1 cell array of valid sampling times.
## If not, it tries to convert tsam accordingly.
-## Empty tsam are filled with default value 1.
+## Empty tsam are filled with default value -1.
## Author: Lukas Reichlin <luk...@gm...>
## Created: February 2012
@@ -27,14 +27,14 @@
function tsam = __adjust_iddata_tsam__ (tsam, e)
if (isempty (tsam))
- tsam = num2cell (ones (e, 1));
+ tsam = num2cell (-ones (e, 1));
elseif (iscell (tsam))
tsam = reshape (tsam, [], 1);
else
tsam = {tsam};
endif
- tmp = cellfun (@issample, tsam);
+ tmp = cellfun (@(x) issample (x, -1), tsam);
if (any (! tmp))
error ("iddata: invalid sampling time");
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <as...@us...> - 2012-05-01 17:16:16
|
Revision: 10349
http://octave.svn.sourceforge.net/octave/?rev=10349&view=rev
Author: asnelt
Date: 2012-05-01 17:16:10 +0000 (Tue, 01 May 2012)
Log Message:
-----------
Fixed a texinfo typo.
Modified Paths:
--------------
trunk/octave-forge/main/statistics/inst/kmeans.m
Modified: trunk/octave-forge/main/statistics/inst/kmeans.m
===================================================================
--- trunk/octave-forge/main/statistics/inst/kmeans.m 2012-05-01 16:58:01 UTC (rev 10348)
+++ trunk/octave-forge/main/statistics/inst/kmeans.m 2012-05-01 17:16:10 UTC (rev 10349)
@@ -15,7 +15,7 @@
## this program; if not, see <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
-## @deftypefn {Function File} {[@var{idx}, @var{centers}] =} kmeans (@var{data}, @var{k}, @var{param1}, @var{value1, @dots{})
+## @deftypefn {Function File} {[@var{idx}, @var{centers}] =} kmeans (@var{data}, @var{k}, @var{param1}, @var{value1}, @dots{})
## K-means clustering.
##
## @seealso{linkage}
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
48 messages has been excluded from this view by a project administrator.
