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-03-29 15:45:36
|
Revision: 10098
http://octave.svn.sourceforge.net/octave/?rev=10098&view=rev
Author: paramaniac
Date: 2012-03-29 15:45:27 +0000 (Thu, 29 Mar 2012)
Log Message:
-----------
control-devel: update debugging stuff for ib01ad (2)
Modified Paths:
--------------
trunk/octave-forge/extra/control-devel/devel/PowerPlantTest.m
trunk/octave-forge/extra/control-devel/src/slitest.cc
Modified: trunk/octave-forge/extra/control-devel/devel/PowerPlantTest.m
===================================================================
--- trunk/octave-forge/extra/control-devel/devel/PowerPlantTest.m 2012-03-29 15:33:41 UTC (rev 10097)
+++ trunk/octave-forge/extra/control-devel/devel/PowerPlantTest.m 2012-03-29 15:45:27 UTC (rev 10098)
@@ -77,6 +77,7 @@
for k = 2 : l
err(k-1) = norm (abs(r{1}) - abs(r{k}), 1);
+ % err(k-1) = norm (r{1} - r{k}, 1);
% err(k-1) = norm (abs(abs(r{1}) - abs(r{k})), 1);
endfor
Modified: trunk/octave-forge/extra/control-devel/src/slitest.cc
===================================================================
--- trunk/octave-forge/extra/control-devel/src/slitest.cc 2012-03-29 15:33:41 UTC (rev 10097)
+++ trunk/octave-forge/extra/control-devel/src/slitest.cc 2012-03-29 15:45:27 UTC (rev 10098)
@@ -106,15 +106,15 @@
////////////////////////////////////////////////////////////////////////////////////
// arguments in
- char meth = 'N'; // <--- not used, use metha
- char alg = 'C';
- char jobd = 'N';
+ char meth = 'M'; // <--- not used, use metha
+ char alg = 'Q';
+ char jobd = 'M';
char batch = 'O';
char conct = 'N';
char ctrl = 'N';
- char metha = 'N';
+ char metha = 'M';
// ??? char jobda ;
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-03-29 15:33:52
|
Revision: 10097
http://octave.svn.sourceforge.net/octave/?rev=10097&view=rev
Author: paramaniac
Date: 2012-03-29 15:33:41 +0000 (Thu, 29 Mar 2012)
Log Message:
-----------
control-devel: update debugging stuff for ib01ad
Modified Paths:
--------------
trunk/octave-forge/extra/control-devel/devel/identtest.m
trunk/octave-forge/extra/control-devel/src/slident.cc
trunk/octave-forge/extra/control-devel/src/slitest.cc
Modified: trunk/octave-forge/extra/control-devel/devel/identtest.m
===================================================================
--- trunk/octave-forge/extra/control-devel/devel/identtest.m 2012-03-29 15:21:47 UTC (rev 10096)
+++ trunk/octave-forge/extra/control-devel/devel/identtest.m 2012-03-29 15:33:41 UTC (rev 10097)
@@ -1,52 +1,5 @@
function [r, sv, n] = identtest (dat, s = [], n = [], ldwork)
+ [r, sv, n] = slitest (dat.y{1}, dat.u{1}, s, ldwork);
-
- %nobr = 15;
- meth = 2;
- alg = 0;
- 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] = slitest (dat.y{1}, dat.u{1}, nobr, n, meth, alg, jobd, batch, conct, ctrl, rcond, tol, ldwork);
-
endfunction
Modified: trunk/octave-forge/extra/control-devel/src/slident.cc
===================================================================
--- trunk/octave-forge/extra/control-devel/src/slident.cc 2012-03-29 15:21:47 UTC (rev 10096)
+++ trunk/octave-forge/extra/control-devel/src/slident.cc 2012-03-29 15:33:41 UTC (rev 10097)
@@ -113,7 +113,7 @@
char conct;
char ctrl;
char metha;
- char jobda;
+ char jobda; // ??? unused
Matrix y = args(0).matrix_value ();
Matrix u = args(1).matrix_value ();
Modified: trunk/octave-forge/extra/control-devel/src/slitest.cc
===================================================================
--- trunk/octave-forge/extra/control-devel/src/slitest.cc 2012-03-29 15:21:47 UTC (rev 10096)
+++ trunk/octave-forge/extra/control-devel/src/slitest.cc 2012-03-29 15:33:41 UTC (rev 10097)
@@ -95,7 +95,7 @@
int nargin = args.length ();
octave_value_list retval;
- if (nargin != 13)
+ if (nargin != 4)
{
print_usage ();
}
@@ -106,127 +106,36 @@
////////////////////////////////////////////////////////////////////////////////////
// arguments in
- char meth;
- char alg;
- char jobd;
- char batch;
- char conct;
- char ctrl;
- char metha;
- char jobda;
+ char meth = 'N'; // <--- not used, use metha
+ char alg = 'C';
+ char jobd = 'N';
+ char batch = 'O';
+ char conct = 'N';
+ char ctrl = 'N';
+
+ char metha = 'N';
+
+ // ??? char jobda ;
+
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 = 0.0;
+ double tol = -1.0;
- double rcond = args(10).double_value ();
- double tol = args(11).double_value ();
- double tolb = args(10).double_value (); // tolb = rcond
-
- int ldwork = args(12).int_value ();
+ int ldwork = args(3).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 ldu = nsmp;
int ldy = nsmp;
// arguments out
@@ -322,19 +231,7 @@
*/
-/*
-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);
@@ -401,13 +298,6 @@
int rs = 2*(m+l)*nobr;
r.resize (rs, rs);
- if (nuser > 0)
- {
- if (nuser < nobr)
- n = nuser;
- else
- error ("ident: 'nuser' invalid");
- }
retval(0) = r;
retval(1) = sv;
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <cde...@us...> - 2012-03-29 15:21:56
|
Revision: 10096
http://octave.svn.sourceforge.net/octave/?rev=10096&view=rev
Author: cdemills
Date: 2012-03-29 15:21:47 +0000 (Thu, 29 Mar 2012)
Log Message:
-----------
- Take more than one field for datefmt iff sep contains space
Modified Paths:
--------------
trunk/octave-forge/extra/dataframe/inst/@dataframe/dataframe.m
Modified: trunk/octave-forge/extra/dataframe/inst/@dataframe/dataframe.m
===================================================================
--- trunk/octave-forge/extra/dataframe/inst/@dataframe/dataframe.m 2012-03-29 15:09:25 UTC (rev 10095)
+++ trunk/octave-forge/extra/dataframe/inst/@dataframe/dataframe.m 2012-03-29 15:21:47 UTC (rev 10096)
@@ -162,7 +162,12 @@
if (~isempty (datefmt))
%# replace consecutive spaces by one
datefmt = regexprep (datefmt, '[ ]+', ' ');
- datefields = 1 + length (regexp (datefmt, ' '));
+ %# is "space" used as separator ? Then we may take more than one field.
+ if (~isempty (regexp (sep, ' ')))
+ datefields = 1 + length (regexp (datefmt, ' '));
+ else
+ datefields = 1;
+ endif
else
datefields = 1;
endif
@@ -323,6 +328,9 @@
for indc = (2:datefields)
datetime = cstrcat(datetime, ' ', dummy{indk+indc-1});
endfor
+ else
+ %# ensure spaces are unique
+ datetime = regexprep (datetime, '[ ]+', ' ');
endif
try
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <cde...@us...> - 2012-03-29 15:09:35
|
Revision: 10095
http://octave.svn.sourceforge.net/octave/?rev=10095&view=rev
Author: cdemills
Date: 2012-03-29 15:09:25 +0000 (Thu, 29 Mar 2012)
Log Message:
-----------
- If the datefmt contains spaces, then concatenate more than one field.
Modified Paths:
--------------
trunk/octave-forge/extra/dataframe/inst/@dataframe/dataframe.m
Modified: trunk/octave-forge/extra/dataframe/inst/@dataframe/dataframe.m
===================================================================
--- trunk/octave-forge/extra/dataframe/inst/@dataframe/dataframe.m 2012-03-29 15:08:03 UTC (rev 10094)
+++ trunk/octave-forge/extra/dataframe/inst/@dataframe/dataframe.m 2012-03-29 15:09:25 UTC (rev 10095)
@@ -159,6 +159,14 @@
endwhile
endif
+if (~isempty (datefmt))
+ %# replace consecutive spaces by one
+ datefmt = regexprep (datefmt, '[ ]+', ' ');
+ datefields = 1 + length (regexp (datefmt, ' '));
+else
+ datefields = 1;
+endif
+
if (~isempty (seeked) && ~isempty (trigger))
error ('seeked and trigger are mutuallly incompatible arguments');
endif
@@ -283,45 +291,62 @@
'UniformOutput', false);
endif
- for indk = (1:size (the_line, 2))
+ indk = 1; indm = 1;
+ while (indk <= size (the_line, 2))
if (isempty (the_line{indk}) || any (size (the_line{indk}) > 1))
%#if indi > 1 && indk > 1, disp('line 117 '); keyboard; %#endif
if (unquot)
try
%# remove quotes and leading space(s)
- x(indj, indk) = regexp (dummy{indk}, '[^'' ].*[^'']', 'match'){1};
+ x(indj, indm) = regexp (dummy{indk}, '[^'' ].*[^'']', 'match'){1};
catch
%# if the previous test fails, try a simpler one
in = regexp (dummy{indk}, '[^'' ]+', 'match');
if (~isempty (in))
- x(indj, indk) = in{1};
+ x(indj, indm) = in{1};
%# else
%# x(indj, indk) = [];
endif
end_try_catch
else
%# no conversion possible, store and remove leading space(s)
- x(indj, indk) = regexp (dummy{indk}, '[^ ].*', 'match');
+ x(indj, indm) = regexp (dummy{indk}, '[^ ].*', 'match');
endif
elseif (~isempty (regexp (dummy{indk}, '[/:-]')) && ...
~isempty (datefmt))
+ %# does it look like a date ?
+ datetime = dummy{indk};
+
+ if (datefields > 1)
+ %# concatenate the required number of fields
+ indc = 1;
+ for indc = (2:datefields)
+ datetime = cstrcat(datetime, ' ', dummy{indk+indc-1});
+ endfor
+ endif
+
try
- datetime = datevec (dummy{indk}, datefmt);
+ datetime = datevec (datetime, datefmt);
timeval = struct ("usec", 0, "sec", floor (datetime (6)),
"min", datetime(5), "hour", datetime(4),
"mday", datetime(3), "mon", datetime(2)-1,
"year", datetime(1)-1900);
- timeval.usec = 1e6*(datetime(6)-timeval.sec);
- x(indj, indk) = str2num (strftime ([char(37) 's'], timeval)) + ...
+ timeval.usec = 1e6*(datetime(6) - timeval.sec);
+ x(indj, indm) = str2num (strftime ([char(37) 's'], timeval)) + ...
timeval.usec * 1e-6;
+ if (datefields > 1)
+ %# skip fields successfully converted
+ indk = indk + (datefields - 1);
+ endif
catch
%# store it as is
- x(indj, indk) = the_line{indk};
+ x(indj, indm) = the_line{indk};
end_try_catch
else
- x(indj, indk) = the_line{indk};
+ x(indj, indm) = the_line{indk};
endif
- endfor
+ indk = indk + 1; indm = indm + 1;
+ endwhile
indl = indl + 1; indj = indj + 1;
endwhile
@@ -336,12 +361,12 @@
endif
clear UTF8_BOM fid in lines indl the_line content empty_lines
- clear timeval timestr nfields idx
+ clear datetime timeval idx
endif
end_try_catch
endif
-
+
%# fallback, avoiding a recursive call
idx.type = '()';
if (~isa (x, 'char'))
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <cde...@us...> - 2012-03-29 15:08:14
|
Revision: 10094
http://octave.svn.sourceforge.net/octave/?rev=10094&view=rev
Author: cdemills
Date: 2012-03-29 15:08:03 +0000 (Thu, 29 Mar 2012)
Log Message:
-----------
- Made the cast to 'char' more usefull
Modified Paths:
--------------
trunk/octave-forge/extra/dataframe/inst/@dataframe/subsasgn.m
Modified: trunk/octave-forge/extra/dataframe/inst/@dataframe/subsasgn.m
===================================================================
--- trunk/octave-forge/extra/dataframe/inst/@dataframe/subsasgn.m 2012-03-29 15:07:17 UTC (rev 10093)
+++ trunk/octave-forge/extra/dataframe/inst/@dataframe/subsasgn.m 2012-03-29 15:08:03 UTC (rev 10094)
@@ -30,7 +30,7 @@
if (isnull (df))
error ('dataframe subsasgn: first argument may not be empty');
endif
-
+
switch (S(1).type)
case '{}'
error ('Invalid dataframe as cell assignement');
@@ -65,6 +65,19 @@
if (isnull(RHS)) error("Types can't be nulled"); endif
if (1 == length (S))
%# perform explicit cast on each column
+ switch (RHS)
+ case {'char'}
+ for indj = (1:df._cnt(2))
+ if (isnumeric (df._data{indj}) || islogical (df._data{indj}))
+ df._data(indj) = cellfun (@(x) cellstr (num2str(x, "%f")), \
+ df._data(indj),
+ "UniformOutput", false);
+ endif
+ endfor
+ otherwise
+ df._data = cellfun (@(x) cast (x, RHS), df._data,
+ "UniformOutput", false);
+ endswitch
df._data = cellfun (@(x) cast (x, RHS), df._data,
"UniformOutput", false);
df._type(1:end) = RHS;
@@ -81,8 +94,17 @@
else
indj = S(2).subs{1}; ncol = length (indj);
endif
- df._data(indj) = cellfun (@(x) cast (x, RHS), df._data(indj),
- "UniformOutput", false);
+ switch (RHS)
+ case {'char'}
+ if (isnumeric (df._data{indj}) || islogical (df._data{indj}))
+ df._data(indj) = cellfun (@(x) cellstr (num2str(x, "%f")), \
+ df._data(indj),
+ "UniformOutput", false);
+ endif
+ otherwise
+ df._data(indj) = cellfun (@(x) cast (x, RHS), df._data(indj),
+ "UniformOutput", false);
+ endswitch
df._type(indj) = {RHS};
endif
return
@@ -180,9 +202,7 @@
if (~isnull (RHS))
S(1).subs{1} = indr; S(1).subs{2} = indc;
endif
-
df = df_matassign (df, S, indc, ncol, RHS);
-
endswitch
%# disp("end of subasgn"); keyboard
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-03-29 15:07:27
|
Revision: 10093
http://octave.svn.sourceforge.net/octave/?rev=10093&view=rev
Author: paramaniac
Date: 2012-03-29 15:07:17 +0000 (Thu, 29 Mar 2012)
Log Message:
-----------
control-devel: debugging stuff for ib01ad
Modified Paths:
--------------
trunk/octave-forge/extra/control-devel/inst/@iddata/set.m
trunk/octave-forge/extra/control-devel/src/devel_slicot_functions.cc
Added Paths:
-----------
trunk/octave-forge/extra/control-devel/devel/PowerPlantTest.m
trunk/octave-forge/extra/control-devel/devel/identtest.m
trunk/octave-forge/extra/control-devel/src/slitest.cc
Added: trunk/octave-forge/extra/control-devel/devel/PowerPlantTest.m
===================================================================
--- trunk/octave-forge/extra/control-devel/devel/PowerPlantTest.m (rev 0)
+++ trunk/octave-forge/extra/control-devel/devel/PowerPlantTest.m 2012-03-29 15:07:17 UTC (rev 10093)
@@ -0,0 +1,83 @@
+%{
+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);
+
+%}
+
+clear all, 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);
+
+
+% ldwork = [401, 802, 1203, 1604]
+% warning: implicit conversion from real matrix to real scalar
+
+ldwork = [802, 1203, 1604]
+
+r = arrayfun (@(x) identtest (dat, 10, 8, x), ldwork, 'uniformoutput', false);
+
+ % s=10, n=8
+
+l = length (ldwork);
+
+err = cell (l-1, 1);
+
+for k = 2 : l
+ err(k-1) = norm (abs(r{1}) - abs(r{k}), 1);
+ % err(k-1) = norm (abs(abs(r{1}) - abs(r{k})), 1);
+endfor
+
+err
Added: trunk/octave-forge/extra/control-devel/devel/identtest.m
===================================================================
--- trunk/octave-forge/extra/control-devel/devel/identtest.m (rev 0)
+++ trunk/octave-forge/extra/control-devel/devel/identtest.m 2012-03-29 15:07:17 UTC (rev 10093)
@@ -0,0 +1,52 @@
+function [r, sv, n] = identtest (dat, s = [], n = [], ldwork)
+
+
+
+ %nobr = 15;
+ meth = 2;
+ alg = 0;
+ 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] = slitest (dat.y{1}, dat.u{1}, nobr, n, meth, alg, jobd, batch, conct, ctrl, rcond, tol, ldwork);
+
+endfunction
Modified: trunk/octave-forge/extra/control-devel/inst/@iddata/set.m
===================================================================
--- trunk/octave-forge/extra/control-devel/inst/@iddata/set.m 2012-03-29 14:33:05 UTC (rev 10092)
+++ trunk/octave-forge/extra/control-devel/inst/@iddata/set.m 2012-03-29 15:07:17 UTC (rev 10093)
@@ -84,8 +84,6 @@
dat.outunit = __adjust_labels__ (val, p);
case {"inunit", "inputunit"}
dat.inunit = __adjust_labels__ (val, m);
- case {"tsam", "ts"}
- dat.tsam;
case {"timeunit"}
dat.timeunit
case {"expname", "experimentname"}
Modified: trunk/octave-forge/extra/control-devel/src/devel_slicot_functions.cc
===================================================================
--- trunk/octave-forge/extra/control-devel/src/devel_slicot_functions.cc 2012-03-29 14:33:05 UTC (rev 10092)
+++ trunk/octave-forge/extra/control-devel/src/devel_slicot_functions.cc 2012-03-29 15:07:17 UTC (rev 10093)
@@ -1,2 +1,3 @@
// #include "slib01ad.cc" // preprocess the input-output data
#include "slident.cc" // system identification
+#include "slitest.cc" // debug system identification
Added: trunk/octave-forge/extra/control-devel/src/slitest.cc
===================================================================
--- trunk/octave-forge/extra/control-devel/src/slitest.cc (rev 0)
+++ trunk/octave-forge/extra/control-devel/src/slitest.cc 2012-03-29 15:07:17 UTC (rev 10093)
@@ -0,0 +1,418 @@
+/*
+
+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 ("slitest", "devel_slicot_functions.oct");
+DEFUN_DLD (slitest, 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 != 13)
+ {
+ 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;
+
+ 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
+
+ int ldwork = args(12).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;
+
+ 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;
+ }
+ }
+////////////////////////////////////////////////////////////////////
+// TO BE REMOVED !!!
+////////////////////////////////////////////////////////////////////
+ldwork *= 2;
+//ldwork = 1000000;
+////////////////////////////////////////////////////////////////////
+
+*/
+
+/*
+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);
+
+ if (nuser > 0)
+ {
+ if (nuser < nobr)
+ n = nuser;
+ else
+ error ("ident: 'nuser' invalid");
+ }
+
+ retval(0) = r;
+ retval(1) = sv;
+ retval(2) = octave_value (n);
+ }
+
+ return retval;
+}
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <cde...@us...> - 2012-03-29 14:33:15
|
Revision: 10092
http://octave.svn.sourceforge.net/octave/?rev=10092&view=rev
Author: cdemills
Date: 2012-03-29 14:33:05 +0000 (Thu, 29 Mar 2012)
Log Message:
-----------
- A few stylistic changes
Modified Paths:
--------------
trunk/octave-forge/extra/dataframe/inst/@dataframe/private/df_matassign.m
Modified: trunk/octave-forge/extra/dataframe/inst/@dataframe/private/df_matassign.m
===================================================================
--- trunk/octave-forge/extra/dataframe/inst/@dataframe/private/df_matassign.m 2012-03-29 07:52:00 UTC (rev 10091)
+++ trunk/octave-forge/extra/dataframe/inst/@dataframe/private/df_matassign.m 2012-03-29 14:33:05 UTC (rev 10092)
@@ -391,7 +391,7 @@
catch
dummy = \
sprintf ("Assignement failed for colum %d, of type %s and length %d,\nwith new content\n%s", \
- indj, df._type{indc(indi)}, length (indr), disp(RHS(:, indj)));
+ indj, df._type{indc(indi)}, length (indr), disp (RHS(:, indj)));
error (dummy);
end_try_catch
endif
@@ -426,14 +426,14 @@
[df, S] = df_cow(df, S, indc(indi));
if (strcmp (df._type(indc(indi)), RHS._type(indi)))
try
- df._data{indc(indi)} = feval(@subsasgn, df._data{indc(indi)}, S, \
- RHS._data{indi}(:, RHS._rep{indi}));
+ df._data{indc(indi)} = feval (@subsasgn, df._data{indc(indi)}, S, \
+ RHS._data{indi}(:, RHS._rep{indi}));
catch
- disp(lasterr()); disp('line 516 ???'); keyboard
+ disp (lasterr ()); disp('line 516 ???'); keyboard
end_try_catch
else
- df._data{indc(indi)} = feval(@subsasgn, df._data{indc(indi)}, S, \
- cast(RHS._data{indi}(:, RHS._rep{indi}),\
+ df._data{indc(indi)} = feval (@subsasgn, df._data{indc(indi)}, S, \
+ cast (RHS._data{indi}(:, RHS._rep{indi}),\
df._type(indc(indi))));
endif
S = Sorig;
@@ -480,8 +480,8 @@
df._data{indc(indi)} = fillfunc (df._data{indc(indi)}, S, indi);
S = Sorig;
catch
- disp(lasterr)
- disp('line 470 '); keyboard
+ disp (lasterr ())
+ disp ('line 470 '); keyboard
end_try_catch
# catch
# if ndims(df._data{indc(indi)}) > 2,
@@ -557,10 +557,10 @@
if (length (dummy) == ncol)
df._name{2}(indc, 1) = dummy;
else
- disp('line 528 '); keyboard
+ disp ('line 528 '); keyboard
endif
else
- disp('line 531 '); keyboard
+ disp ('line 531 '); keyboard
endif
end_try_catch
df._over{2}(1, indc) = false;
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <cde...@us...> - 2012-03-29 07:52:07
|
Revision: 10091
http://octave.svn.sourceforge.net/octave/?rev=10091&view=rev
Author: cdemills
Date: 2012-03-29 07:52:00 +0000 (Thu, 29 Mar 2012)
Log Message:
-----------
- Updated revision number
Modified Paths:
--------------
trunk/octave-forge/extra/dataframe/DESCRIPTION
Modified: trunk/octave-forge/extra/dataframe/DESCRIPTION
===================================================================
--- trunk/octave-forge/extra/dataframe/DESCRIPTION 2012-03-29 07:38:06 UTC (rev 10090)
+++ trunk/octave-forge/extra/dataframe/DESCRIPTION 2012-03-29 07:52:00 UTC (rev 10091)
@@ -1,6 +1,6 @@
Name: dataframe
-Version: 0.9.1
-Date: 2012-02-10
+Version: 0.9.2
+Date: 2012-03-28
Author: Pascal Dupuis
Maintainer: The Octave Community
Title: Data Frame
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <cde...@us...> - 2012-03-29 07:38:13
|
Revision: 10090
http://octave.svn.sourceforge.net/octave/?rev=10090&view=rev
Author: cdemills
Date: 2012-03-29 07:38:06 +0000 (Thu, 29 Mar 2012)
Log Message:
-----------
- Style fixup
Modified Paths:
--------------
trunk/octave-forge/extra/dataframe/inst/@dataframe/ldivide.m
Modified: trunk/octave-forge/extra/dataframe/inst/@dataframe/ldivide.m
===================================================================
--- trunk/octave-forge/extra/dataframe/inst/@dataframe/ldivide.m 2012-03-29 07:37:44 UTC (rev 10089)
+++ trunk/octave-forge/extra/dataframe/inst/@dataframe/ldivide.m 2012-03-29 07:38:06 UTC (rev 10090)
@@ -28,10 +28,10 @@
%#
try
- resu = df_func(@ldivide, A, B);
+ resu = df_func (@ldivide, A, B);
catch
- disp(lasterr());
- error("Operator .\\ problem for %s vs. %s", class(A), class(B));
+ disp (lasterr ());
+ error ("Operator .\\ problem for %s vs. %s", class(A), class(B));
end_try_catch
endfunction
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <cde...@us...> - 2012-03-29 07:37:56
|
Revision: 10089
http://octave.svn.sourceforge.net/octave/?rev=10089&view=rev
Author: cdemills
Date: 2012-03-29 07:37:44 +0000 (Thu, 29 Mar 2012)
Log Message:
-----------
- Style fixup
Modified Paths:
--------------
trunk/octave-forge/extra/dataframe/inst/@dataframe/mldivide.m
Modified: trunk/octave-forge/extra/dataframe/inst/@dataframe/mldivide.m
===================================================================
--- trunk/octave-forge/extra/dataframe/inst/@dataframe/mldivide.m 2012-03-29 07:34:43 UTC (rev 10088)
+++ trunk/octave-forge/extra/dataframe/inst/@dataframe/mldivide.m 2012-03-29 07:37:44 UTC (rev 10089)
@@ -28,10 +28,10 @@
%#
try
- resu = df_func(@mldivide, A, B, false, [true false]);
+ resu = df_func (@mldivide, A, B, false, [true false]);
catch
- disp(lasterr());
- error("Operator \\ problem for %s vs. %s", class(A), class(B));
+ disp (lasterr ());
+ error ("Operator \\ problem for %s vs. %s", class (A), class (B));
end_try_catch
endfunction
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <cde...@us...> - 2012-03-29 07:34:49
|
Revision: 10088
http://octave.svn.sourceforge.net/octave/?rev=10088&view=rev
Author: cdemills
Date: 2012-03-29 07:34:43 +0000 (Thu, 29 Mar 2012)
Log Message:
-----------
- Replaced a few else if by elseif
Modified Paths:
--------------
trunk/octave-forge/extra/dataframe/inst/@dataframe/dataframe.m
Modified: trunk/octave-forge/extra/dataframe/inst/@dataframe/dataframe.m
===================================================================
--- trunk/octave-forge/extra/dataframe/inst/@dataframe/dataframe.m 2012-03-28 20:30:21 UTC (rev 10087)
+++ trunk/octave-forge/extra/dataframe/inst/@dataframe/dataframe.m 2012-03-29 07:34:43 UTC (rev 10088)
@@ -90,7 +90,7 @@
endif
%# default values
-seeked = []; trigger =[]; unquot = true; sep = "\t,"; cmt_lines = [];
+seeked = []; trigger = []; unquot = true; sep = "\t,"; cmt_lines = [];
locales = "C"; datefmt = '';
if (length (varargin) > 0)
@@ -125,8 +125,8 @@
%# detect assignment - functions calls - ranges
dummy = cellfun ('size', cellfun (@(x) strsplit (x, ":=("), df._name{2}, \
"UniformOutput", false), 2);
- if (any(dummy > 1))
- warning('dataframe colnames taken literally and not interpreted');
+ if (any (dummy > 1))
+ warning ('dataframe colnames taken literally and not interpreted');
endif
df._name{2} = genvarname (df._name{2});
df._over{2}(1, 1:length (df._name{2})) = false;
@@ -276,10 +276,11 @@
the_line = cellfun (@(x) sscanf (x, "%f", locales), dummy, \
'UniformOutput', false);
else
- %# this code require a patch to src/file-io.cc in main
- %# Octave tree
+ %# this faster code requires a patch to src/file-io.cc in
+ %# main Octave tree
the_line = sscanf (dummy, "%f", locales);
- the_line = cellfun (@(x) x{1}, the_line, 'UniformOutput', false);
+ the_line = cellfun (@(x) x{1}, the_line, \
+ 'UniformOutput', false);
endif
for indk = (1:size (the_line, 2))
@@ -302,40 +303,41 @@
%# no conversion possible, store and remove leading space(s)
x(indj, indk) = regexp (dummy{indk}, '[^ ].*', 'match');
endif
+ elseif (~isempty (regexp (dummy{indk}, '[/:-]')) && ...
+ ~isempty (datefmt))
+ try
+ datetime = datevec (dummy{indk}, datefmt);
+ timeval = struct ("usec", 0, "sec", floor (datetime (6)),
+ "min", datetime(5), "hour", datetime(4),
+ "mday", datetime(3), "mon", datetime(2)-1,
+ "year", datetime(1)-1900);
+ timeval.usec = 1e6*(datetime(6)-timeval.sec);
+ x(indj, indk) = str2num (strftime ([char(37) 's'], timeval)) + ...
+ timeval.usec * 1e-6;
+ catch
+ %# store it as is
+ x(indj, indk) = the_line{indk};
+ end_try_catch
else
- if (~isempty (regexp (dummy{indk}, '[/:-]')) && ...
- ~isempty (datefmt))
-
- try
- datetime = datevec (dummy{indk}, datefmt);
- timeval = struct ("usec", 0, "sec", floor (datetime (6)),
- "min", datetime(5), "hour", datetime(4),
- "mday", datetime(3), "mon", datetime(2)-1,
- "year", datetime(1)-1900);
- timeval.usec = 1e6*(datetime(6)-timeval.sec);
- x(indj, indk) = str2num (strftime ([char(37) 's'], timeval)) + ...
- timeval.usec * 1e-6;
- catch
- %# store it as is
- x(indj, indk) = the_line{indk};
- end_try_catch
- else
- x(indj, indk) = the_line{indk};
- endif
+ x(indj, indk) = the_line{indk};
endif
endfor
indl = indl + 1; indj = indj + 1;
endwhile
+
if (~isempty (empty_lines))
x(empty_lines, :) = [];
endif
+
%# detect empty columns
empty_lines = find (0 == sum (cellfun ('size', x, 2)));
if (~isempty (empty_lines))
x(:, empty_lines) = [];
endif
+
clear UTF8_BOM fid in lines indl the_line content empty_lines
clear timeval timestr nfields idx
+
endif
end_try_catch
endif
@@ -365,11 +367,10 @@
df = df_pad (df, 2, [length(dummy) indc], dummy);
x = x{2};
indj = indc + (1:size (x, 2)); %# redefine target range
- else
- if (isa (x{1}, 'cell'))
- x = x{1}; %# remove one cell level
- endif
+ elseif (isa (x{1}, 'cell'))
+ x = x{1}; %# remove one cell level
endif
+
if (length (df._name{2}) < indj(1) || isempty (df._name{2}(indj)))
[df._name{2}(indj, 1), df._over{2}(1, indj)] ...
= df_colnames (inputname(indi), indj);
@@ -377,19 +378,19 @@
endif
%# allow overwriting of column names
df._over{2}(1, indj) = true;
- else
- if (~isempty (indj))
- if (1 == length (df._name{2}) && length (df._name{2}) < \
- length (indj))
- [df._name{2}(indj, 1), df._over{2}(1, indj)] ...
- = df_colnames (char (df._name{2}), indj);
- elseif (length (df._name{2}) < indj(1) || isempty (df._name{2}(indj)))
- [df._name{2}(indj, 1), df._over{2}(1, indj)] ...
- = df_colnames (inputname(indi), indj);
- endif
- df._name{2} = genvarname (df._name{2});
+
+ elseif (~isempty (indj))
+ if (1 == length (df._name{2}) && length (df._name{2}) < \
+ length (indj))
+ [df._name{2}(indj, 1), df._over{2}(1, indj)] ...
+ = df_colnames (char (df._name{2}), indj);
+ elseif (length (df._name{2}) < indj(1) || isempty (df._name{2}(indj)))
+ [df._name{2}(indj, 1), df._over{2}(1, indj)] ...
+ = df_colnames (inputname(indi), indj);
endif
+ df._name{2} = genvarname (df._name{2});
endif
+
if (~isempty (indj))
%# the exact row size will be determined latter
idx.subs = {'', indj};
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <jpi...@us...> - 2012-03-28 20:30:27
|
Revision: 10087
http://octave.svn.sourceforge.net/octave/?rev=10087&view=rev
Author: jpicarbajal
Date: 2012-03-28 20:30:21 +0000 (Wed, 28 Mar 2012)
Log Message:
-----------
system-identification: TISEAN d2 wrapper
Modified Paths:
--------------
trunk/octave-forge/main/system-identification/inst/tisean/corrdim.m
Modified: trunk/octave-forge/main/system-identification/inst/tisean/corrdim.m
===================================================================
--- trunk/octave-forge/main/system-identification/inst/tisean/corrdim.m 2012-03-28 20:21:17 UTC (rev 10086)
+++ trunk/octave-forge/main/system-identification/inst/tisean/corrdim.m 2012-03-28 20:30:21 UTC (rev 10087)
@@ -20,13 +20,15 @@
%%
%% @end deftypefn
-function [c d h] = corrdim (data, comp=1, maxedim=10, varargin)
+function [c d h] = corrdim (data, varargin)
[nT M] = size (data);
# --- Parse arguments --- #
parser = inputParser ();
parser.FunctionName = "corrdim";
+ parser = addParamValue (parser,'Components', 1, @(x)x>0);
+ parser = addParamValue (parser,'MaxEdim', 10, @(x)x>0);
parser = addParamValue (parser,'Delay', 1, @(x)x>0);
parser = addParamValue (parser,'TheilerWindow', 0, @(x)x>=0);
parser = addParamValue (parser,'ScaleSpan', nT*[1e-3 1], @(x)all(x>0));
@@ -50,7 +52,8 @@
%% Prepare format of the embedding vector
syscmd = sprintf ("d2 -d%d -M%d,%d -t%d -r%d -R%d -#%d -N%d %s -o%s -V0 %s", ...
parser.Results.Delay, ...
- comp, maxedim, ...
+ parser.Results.Components, ...
+ parser.Results.MaxEdim, ...
parser.Results.TheilerWindow, ...
parser.Results.ScaleSpan, ...
parser.Results.EpsilonCount, ...
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <jpi...@us...> - 2012-03-28 20:21:23
|
Revision: 10086
http://octave.svn.sourceforge.net/octave/?rev=10086&view=rev
Author: jpicarbajal
Date: 2012-03-28 20:21:17 +0000 (Wed, 28 Mar 2012)
Log Message:
-----------
system-identification: Adding wrapper to TISEAN function d2
Added Paths:
-----------
trunk/octave-forge/main/system-identification/inst/tisean/corrdim.m
Added: trunk/octave-forge/main/system-identification/inst/tisean/corrdim.m
===================================================================
--- trunk/octave-forge/main/system-identification/inst/tisean/corrdim.m (rev 0)
+++ trunk/octave-forge/main/system-identification/inst/tisean/corrdim.m 2012-03-28 20:21:17 UTC (rev 10086)
@@ -0,0 +1,68 @@
+%% Copyright (c) 2012 Juan Pablo Carbajal <car...@if...>
+%%
+%% 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
+%% 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{c} @var{d} @var{h} @var{stat}] = } corrdim (@var{data}, @var{edim})
+%% Correlation dimension from @var{data}. This function calls @code{d2} @
+%% from the TISEAN package.
+%%
+%% @end deftypefn
+
+function [c d h] = corrdim (data, comp=1, maxedim=10, varargin)
+
+
+ [nT M] = size (data);
+ # --- Parse arguments --- #
+ parser = inputParser ();
+ parser.FunctionName = "corrdim";
+ parser = addParamValue (parser,'Delay', 1, @(x)x>0);
+ parser = addParamValue (parser,'TheilerWindow', 0, @(x)x>=0);
+ parser = addParamValue (parser,'ScaleSpan', nT*[1e-3 1], @(x)all(x>0));
+ parser = addParamValue (parser,'EpsilonCount', 100, @(x)x>0);
+ parser = addParamValue (parser,'PairCount', 1000, @(x)x>=0);
+ parser = addParamValue (parser,'Normalize', false);
+ parser = addParamValue (parser,'Verbose', false);
+ parser = parse(parser,varargin{:});
+
+ flag.E = "";
+ if parser.Results.Normalize
+ flag.E = "-E";
+ end
+
+ infile = tmpnam ();
+ outfile = tmpnam ();
+
+ %% Write data to file
+ save ('-ascii',infile, 'data');
+
+ %% Prepare format of the embedding vector
+ syscmd = sprintf ("d2 -d%d -M%d,%d -t%d -r%d -R%d -#%d -N%d %s -o%s -V0 %s", ...
+ parser.Results.Delay, ...
+ comp, maxedim, ...
+ parser.Results.TheilerWindow, ...
+ parser.Results.ScaleSpan, ...
+ parser.Results.EpsilonCount, ...
+ parser.Results.PairCount, ...
+ flag.E, outfile, infile);
+
+ %% Function call
+ system (syscmd);
+
+ c = load ([outfile ".c2"]);
+ d = load ([outfile ".d2"]);
+ h = load ([outfile ".h2"]);
+
+
+endfunction
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <jpi...@us...> - 2012-03-28 17:15:24
|
Revision: 10085
http://octave.svn.sourceforge.net/octave/?rev=10085&view=rev
Author: jpicarbajal
Date: 2012-03-28 17:15:15 +0000 (Wed, 28 Mar 2012)
Log Message:
-----------
system-identification: adding first TISEAN wrapper
Added Paths:
-----------
trunk/octave-forge/main/system-identification/inst/tisean/delayvec.m
Added: trunk/octave-forge/main/system-identification/inst/tisean/delayvec.m
===================================================================
--- trunk/octave-forge/main/system-identification/inst/tisean/delayvec.m (rev 0)
+++ trunk/octave-forge/main/system-identification/inst/tisean/delayvec.m 2012-03-28 17:15:15 UTC (rev 10085)
@@ -0,0 +1,64 @@
+%% Copyright (c) 2012 Juan Pablo Carbajal <car...@if...>
+%%
+%% 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
+%% 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{delayed} = } delayvec (@var{data}, @var{edim})
+%% Produces delay vectors from @var{data}. This function calls @code{delay} @
+%% from the TISEAN package.
+%%
+%% @end deftypefn
+
+function delayed = delayvec (data, edim=1, dformat=[], delays=[])
+
+
+ [nT M] = size (data);
+ flag.F = "";
+ flag.D = "";
+ %% Check arguments
+ if (edim > M && mod (edim, M) && isempty (dformat) )
+
+ print_usage ();
+
+ elseif edim < M && isempty (dformat)
+
+ print_usage ();
+
+ elseif !isempty (dformat)
+
+ flag.F = sprintf(["-F%d" repmat(',%d',1,length(dformat)-1)],dformat);
+
+ end
+
+ if !isempty (delays)
+
+ flag.D = sprintf(["-D%d" repmat(',%d',1,length(delays)-1)],delays);
+
+ end
+ infile = tmpnam ();
+ outfile = tmpnam ();
+
+ %% Write data to file
+ save ('-ascii',infile, 'data');
+
+ %% Prepare format of the embedding vector
+ syscmd = sprintf("delay -M%d -m%d %s %s -o%s -V0 %s", ...
+ M, edim, flag.F, flag.D, outfile, infile);
+
+ %% Function call
+ system (syscmd);
+
+ delayed = load (outfile);
+
+endfunction
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <cde...@us...> - 2012-03-28 15:46:06
|
Revision: 10084
http://octave.svn.sourceforge.net/octave/?rev=10084&view=rev
Author: cdemills
Date: 2012-03-28 15:46:00 +0000 (Wed, 28 Mar 2012)
Log Message:
-----------
- new parameter, datefmt, and MUCH simpler date processing
Modified Paths:
--------------
trunk/octave-forge/extra/dataframe/inst/@dataframe/dataframe.m
Modified: trunk/octave-forge/extra/dataframe/inst/@dataframe/dataframe.m
===================================================================
--- trunk/octave-forge/extra/dataframe/inst/@dataframe/dataframe.m 2012-03-28 14:39:57 UTC (rev 10083)
+++ trunk/octave-forge/extra/dataframe/inst/@dataframe/dataframe.m 2012-03-28 15:46:00 UTC (rev 10084)
@@ -27,6 +27,7 @@
%# @item colnames : take the values as initialiser for column names
%# @item seeked : a (kept) field value which triggers start of processing.
%# @item trigger : a (unkept) field value which triggers start of processing.
+ %# @item datefmt: date format, see datestr help
%# Each preceeding line is silently skipped. Default: none
%# @item unquot: a logical switch telling wheter or not strings should
%# be unquoted before storage, default = true;
@@ -90,7 +91,7 @@
%# default values
seeked = []; trigger =[]; unquot = true; sep = "\t,"; cmt_lines = [];
-locales = "C";
+locales = "C"; datefmt = '';
if (length (varargin) > 0)
indi = 1;
@@ -145,6 +146,9 @@
case 'locales'
locales = varargin{indi + 1};
varargin(indi:indi+1) = [];
+ case 'datefmt'
+ datefmt = varargin{indi + 1};
+ varargin(indi:indi+1) = [];
otherwise %# FIXME: just skip it for now
disp (sprintf ("Ignoring unkown argument %s", varargin{indi}));
indi = indi + 1;
@@ -160,10 +164,10 @@
endif
indi = 0;
-while (indi <= size(varargin, 2))
+while (indi <= size (varargin, 2))
indi = indi + 1;
if (~isa (x, 'dataframe'))
- if (isa(x, 'char') && size(x, 1) < 2)
+ if (isa (x, 'char') && size (x, 1) < 2)
%# read the data frame from a file
try
dummy = tilde_expand (x);
@@ -171,7 +175,7 @@
df._src{end+1, 1} = dummy;
catch
%# try our own method
- UTF8_BOM = char([0xEF 0xBB 0xBF]);
+ UTF8_BOM = char ([0xEF 0xBB 0xBF]);
unwind_protect
dummy = tilde_expand (x);
fid = fopen (dummy);
@@ -250,7 +254,7 @@
endif
endwhile
endif
- x = cell (1+length (lines)-indl, size(dummy, 2));
+ x = cell (1+length (lines)-indl, size (dummy, 2));
empty_lines = []; cmt_lines = [];
while (indl <= length (lines))
dummy = content{indl};
@@ -268,9 +272,9 @@
endif
%# try to convert to float
- if (1)
- the_line = cellfun (@(x) sscanf (x, "%f", locales), dummy, \
- 'UniformOutput', false);
+ if (1)
+ the_line = cellfun (@(x) sscanf (x, "%f", locales), dummy, \
+ 'UniformOutput', false);
else
%# this code require a patch to src/file-io.cc in main
%# Octave tree
@@ -299,29 +303,22 @@
x(indj, indk) = regexp (dummy{indk}, '[^ ].*', 'match');
endif
else
- if (~isempty (regexp (dummy{indk}, '[/:-]')))
- %# try to convert to a date
- [timeval, nfields] = strptime( dummy{indk},
- [char(37) 'd/' char(37) 'm/' char(37) 'Y ' char(37) 'T']);
- if (nfields > 0) %# at least a few fields are OK
- keyboard
- timestr = strftime ([char(37) 'H:' char(37) 'M:' char(37) 'S'],
- timeval);
- %# try to extract the usec field, if any
- idx = regexp (dummy{indk}, timestr, 'end');
- if (~isempty (idx))
- idx = idx + 1;
- if (ispunct (dummy{indk}(idx)))
- idx = idx + 1;
- endif
- timeval.usec = str2num(dummy{indk}(idx:end));
- endif
+ if (~isempty (regexp (dummy{indk}, '[/:-]')) && ...
+ ~isempty (datefmt))
+
+ try
+ datetime = datevec (dummy{indk}, datefmt);
+ timeval = struct ("usec", 0, "sec", floor (datetime (6)),
+ "min", datetime(5), "hour", datetime(4),
+ "mday", datetime(3), "mon", datetime(2)-1,
+ "year", datetime(1)-1900);
+ timeval.usec = 1e6*(datetime(6)-timeval.sec);
x(indj, indk) = str2num (strftime ([char(37) 's'], timeval)) + ...
timeval.usec * 1e-6;
- else
+ catch
%# store it as is
x(indj, indk) = the_line{indk};
- endif
+ end_try_catch
else
x(indj, indk) = the_line{indk};
endif
@@ -349,10 +346,10 @@
indj = df._cnt(2)+(1:size (x, 2));
else
%# at this point, reading some filename failed
- error("dataframe: can't open '%s' for reading data", x);
+ error ("dataframe: can't open '%s' for reading data", x);
endif;
- if (iscell(x))
+ if (iscell (x))
if (2 == length (x))
%# use the intermediate value as destination column
[indc, ncol] = df_name2idx (df._name{2}, x{1}, df._cnt(2), "column");
@@ -367,7 +364,7 @@
end_try_catch
df = df_pad (df, 2, [length(dummy) indc], dummy);
x = x{2};
- indj = indc + (1:size(x, 2)); %# redefine target range
+ indj = indc + (1:size (x, 2)); %# redefine target range
else
if (isa (x{1}, 'cell'))
x = x{1}; %# remove one cell level
@@ -385,7 +382,7 @@
if (1 == length (df._name{2}) && length (df._name{2}) < \
length (indj))
[df._name{2}(indj, 1), df._over{2}(1, indj)] ...
- = df_colnames (char(df._name{2}), indj);
+ = df_colnames (char (df._name{2}), indj);
elseif (length (df._name{2}) < indj(1) || isempty (df._name{2}(indj)))
[df._name{2}(indj, 1), df._over{2}(1, indj)] ...
= df_colnames (inputname(indi), indj);
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <jpi...@us...> - 2012-03-28 14:40:08
|
Revision: 10083
http://octave.svn.sourceforge.net/octave/?rev=10083&view=rev
Author: jpicarbajal
Date: 2012-03-28 14:39:57 +0000 (Wed, 28 Mar 2012)
Log Message:
-----------
system-identification: testing the port procedure
Modified Paths:
--------------
trunk/octave-forge/main/system-identification/devel/tisean/src/delay.c
trunk/octave-forge/main/system-identification/devel/tisean/src/tisean_functions.cc
Modified: trunk/octave-forge/main/system-identification/devel/tisean/src/delay.c
===================================================================
--- trunk/octave-forge/main/system-identification/devel/tisean/src/delay.c 2012-03-28 13:57:33 UTC (rev 10082)
+++ trunk/octave-forge/main/system-identification/devel/tisean/src/delay.c 2012-03-28 14:39:57 UTC (rev 10083)
@@ -1,33 +1,108 @@
/*
- * This file is part of TISEAN
- *
- * Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
- *
- * TISEAN 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.
- *
- * TISEAN 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 TISEAN; if not, write to the Free Software
- * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
- */
-/*Author: Rainer Hegger. Last modified (rewritten in C) Aug 22, 2004*/
+Copyright (c) 2012 Juan Pablo Carbajal <car...@if...>
+
+ 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
+ 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/>.
+
+Delay vectors either from time series.
+Uses TISEAN 3.0.1 delay.c by courtesy of Rainer Hegger, Holger Kantz and Thomas Schreiber. Original-Author: Rainer Hegger. Last modified (rewritten in C) Aug 22, 2004
+<http://www.mpipks-dresden.mpg.de/~tisean/Tisean_3.0.1/>
+
+
+Author: Juan Pablo Carbajal <car...@if...>
+Created: Wednesday, March 28 2012
+Version: 0.1
+*/
+
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <limits.h>
#include <ctype.h>
#include "tsa.h"
+#include <octave/oct.h>
#define WID_STR "Produces delay vectors"
+DEFUN_DLD (delay, args, nargout,
+ "-*- texinfo -*-\n\
+Not documented.")
+{
+ int nargin = args.length ();
+ octave_value_list retval;
+ if (nargin != 3)
+ {
+ print_usage ();
+ }
+ else
+ {
+ // arguments in
+ // TODO Prepare input arguments
+ // arguments out
+ // TODO Prepare output arguments
+ // workspace
+ // TODO what is this?
+ // error indicators
+ // TODO Prepare error indicators
+
+ // Call TISEAN SLICOT routine
+ // TODO call function
+
+ //Check errors/warnings
+
+ // return values
+ // TODO Prepare return values
+ }
+
+ return retval;
+}
+
+void create_delay_list(void)
+{
+ unsigned int i=0,num=0;
+
+ while (multidelay[i]) {
+ if (!(isdigit(multidelay[i])) && !(multidelay[i] == ',')) {
+ fprintf(stderr,"Wrong format of -D parameter. Exiting!\n");
+ exit(DELAY_WRONG_FORMAT_D);
+ }
+ i++;
+ }
+
+ i=0;
+ while (multidelay[i]) {
+ if (multidelay[i++] == ',')
+ num++;
+ }
+
+ check_alloc(delaylist=(unsigned int*)malloc(sizeof(int)*(num+1)));
+ for (i=0;i<=num;i++) {
+ sscanf(multidelay,"%d",&delaylist[i]);
+ if (i<num) {
+ while ((*multidelay) != ',')
+ multidelay++;
+ }
+ multidelay++;
+ }
+
+ if ((num+1) != (embdim-indim)) {
+ fprintf(stderr,"Wrong number of delays. See man page. Exiting!\n");
+ exit(DELAY_WRONG_NUM_D);
+ }
+}
+
+/*
unsigned long length=ULONG_MAX;
unsigned long exclude=0;
unsigned int verbosity=0xff;
@@ -152,40 +227,7 @@
embdim=sum;
}
-void create_delay_list(void)
-{
- unsigned int i=0,num=0;
- while (multidelay[i]) {
- if (!(isdigit(multidelay[i])) && !(multidelay[i] == ',')) {
- fprintf(stderr,"Wrong format of -D parameter. Exiting!\n");
- exit(DELAY_WRONG_FORMAT_D);
- }
- i++;
- }
-
- i=0;
- while (multidelay[i]) {
- if (multidelay[i++] == ',')
- num++;
- }
-
- check_alloc(delaylist=(unsigned int*)malloc(sizeof(int)*(num+1)));
- for (i=0;i<=num;i++) {
- sscanf(multidelay,"%d",&delaylist[i]);
- if (i<num) {
- while ((*multidelay) != ',')
- multidelay++;
- }
- multidelay++;
- }
-
- if ((num+1) != (embdim-indim)) {
- fprintf(stderr,"Wrong number of delays. See man page. Exiting!\n");
- exit(DELAY_WRONG_NUM_D);
- }
-}
-
int main(int argc,char **argv)
{
char stin=0;
@@ -329,3 +371,4 @@
return 0;
}
+*/
Modified: trunk/octave-forge/main/system-identification/devel/tisean/src/tisean_functions.cc
===================================================================
--- trunk/octave-forge/main/system-identification/devel/tisean/src/tisean_functions.cc 2012-03-28 13:57:33 UTC (rev 10082)
+++ trunk/octave-forge/main/system-identification/devel/tisean/src/tisean_functions.cc 2012-03-28 14:39:57 UTC (rev 10083)
@@ -1 +1 @@
-#include "dealy.c" // produces delay vectors either from a scalar or from a multivariate time series.
+#include "delay.c" // produces delay vectors either from a scalar or from a multivariate time series.
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <jpi...@us...> - 2012-03-28 13:57:43
|
Revision: 10082
http://octave.svn.sourceforge.net/octave/?rev=10082&view=rev
Author: jpicarbajal
Date: 2012-03-28 13:57:33 +0000 (Wed, 28 Mar 2012)
Log Message:
-----------
system-identification: removing big collection and adding only under development
Modified Paths:
--------------
trunk/octave-forge/main/system-identification/devel/tisean/src/delay.c
Added Paths:
-----------
trunk/octave-forge/main/system-identification/devel/tisean/src/tisean_cec.h
trunk/octave-forge/main/system-identification/devel/tisean/src/tsa.h
Modified: trunk/octave-forge/main/system-identification/devel/tisean/src/delay.c
===================================================================
--- trunk/octave-forge/main/system-identification/devel/tisean/src/delay.c 2012-03-28 13:55:12 UTC (rev 10081)
+++ trunk/octave-forge/main/system-identification/devel/tisean/src/delay.c 2012-03-28 13:57:33 UTC (rev 10082)
@@ -23,7 +23,7 @@
#include <string.h>
#include <limits.h>
#include <ctype.h>
-#include "routines/tsa.h"
+#include "tsa.h"
#define WID_STR "Produces delay vectors"
Added: trunk/octave-forge/main/system-identification/devel/tisean/src/tisean_cec.h
===================================================================
--- trunk/octave-forge/main/system-identification/devel/tisean/src/tisean_cec.h (rev 0)
+++ trunk/octave-forge/main/system-identification/devel/tisean/src/tisean_cec.h 2012-03-28 13:57:33 UTC (rev 10082)
@@ -0,0 +1,85 @@
+/*
+ * This file is part of TISEAN
+ *
+ * Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
+ *
+ * TISEAN 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.
+ *
+ * TISEAN 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 TISEAN; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+/*Author: Rainer Hegger Last modified: May 26, 2000*/
+
+/* These definitions give the exit codes for the C part of the Tisean package.
+ Typically the name is build up of, first, the name of the routine creating
+ the exception, secondly, sort of an description of the exception.
+ */
+
+#ifndef _TISEAN_CEC_H
+#define _TISEAN_CEC_H
+
+/* These are the codes for the routines subtree */
+#define RESCALE_DATA_ZERO_INTERVAL 11
+#define CHECK_ALLOC_NOT_ENOUGH_MEMORY 12
+#define CHECK_OPTION_NOT_UNSIGNED 13
+#define CHECK_OPTION_NOT_INTEGER 14
+#define CHECK_OPTION_NOT_FLOAT 15
+#define CHECK_OPTION_NOT_TWO 16
+#define CHECK_OPTION_C_NO_VALUE 17
+#define TEST_OUTFILE_NO_WRITE_ACCESS 18
+#define SOLVELE_SINGULAR_MATRIX 19
+#define GET_SERIES_NO_LINES 20
+#define GET_MULTI_SERIES_WRONG_TYPE_OF_C 21
+#define GET_MULTI_SERIES_NO_LINES 22
+#define VARIANCE_VAR_EQ_ZERO 23
+#define EIG2_TOO_MANY_ITERATIONS 24
+#define CHECK_OPTION_NOT_THREE 25
+
+/* These are the codes for the main routines */
+#define LYAP_SPEC_NOT_ENOUGH_NEIGHBORS 50
+#define LYAP_SPEC_DATA_TOO_SHORT 51
+#define AR_MODEL_TOO_MANY_POLES 52
+#define EXTREMA_STRANGE_COMPONENT 53
+#define FALSE_NEAREST_NOT_ENOUGH_POINTS 54
+#define FSLE__TOO_LARGE_MINEPS 55
+#define GHKSS__TOO_MANY_NEIGHBORS 56
+#define NSTAT_Z__INVALID_STRING_FOR_OPTION 57
+#define NSTAT_Z__NOT_UNSIGNED_FOR_OPTION 58
+#define NSTAT_Z__TOO_LARGE_FOR_OPTION 59
+#define NSTAT_Z__OPTION_NOT_SET 60
+#define NSTAT_Z__TOO_MANY_PIECES 61
+#define NSTEP__ESCAPE_REGION 62
+#define POINCARE__WRONG_COMPONENT 63
+#define POINCARE__OUTSIDE_REGION 64
+#define POLYBACK__WRONG_PARAMETER_FILE 65
+#define POLYNOMP__WRONG_PARAMETER_FILE 66
+#define RESCALE__WRONG_INTERVAL 67
+#define SAV_GOL__UNDERDETERMINED 68
+#define SAV_GOL__TOO_LARGE_DERIVATIVE 69
+#define MAKENOISE__FLAGS_REQUIRED 70
+#define ZEROTH__STEP_TOO_LARGE 71
+#define LYAP_K__MAXITER_TOO_LARGE 72
+#define DELAY_WRONG_FORMAT_F 73
+#define DELAY_DIM_NOT_EQUAL_F_M 74
+#define DELAY_DIM_NOT_EQUAL_F_m 75
+#define DELAY_WRONG_FORMAT_D 76
+#define DELAY_WRONG_NUM_D 77
+#define DELAY_INCONS_d_D 78
+#define DELAY_SMALL_ZERO 79
+#define DELAY_INCONS_m_M 80
+#define ONESTEP_TOO_FEW_POINTS 81
+#define MEM_SPEC_TOO_MANY_POLES 82
+
+/* Global stuff */
+#define VECTOR_TOO_LARGE_FOR_LENGTH 100
+
+#endif
Added: trunk/octave-forge/main/system-identification/devel/tisean/src/tsa.h
===================================================================
--- trunk/octave-forge/main/system-identification/devel/tisean/src/tsa.h (rev 0)
+++ trunk/octave-forge/main/system-identification/devel/tisean/src/tsa.h 2012-03-28 13:57:33 UTC (rev 10082)
@@ -0,0 +1,105 @@
+/*
+ * This file is part of TISEAN
+ *
+ * Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
+ *
+ * TISEAN 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.
+ *
+ * TISEAN 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 TISEAN; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+/*Author: Rainer Hegger Last modified: Sep 3, 1999 */
+
+#ifndef _TSA_ROUTINES_H
+#define _TSA_ROUTINES_H
+
+#ifndef _TISEAN_CEC_H
+#include "tisean_cec.h"
+#endif
+
+/* size of the string which reads the input data
+ if your lines are longer than some 500 reals, increase the value
+ */
+#define INPUT_SIZE 1024
+
+/* The possible names of the verbosity levels */
+#define VER_INPUT 0x1
+#define VER_USR1 0x2
+#define VER_USR2 0x4
+#define VER_USR3 0x8
+#define VER_USR4 0x10
+#define VER_USR5 0x20
+#define VER_USR6 0x40
+#define VER_FIRST_LINE 0x80
+
+/* Uncomment the variable to get rid of the initial Version message */
+/*#define OMIT_WHAT_I_DO*/
+
+#define sqr(x) ((x)*(x))
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+extern int scan_help(int,char**);
+extern double *get_series(char *,unsigned long *,unsigned long,
+ unsigned int,unsigned int);
+extern double **get_multi_series(char *,unsigned long *,unsigned long,
+ unsigned int *,char *,char,unsigned int);
+extern void rescale_data(double *,unsigned long,double *,double *);
+extern void variance(double *,unsigned long,double *,double *);
+extern void make_box(double *,long **,long *,unsigned long,
+ unsigned int,unsigned int,unsigned int,double);
+extern unsigned long find_neighbors(double *,long **,long *,double *,
+ unsigned long,unsigned int,unsigned int,
+ unsigned int,double,unsigned long *);
+extern char* search_datafile(int, char**,unsigned int*,unsigned int);
+extern char* check_option(char**,int,int,int);
+extern void solvele(double**,double *,unsigned int);
+extern void test_outfile(char*);
+extern double** invert_matrix(double**,unsigned int);
+extern unsigned long exclude_interval(unsigned long,long,long,
+ unsigned long*,unsigned long*);
+extern void make_multi_box(double **,long **,long *,unsigned long,
+ unsigned int,unsigned int,unsigned int,
+ unsigned int,double);
+ /*only used for nrlazy. Will be removed with nrlazy */
+extern void make_multi_box2(double **,long **,long *,unsigned long,
+ unsigned int,unsigned int,unsigned int,
+ unsigned int,double);
+extern unsigned long find_multi_neighbors(double **,long **,long *,double **,
+ unsigned long,unsigned int,
+ unsigned int,unsigned int,
+ unsigned int,double,unsigned long *);
+extern unsigned int** make_multi_index(unsigned int,unsigned int,unsigned int);
+
+extern void check_alloc(void *);
+extern char* myfgets(char *,int *,FILE *,unsigned int);
+extern void what_i_do(char *, char *);
+extern double* rand_arb_dist(double *,unsigned long,unsigned long,
+ unsigned int,unsigned long);
+
+/* routines from rand.c */
+extern void rnd_init(unsigned long);
+extern unsigned long rnd_long();
+extern unsigned long rnd_1279();
+extern unsigned long rnd69069();
+extern double gaussian(double);
+
+/* routines from eigen.c */
+extern void eigen(double**,unsigned long,double*);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <jpi...@us...> - 2012-03-28 13:55:23
|
Revision: 10081
http://octave.svn.sourceforge.net/octave/?rev=10081&view=rev
Author: jpicarbajal
Date: 2012-03-28 13:55:12 +0000 (Wed, 28 Mar 2012)
Log Message:
-----------
system-identification: removing big collection and adding only under development
Added Paths:
-----------
trunk/octave-forge/main/system-identification/devel/tisean/src/
trunk/octave-forge/main/system-identification/devel/tisean/src/delay.c
trunk/octave-forge/main/system-identification/devel/tisean/src/tisean_functions.cc
Removed Paths:
-------------
trunk/octave-forge/main/system-identification/devel/tisean/source_c/
trunk/octave-forge/main/system-identification/devel/tisean/source_f/
Added: trunk/octave-forge/main/system-identification/devel/tisean/src/delay.c
===================================================================
--- trunk/octave-forge/main/system-identification/devel/tisean/src/delay.c (rev 0)
+++ trunk/octave-forge/main/system-identification/devel/tisean/src/delay.c 2012-03-28 13:55:12 UTC (rev 10081)
@@ -0,0 +1,331 @@
+/*
+ * This file is part of TISEAN
+ *
+ * Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
+ *
+ * TISEAN 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.
+ *
+ * TISEAN 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 TISEAN; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+/*Author: Rainer Hegger. Last modified (rewritten in C) Aug 22, 2004*/
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <limits.h>
+#include <ctype.h>
+#include "routines/tsa.h"
+
+#define WID_STR "Produces delay vectors"
+
+
+unsigned long length=ULONG_MAX;
+unsigned long exclude=0;
+unsigned int verbosity=0xff;
+int delay=1;
+unsigned int indim=1,embdim=2;
+char *column=NULL,*format=NULL,*multidelay=NULL;
+char *outfile=NULL;
+char *infile=NULL;
+char dimset=0,formatset=0,embset=0,mdelayset=0,delayset=0;
+char stdo=1;
+
+double **series;
+unsigned int *formatlist,*delaylist;
+
+void show_options(char *progname)
+{
+ what_i_do(progname,WID_STR);
+ fprintf(stderr,"\nUsage: %s [options]\n",progname);
+ fprintf(stderr,"Options:\n");
+ fprintf(stderr,"Everything not being a valid option will be interpreted as a"
+ " possible datafile.\nIf no datafile is given stdin is read."
+ " Just - also means stdin\n");
+ fprintf(stderr,"\t-l # of data [default: whole file]\n");
+ fprintf(stderr,"\t-x # of rows to ignore [default: 0]\n");
+ fprintf(stderr,"\t-M num. of columns to read [default: %u]\n",indim);
+ fprintf(stderr,"\t-c columns to read [default: 1,...,M]\n");
+ fprintf(stderr,"\t-m dimension [default: 2]\n");
+ fprintf(stderr,"\t-F format of the delay vector (see man page)\n");
+ fprintf(stderr,"\t-d delay [default: 1]\n");
+ fprintf(stderr,"\t-D multi delay list (see man page)\n");
+ fprintf(stderr,"\t-V verbosity level [default: 1]\n\t\t"
+ "0='only panic messages'\n\t\t"
+ "1='+ input/output messages'\n");
+ fprintf(stderr,"\t-o output file [default: 'datafile'.del, "
+ "without -o: stdout]\n");
+ fprintf(stderr,"\t-h show these options\n");
+ exit(0);
+}
+
+void scan_options(int n,char **str)
+{
+ char *out;
+
+ if ((out=check_option(str,n,'l','u')) != NULL)
+ sscanf(out,"%lu",&length);
+ if ((out=check_option(str,n,'x','u')) != NULL)
+ sscanf(out,"%lu",&exclude);
+ if ((out=check_option(str,n,'c','s')) != NULL)
+ column=out;
+ if ((out=check_option(str,n,'M','u')) != NULL) {
+ sscanf(out,"%u",&indim);
+ dimset=1;
+ }
+ if ((out=check_option(str,n,'F','s')) != NULL) {
+ format=out;
+ formatset=1;
+ }
+ if ((out=check_option(str,n,'m','u')) != NULL) {
+ sscanf(out,"%u",&embdim);
+ embset=1;
+ }
+ if ((out=check_option(str,n,'d','u')) != NULL) {
+ sscanf(out,"%u",&delay);
+ delayset=1;
+ }
+ if ((out=check_option(str,n,'D','s')) != NULL) {
+ multidelay=out;
+ mdelayset=1;
+ }
+ if ((out=check_option(str,n,'V','u')) != NULL)
+ sscanf(out,"%u",&verbosity);
+ if ((out=check_option(str,n,'o','o')) != NULL) {
+ stdo=0;
+ if (strlen(out) > 0)
+ outfile=out;
+ }
+}
+
+void create_format_list(void)
+{
+ unsigned int i=0,num=0,sum=0;
+
+ while (format[i]) {
+ if (!(isdigit(format[i])) && !(format[i] == ',')) {
+ fprintf(stderr,"Wrong format of -F parameter. Exiting!\n");
+ exit(DELAY_WRONG_FORMAT_F);
+ }
+ i++;
+ }
+
+ i=0;
+ while (format[i]) {
+ if (format[i++] == ',')
+ num++;
+ }
+
+ check_alloc(formatlist=(unsigned int*)malloc(sizeof(int)*(num+1)));
+ for (i=0;i<=num;i++) {
+ sscanf(format,"%d",&formatlist[i]);
+ if (i<num) {
+ while ((*format) != ',')
+ format++;
+ }
+ format++;
+ }
+
+ if (dimset && ((num+1) != indim)) {
+ fprintf(stderr,"Number of dimensions in -F is not equal to -M. Exiting!\n");
+ exit(DELAY_DIM_NOT_EQUAL_F_M);
+ }
+
+ for (i=0;i<=num;i++)
+ sum += formatlist[i];
+ if (embset && (sum != embdim)) {
+ fprintf(stderr,"The dimensions given in -m and -F are not equal!"
+ " Exiting\n");
+ exit(DELAY_DIM_NOT_EQUAL_F_m);
+ }
+ if (!dimset)
+ indim=num+1;
+ if (!embset)
+ embdim=sum;
+}
+
+void create_delay_list(void)
+{
+ unsigned int i=0,num=0;
+
+ while (multidelay[i]) {
+ if (!(isdigit(multidelay[i])) && !(multidelay[i] == ',')) {
+ fprintf(stderr,"Wrong format of -D parameter. Exiting!\n");
+ exit(DELAY_WRONG_FORMAT_D);
+ }
+ i++;
+ }
+
+ i=0;
+ while (multidelay[i]) {
+ if (multidelay[i++] == ',')
+ num++;
+ }
+
+ check_alloc(delaylist=(unsigned int*)malloc(sizeof(int)*(num+1)));
+ for (i=0;i<=num;i++) {
+ sscanf(multidelay,"%d",&delaylist[i]);
+ if (i<num) {
+ while ((*multidelay) != ',')
+ multidelay++;
+ }
+ multidelay++;
+ }
+
+ if ((num+1) != (embdim-indim)) {
+ fprintf(stderr,"Wrong number of delays. See man page. Exiting!\n");
+ exit(DELAY_WRONG_NUM_D);
+ }
+}
+
+int main(int argc,char **argv)
+{
+ char stin=0;
+ unsigned long i;
+ int j,k;
+ unsigned int alldim,maxemb,emb,rundel,delsum,runmdel;
+ unsigned int *inddelay;
+ FILE *fout;
+
+ if (scan_help(argc,argv))
+ show_options(argv[0]);
+
+ scan_options(argc,argv);
+#ifndef OMIT_WHAT_I_DO
+ if (verbosity&VER_INPUT)
+ what_i_do(argv[0],WID_STR);
+#endif
+
+
+ infile=search_datafile(argc,argv,NULL,verbosity);
+ if (infile == NULL)
+ stin=1;
+
+ if (outfile == NULL) {
+ if (!stin) {
+ check_alloc(outfile=(char*)calloc(strlen(infile)+5,1));
+ strcpy(outfile,infile);
+ strcat(outfile,".del");
+ }
+ else {
+ check_alloc(outfile=(char*)calloc(10,1));
+ strcpy(outfile,"stdin.del");
+ }
+ }
+ if (!stdo)
+ test_outfile(outfile);
+
+ if (delayset && mdelayset) {
+ fprintf(stderr,"-d and -D can't be used simultaneously. Exiting!\n");
+ exit(DELAY_INCONS_d_D);
+ }
+
+ if (delay < 1) {
+ fprintf(stderr,"Delay has to be larger than 0. Exiting!\n");
+ exit(DELAY_SMALL_ZERO);
+ }
+
+ if (!formatset && (embdim%indim)) {
+ fprintf(stderr,"Inconsistent -m and -M. Please set -F\n");
+ exit(DELAY_INCONS_m_M);
+ }
+ if (formatset) {
+ create_format_list();
+ }
+ else {
+ check_alloc(formatlist=(unsigned int*)malloc(sizeof(int)*indim));
+ for (i=0;i<indim;i++) {
+ formatlist[i]=embdim/indim;
+ }
+ }
+
+ alldim=0;
+ for (i=0;i<indim;i++)
+ alldim += formatlist[i];
+
+ if (mdelayset) {
+ create_delay_list();
+ }
+
+ check_alloc(inddelay=(unsigned int*)malloc(sizeof(int)*alldim));
+
+ rundel=0;
+ if (!mdelayset) {
+ for (i=0;i<indim;i++) {
+ delsum=0;
+ inddelay[rundel++]=delsum;
+ for (j=1;j<formatlist[i];j++) {
+ delsum += delay;
+ inddelay[rundel++]=delsum;
+ }
+ }
+ }
+ else {
+ runmdel=0;
+ for (i=0;i<indim;i++) {
+ delsum=0;
+ inddelay[rundel++]=delsum;
+ for (j=1;j<formatlist[i];j++) {
+ delsum += delaylist[runmdel++];
+ inddelay[rundel++]=delsum;
+ }
+ }
+ }
+
+ maxemb=0;
+ for (i=0;i<alldim;i++)
+ maxemb=(maxemb<inddelay[i])?inddelay[i]:maxemb;
+
+ if (column == NULL) {
+ series=get_multi_series(infile,&length,exclude,&indim,"",dimset,verbosity);
+ }
+ else {
+ series=get_multi_series(infile,&length,exclude,&indim,column,dimset,
+ verbosity);
+ }
+
+ if (stdo) {
+ if (verbosity)
+ fprintf(stderr,"Writing to stdout\n");
+ for (i=maxemb;i<length;i++) {
+ rundel=0;
+ for (j=0;j<indim;j++) {
+ emb=formatlist[j];
+ for (k=0;k<emb;k++)
+ fprintf(stdout,"%e ",series[j][i-inddelay[rundel++]]);
+ }
+ fprintf(stdout,"\n");
+ }
+ }
+ else {
+ fout=fopen(outfile,"w");
+ if (verbosity)
+ fprintf(stderr,"Opened %s for writing\n",outfile);
+ for (i=maxemb;i<length;i++) {
+ for (j=0;j<indim;j++) {
+ rundel=0;
+ emb=formatlist[j];
+ for (k=0;k<emb;k++)
+ fprintf(fout,"%e ",series[j][i-inddelay[rundel++]]);
+ }
+ fprintf(fout,"\n");
+ }
+ fclose(fout);
+ }
+
+ if (formatlist != NULL)
+ free(formatlist);
+ if (delaylist != NULL)
+ free(delaylist);
+ free(inddelay);
+
+ return 0;
+}
Added: trunk/octave-forge/main/system-identification/devel/tisean/src/tisean_functions.cc
===================================================================
--- trunk/octave-forge/main/system-identification/devel/tisean/src/tisean_functions.cc (rev 0)
+++ trunk/octave-forge/main/system-identification/devel/tisean/src/tisean_functions.cc 2012-03-28 13:55:12 UTC (rev 10081)
@@ -0,0 +1 @@
+#include "dealy.c" // produces delay vectors either from a scalar or from a multivariate time series.
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <jpi...@us...> - 2012-03-28 13:32:52
|
Revision: 10080
http://octave.svn.sourceforge.net/octave/?rev=10080&view=rev
Author: jpicarbajal
Date: 2012-03-28 13:32:37 +0000 (Wed, 28 Mar 2012)
Log Message:
-----------
system-identitifaction: Adding devel TISEAN files
Added Paths:
-----------
trunk/octave-forge/main/system-identification/devel/tisean/
trunk/octave-forge/main/system-identification/devel/tisean/source_c/
trunk/octave-forge/main/system-identification/devel/tisean/source_c/Makefile.in
trunk/octave-forge/main/system-identification/devel/tisean/source_c/ar-model.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/arima-model.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/av-d2.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/boxcount.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/corr.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/d2.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/delay.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/extrema.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/false_nearest.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/fsle.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/ghkss.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/histogram.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/lfo-ar.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/lfo-run.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/lfo-test.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/low121.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/lyap_k.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/lyap_r.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/lyap_spec.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/lzo-gm.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/lzo-run.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/lzo-test.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/makenoise.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/mem_spec.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/mutual.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/new.tgz
trunk/octave-forge/main/system-identification/devel/tisean/source_c/nrlazy.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/nstat_z.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/pca.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/poincare.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/polyback.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/polynom.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/polynomp.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/polypar.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/rbf.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/recurr.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/resample.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/rescale.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/routines/
trunk/octave-forge/main/system-identification/devel/tisean/source_c/routines/Makefile.in
trunk/octave-forge/main/system-identification/devel/tisean/source_c/routines/arima.tgz
trunk/octave-forge/main/system-identification/devel/tisean/source_c/routines/check_alloc.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/routines/check_option.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/routines/diffc.log
trunk/octave-forge/main/system-identification/devel/tisean/source_c/routines/diffh.log
trunk/octave-forge/main/system-identification/devel/tisean/source_c/routines/eigen.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/routines/exclude_interval.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/routines/find_multi_neighbors.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/routines/find_neighbors.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/routines/get_multi_series.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/routines/get_series.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/routines/invert_matrix.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/routines/make_box.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/routines/make_multi_box.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/routines/make_multi_box2.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/routines/make_multi_index.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/routines/myfgets.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/routines/rand.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/routines/rand_arb_dist.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/routines/rescale_data.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/routines/scan_help.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/routines/search_datafile.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/routines/solvele.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/routines/test_outfile.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/routines/tisean_cec.h
trunk/octave-forge/main/system-identification/devel/tisean/source_c/routines/tsa.h
trunk/octave-forge/main/system-identification/devel/tisean/source_c/routines/variance.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/routines/what_i_do.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/sav_gol.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/xcor.c
trunk/octave-forge/main/system-identification/devel/tisean/source_c/xzero.c
trunk/octave-forge/main/system-identification/devel/tisean/source_f/
trunk/octave-forge/main/system-identification/devel/tisean/source_f/Makefile.in
trunk/octave-forge/main/system-identification/devel/tisean/source_f/addnoise.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/any_s.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/ar-run.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/arguments.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/autocor.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/c1.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/c2d.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/c2g.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/c2naive.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/c2t.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/choose.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/cluster.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/commandline.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/compare.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/d1.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/endtoend.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/events.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/gpl.txt
trunk/octave-forge/main/system-identification/devel/tisean/source_f/help.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/henon.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/ikeda.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/intervals.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/istdio_temp.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/lazy.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/lorenz.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/neigh.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/nmore.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/normal.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/notch.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/others/
trunk/octave-forge/main/system-identification/devel/tisean/source_f/pc.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/predict.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/project.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/randomize/
trunk/octave-forge/main/system-identification/devel/tisean/source_f/randomize/Makefile.in
trunk/octave-forge/main/system-identification/devel/tisean/source_f/randomize/cool/
trunk/octave-forge/main/system-identification/devel/tisean/source_f/randomize/cool/exp.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/randomize/cost/
trunk/octave-forge/main/system-identification/devel/tisean/source_f/randomize/cost/auto.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/randomize/cost/autop.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/randomize/cost/spikeauto.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/randomize/cost/spikespec.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/randomize/cost/uneven.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/randomize/perm/
trunk/octave-forge/main/system-identification/devel/tisean/source_f/randomize/perm/event.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/randomize/perm/random.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/randomize/randomize.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/rank.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/readfile.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/rms.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/Makefile.in
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/chkder.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/d1mach.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/dqk15.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/enorm.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/fdjac3.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/fdump.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/i1mach.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/j4save.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/lmpar.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/pythag.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/qrfac.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/qrsolv.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/r1mach.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/radb2.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/radb3.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/radb4.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/radb5.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/radbg.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/radf2.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/radf3.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/radf4.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/radf5.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/radfg.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/rand.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/rfftb1.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/rfftf1.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/rffti1.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/rgauss.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/rs.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/rwupdt.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/snls1.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/tql2.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/tqlrat.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/tred1.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/tred2.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/xercnt.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/xerhlt.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/xermsg.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/xerprn.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/xersve.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/slatec/xgetua.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/spectrum.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/spikeauto.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/spikespec.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/store_spec.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/stp.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/surrogates.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/timerev.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/tospec.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/totospec.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/upo.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/upoembed.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/verbose.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/wiener1.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/wiener2.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/xc2.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/xreadfile.f
trunk/octave-forge/main/system-identification/devel/tisean/source_f/xrecur.f
Added: trunk/octave-forge/main/system-identification/devel/tisean/source_c/Makefile.in
===================================================================
--- trunk/octave-forge/main/system-identification/devel/tisean/source_c/Makefile.in (rev 0)
+++ trunk/octave-forge/main/system-identification/devel/tisean/source_c/Makefile.in 2012-03-28 13:32:37 UTC (rev 10080)
@@ -0,0 +1,47 @@
+SHELL = /bin/sh
+
+prefix = @prefix@
+exec_prefix = @exec_prefix@
+BINDIR = ${exec_prefix}/@bindir@
+
+CC = @CC@
+CFLAGS = @CFLAGS@
+AR = @AR@
+ARFLAGS = @ARFLAGS@
+INSTALL = @INSTALL@
+
+LOADLIBS = routines/libddtsa.a -lm
+
+# list of executables we want to produce
+ ALL = poincare extrema rescale recurr corr mutual false_nearest \
+ lyap_r lyap_k lyap_spec d2 av-d2 makenoise nrlazy low121 \
+ lzo-test lfo-run lfo-test rbf polynom polyback polynomp polypar \
+ ar-model mem_spec pca ghkss lfo-ar xzero xcor boxcount fsle \
+ resample histogram nstat_z sav_gol delay lzo-gm arima-model \
+ lzo-run
+
+all: $(ALL)
+
+routines/libddtsa.a:
+ (cd routines && $(MAKE))
+
+$(ALL): routines/libddtsa.a *.c
+ -$(CC) $(CFLAGS) $(COPTS) -o $@ $@.c $(LOADLIBS)
+
+install: all
+ -for bin in $(ALL); do $(INSTALL) $$bin $(BINDIR); done
+
+clean:
+ @rm -f *.o *~ #*#
+ @rm -f $(ALL)
+ -(cd routines && $(MAKE) clean)
+
+missing:
+ -@for bin in $(ALL); do \
+ test -z "`$$bin -h 2>&1 | grep Usage`" \
+ && echo $$bin "(Dresden C)" >> ../missing.log; \
+ $$bin -h 2>&1 | cat >> ../install.log; \
+ done; :
+
+uninstall:
+ -@for bin in $(ALL); do rm -f $(BINDIR)/$$bin; done
Added: trunk/octave-forge/main/system-identification/devel/tisean/source_c/ar-model.c
===================================================================
--- trunk/octave-forge/main/system-identification/devel/tisean/source_c/ar-model.c (rev 0)
+++ trunk/octave-forge/main/system-identification/devel/tisean/source_c/ar-model.c 2012-03-28 13:32:37 UTC (rev 10080)
@@ -0,0 +1,395 @@
+/*
+ * This file is part of TISEAN
+ *
+ * Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
+ *
+ * TISEAN 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.
+ *
+ * TISEAN 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 TISEAN; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+/*Author: Rainer Hegger*/
+/*Changes:
+ Jun 24, 2005: Output average error for multivariate data
+ Nov 25, 2005: Handle model order = 0
+ Jan 31, 2006: Add verbosity 4 to print data+residuals
+ */
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <limits.h>
+#include <math.h>
+#include "routines/tsa.h"
+
+#define WID_STR "Fits an multivariate AR model to the data and gives\
+ the coefficients\n\tand the residues (or an iterated model)"
+
+unsigned long length=ULONG_MAX,exclude=0;
+unsigned int dim=1,poles=1,ilength;
+unsigned int verbosity=1;
+char *outfile=NULL,*column=NULL,stdo=1,dimset=0,run_model=0;
+char *infile=NULL;
+double **series,*my_average;
+
+void show_options(char *progname)
+{
+ what_i_do(progname,WID_STR);
+ fprintf(stderr," Usage: %s [options]\n",progname);
+ fprintf(stderr," Options:\n");
+ fprintf(stderr,"Everything not being a valid option will be interpreted"
+ " as a possible"
+ " datafile.\nIf no datafile is given stdin is read. Just - also"
+ " means stdin\n");
+ fprintf(stderr,"\t-l length of file [default is whole file]\n");
+ fprintf(stderr,"\t-x # of lines to be ignored [default is 0]\n");
+ fprintf(stderr,"\t-m dimension [default is 1]\n");
+ fprintf(stderr,"\t-c columns to read [default is 1,...,dimension]\n");
+ fprintf(stderr,"\t-p #order of AR-Fit [default is 1]\n");
+ fprintf(stderr,"\t-s length of iterated model [default no iteration]\n");
+ fprintf(stderr,"\t-o output file name [default is 'datafile'.ar]\n");
+ fprintf(stderr,"\t-V verbosity level [default is 1]\n\t\t"
+ "0='only panic messages'\n\t\t"
+ "1='+ input/output messages'\n\t\t"
+ "2='+ print residuals though iterating a model'\n\t\t"
+ "4='+ print original data plus residuals'\n");
+ fprintf(stderr,"\t-h show these options\n\n");
+ exit(0);
+}
+
+void scan_options(int argc,char **argv)
+{
+ char *out;
+
+ if ((out=check_option(argv,argc,'p','u')) != NULL) {
+ sscanf(out,"%u",&poles);
+ if (poles < 1) {
+ fprintf(stderr,"The order should at least be one!\n");
+ exit(127);
+ }
+ }
+ if ((out=check_option(argv,argc,'l','u')) != NULL)
+ sscanf(out,"%lu",&length);
+ if ((out=check_option(argv,argc,'x','u')) != NULL)
+ sscanf(out,"%lu",&exclude);
+ if ((out=check_option(argv,argc,'m','u')) != NULL) {
+ sscanf(out,"%u",&dim);
+ dimset=1;
+ }
+ if ((out=check_option(argv,argc,'c','u')) != NULL)
+ column=out;
+ if ((out=check_option(argv,argc,'V','u')) != NULL)
+ sscanf(out,"%u",&verbosity);
+ if ((out=check_option(argv,argc,'s','u')) != NULL) {
+ sscanf(out,"%u",&ilength);
+ run_model=1;
+ }
+ if ((out=check_option(argv,argc,'o','o')) != NULL) {
+ stdo=0;
+ if (strlen(out) > 0)
+ outfile=out;
+ }
+}
+
+void set_averages_to_zero(void)
+{
+ double var;
+ long i,j;
+
+ for (i=0;i<dim;i++) {
+ variance(series[i],length,&my_average[i],&var);
+ for (j=0;j<length;j++)
+ series[i][j] -= my_average[i];
+ }
+}
+
+double** build_matrix(double **mat)
+{
+ long n,i1,j1,i2,j2,hi,hj;
+ double norm;
+
+ norm=1./((double)length-(double)poles);
+
+ for (i1=0;i1<dim;i1++)
+ for (i2=0;i2<poles;i2++) {
+ hi=i1*poles+i2;
+ for (j1=0;j1<dim;j1++)
+ for (j2=0;j2<poles;j2++) {
+ hj=j1*poles+j2;
+ mat[hi][hj]=0.0;
+ for (n=poles-1;n<length-1;n++)
+ mat[hi][hj] += series[i1][n-i2]*series[j1][n-j2];
+ mat[hi][hj] *= norm;
+ }
+ }
+
+ return invert_matrix(mat,(unsigned int)(dim*poles));
+}
+
+void build_vector(double *vec,long comp)
+{
+ long i1,i2,hi,n;
+ double norm;
+
+ norm=1./((double)length-(double)poles);
+
+ for (i1=0;i1<poles*dim;i1++)
+ vec[i1]=0.0;
+
+ for (i1=0;i1<dim;i1++)
+ for (i2=0;i2<poles;i2++) {
+ hi=i1*poles+i2;
+ for (n=poles-1;n<length-1;n++)
+ vec[hi] += series[comp][n+1]*series[i1][n-i2];
+ vec[hi] *= norm;
+ }
+}
+
+double* multiply_matrix_vector(double **mat,double *vec)
+{
+ long i,j;
+ double *new_vec;
+
+ check_alloc(new_vec=(double*)malloc(sizeof(double)*poles*dim));
+
+ for (i=0;i<poles*dim;i++) {
+ new_vec[i]=0.0;
+ for (j=0;j<poles*dim;j++)
+ new_vec[i] += mat[i][j]*vec[j];
+ }
+ return new_vec;
+}
+
+double* make_residuals(double **diff,double **coeff)
+{
+ long n,d,i,j;
+ double *resi;
+
+ check_alloc(resi=(double*)malloc(sizeof(double)*dim));
+ for (i=0;i<dim;i++)
+ resi[i]=0.0;
+
+ for (n=poles-1;n<length-1;n++) {
+ for (d=0;d<dim;d++) {
+ diff[d][n+1]=series[d][n+1];
+ for (i=0;i<dim;i++)
+ for (j=0;j<poles;j++)
+ diff[d][n+1] -= coeff[d][i*poles+j]*series[i][n-j];
+ resi[d] += sqr(diff[d][n+1]);
+ }
+ }
+ for (i=0;i<dim;i++)
+ resi[i]=sqrt(resi[i]/((double)length-(double)poles));
+
+ return resi;
+}
+
+void iterate_model(double **coeff,double *sigma,FILE *file)
+{
+ long i,j,i1,i2,n,d;
+ double **iterate,*swap;
+
+ check_alloc(iterate=(double**)malloc(sizeof(double*)*(poles+1)));
+ for (i=0;i<=poles;i++)
+ check_alloc(iterate[i]=(double*)malloc(sizeof(double)*dim));
+ rnd_init(0x44325);
+ for (i=0;i<1000;i++)
+ gaussian(1.0);
+ for (i=0;i<dim;i++)
+ for (j=0;j<poles;j++)
+ iterate[j][i]=gaussian(sigma[i]);
+
+ for (n=0;n<ilength;n++) {
+ for (d=0;d<dim;d++) {
+ iterate[poles][d]=gaussian(sigma[d]);
+ for (i1=0;i1<dim;i1++)
+ for (i2=0;i2<poles;i2++)
+ iterate[poles][d] += coeff[d][i1*poles+i2]*iterate[poles-1-i2][i1];
+ }
+ if (file != NULL) {
+ for (d=0;d<dim;d++)
+ fprintf(file,"%e ",iterate[poles][d]);
+ fprintf(file,"\n");
+ }
+ else {
+ for (d=0;d<dim;d++)
+ printf("%e ",iterate[poles][d]);
+ printf("\n");
+ }
+
+ swap=iterate[0];
+ for (i=0;i<poles;i++)
+ iterate[i]=iterate[i+1];
+ iterate[poles]=swap;
+ }
+
+ for (i=0;i<=poles;i++)
+ free(iterate[i]);
+ free(iterate);
+}
+
+int main(int argc,char **argv)
+{
+ char stdi=0;
+ double *pm;
+ long i,j;
+ FILE *file;
+ double **mat,**inverse,*vec,**coeff,**diff,avpm;
+
+ if (scan_help(argc,argv))
+ show_options(argv[0]);
+
+ scan_options(argc,argv);
+#ifndef OMIT_WHAT_I_DO
+ if (verbosity&VER_INPUT)
+ what_i_do(argv[0],WID_STR);
+#endif
+
+ infile=search_datafile(argc,argv,NULL,verbosity);
+ if (infile == NULL)
+ stdi=1;
+
+ if (outfile == NULL) {
+ if (!stdi) {
+ check_alloc(outfile=(char*)calloc(strlen(infile)+4,(size_t)1));
+ strcpy(outfile,infile);
+ strcat(outfile,".ar");
+ }
+ else {
+ check_alloc(outfile=(char*)calloc((size_t)9,(size_t)1));
+ strcpy(outfile,"stdin.ar");
+ }
+ }
+ if (!stdo)
+ test_outfile(outfile);
+
+ if (column == NULL)
+ series=(double**)get_multi_series(infile,&length,exclude,&dim,"",dimset,
+ verbosity);
+ else
+ series=(double**)get_multi_series(infile,&length,exclude,&dim,column,
+ dimset,verbosity);
+
+ check_alloc(my_average=(double*)malloc(sizeof(double)*dim));
+ set_averages_to_zero();
+
+ if (poles >= length) {
+ fprintf(stderr,"It makes no sense to have more poles than data! Exiting\n");
+ exit(AR_MODEL_TOO_MANY_POLES);
+ }
+
+
+ check_alloc(vec=(double*)malloc(sizeof(double)*poles*dim));
+ check_alloc(mat=(double**)malloc(sizeof(double*)*poles*dim));
+ for (i=0;i<poles*dim;i++)
+ check_alloc(mat[i]=(double*)malloc(sizeof(double)*poles*dim));
+
+ check_alloc(coeff=(double**)malloc(sizeof(double*)*dim));
+ inverse=build_matrix(mat);
+ for (i=0;i<dim;i++) {
+ build_vector(vec,i);
+ coeff[i]=multiply_matrix_vector(inverse,vec);
+ }
+
+ check_alloc(diff=(double**)malloc(sizeof(double*)*dim));
+ for (i=0;i<dim;i++)
+ check_alloc(diff[i]=(double*)malloc(sizeof(double)*length));
+
+ pm=make_residuals(diff,coeff);
+
+ if (stdo) {
+ avpm=pm[0]*pm[0];
+ for (i=1;i<dim;i++)
+ avpm += pm[i]*pm[i];
+ avpm=sqrt(avpm/dim);
+ printf("#average forcast error= %e\n",avpm);
+ printf("#individual forecast errors: ");
+ for (i=0;i<dim;i++)
+ printf("%e ",pm[i]);
+ printf("\n");
+ for (i=0;i<dim*poles;i++) {
+ printf("# ");
+ for (j=0;j<dim;j++)
+ printf("%e ",coeff[j][i]);
+ printf("\n");
+ }
+ if (!run_model || (verbosity&VER_USR1)) {
+ for (i=poles;i<length;i++) {
+ if (run_model)
+ printf("#");
+ for (j=0;j<dim;j++)
+ if (verbosity&VER_USR2)
+ printf("%e %e ",series[j][i]+my_average[j],diff[j][i]);
+ else
+ printf("%e ",diff[j][i]);
+ printf("\n");
+ }
+ }
+ if (run_model && (ilength > 0))
+ iterate_model(coeff,pm,NULL);
+ }
+ else {
+ file=fopen(outfile,"w");
+ if (verbosity&VER_INPUT)
+ fprintf(stderr,"Opened %s for output\n",outfile);
+ avpm=pm[0]*pm[0];
+ for (i=1;i<dim;i++)
+ avpm += pm[i]*pm[i];
+ avpm=sqrt(avpm/dim);
+ fprintf(file,"#average forcast error= %e\n",avpm);
+ fprintf(file,"#individual forecast errors: ");
+ for (i=0;i<dim;i++)
+ fprintf(file,"%e ",pm[i]);
+ fprintf(file,"\n");
+ for (i=0;i<dim*poles;i++) {
+ fprintf(file,"# ");
+ for (j=0;j<dim;j++)
+ fprintf(file,"%e ",coeff[j][i]);
+ fprintf(file,"\n");
+ }
+ if (!run_model || (verbosity&VER_USR1)) {
+ for (i=poles;i<length;i++) {
+ if (run_model)
+ fprintf(file,"#");
+ for (j=0;j<dim;j++)
+ if (verbosity&VER_USR2)
+ fprintf(file,"%e %e ",series[j][i]+my_average[j],diff[j][i]);
+ else
+ fprintf(file,"%e ",diff[j][i]);
+ fprintf(file,"\n");
+ }
+ }
+ if (run_model && (ilength > 0))
+ iterate_model(coeff,pm,file);
+ fclose(file);
+ }
+
+ if (outfile != NULL)
+ free(outfile);
+ if (infile != NULL)
+ free(infile);
+ free(vec);
+ for (i=0;i<poles*dim;i++) {
+ free(mat[i]);
+ free(inverse[i]);
+ }
+ free(mat);
+ free(inverse);
+ for (i=0;i<dim;i++) {
+ free(coeff[i]);
+ free(diff[i]);
+ }
+ free(coeff);
+ free(diff);
+ free(pm);
+
+ return 0;
+}
Added: trunk/octave-forge/main/system-identification/devel/tisean/source_c/arima-model.c
===================================================================
--- trunk/octave-forge/main/system-identification/devel/tisean/source_c/arima-model.c (rev 0)
+++ trunk/octave-forge/main/system-identification/devel/tisean/source_c/arima-model.c 2012-03-28 13:32:37 UTC (rev 10080)
@@ -0,0 +1,721 @@
+/*
+ * This file is part of TISEAN
+ *
+ * Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
+ *
+ * TISEAN 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.
+ *
+ * TISEAN 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 TISEAN; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+/*Author: Rainer Hegger, Last modified: Feb 6, 2006 */
+/*Changes:
+ Feb 4, 2006: First version
+ Feb 6, 2006: Find and remove bugs (1)
+ Feb 11, 2006: Add rand_arb_dist to iterate_***_model
+ */
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <limits.h>
+#include <math.h>
+#include "routines/tsa.h"
+
+#define WID_STR "Fits an multivariate ARIMA model to the data and gives\
+ the coefficients\n\tand the residues (or an iterated model)"
+
+unsigned long length=ULONG_MAX,exclude=0;
+unsigned int dim=1,poles=10,ilength,ITER=50;
+unsigned int arpoles=0,ipoles=0,mapoles=0,offset;
+unsigned int verbosity=1;
+char *outfile=NULL,*column=NULL,stdo=1,dimset=0,run_model=0,arimaset=0;
+char *infile=NULL;
+double **series,convergence=1.0e-3;
+
+double *my_average;
+unsigned long ardim,armadim;
+unsigned int **aindex;
+
+void show_options(char *progname)
+{
+ what_i_do(progname,WID_STR);
+ fprintf(stderr," Usage: %s [options]\n",progname);
+ fprintf(stderr," Options:\n");
+ fprintf(stderr,"Everything not being a valid option will be interpreted"
+ " as a possible"
+ " datafile.\nIf no datafile is given stdin is read. Just - also"
+ " means stdin\n");
+ fprintf(stderr,"\t-l length of file [default is whole file]\n");
+ fprintf(stderr,"\t-x # of lines to be ignored [default is 0]\n");
+ fprintf(stderr,"\t-m dimension [default is 1]\n");
+ fprintf(stderr,"\t-c columns to read [default is 1,...,dimension]\n");
+ fprintf(stderr,"\t-p order of initial AR-Fit [default is %u]\n",poles);
+ fprintf(stderr,"\t-P order of AR,I,MA-Fit [default is %u,%u,%u]\n",
+ arpoles,ipoles,mapoles);
+ fprintf(stderr,"\t-I # of arima iterations [default is %u]\n",ITER);
+ fprintf(stderr,"\t-e accuracy of convergence [default is %lf]\n",convergence);
+ fprintf(stderr,"\t-s length of iterated model [default no iteration]\n");
+ fprintf(stderr,"\t-o output file name [default is 'datafile'.ari]\n");
+ fprintf(stderr,"\t-V verbosity level [default is 1]\n\t\t"
+ "0='only panic messages'\n\t\t"
+ "1='+ input/output messages'\n\t\t"
+ "2='+ print residuals though iterating a model'\n\t\t"
+ "4='+ print original data plus residuals'\n");
+ fprintf(stderr,"\t-h show these options\n\n");
+ exit(0);
+}
+
+void scan_options(int argc,char **argv)
+{
+ char *out;
+
+ if ((out=check_option(argv,argc,'p','u')) != NULL) {
+ sscanf(out,"%u",&poles);
+ if (poles < 1) {
+ fprintf(stderr,"The order should at least be one!\n");
+ exit(127);
+ }
+ }
+ if ((out=check_option(argv,argc,'l','u')) != NULL)
+ sscanf(out,"%lu",&length);
+ if ((out=check_option(argv,argc,'x','u')) != NULL)
+ sscanf(out,"%lu",&exclude);
+ if ((out=check_option(argv,argc,'m','u')) != NULL) {
+ sscanf(out,"%u",&dim);
+ dimset=1;
+ }
+ if ((out=check_option(argv,argc,'P','3')) != NULL) {
+ sscanf(out,"%u,%u,%u",&arpoles,&ipoles,&mapoles);
+ if ((arpoles+ipoles+mapoles)>0)
+ arimaset=1;
+ }
+ if ((out=check_option(argv,argc,'I','u')) != NULL)
+ sscanf(out,"%u",&ITER);
+ if ((out=check_option(argv,argc,'e','f')) != NULL)
+ sscanf(out,"%lf",&convergence);
+ if ((out=check_option(argv,argc,'c','u')) != NULL)
+ column=out;
+ if ((out=check_option(argv,argc,'V','u')) != NULL)
+ sscanf(out,"%u",&verbosity);
+ if ((out=check_option(argv,argc,'s','u')) != NULL) {
+ sscanf(out,"%u",&ilength);
+ run_model=1;
+ }
+ if ((out=check_option(argv,argc,'o','o')) != NULL) {
+ stdo=0;
+ if (strlen(out) > 0)
+ outfile=out;
+ }
+}
+
+void make_difference(void)
+{
+ unsigned long i,d;
+
+ for (i=length-1;i>0;i--)
+ for (d=0;d<dim;d++)
+ series[d][i]=series[d][i]-series[d][i-1];
+}
+
+unsigned int** make_ar_index(void)
+{
+ unsigned int** ar_index;
+ unsigned long i;
+
+ check_alloc(ar_index=(unsigned int**)malloc(sizeof(unsigned int*)*2));
+ for (i=0;i<2;i++)
+ check_alloc(ar_index[i]=(unsigned int*)
+ malloc(sizeof(unsigned int)*ardim));
+ for (i=0;i<ardim;i++) {
+ ar_index[0][i]=i/poles;
+ ar_index[1][i]=i%poles;
+ }
+ return ar_index;
+}
+
+unsigned int** make_arima_index(unsigned int ars,unsigned int mas)
+{
+ unsigned int** arima_index;
+ unsigned int armad;
+ unsigned long i,i0;
+
+ armad=(ars+mas)*dim;
+ check_alloc(arima_index=(unsigned int**)malloc(sizeof(unsigned int*)*2));
+ for (i=0;i<2;i++)
+ check_alloc(arima_index[i]=(unsigned int*)
+ malloc(sizeof(unsigned int)*armad));
+ for (i=0;i<ars*dim;i++) {
+ arima_index[0][i]=i/ars;
+ arima_index[1][i]=i%ars;
+ }
+ i0=ars*dim;
+ for (i=0;i<mas*dim;i++) {
+ arima_index[0][i+i0]=dim+i/mas;
+ arima_index[1][i+i0]=i%mas;
+ }
+
+ return arima_index;
+}
+
+void set_averages_to_zero(void)
+{
+ double var;
+ long i,j;
+
+ for (i=0;i<dim;i++) {
+ variance(series[i],length,&my_average[i],&var);
+ for (j=0;j<length;j++)
+ series[i][j] -= my_average[i];
+ }
+}
+
+double** build_matrix(double **mat,unsigned int size)
+{
+ long n,i,j,is,id,js,jd;
+ double norm;
+
+ norm=1./((double)length-1.0-(double)poles-(double)offset);
+
+ for (i=0;i<size;i++) {
+ id=aindex[0][i];
+ is=aindex[1][i];
+ for (j=i;j<size;j++) {
+ jd=aindex[0][j];
+ js=aindex[1][j];
+ mat[i][j]=0.0;
+ for (n=offset+poles-1;n<length-1;n++)
+ mat[i][j] += series[id][n-is]*series[jd][n-js];
+ mat[i][j] *= norm;
+ mat[j][i]=mat[i][j];
+ }
+ }
+
+ return invert_matrix(mat,size);
+}
+
+void build_vector(double *vec,unsigned int size,long comp)
+{
+ long i,is,id,n;
+ double norm;
+
+ norm=1./((double)length-1.0-(double)poles-(double)offset);
+
+ for (i=0;i<size;i++) {
+ id=aindex[0][i];
+ is=aindex[1][i];
+ vec[i]=0.0;
+ for (n=offset+poles-1;n<length-1;n++)
+ vec[i] += series[comp][n+1]*series[id][n-is];
+ vec[i] *= norm;
+ }
+}
+
+double* multiply_matrix_vector(double **mat,double *vec,unsigned int size)
+{
+ long i,j;
+ double *new_vec;
+
+ check_alloc(new_vec=(double*)malloc(sizeof(double)*size));
+
+ for (i=0;i<size;i++) {
+ new_vec[i]=0.0;
+ for (j=0;j<size;j++)
+ new_vec[i] += mat[i][j]*vec[j];
+ }
+
+ return new_vec;
+}
+
+double* make_residuals(double **diff,double **coeff,unsigned int size)
+{
+ long n,n1,d,i,is,id;
+ double *resi;
+
+ check_alloc(resi=(double*)malloc(sizeof(double)*dim));
+ for (i=0;i<dim;i++)
+ resi[i]=0.0;
+
+ for (n=poles-1;n<length-1;n++) {
+ n1=n+1;
+ for (d=0;d<dim;d++) {
+ diff[d][n1]=series[d][n1];
+ for (i=0;i<size;i++) {
+ id=aindex[0][i];
+ is=aindex[1][i];
+ diff[d][n1] -= coeff[d][i]*series[id][n-is];
+ }
+ resi[d] += sqr(diff[d][n1]);
+ }
+ }
+
+ for (i=0;i<dim;i++)
+ resi[i]=sqrt(resi[i]/((double)length-(double)poles));
+
+ return resi;
+}
+
+void iterate_model(double **coeff,double *sigma,double **diff,FILE *file)
+{
+ long i,j,i1,i2,n,d;
+ double **iterate,*swap,**myrand;
+
+ check_alloc(iterate=(double**)malloc(sizeof(double*)*(poles+1)));
+ for (i=0;i<=poles;i++)
+ check_alloc(iterate[i]=(double*)malloc(sizeof(double)*dim));
+
+ check_alloc(myrand=(double**)malloc(sizeof(double*)*dim));
+ for (i=0;i<dim;i++)
+ myrand[i]=rand_arb_dist(diff[i],length,ilength+poles,100,0x44325);
+
+ rnd_init(0x44325);
+ for (i=0;i<1000;i++)
+ rnd_long();
+ for (i=0;i<dim;i++)
+ for (j=0;j<poles;j++)
+ iterate[j][i]=myrand[i][j];
+
+ for (n=0;n<ilength;n++) {
+ for (d=0;d<dim;d++) {
+ iterate[poles][d]=myrand[d][n+poles];
+ for (i1=0;i1<dim;i1++)
+ for (i2=0;i2<poles;i2++)
+ iterate[poles][d] += coeff[d][i1*poles+i2]*iterate[poles-1-i2][i1];
+ }
+ if (file != NULL) {
+ for (d=0;d<dim;d++)
+ fprintf(file,"%e ",iterate[poles][d]);
+ fprintf(file,"\n");
+ }
+ else {
+ for (d=0;d<dim;d++)
+ printf("%e ",iterate[poles][d]);
+ printf("\n");
+ }
+
+ swap=iterate[0];
+ for (i=0;i<poles;i++)
+ iterate[i]=iterate[i+1];
+ iterate[poles]=swap;
+ }
+
+ for (i=0;i<=poles;i++)
+ free(iterate[i]);
+ free(iterate);
+
+ for (i=0;i<dim;i++)
+ free(myrand[i]);
+ free(myrand);
+}
+
+void iterate_arima_model(double **coeff,double *sigma,double **diff,FILE *file)
+{
+ double **iterate,*swap,**myrand;
+ unsigned long i,j,n,is,id;
+
+ check_alloc(iterate=(double**)malloc(sizeof(double*)*(poles+1)));
+ for (i=0;i<=poles;i++)
+ check_alloc(iterate[i]=(double*)malloc(sizeof(double)*2*dim));
+
+ check_alloc(myrand=(double**)malloc(sizeof(double*)*dim));
+ for (i=0;i<dim;i++)
+ myrand[i]=rand_arb_dist(diff[i],length,ilength+poles,100,0x44325);
+
+ rnd_init(0x44325);
+ for (i=0;i<1000;i++)
+ rnd_long();
+ for (i=0;i<dim;i++)
+ for (j=0;j<poles;j++)
+ iterate[j][i]=iterate[j][dim+i]=myrand[i][j];
+
+ for (n=0;n<ilength;n++) {
+ for (i=0;i<dim;i++)
+ iterate[poles][i]=iterate[poles][i+dim]=myrand[i][n+poles];
+
+ for (j=0;j<dim;j++) {
+ for (i=0;i<armadim;i++) {
+ id=aindex[0][i];
+ is=aindex[1][i];
+ iterate[poles][j] += coeff[j][i]*iterate[poles-1-is][id];
+ }
+ }
+
+ if (file != NULL) {
+ for (i=0;i<dim;i++)
+ fprintf(file,"%e ",iterate[poles][i]);
+ fprintf(file,"\n");
+ }
+ else {
+ for (i=0;i<dim;i++)
+ printf("%e ",iterate[poles][i]);
+ printf("\n");
+ }
+
+ swap=iterate[0];
+ for (i=0;i<poles;i++)
+ iterate[i]=iterate[i+1];
+ iterate[poles]=swap;
+ }
+
+ for (i=0;i<=poles;i++)
+ free(iterate[i]);
+ free(iterate);
+ for (i=0;i<dim;i++)
+ free(myrand[i]);
+ free(myrand);
+}
+
+int main(int argc,char **argv)
+{
+ char stdi=0;
+ double *pm;
+ long i,j,iter,hj,realiter=0;
+ unsigned int size,is,id;
+ FILE *file;
+ double **mat,**inverse,*vec,**coeff,**diff,**hseries;
+ double **oldcoeff,*diffcoeff=NULL;
+ double hdiff,**xdiff=NULL,avpm;
+ double loglikelihood,aic,alldiff;
+
+ if (scan_help(argc,argv))
+ show_options(argv[0]);
+
+ scan_options(argc,argv);
+#ifndef OMIT_WHAT_I_DO
+ if (verbosity&VER_INPUT)
+ what_i_do(argv[0],WID_STR);
+#endif
+
+ infile=search_datafile(argc,argv,NULL,verbosity);
+ if (infile == NULL)
+ stdi=1;
+
+ if (outfile == NULL) {
+ if (!stdi) {
+ check_alloc(outfile=(char*)calloc(strlen(infile)+5,(size_t)1));
+ strcpy(outfile,infile);
+ strcat(outfile,".ari");
+ }
+ else {
+ check_alloc(outfile=(char*)calloc((size_t)10,(size_t)1));
+ strcpy(outfile,"stdin.ari");
+ }
+ }
+ if (!stdo)
+ test_outfile(outfile);
+
+ if (column == NULL)
+ series=(double**)get_multi_series(infile,&length,exclude,&dim,"",dimset,
+ verbosity);
+ else
+ series=(double**)get_multi_series(infile,&length,exclude,&dim,column,
+ dimset,verbosity);
+
+ check_alloc(my_average=(double*)malloc(sizeof(double)*dim));
+
+ for (i=0;i<ipoles;i++)
+ make_difference();
+
+ for (i=0;i<dim;i++)
+ series[i] += ipoles;
+ length -= ipoles;
+
+ set_averages_to_zero();
+
+ if (poles >= length) {
+ fprintf(stderr,"It makes no sense to have more poles than data! Exiting\n");
+ exit(AR_MODEL_TOO_MANY_POLES);
+ }
+ if (arimaset) {
+ if ((arpoles >= length) || (mapoles >= length)) {
+ fprintf(stderr,"It makes no sense to have more poles than data! Exiting\n");
+ exit(AR_MODEL_TOO_MANY_POLES);
+ }
+ }
+
+ ardim=poles*dim;
+ aindex=make_ar_index();
+
+ check_alloc(vec=(double*)malloc(sizeof(double)*ardim));
+ check_alloc(mat=(double**)malloc(sizeof(double*)*ardim));
+ for (i=0;i<ardim;i++)
+ check_alloc(mat[i]=(double*)malloc(sizeof(double)*ardim));
+
+ check_alloc(coeff=(double**)malloc(sizeof(double*)*dim));
+ inverse=build_matrix(mat,ardim);
+ for (i=0;i<dim;i++) {
+ build_vector(vec,ardim,i);
+ coeff[i]=multiply_matrix_vector(inverse,vec,ardim);
+ }
+
+ check_alloc(diff=(double**)malloc(sizeof(double*)*dim));
+ for (i=0;i<dim;i++)
+ check_alloc(diff[i]=(double*)malloc(sizeof(double)*length));
+
+ pm=make_residuals(diff,coeff,ardim);
+
+ free(vec);
+ for (i=0;i<ardim;i++) {
+ free(mat[i]);
+ free(inverse[i]);
+ }
+ free(mat);
+ free(inverse);
+ size=ardim;
+
+ if (arimaset) {
+ offset=poles;
+ for (i=0;i<2;i++)
+ free(aindex[i]);
+ free(aindex);
+
+ for (i=0;i<dim;i++)
+ free(coeff[i]);
+ free(coeff);
+ check_alloc(xdiff=(double**)malloc(sizeof(double*)*ITER));
+ for (i=0;i<ITER;i++)
+ check_alloc(xdiff[i]=(double*)malloc(sizeof(double)*dim));
+
+ armadim=(arpoles+mapoles)*dim;
+ aindex=make_arima_index(arpoles,mapoles);
+ size=armadim;
+
+ check_alloc(hseries=(double**)malloc(sizeof(double*)*2*dim));
+ for (i=0;i<dim;i++) {
+ check_alloc(hseries[i]=(double*)malloc(sizeof(double)*length));
+ check_alloc(hseries[i+dim]=(double*)malloc(sizeof(double)*length));
+ for (j=0;j<length;j++) {
+ hseries[i][j]=series[i][j];
+ hseries[i+dim][j]=diff[i][j];
+ }
+ }
+
+ for (i=0;i<dim;i++)
+ free(series[i]-ipoles);
+ free(series);
+
+ series=hseries;
+
+ check_alloc(oldcoeff=(double**)malloc(sizeof(double*)*dim));
+ for (i=0;i<dim;i++) {
+ check_alloc(oldcoeff[i]=(double*)malloc(sizeof(double)*armadim));
+ for (j=0;j<armadim;j++)
+ oldcoeff[i][j]=0.0;
+ }
+ check_alloc(diffcoeff=(double*)malloc(sizeof(double)*ITER));
+
+ for (iter=1;iter<=ITER;iter++) {
+ check_alloc(vec=(double*)malloc(sizeof(double)*armadim));
+ check_alloc(mat=(double**)malloc(sizeof(double*)*armadim));
+ for (i=0;i<armadim;i++)
+ check_alloc(mat[i]=(double*)malloc(sizeof(double)*armadim));
+
+ check_alloc(coeff=(double**)malloc(sizeof(double*)*dim));
+
+ poles=(arpoles > mapoles)? arpoles:mapoles;
+
+ offset += poles;
+ inverse=build_matrix(mat,armadim);
+
+ for (i=0;i<dim;i++) {
+ build_vector(vec,armadim,i);
+ coeff[i]=multiply_matrix_vector(inverse,vec,armadim);
+ }
+
+ pm=make_residuals(diff,coeff,armadim);
+
+ for (j=0;j<dim;j++) {
+ hdiff=0.0;
+ hj=j+dim;
+ for (i=offset;i<length;i++)
+ hdiff += sqr(series[hj][i]-diff[j][i]);
+ for (i=0;i<length;i++) {
+ series[hj][i]=diff[j][i];
+ }
+ xdiff[iter-1][j]=sqrt(hdiff/(double)(length-offset));
+ }
+
+ free(vec);
+ for (i=0;i<armadim;i++) {
+ free(mat[i]);
+ free(inverse[i]);
+ }
+ free(mat);
+ free(inverse);
+
+ diffcoeff[iter-1]=0.0;
+ for (i=0;i<dim;i++)
+ for (j=0;j<dim;j++) {
+ diffcoeff[iter-1] += sqr(coeff[i][j]-oldcoeff[i][j]);
+ oldcoeff[i][j]=coeff[i][j];
+ }
+ diffcoeff[iter-1]=sqrt(diffcoeff[iter-1]/(double)armadim);
+ alldiff=xdiff[iter-1][0];
+ for (i=1;i<dim;i++)
+ if (xdiff[iter-1][i] > alldiff)
+ alldiff=xdiff[iter-1][i];
+ realiter=iter;
+ if (alldiff < convergence)
+ iter=ITER;
+
+ if (iter < ITER) {
+ for (i=0;i<dim;i++)
+ free(coeff[i]);
+ free(coeff);
+ }
+ }
+ }
+
+ if (stdo) {
+ if (arimaset) {
+ printf("#convergence of residuals in arima fit\n");
+ for (i=0;i<realiter;i++) {
+ printf("#iteration %ld ",i+1);
+ for (j=0;j<dim;j++)
+ printf("%e ",xdiff[i][j]);
+ printf("%e",diffcoeff[i]);
+ printf("\n");
+ }
+ }
+ avpm=pm[0]*pm[0];
+ loglikelihood= -log(pm[0]);
+ for (i=1;i<dim;i++) {
+ avpm += pm[i]*pm[i];
+ loglikelihood -= log(pm[i]);
+ }
+ loglikelihood *= ((double)length);
+ loglikelihood += -((double)length)*
+ ((1.0+log(2.*M_PI))*dim)/2.0;
+ avpm=sqrt(avpm/dim);
+ printf("#average forcast error= %e\n",avpm);
+ printf("#individual forecast errors: ");
+ for (i=0;i<dim;i++)
+ printf("%e ",pm[i]);
+ printf("\n");
+ if (arimaset)
+ aic=2.0*(arpoles+mapoles)-2.0*loglikelihood;
+ else
+ aic=2.0*poles-2.0*loglikelihood;
+ printf("#Log-Likelihood= %e\t AIC= %e\n",loglikelihood,aic);
+ for (i=0;i<size;i++) {
+ id=aindex[0][i];
+ is=aindex[1][i];
+ if (id < dim)
+ printf("#x_%u(n-%u) ",id+1,is);
+ else
+ printf("#e_%u(n-%u) ",id+1-dim,is);
+ for (j=0;j<dim;j++)
+ printf("%e ",coeff[j][i]);
+ printf("\n");
+ }
+ if (!run_model || (verbosity&VER_USR1)) {
+ for (i=poles;i<length;i++) {
+ if (run_model)
+ printf("#");
+ for (j=0;j<dim;j++)
+ if (verbosity&VER_USR2)
+ printf("%e %e ",series[j][i]+my_average[j],diff[j][i]);
+ else
+ printf("%e ",diff[j][i]);
+ printf("\n");
+ }
+ }
+ if (run_model && (ilength > 0)) {
+ if (!arimaset)
+ iterate_model(coeff,pm,diff,NULL);
+ else
+ iterate_arima_model(coeff,pm,diff,NULL);
+ }
+ }
+ else {
+ file=fopen(outfile,"w");
+ if (verbosity&VER_INPUT)
+ fprintf(stderr,"Opened %s for output\n",outfile);
+ if (arimaset) {
+ fprintf(file,"#convergence of residuals in arima fit\n");
+ for (i=0;i<realiter;i++) {
+ fprintf(file,"#iteration %ld ",i+1);
+ for (j=0;j<dim;j++)
+ fprintf(file,"%e ",xdiff[i][j]);
+ fprintf(file,"%e",diffcoeff[i]);
+ fprintf(file,"\n");
+ }
+ }
+ avpm=pm[0]*pm[0];
+ loglikelihood= -log(pm[0]);
+ for (i=1;i<dim;i++) {
+ avpm += pm[i]*pm[i];
+ loglikelihood -= log(pm[i]);
+ }
+ loglikelihood *= ((double)length);
+ loglikelihood += -((double)length)*
+ ((1.0+log(2.*M_PI))*dim)/2.0;
+ avpm=sqrt(avpm/dim);
+ fprintf(file,"#average forcast error= %e\n",avpm);
+ fprintf(file,"#individual forecast errors: ");
+ for (i=0;i<dim;i++)
+ fprintf(file,"%e ",pm[i]);
+ fprintf(file,"\n");
+ if (arimaset)
+ aic=2.0*(arpoles+mapoles)-2.0*loglikelihood;
+ else
+ aic=2.0*poles-2.0*loglikelihood;
+ fprintf(file,"#Log-Likelihood= %e\t AIC= %e\n",loglikelihood,aic);
+ for (i=0;i<size;i++) {
+ id=aindex[0][i];
+ is=aindex[1][i];
+ if (id < dim)
+ fprintf(file,"#x_%u(n-%u) ",id+1,is);
+ else
+ fprintf(file,"#e_%u(n-%u) ",id+1-dim,is);
+ for (j=0;j<dim;j++)
+ fprintf(file,"%e ",coeff[j][i]);
+ fprintf(file,"\n");
+ }
+ if (!run_model || (verbosity&VER_USR1)) {
+ for (i=poles;i<length;i++) {
+ if (run_model)
+ fprintf(file,"#");
+ for (j=0;j<dim;j++)
+ if (verbosity&VER_USR2)
+ fprintf(file,"%e %e ",series[j][i]+my_average[j],diff[j][i]);
+ else
+ fprintf(file,"%e ",diff[j][i]);
+ fprintf(file,"\n");
+ }
+ }
+ if (run_model && (ilength > 0)) {
+ if (!arimaset)
+ iterate_model(coeff,pm,diff,file);
+ else
+ iterate_arima_model(coeff,pm,diff,file);
+ }
+ fclose(file);
+ }
+ if (outfile != NULL)
+ free(outfile);
+ if (infile != NULL)
+ free(infile);
+ for (i=0;i<dim;i++) {
+ free(coeff[i]);
+ free(diff[i]);
+ free(series[i]);
+ if (arimaset)
+ free(series[i+dim]);
+ }
+ free(coeff);
+ free(diff);
+ free(series);
+
+ free(pm);
+
+ return 0;
+}
Added: trunk/octave-forge/main/system-identification/devel/tisean/source_c/av-d2.c
===================================================================
--- trunk/octave-forge/main/system-identification/devel/tisean/source_c/av-d2.c (rev 0)
+++ trunk/octave-forge/main/system-identification/devel/tisean/source_c/av-d2.c 2012-03-28 13:32:37 UTC (rev 10080)
@@ -0,0 +1,188 @@
+/*
+ * This file is part of TISEAN
+ *
+ * Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
+ *
+ * TISEAN 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.
+ *
+ * TISEAN 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 TISEAN; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+/*Author: Rainer Hegger Last modified: Sep 3, 1999 */
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <limits.h>
+#include "routines/tsa.h"
+
+#define WID_STR "Smoothes the output of the d2 program"
+
+#define MAXLENGTH 1000
+
+unsigned int maxdim=UINT_MAX,mindim=1;
+unsigned int verbosity=0xff;
+int aver=1;
+char rescaled=0;
+char stout=1;
+char *outfile=NULL;
+char *infile=NULL;
+
+void show_options(char *progname)
+{
+ what_i_do(progname,WID_STR);
+ fprintf(stderr,"Usage: %s [options]\n",progname);
+ fprintf(stderr,"Options:\n");
+ fprintf(stderr,"Everything not being a valid option will be interpreted"
+ " as a possible datafile.\nStdin does NOT work.\n");
+ fprintf(stderr,"\t-m dimension to start with [Default: 1]\n");
+ fprintf(stderr,"\t-M dimension to end with [Default: whole file]\n");
+ fprintf(stderr,"\t-a n average over (2n+1) values [Default: 1]\n");
+ fprintf(stderr,"\t-E use rescaled data for the length scales\n\t\t"
+ "[Default: use units of data]\n");
+ fprintf(stderr,"\t-o name of output file [Default: stdout,\n\t\t"
+ "-o without value means 'datafile'.av]\n");
+ fprintf(stderr,"\t-V verbosity level [Default: 1]\n\t\t"
+ "0='only panic messages'\n\t\t"
+ "1='+ input/output messages'\n");
+ fprintf(stderr,"\t-h show these options\n");
+ exit(0);
+}
+
+void scan_options(int n,char **in)
+{
+ char *out;
+
+ if ((out=check_option(in,n,'m','u')) != NULL)
+ sscanf(out,"%u",&mindim);
+ if ((out=check_option(in,n,'M','u')) != NULL)
+ sscanf(out,"%u",&maxdim);
+ if ((out=check_option(in,n,'a','u')) != NULL)
+ sscanf(out,"%u",&aver);
+ if ((out=check_option(in,n,'V','u')) != NULL)
+ sscanf(out,"%u",&verbosity);
+ if ((out=check_option(in,n,'E','n')) != NULL)
+ rescaled=1;
+ if ((out=check_option(in,n,'o','o')) != NULL) {
+ stout=0;
+ if (strlen(out) > 0)
+ outfile=out;
+ }
+}
+
+int main(int argc,char **argv)
+{
+ char instr[1024];
+ char *form1="%lf%lf",*form2="%*lf%lf%lf";
+ char empty=0;
+ unsigned int howmany,size=1;
+ int j,k;
+ long dim;
+ double *eps,*y;
+ double avy,aveps,norm;
+ FILE *file,*fout=NULL;
+
+ if ((argc < 2) || scan_help(argc,argv))
+ show_options(argv[0]);
+
+ scan_options(argc,argv);
+#ifndef OMIT_WHAT_I_DO
+ if (verbosity&VER_INPUT)
+ what_i_do(argv[0],WID_STR);
+#endif
+
+ infile=search_datafile(argc,argv,0L,verbosity);
+ if (infile == NULL) {
+ fprintf(stderr,"You have to give a datafile. Exiting!\n");
+ exit(127);
+ }
+ if (outfile == NULL) {
+ check_alloc(outfile=(char*)calloc(strlen(infile)+4,(size_t)1));
+ sprintf(outfile,"%s.av",infile);
+ }
+
+ check_alloc(eps=(double*)malloc(sizeof(double)*MAXLENGTH));
+ check_alloc(y=(double*)malloc(sizeof(double)*MAXLENGTH));
+
+ file=fopen(infile,"r");
+
+ if (!stout) {
+ test_outfile(outfile);
+ fout=fopen(outfile,"w");
+ if (verbosity&VER_INPUT)
+ fprintf(stderr,"Opened %s for writing\n",outfile);
+ }
+ else {
+ if (verbosity&VER_INPUT)
+ fprintf(stderr,"Writing to stdout\n");
+ }
+
+ if (mindim > maxdim)
+ mindim=maxdim;
+ norm=2.0*aver+1.0;
+
+ while (fgets(instr,1024,file) != NULL) {
+ if (strlen(instr) != 1) {
+ if (instr[0] == '#') {
+ if (strstr(instr,"m= ") != NULL) {
+ sscanf(instr,"%*s %ld",&dim);
+ if ((dim >= mindim) && (dim <= maxdim)) {
+ howmany=0;
+ empty=0;
+ do {
+ if (fgets(instr,1024,file) == NULL)
+ exit(127);
+ if (strlen(instr) == 1)
+ empty=1;
+ if (!empty && (instr[0] != '#')) {
+ if (!rescaled)
+ sscanf(instr,form1,&eps[howmany],&y[howmany]);
+ else
+ sscanf(instr,form2,&y[howmany],&eps[howmany]);
+ howmany++;
+ if (!(howmany%MAXLENGTH)) {
+ check_alloc(realloc(eps,size*MAXLENGTH*sizeof(double)));
+ check_alloc(realloc(y,size*MAXLENGTH*sizeof(double)));
+ size++;
+ }
+ }
+ } while (!empty);
+ for (k=aver;k<howmany-aver;k++) {
+ avy=aveps=0.0;
+ for (j= -aver;j<=aver;j++) {
+ avy += y[k+j];
+ aveps += eps[k+j];
+ }
+ if (!stout)
+ fprintf(fout,"%e %e\n",aveps/norm,avy/norm);
+ else
+ fprintf(stdout,"%e %e\n",aveps/norm,avy/norm);
+ }
+ if (!stout)
+ fprintf(fout,"\n");
+ else
+ fprintf(stdout,"\n");
+ }
+ }
+ }
+ }
+ }
+
+ if (outfile != NULL)
+ free(outfile);
+ if (infile != NULL)
+ free(infile);
+ free(eps);
+ free(y);
+
+ return 0;
+}
Added: trunk/octave-forge/main/system-identification/devel/tisean/source_c/boxcount.c
===================================================================
--- trunk/octave-forge/main/system-identification/devel/tisean/source_c/boxcount.c (rev 0)
+++ trunk/octave-forge/main/system-identification/devel/tisean/source_c/boxcount.c 2012-03-28 13:32:37 UTC (rev 10080)
@@ -0,0 +1,369 @@
+/*
+ * This file is part of TISEAN
+ *
+ * Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
+ *
+ * TISEAN 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.
+ *
+ * TISEAN 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 TISEAN; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+/* Author: Rainer Hegger Last modified: Feb 22, 2006 */
+/* Changes:
+ 02/22/06: Remove this strange else in start_box that
+ did not compile anyways
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+#include <limits.h>
+#include "routines/tsa.h"
+
+#define WID_STR "Estimates the Renyi entropy of Qth order\n\t\
+using a partition instead of a covering."
+
+typedef struct {
+ double *hist;
+ void *ptr;
+} hliste;
+
+unsigned long LENGTH=ULONG_MAX,exclude=0;
+unsigned int maxembed=10,dimension=1,DELAY=1,EPSCOUNT=20;
+unsigned int verbosity=0xff;
+double Q=2.0,EPSMIN=1.e-3,EPSMAX=1.0;
+char dimset=0,epsminset=0,epsmaxset=0;
+char *outfile=NULL;
+char *column=NULL;
+
+int epsi;
+unsigned long length;
+double EPSFAKTOR;
+unsigned int **which_dims;
+double *histo;
+double **series;
+
+void show_options(char *progname)
+{
+ what_i_do(progname,WID_STR);
+ fprintf(stderr,"Usage: %s [Options]\n",progname);
+ fprintf(stderr,"Options:\n");
+ fprintf(stderr,"\t-l # of datapoints [Default: whole file]\n");
+ fprintf(stderr,"\t-x # of lines to ignore [Default: %lu]\n",exclude);
+ fprintf(stderr,"\t-M # of columns,maximal embedding dimension "
+ "[Default: %u,%u]\n",dimension,maxembed);
+ fprintf(stderr,"\t-c columns to read [Default: 1,...,#of compon.]\n");
+ fprintf(stderr,"\t-d delay [Default: %u]\n",DELAY);
+ fprintf(stderr,"\t-Q order of the Renyi entropy [Default: %.1f]\n",Q);
+ fprintf(stderr,"\t-r minimal epsilon [Default: (data interval)/1000]\n");
+ fprintf(stderr,"\t-R maximal epsilon [Default: data interval]\n");
+ fprintf(stderr,"\t-# # of epsilons to use [Default: %u]\n",EPSCOUNT);
+ fprintf(stderr,"\t-o output file name [Default: 'datafile'.box]\n");
+ fprintf(stderr,"\t-V verbosity level [Default: 1]\n\t\t"
+ "0='only panic messages'\n\t\t"
+ "1='+ input/output messages'\n");
+ fprintf(stderr,"\t-h show these options\n\n");
+ exit(0);
+}
+
+void scan_options(int n,char **in)
+{
+ char *out;
+
+ if ((out=check_option(in,n,'l','u')) != NULL)
+ sscanf(out,"%lu",&LENGTH);
+ if ((out=check_option(in,n,'c','s')) != NULL)
+ column=out;
+ if ((out=check_option(in,n,'x','u')) != NULL)
+ sscanf(out,"%lu",&exclude);
+ if ((out=check_option(in,n,'M','2')) != NULL) {
+ sscanf(out,"%u,%u",&dimension,&maxembed);
+ dimset=1;
+ }
+ if ((out=check_option(in,n,'d','u')) != NULL)
+ sscanf(out,"%u",&DELAY);
+ if ((out=check_option(in,n,'Q','f')) != NULL)
+ sscanf(out,"%lf",&Q);
+ if ((out=check_option(in,n,'r','f')) != NULL) {
+ sscanf(out,"%lf",&EPSMIN);
+ epsminset=1;
+ }
+ if ((out=check_option(in,n,'R','f')) != NULL) {
+ sscanf(out,"%lf",&EPSMAX);
+ epsmaxset=1;
+ }
+ if ((out=check_option(in,n,'#','u')) != NULL)
+ sscanf(out,"%u",&EPSCOUNT);
+ if ((out=check_option(in,n,'V','u')) != NULL)
+ sscanf(out,"%u",&verbosity);
+ if ((out=check_option(in,n,'o','s')) != NULL)
+ outfile=out;
+}
+
+hliste *make_histo(void)
+{
+ int i;
+ hliste *element;
+
+ check_alloc(element=(hliste*)malloc(sizeof(hliste)));
+ element->ptr=NULL;
+ check_alloc(element->hist=(double*)malloc(sizeof(double)*maxembed*dimension));
+ for (i=0;i<maxembed*dimension;i++)
+ element->hist[i]=0.0;
+
+ return element;
+}
+
+void next_dim(int wd,int n,unsigned int *first)
+{
+ int i,which,d1,comp;
+ double epsinv,norm,p;
+ unsigned int **act;
+ int *found,hf;
+
+ comp=which_dims[wd][0];
+ d1=which_dims[wd][1]*DELAY;
+
+ epsinv=(double)epsi;
+ norm=(double)length;
+
+ check_alloc(act=(unsigned int**)malloc(epsi*sizeof(int*)));
+ check_alloc(found=(int*)malloc(epsi*sizeof(int)));
+
+ for (i=0;i<epsi;i++) {
+ found[i]=0;
+ act[i]=NULL;
+ }
+
+ for (i=0;i<n;i++) {
+ which=(int)(series[comp][first[i]+d1]*epsinv);
+ hf= ++found[which];
+ check_alloc(act[which]=
+ realloc((unsigned int*)act[which],hf*sizeof(unsigned int)));
+ act[which][hf-1]=first[i];
+ }
+
+ for (i=0;i<epsi;i++)
+ if (found[i]) {
+ p=(double)(found[i])/(norm);
+ if (Q == 1.0)
+ histo[wd] -= p*log(p);
+ else
+ histo[wd] += pow(p,Q);
+ }
+
+ if (wd<(maxembed*dimension-1))
+ for (i=0;i<epsi;i++)
+ if (found[i])
+ next_dim(wd+1,found[i],act[i]);
+
+ for (i=0;i<epsi;i++)
+ if (found[i])
+ free(act[i]);
+
+ free(act);
+ free(found);
+}
+
+void start_box(void)
+{
+ int i,which;
+ double epsinv,norm,p;
+ unsigned int **act;
+ int *found,hf;
+ void next_dim();
+
+ epsinv=(double)epsi;
+ norm=(double)length;
+
+ check_alloc(act=(unsigned int**)malloc(epsi*sizeof(int*)));
+ check_alloc(found=(int*)malloc(epsi*sizeof(int)));
+
+ for (i=0;i<epsi;i++) {
+ found[i]=0;
+ act[i]=NULL;
+ }
+
+ for (i=0;i<length;i++) {
+ which=(int)(series[0][i]*epsinv);
+ hf= ++found[w...
[truncated message content] |
|
From: <jpi...@us...> - 2012-03-28 13:27:25
|
Revision: 10079
http://octave.svn.sourceforge.net/octave/?rev=10079&view=rev
Author: jpicarbajal
Date: 2012-03-28 13:27:12 +0000 (Wed, 28 Mar 2012)
Log Message:
-----------
system-identification: adding structure
Modified Paths:
--------------
trunk/octave-forge/main/geometry/PKG_DEL
trunk/octave-forge/main/mechanics/inst/core/masscenter.m
Added Paths:
-----------
trunk/octave-forge/main/system-identification/
trunk/octave-forge/main/system-identification/COPYING
trunk/octave-forge/main/system-identification/DESCRIPTION
trunk/octave-forge/main/system-identification/INDEX
trunk/octave-forge/main/system-identification/NEWS
trunk/octave-forge/main/system-identification/PKG_ADD
trunk/octave-forge/main/system-identification/PKG_DEL
trunk/octave-forge/main/system-identification/devel/
trunk/octave-forge/main/system-identification/inst/
trunk/octave-forge/main/system-identification/inst/tisean/
trunk/octave-forge/main/system-identification/inst/tisean/src/
trunk/octave-forge/main/system-identification/post_install.m
trunk/octave-forge/main/system-identification/pre_install.m
Modified: trunk/octave-forge/main/geometry/PKG_DEL
===================================================================
--- trunk/octave-forge/main/geometry/PKG_DEL 2012-03-28 03:16:58 UTC (rev 10078)
+++ trunk/octave-forge/main/geometry/PKG_DEL 2012-03-28 13:27:12 UTC (rev 10079)
@@ -3,11 +3,16 @@
dirname = fileparts (canonicalize_file_name (mfilename ("fullpath")));
pp = strsplit (dirname,filesep (), true);
-%% Check if prefix was used
-[pkg_folder dep_folder] = pkg ("prefix");
-pkg_folder = [pkg_folder filesep() strcat(filesep(),{pp{end-1:end}}){:} ];
-dep_folder = [dep_folder filesep() strcat(filesep(),{pp{end-1:end}}){:} ];
+%% Get the correct path
+p1 = pkg_info ("miscellaneous", "archprefix");
+p2 = octave_config_info ("canonical_host_type");
+p2 = octave_config_info ("api_version");
+arch_dep_fldr = [p1 filesep "-" p3];
+%[pkg_folder dep_folder] = pkg ("prefix");
+%pkg_folder = [pkg_folder filesep() strcat(filesep(),{pp{end-1:end}}){:} ];
+%dep_folder = [dep_folder filesep() strcat(filesep(),{pp{end-1:end}}){:} ];
+
%% If we are not in Architecture dependent folder
arch = cstrcat (octave_config_info ("canonical_host_type"),
"-", octave_config_info ("api_version"));
Modified: trunk/octave-forge/main/mechanics/inst/core/masscenter.m
===================================================================
--- trunk/octave-forge/main/mechanics/inst/core/masscenter.m 2012-03-28 03:16:58 UTC (rev 10078)
+++ trunk/octave-forge/main/mechanics/inst/core/masscenter.m 2012-03-28 13:27:12 UTC (rev 10079)
@@ -36,8 +36,8 @@
px = x(1,:);
py = x(2,:);
- Px = polyint (conv(conv (px , px)/2 , polyder (py)));
- Py = polyint (conv(-conv (py , py)/2 , polyder (px)));
+ Px = polyint (conv(conv (-px , py) , polyder (px)));
+ Py = polyint (conv(conv (px , py) , polyder (py)));
dcm = zeros (1,2);
dcm(1) = diff(polyval(Px,[0 1]));
Added: trunk/octave-forge/main/system-identification/COPYING
===================================================================
--- trunk/octave-forge/main/system-identification/COPYING (rev 0)
+++ trunk/octave-forge/main/system-identification/COPYING 2012-03-28 13:27:12 UTC (rev 10079)
@@ -0,0 +1,674 @@
+ GNU GENERAL PUBLIC LICENSE
+ Version 3, 29 June 2007
+
+ Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+ Preamble
+
+ The GNU General Public License is a free, copyleft license for
+software and other kinds of works.
+
+ The licenses for most software and other practical works are designed
+to take away your freedom to share and change the works. By contrast,
+the GNU General Public License is intended to guarantee your freedom to
+share and change all versions of a program--to make sure it remains free
+software for all its users. We, the Free Software Foundation, use the
+GNU General Public License for most of our software; it applies also to
+any other work released this way by its authors. You can apply it to
+your programs, too.
+
+ When we speak of free software, we are referring to freedom, not
+price. Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+them if you wish), that you receive source code or can get it if you
+want it, that you can change the software or use pieces of it in new
+free programs, and that you know you can do these things.
+
+ To protect your rights, we need to prevent others from denying you
+these rights or asking you to surrender the rights. Therefore, you have
+certain responsibilities if you distribute copies of the software, or if
+you modify it: responsibilities to respect the freedom of others.
+
+ For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must pass on to the recipients the same
+freedoms that you received. You must make sure that they, too, receive
+or can get the source code. And you must show them these terms so they
+know their rights.
+
+ Developers that use the GNU GPL protect your rights with two steps:
+(1) assert copyright on the software, and (2) offer you this License
+giving you legal permission to copy, distribute and/or modify it.
+
+ For the developers' and authors' protection, the GPL clearly explains
+that there is no warranty for this free software. For both users' and
+authors' sake, the GPL requires that modified versions be marked as
+changed, so that their problems will not be attributed erroneously to
+authors of previous versions.
+
+ Some devices are designed to deny users access to install or run
+modified versions of the software inside them, although the manufacturer
+can do so. This is fundamentally incompatible with the aim of
+protecting users' freedom to change the software. The systematic
+pattern of such abuse occurs in the area of products for individuals to
+use, which is precisely where it is most unacceptable. Therefore, we
+have designed this version of the GPL to prohibit the practice for those
+products. If such problems arise substantially in other domains, we
+stand ready to extend this provision to those domains in future versions
+of the GPL, as needed to protect the freedom of users.
+
+ Finally, every program is threatened constantly by software patents.
+States should not allow patents to restrict development and use of
+software on general-purpose computers, but in those that do, we wish to
+avoid the special danger that patents applied to a free program could
+make it effectively proprietary. To prevent this, the GPL assures that
+patents cannot be used to render the program non-free.
+
+ The precise terms and conditions for copying, distribution and
+modification follow.
+
+ TERMS AND CONDITIONS
+
+ 0. Definitions.
+
+ "This License" refers to version 3 of the GNU General Public License.
+
+ "Copyright" also means copyright-like laws that apply to other kinds of
+works, such as semiconductor masks.
+
+ "The Program" refers to any copyrightable work licensed under this
+License. Each licensee is addressed as "you". "Licensees" and
+"recipients" may be individuals or organizations.
+
+ To "modify" a work means to copy from or adapt all or part of the work
+in a fashion requiring copyright permission, other than the making of an
+exact copy. The resulting work is called a "modified version" of the
+earlier work or a work "based on" the earlier work.
+
+ A "covered work" means either the unmodified Program or a work based
+on the Program.
+
+ To "propagate" a work means to do anything with it that, without
+permission, would make you directly or secondarily liable for
+infringement under applicable copyright law, except executing it on a
+computer or modifying a private copy. Propagation includes copying,
+distribution (with or without modification), making available to the
+public, and in some countries other activities as well.
+
+ To "convey" a work means any kind of propagation that enables other
+parties to make or receive copies. Mere interaction with a user through
+a computer network, with no transfer of a copy, is not conveying.
+
+ An interactive user interface displays "Appropriate Legal Notices"
+to the extent that it includes a convenient and prominently visible
+feature that (1) displays an appropriate copyright notice, and (2)
+tells the user that there is no warranty for the work (except to the
+extent that warranties are provided), that licensees may convey the
+work under this License, and how to view a copy of this License. If
+the interface presents a list of user commands or options, such as a
+menu, a prominent item in the list meets this criterion.
+
+ 1. Source Code.
+
+ The "source code" for a work means the preferred form of the work
+for making modifications to it. "Object code" means any non-source
+form of a work.
+
+ A "Standard Interface" means an interface that either is an official
+standard defined by a recognized standards body, or, in the case of
+interfaces specified for a particular programming language, one that
+is widely used among developers working in that language.
+
+ The "System Libraries" of an executable work include anything, other
+than the work as a whole, that (a) is included in the normal form of
+packaging a Major Component, but which is not part of that Major
+Component, and (b) serves only to enable use of the work with that
+Major Component, or to implement a Standard Interface for which an
+implementation is available to the public in source code form. A
+"Major Component", in this context, means a major essential component
+(kernel, window system, and so on) of the specific operating system
+(if any) on which the executable work runs, or a compiler used to
+produce the work, or an object code interpreter used to run it.
+
+ The "Corresponding Source" for a work in object code form means all
+the source code needed to generate, install, and (for an executable
+work) run the object code and to modify the work, including scripts to
+control those activities. However, it does not include the work's
+System Libraries, or general-purpose tools or generally available free
+programs which are used unmodified in performing those activities but
+which are not part of the work. For example, Corresponding Source
+includes interface definition files associated with source files for
+the work, and the source code for shared libraries and dynamically
+linked subprograms that the work is specifically designed to require,
+such as by intimate data communication or control flow between those
+subprograms and other parts of the work.
+
+ The Corresponding Source need not include anything that users
+can regenerate automatically from other parts of the Corresponding
+Source.
+
+ The Corresponding Source for a work in source code form is that
+same work.
+
+ 2. Basic Permissions.
+
+ All rights granted under this License are granted for the term of
+copyright on the Program, and are irrevocable provided the stated
+conditions are met. This License explicitly affirms your unlimited
+permission to run the unmodified Program. The output from running a
+covered work is covered by this License only if the output, given its
+content, constitutes a covered work. This License acknowledges your
+rights of fair use or other equivalent, as provided by copyright law.
+
+ You may make, run and propagate covered works that you do not
+convey, without conditions so long as your license otherwise remains
+in force. You may convey covered works to others for the sole purpose
+of having them make modifications exclusively for you, or provide you
+with facilities for running those works, provided that you comply with
+the terms of this License in conveying all material for which you do
+not control copyright. Those thus making or running the covered works
+for you must do so exclusively on your behalf, under your direction
+and control, on terms that prohibit them from making any copies of
+your copyrighted material outside their relationship with you.
+
+ Conveying under any other circumstances is permitted solely under
+the conditions stated below. Sublicensing is not allowed; section 10
+makes it unnecessary.
+
+ 3. Protecting Users' Legal Rights From Anti-Circumvention Law.
+
+ No covered work shall be deemed part of an effective technological
+measure under any applicable law fulfilling obligations under article
+11 of the WIPO copyright treaty adopted on 20 December 1996, or
+similar laws prohibiting or restricting circumvention of such
+measures.
+
+ When you convey a covered work, you waive any legal power to forbid
+circumvention of technological measures to the extent such circumvention
+is effected by exercising rights under this License with respect to
+the covered work, and you disclaim any intention to limit operation or
+modification of the work as a means of enforcing, against the work's
+users, your or third parties' legal rights to forbid circumvention of
+technological measures.
+
+ 4. Conveying Verbatim Copies.
+
+ You may convey verbatim copies of the Program's source code as you
+receive it, in any medium, provided that you conspicuously and
+appropriately publish on each copy an appropriate copyright notice;
+keep intact all notices stating that this License and any
+non-permissive terms added in accord with section 7 apply to the code;
+keep intact all notices of the absence of any warranty; and give all
+recipients a copy of this License along with the Program.
+
+ You may charge any price or no price for each copy that you convey,
+and you may offer support or warranty protection for a fee.
+
+ 5. Conveying Modified Source Versions.
+
+ You may convey a work based on the Program, or the modifications to
+produce it from the Program, in the form of source code under the
+terms of section 4, provided that you also meet all of these conditions:
+
+ a) The work must carry prominent notices stating that you modified
+ it, and giving a relevant date.
+
+ b) The work must carry prominent notices stating that it is
+ released under this License and any conditions added under section
+ 7. This requirement modifies the requirement in section 4 to
+ "keep intact all notices".
+
+ c) You must license the entire work, as a whole, under this
+ License to anyone who comes into possession of a copy. This
+ License will therefore apply, along with any applicable section 7
+ additional terms, to the whole of the work, and all its parts,
+ regardless of how they are packaged. This License gives no
+ permission to license the work in any other way, but it does not
+ invalidate such permission if you have separately received it.
+
+ d) If the work has interactive user interfaces, each must display
+ Appropriate Legal Notices; however, if the Program has interactive
+ interfaces that do not display Appropriate Legal Notices, your
+ work need not make them do so.
+
+ A compilation of a covered work with other separate and independent
+works, which are not by their nature extensions of the covered work,
+and which are not combined with it such as to form a larger program,
+in or on a volume of a storage or distribution medium, is called an
+"aggregate" if the compilation and its resulting copyright are not
+used to limit the access or legal rights of the compilation's users
+beyond what the individual works permit. Inclusion of a covered work
+in an aggregate does not cause this License to apply to the other
+parts of the aggregate.
+
+ 6. Conveying Non-Source Forms.
+
+ You may convey a covered work in object code form under the terms
+of sections 4 and 5, provided that you also convey the
+machine-readable Corresponding Source under the terms of this License,
+in one of these ways:
+
+ a) Convey the object code in, or embodied in, a physical product
+ (including a physical distribution medium), accompanied by the
+ Corresponding Source fixed on a durable physical medium
+ customarily used for software interchange.
+
+ b) Convey the object code in, or embodied in, a physical product
+ (including a physical distribution medium), accompanied by a
+ written offer, valid for at least three years and valid for as
+ long as you offer spare parts or customer support for that product
+ model, to give anyone who possesses the object code either (1) a
+ copy of the Corresponding Source for all the software in the
+ product that is covered by this License, on a durable physical
+ medium customarily used for software interchange, for a price no
+ more than your reasonable cost of physically performing this
+ conveying of source, or (2) access to copy the
+ Corresponding Source from a network server at no charge.
+
+ c) Convey individual copies of the object code with a copy of the
+ written offer to provide the Corresponding Source. This
+ alternative is allowed only occasionally and noncommercially, and
+ only if you received the object code with such an offer, in accord
+ with subsection 6b.
+
+ d) Convey the object code by offering access from a designated
+ place (gratis or for a charge), and offer equivalent access to the
+ Corresponding Source in the same way through the same place at no
+ further charge. You need not require recipients to copy the
+ Corresponding Source along with the object code. If the place to
+ copy the object code is a network server, the Corresponding Source
+ may be on a different server (operated by you or a third party)
+ that supports equivalent copying facilities, provided you maintain
+ clear directions next to the object code saying where to find the
+ Corresponding Source. Regardless of what server hosts the
+ Corresponding Source, you remain obligated to ensure that it is
+ available for as long as needed to satisfy these requirements.
+
+ e) Convey the object code using peer-to-peer transmission, provided
+ you inform other peers where the object code and Corresponding
+ Source of the work are being offered to the general public at no
+ charge under subsection 6d.
+
+ A separable portion of the object code, whose source code is excluded
+from the Corresponding Source as a System Library, need not be
+included in conveying the object code work.
+
+ A "User Product" is either (1) a "consumer product", which means any
+tangible personal property which is normally used for personal, family,
+or household purposes, or (2) anything designed or sold for incorporation
+into a dwelling. In determining whether a product is a consumer product,
+doubtful cases shall be resolved in favor of coverage. For a particular
+product received by a particular user, "normally used" refers to a
+typical or common use of that class of product, regardless of the status
+of the particular user or of the way in which the particular user
+actually uses, or expects or is expected to use, the product. A product
+is a consumer product regardless of whether the product has substantial
+commercial, industrial or non-consumer uses, unless such uses represent
+the only significant mode of use of the product.
+
+ "Installation Information" for a User Product means any methods,
+procedures, authorization keys, or other information required to install
+and execute modified versions of a covered work in that User Product from
+a modified version of its Corresponding Source. The information must
+suffice to ensure that the continued functioning of the modified object
+code is in no case prevented or interfered with solely because
+modification has been made.
+
+ If you convey an object code work under this section in, or with, or
+specifically for use in, a User Product, and the conveying occurs as
+part of a transaction in which the right of possession and use of the
+User Product is transferred to the recipient in perpetuity or for a
+fixed term (regardless of how the transaction is characterized), the
+Corresponding Source conveyed under this section must be accompanied
+by the Installation Information. But this requirement does not apply
+if neither you nor any third party retains the ability to install
+modified object code on the User Product (for example, the work has
+been installed in ROM).
+
+ The requirement to provide Installation Information does not include a
+requirement to continue to provide support service, warranty, or updates
+for a work that has been modified or installed by the recipient, or for
+the User Product in which it has been modified or installed. Access to a
+network may be denied when the modification itself materially and
+adversely affects the operation of the network or violates the rules and
+protocols for communication across the network.
+
+ Corresponding Source conveyed, and Installation Information provided,
+in accord with this section must be in a format that is publicly
+documented (and with an implementation available to the public in
+source code form), and must require no special password or key for
+unpacking, reading or copying.
+
+ 7. Additional Terms.
+
+ "Additional permissions" are terms that supplement the terms of this
+License by making exceptions from one or more of its conditions.
+Additional permissions that are applicable to the entire Program shall
+be treated as though they were included in this License, to the extent
+that they are valid under applicable law. If additional permissions
+apply only to part of the Program, that part may be used separately
+under those permissions, but the entire Program remains governed by
+this License without regard to the additional permissions.
+
+ When you convey a copy of a covered work, you may at your option
+remove any additional permissions from that copy, or from any part of
+it. (Additional permissions may be written to require their own
+removal in certain cases when you modify the work.) You may place
+additional permissions on material, added by you to a covered work,
+for which you have or can give appropriate copyright permission.
+
+ Notwithstanding any other provision of this License, for material you
+add to a covered work, you may (if authorized by the copyright holders of
+that material) supplement the terms of this License with terms:
+
+ a) Disclaiming warranty or limiting liability differently from the
+ terms of sections 15 and 16 of this License; or
+
+ b) Requiring preservation of specified reasonable legal notices or
+ author attributions in that material or in the Appropriate Legal
+ Notices displayed by works containing it; or
+
+ c) Prohibiting misrepresentation of the origin of that material, or
+ requiring that modified versions of such material be marked in
+ reasonable ways as different from the original version; or
+
+ d) Limiting the use for publicity purposes of names of licensors or
+ authors of the material; or
+
+ e) Declining to grant rights under trademark law for use of some
+ trade names, trademarks, or service marks; or
+
+ f) Requiring indemnification of licensors and authors of that
+ material by anyone who conveys the material (or modified versions of
+ it) with contractual assumptions of liability to the recipient, for
+ any liability that these contractual assumptions directly impose on
+ those licensors and authors.
+
+ All other non-permissive additional terms are considered "further
+restrictions" within the meaning of section 10. If the Program as you
+received it, or any part of it, contains a notice stating that it is
+governed by this License along with a term that is a further
+restriction, you may remove that term. If a license document contains
+a further restriction but permits relicensing or conveying under this
+License, you may add to a covered work material governed by the terms
+of that license document, provided that the further restriction does
+not survive such relicensing or conveying.
+
+ If you add terms to a covered work in accord with this section, you
+must place, in the relevant source files, a statement of the
+additional terms that apply to those files, or a notice indicating
+where to find the applicable terms.
+
+ Additional terms, permissive or non-permissive, may be stated in the
+form of a separately written license, or stated as exceptions;
+the above requirements apply either way.
+
+ 8. Termination.
+
+ You may not propagate or modify a covered work except as expressly
+provided under this License. Any attempt otherwise to propagate or
+modify it is void, and will automatically terminate your rights under
+this License (including any patent licenses granted under the third
+paragraph of section 11).
+
+ However, if you cease all violation of this License, then your
+license from a particular copyright holder is reinstated (a)
+provisionally, unless and until the copyright holder explicitly and
+finally terminates your license, and (b) permanently, if the copyright
+holder fails to notify you of the violation by some reasonable means
+prior to 60 days after the cessation.
+
+ Moreover, your license from a particular copyright holder is
+reinstated permanently if the copyright holder notifies you of the
+violation by some reasonable means, this is the first time you have
+received notice of violation of this License (for any work) from that
+copyright holder, and you cure the violation prior to 30 days after
+your receipt of the notice.
+
+ Termination of your rights under this section does not terminate the
+licenses of parties who have received copies or rights from you under
+this License. If your rights have been terminated and not permanently
+reinstated, you do not qualify to receive new licenses for the same
+material under section 10.
+
+ 9. Acceptance Not Required for Having Copies.
+
+ You are not required to accept this License in order to receive or
+run a copy of the Program. Ancillary propagation of a covered work
+occurring solely as a consequence of using peer-to-peer transmission
+to receive a copy likewise does not require acceptance. However,
+nothing other than this License grants you permission to propagate or
+modify any covered work. These actions infringe copyright if you do
+not accept this License. Therefore, by modifying or propagating a
+covered work, you indicate your acceptance of this License to do so.
+
+ 10. Automatic Licensing of Downstream Recipients.
+
+ Each time you convey a covered work, the recipient automatically
+receives a license from the original licensors, to run, modify and
+propagate that work, subject to this License. You are not responsible
+for enforcing compliance by third parties with this License.
+
+ An "entity transaction" is a transaction transferring control of an
+organization, or substantially all assets of one, or subdividing an
+organization, or merging organizations. If propagation of a covered
+work results from an entity transaction, each party to that
+transaction who receives a copy of the work also receives whatever
+licenses to the work the party's predecessor in interest had or could
+give under the previous paragraph, plus a right to possession of the
+Corresponding Source of the work from the predecessor in interest, if
+the predecessor has it or can get it with reasonable efforts.
+
+ You may not impose any further restrictions on the exercise of the
+rights granted or affirmed under this License. For example, you may
+not impose a license fee, royalty, or other charge for exercise of
+rights granted under this License, and you may not initiate litigation
+(including a cross-claim or counterclaim in a lawsuit) alleging that
+any patent claim is infringed by making, using, selling, offering for
+sale, or importing the Program or any portion of it.
+
+ 11. Patents.
+
+ A "contributor" is a copyright holder who authorizes use under this
+License of the Program or a work on which the Program is based. The
+work thus licensed is called the contributor's "contributor version".
+
+ A contributor's "essential patent claims" are all patent claims
+owned or controlled by the contributor, whether already acquired or
+hereafter acquired, that would be infringed by some manner, permitted
+by this License, of making, using, or selling its contributor version,
+but do not include claims that would be infringed only as a
+consequence of further modification of the contributor version. For
+purposes of this definition, "control" includes the right to grant
+patent sublicenses in a manner consistent with the requirements of
+this License.
+
+ Each contributor grants you a non-exclusive, worldwide, royalty-free
+patent license under the contributor's essential patent claims, to
+make, use, sell, offer for sale, import and otherwise run, modify and
+propagate the contents of its contributor version.
+
+ In the following three paragraphs, a "patent license" is any express
+agreement or commitment, however denominated, not to enforce a patent
+(such as an express permission to practice a patent or covenant not to
+sue for patent infringement). To "grant" such a patent license to a
+party means to make such an agreement or commitment not to enforce a
+patent against the party.
+
+ If you convey a covered work, knowingly relying on a patent license,
+and the Corresponding Source of the work is not available for anyone
+to copy, free of charge and under the terms of this License, through a
+publicly available network server or other readily accessible means,
+then you must either (1) cause the Corresponding Source to be so
+available, or (2) arrange to deprive yourself of the benefit of the
+patent license for this particular work, or (3) arrange, in a manner
+consistent with the requirements of this License, to extend the patent
+license to downstream recipients. "Knowingly relying" means you have
+actual knowledge that, but for the patent license, your conveying the
+covered work in a country, or your recipient's use of the covered work
+in a country, would infringe one or more identifiable patents in that
+country that you have reason to believe are valid.
+
+ If, pursuant to or in connection with a single transaction or
+arrangement, you convey, or propagate by procuring conveyance of, a
+covered work, and grant a patent license to some of the parties
+receiving the covered work authorizing them to use, propagate, modify
+or convey a specific copy of the covered work, then the patent license
+you grant is automatically extended to all recipients of the covered
+work and works based on it.
+
+ A patent license is "discriminatory" if it does not include within
+the scope of its coverage, prohibits the exercise of, or is
+conditioned on the non-exercise of one or more of the rights that are
+specifically granted under this License. You may not convey a covered
+work if you are a party to an arrangement with a third party that is
+in the business of distributing software, under which you make payment
+to the third party based on the extent of your activity of conveying
+the work, and under which the third party grants, to any of the
+parties who would receive the covered work from you, a discriminatory
+patent license (a) in connection with copies of the covered work
+conveyed by you (or copies made from those copies), or (b) primarily
+for and in connection with specific products or compilations that
+contain the covered work, unless you entered into that arrangement,
+or that patent license was granted, prior to 28 March 2007.
+
+ Nothing in this License shall be construed as excluding or limiting
+any implied license or other defenses to infringement that may
+otherwise be available to you under applicable patent law.
+
+ 12. No Surrender of Others' Freedom.
+
+ If conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License. If you cannot convey a
+covered work so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you may
+not convey it at all. For example, if you agree to terms that obligate you
+to collect a royalty for further conveying from those to whom you convey
+the Program, the only way you could satisfy both those terms and this
+License would be to refrain entirely from conveying the Program.
+
+ 13. Use with the GNU Affero General Public License.
+
+ Notwithstanding any other provision of this License, you have
+permission to link or combine any covered work with a work licensed
+under version 3 of the GNU Affero General Public License into a single
+combined work, and to convey the resulting work. The terms of this
+License will continue to apply to the part which is the covered work,
+but the special requirements of the GNU Affero General Public License,
+section 13, concerning interaction through a network will apply to the
+combination as such.
+
+ 14. Revised Versions of this License.
+
+ The Free Software Foundation may publish revised and/or new versions of
+the GNU General Public License from time to time. Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+ Each version is given a distinguishing version number. If the
+Program specifies that a certain numbered version of the GNU General
+Public License "or any later version" applies to it, you have the
+option of following the terms and conditions either of that numbered
+version or of any later version published by the Free Software
+Foundation. If the Program does not specify a version number of the
+GNU General Public License, you may choose any version ever published
+by the Free Software Foundation.
+
+ If the Program specifies that a proxy can decide which future
+versions of the GNU General Public License can be used, that proxy's
+public statement of acceptance of a version permanently authorizes you
+to choose that version for the Program.
+
+ Later license versions may give you additional or different
+permissions. However, no additional obligations are imposed on any
+author or copyright holder as a result of your choosing to follow a
+later version.
+
+ 15. Disclaimer of Warranty.
+
+ THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY
+APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT
+HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY
+OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO,
+THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM
+IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF
+ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
+
+ 16. Limitation of Liability.
+
+ IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS
+THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY
+GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE
+USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF
+DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD
+PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS),
+EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF
+SUCH DAMAGES.
+
+ 17. Interpretation of Sections 15 and 16.
+
+ If the disclaimer of warranty and limitation of liability provided
+above cannot be given local legal effect according to their terms,
+reviewing courts shall apply local law that most closely approximates
+an absolute waiver of all civil liability in connection with the
+Program, unless a warranty or assumption of liability accompanies a
+copy of the Program in return for a fee.
+
+ END OF TERMS AND CONDITIONS
+
+ How to Apply These Terms to Your New Programs
+
+ If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+ To do so, attach the following notices to the program. It is safest
+to attach them to the start of each source file to most effectively
+state the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+ <one line to give the program's name and a brief idea of what it does.>
+ Copyright (C) <year> <name of author>
+
+ 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/>.
+
+Also add information on how to contact you by electronic and paper mail.
+
+ If the program does terminal interaction, make it output a short
+notice like this when it starts in an interactive mode:
+
+ <program> Copyright (C) <year> <name of author>
+ This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+ This is free software, and you are welcome to redistribute it
+ under certain conditions; type `show c' for details.
+
+The hypothetical commands `show w' and `show c' should show the appropriate
+parts of the General Public License. Of course, your program's commands
+might be different; for a GUI interface, you would use an "about box".
+
+ You should also get your employer (if you work as a programmer) or school,
+if any, to sign a "copyright disclaimer" for the program, if necessary.
+For more information on this, and how to apply and follow the GNU GPL, see
+<http://www.gnu.org/licenses/>.
+
+ The GNU General Public License does not permit incorporating your program
+into proprietary programs. If your program is a subroutine library, you
+may consider it more useful to permit linking proprietary applications with
+the library. If this is what you want to do, use the GNU Lesser General
+Public License instead of this License. But first, please read
+<http://www.gnu.org/philosophy/why-not-lgpl.html>.
Added: trunk/octave-forge/main/system-identification/DESCRIPTION
===================================================================
--- trunk/octave-forge/main/system-identification/DESCRIPTION (rev 0)
+++ trunk/octave-forge/main/system-identification/DESCRIPTION 2012-03-28 13:27:12 UTC (rev 10079)
@@ -0,0 +1,11 @@
+Name: System Identification
+Version: 0.1.0
+Date: 2012-03-28
+Author: Lukas Reichlin <luk...@gm...>, Juan Pablo Carbajal <car...@if...>
+Maintainer: Lukas Reichlin <luk...@gm...>, Juan Pablo Carbajal <car...@if...>
+Title: System Identification
+Description: Library for the identification of dynamical systems from data (e.g. time series). It contains the TISEAN package.
+Depends: octave (>= 3.6.1)
+Autoload: no
+License: GPLv3+
+Url: http://octave.sf.net, http://www.mpipks-dresden.mpg.de/~tisean/
Added: trunk/octave-forge/main/system-identification/INDEX
===================================================================
--- trunk/octave-forge/main/system-identification/INDEX (rev 0)
+++ trunk/octave-forge/main/system-identification/INDEX 2012-03-28 13:27:12 UTC (rev 10079)
@@ -0,0 +1,4 @@
+system-identification >> System Identification
+TISEAN
+Data types
+
Added: trunk/octave-forge/main/system-identification/NEWS
===================================================================
--- trunk/octave-forge/main/system-identification/NEWS (rev 0)
+++ trunk/octave-forge/main/system-identification/NEWS 2012-03-28 13:27:12 UTC (rev 10079)
@@ -0,0 +1,9 @@
+Summary of important user-visible changes for releases of the system-identification package
+
+===============================================================================
+geometry-0.1.0 Release Date: XX Release Manager: Juan Pablo Carbajal
+===============================================================================
+
+** First official release.
+
+===============================================================================
Added: trunk/octave-forge/main/system-identification/PKG_ADD
===================================================================
--- trunk/octave-forge/main/system-identification/PKG_ADD (rev 0)
+++ trunk/octave-forge/main/system-identification/PKG_ADD 2012-03-28 13:27:12 UTC (rev 10079)
@@ -0,0 +1,35 @@
+%1
+dirlist = {"tisean"};
+dirname = fileparts (canonicalize_file_name (mfilename ("fullpath")));
+pp = strsplit (dirname,filesep (), true);
+
+%% Check if prefix was used
+[pkg_folder dep_folder] = pkg ("prefix");
+pkg_folder = [pkg_folder filesep() strcat(filesep(),{pp{end-1:end}}){:} ];
+dep_folder = [dep_folder filesep() strcat(filesep(),{pp{end-1:end}}){:} ];
+
+%% If we are in Architecture dependent folder add from outside
+arch = cstrcat (octave_config_info ("canonical_host_type"),
+ "-", octave_config_info ("api_version"));
+if strcmp (arch , pp{end})
+ dirname = [strcat(filesep(),{pp{1:end-1}}){:}];
+ pkg_folder = strsplit (pkg_folder,filesep (), true);
+ pkg_folder = [strcat(filesep(),{pkg_folder{1:end-1}}){:}];
+end
+
+if (! exist (fullfile (dirname, "inst"), "dir"))
+%% Installing
+ for ii=1:length (dirlist)
+ addpath ( [ pkg_folder filesep() dirlist{ii}],"-end")
+ endfor
+
+else
+%% Testing
+ warning("system-identitifaction:Devel","Adding path for testing.");
+ for ii=1:length(dirlist)
+ addpath ([ dirname "/inst/" dirlist{ii}])
+ endfor
+endif
+
+warning('off', 'Octave:fopen-file-in-path');
+clear dirlist dirname pp arch pkg_folder dep_folder
Added: trunk/octave-forge/main/system-identification/PKG_DEL
===================================================================
--- trunk/octave-forge/main/system-identification/PKG_DEL (rev 0)
+++ trunk/octave-forge/main/system-identification/PKG_DEL 2012-03-28 13:27:12 UTC (rev 10079)
@@ -0,0 +1,38 @@
+%1
+dirlist = {"tisean"};
+dirname = fileparts (canonicalize_file_name (mfilename ("fullpath")));
+pp = strsplit (dirname,filesep (), true);
+
+%% Get the correct path
+p1 = pkg_info ("miscellaneous", "archprefix");
+p2 = octave_config_info ("canonical_host_type");
+p2 = octave_config_info ("api_version");
+arch_dep_fldr = [p1 filesep "-" p3];
+
+%[pkg_folder dep_folder] = pkg ("prefix");
+%pkg_folder = [pkg_folder filesep() strcat(filesep(),{pp{end-1:end}}){:} ];
+%dep_folder = [dep_folder filesep() strcat(filesep(),{pp{end-1:end}}){:} ];
+
+%% If we are not in Architecture dependent folder
+arch = cstrcat (octave_config_info ("canonical_host_type"),
+ "-", octave_config_info ("api_version"));
+pp = strsplit (dirname,filesep (), true);
+if strcmp(arch , pp{end})
+ dirname = [pkg("prefix") filesep() pp{end-1}];
+ pkg_folder = strsplit (pkg_folder,filesep (), true);
+ pkg_folder = [strcat(filesep(),{pkg_folder{1:end-1}}){:}];
+end
+
+if (! exist (fullfile (dirname, "inst"), "dir"))
+## Run this if the package is installed
+ for ii=1:length (dirlist)
+ rmpath ( [ pkg_folder filesep() dirlist{ii}])
+ endfor
+else
+ warning("system-identification:Devel","Removing path for testing.");
+ for ii=1:length(dirlist)
+ rmpath ([ dirname "/inst/" dirlist{ii}])
+ endfor
+endif
+
+clear dirlist dirname pp arch pkg_folder dep_folder
Added: trunk/octave-forge/main/system-identification/post_install.m
===================================================================
--- trunk/octave-forge/main/system-identification/post_install.m (rev 0)
+++ trunk/octave-forge/main/system-identification/post_install.m 2012-03-28 13:27:12 UTC (rev 10079)
@@ -0,0 +1,7 @@
+function post_install (desc)
+%% Prepares for installation a package that is organized in subfolders
+%% Since src is compiled only in the package main dir
+%% I need to remove the PKG_ADD and PKG_DEL from the architecture dependent folder
+
+
+end
Added: trunk/octave-forge/main/system-identification/pre_install.m
===================================================================
--- trunk/octave-forge/main/system-identification/pre_install.m (rev 0)
+++ trunk/octave-forge/main/system-identification/pre_install.m 2012-03-28 13:27:12 UTC (rev 10079)
@@ -0,0 +1,28 @@
+function pre_install (desc)
+%% Prepares for installation a package that is organized in subfolders
+
+ %% List of folders with src subfolder
+ subfld = {"tisean"};
+
+ %% Create correct strings
+ subfld_ready = strcat ({[pwd() filesep() "inst" filesep()]},
+ subfld,[filesep() "src" filesep() "*"]);
+
+ %% Destination folder
+ to_fld = strcat (pwd (),[filesep() "src"]);
+
+
+ %% Copy files to package/src folder
+ %% TODO handle merging of Makefiles
+ warning ("Copying subfolder src to package main dir, but multiple Makefiles are not handled")
+
+ if !exist("src","dir")
+ system(["mkdir " to_fld]);
+ end
+
+ for from_fld = subfld_ready
+ system (["mv " from_fld{1} " " to_fld]);
+ system (["rm -R " from_fld{1}(1:end-2)]);
+ end
+
+end
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-03-28 03:17:04
|
Revision: 10078
http://octave.svn.sourceforge.net/octave/?rev=10078&view=rev
Author: paramaniac
Date: 2012-03-28 03:16:58 +0000 (Wed, 28 Mar 2012)
Log Message:
-----------
main/quaternion --> extra/quaternion-legacy, extra/quaternion_oo --> main/quaternion
Added Paths:
-----------
trunk/octave-forge/extra/quaternion-legacy/
trunk/octave-forge/main/quaternion/
Removed Paths:
-------------
trunk/octave-forge/extra/quaternion_oo/
trunk/octave-forge/main/quaternion/
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-03-28 02:57:46
|
Revision: 10077
http://octave.svn.sourceforge.net/octave/?rev=10077&view=rev
Author: paramaniac
Date: 2012-03-28 02:57:39 +0000 (Wed, 28 Mar 2012)
Log Message:
-----------
quaternion_oo: prepare release of quaternion-2.0.0
Modified Paths:
--------------
trunk/octave-forge/extra/quaternion_oo/DESCRIPTION
trunk/octave-forge/extra/quaternion_oo/NEWS
Added Paths:
-----------
trunk/octave-forge/extra/quaternion_oo/doc/quaternion.pdf
Modified: trunk/octave-forge/extra/quaternion_oo/DESCRIPTION
===================================================================
--- trunk/octave-forge/extra/quaternion_oo/DESCRIPTION 2012-03-27 18:23:55 UTC (rev 10076)
+++ trunk/octave-forge/extra/quaternion_oo/DESCRIPTION 2012-03-28 02:57:39 UTC (rev 10077)
@@ -1,6 +1,6 @@
Name: quaternion
-Version: 1.9.52
-Date: 2012-01-18
+Version: 2.0.0
+Date: 2012-03-28
Author: Lukas Reichlin <luk...@gm...>
Maintainer: Lukas Reichlin <luk...@gm...>
Title: Quaternion
Modified: trunk/octave-forge/extra/quaternion_oo/NEWS
===================================================================
--- trunk/octave-forge/extra/quaternion_oo/NEWS 2012-03-27 18:23:55 UTC (rev 10076)
+++ trunk/octave-forge/extra/quaternion_oo/NEWS 2012-03-28 02:57:39 UTC (rev 10077)
@@ -1,7 +1,7 @@
Summary of important user-visible changes for releases of the quaternion package
===============================================================================
-quaternion-2.0.0 Release Date: 2012-xx-yy Release Manager: Lukas Reichlin
+quaternion-2.0.0 Release Date: 2012-03-28 Release Manager: Lukas Reichlin
===============================================================================
** First official release. Its main features are:
Added: trunk/octave-forge/extra/quaternion_oo/doc/quaternion.pdf
===================================================================
(Binary files differ)
Property changes on: trunk/octave-forge/extra/quaternion_oo/doc/quaternion.pdf
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <nir...@us...> - 2012-03-27 18:24:01
|
Revision: 10076
http://octave.svn.sourceforge.net/octave/?rev=10076&view=rev
Author: nir-krakauer
Date: 2012-03-27 18:23:55 +0000 (Tue, 27 Mar 2012)
Log Message:
-----------
Added a reference
Modified Paths:
--------------
trunk/octave-forge/main/splines/inst/csaps.m
Modified: trunk/octave-forge/main/splines/inst/csaps.m
===================================================================
--- trunk/octave-forge/main/splines/inst/csaps.m 2012-03-27 17:50:14 UTC (rev 10075)
+++ trunk/octave-forge/main/splines/inst/csaps.m 2012-03-27 18:23:55 UTC (rev 10076)
@@ -33,6 +33,8 @@
## an intermediate amount of smoothing is chosen (the smoothing term and the interpolation term are scaled to be of the same magnitude)
## @end table
##
+## Reference: Carl de Boor (1978), A Practical Guide to Splines, Springer, Chapter XIV
+##
## @end deftypefn
## @seealso{spline, csapi, ppval, gcvspl}
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <nir...@us...> - 2012-03-27 17:50:21
|
Revision: 10075
http://octave.svn.sourceforge.net/octave/?rev=10075&view=rev
Author: nir-krakauer
Date: 2012-03-27 17:50:14 +0000 (Tue, 27 Mar 2012)
Log Message:
-----------
Added new self-contained csaps function in splines for producing cubic smoothing splines
Added Paths:
-----------
trunk/octave-forge/main/splines/inst/csaps.m
Added: trunk/octave-forge/main/splines/inst/csaps.m
===================================================================
--- trunk/octave-forge/main/splines/inst/csaps.m (rev 0)
+++ trunk/octave-forge/main/splines/inst/csaps.m 2012-03-27 17:50:14 UTC (rev 10075)
@@ -0,0 +1,108 @@
+## Copyright (C) 2012 Nir Krakauer
+##
+## 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 2 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{yi} @var{p}] =} csaps(@var{x}, @var{y}, @var{p}, @var{xi}, @var{w}=[])
+## @deftypefnx{Function File}{[@var{pp} @var{p}] =} csaps(@var{x}, @var{y}, @var{p}, [], @var{w}=[])
+##
+## Cubic spline approximation (smoothing)@*
+## approximate [@var{x},@var{y}], weighted by @var{w} (inverse variance; if not given, equal weighting is assumed), at @var{xi}
+##
+## The chosen cubic spline with natural boundary conditions @var{pp}(@var{x}) minimizes @var{p} Sum_i @var{w}_i*(@var{y}_i - @var{pp}(@var{x}_i))^2 + (1-@var{p}) Int @var{pp}''(@var{x}) d@var{x}
+##
+## @var{x} and @var{w} should be n by 1 in size; @var{y} should be n by m; @var{xi} should be k by 1; the values in @var{x} should be distinct; the values in @var{w} should be nonzero
+##
+## @table @asis
+## @item @var{p}=0
+## maximum smoothing: straight line
+## @item @var{p}=1
+## no smoothing: interpolation
+## @item @var{p}<0 or not given
+## an intermediate amount of smoothing is chosen (the smoothing term and the interpolation term are scaled to be of the same magnitude)
+## @end table
+##
+## @end deftypefn
+## @seealso{spline, csapi, ppval, gcvspl}
+
+## Author: Nir Krakauer <nkr...@cc...>
+
+function [ret,p]=csaps(x,y,p,xi,w)
+
+ if(nargin < 5)
+ w = [];
+ if(nargin < 4)
+ xi = [];
+ if(nargin < 3)
+ p = [];
+ endif
+ endif
+ endif
+
+ if(columns(x) > 1)
+ x = x.';
+ y = y.';
+ w = w.';
+ endif
+
+ [x,i] = sort(x);
+ y = y(i, :);
+
+ n = numel(x);
+
+ if isempty(w)
+ w = ones(n, 1);
+ end
+
+ h = diff(x);
+
+ R = spdiags([h(1:end-1) 2*(h(1:end-1) + h(2:end)) h(2:end)], [-1 0 1], n-2, n-2);
+
+ QT = spdiags([1 ./ h(1:end-1) -(1 ./ h(1:end-1) + 1 ./ h(2:end)) 1 ./ h(2:end)], [0 1 2], n-2, n);
+
+## if not given, choose p so that trace(6*(1-p)*QT*diag(1 ./ w)*QT') = trace(pR)
+ if isempty(p) || (p < 0)
+ r = 6*trace(QT*diag(1 ./ w)*QT') / trace(R);
+ p = 1 ./ (1 + r);
+ endif
+
+## solve for the scaled second derivatives u and for the function values a at the knots (if p = 1, a = y)
+ u = (6*(1-p)*QT*diag(1 ./ w)*QT' + p*R) \ (QT*y);
+ a = y - 6*(1-p)*diag(1 ./ w)*QT'*u;
+
+## derivatives at all but the last knot for the piecewise cubic spline
+ aa = a(1:(end-1), :);
+ cc = zeros(size(y));
+ cc(2:(n-1), :) = 6*p*u; #cc([1 n], :) = 0 [natural spline]
+ dd = diff(cc) ./ h;
+ cc = cc(1:(end-1), :);
+ bb = diff(a) ./ h - (cc/2).*h - (dd/6).*(h.^2);
+
+ ret = mkpp (x, cat (2, dd'(:)/6, cc'(:)/2, bb'(:), aa'(:)), size(y, 2));
+
+ if ~isempty(xi)
+ ret = ppval (ret, xi);
+ endif
+
+endfunction
+
+%!shared x,y
+%! x = ([1:10 10.5 11.3])'; y = sin(x);
+
+%!assert (csaps(x,y,1,x), y, 10*eps);
+%!assert (csaps(x,y,1,x'), y', 10*eps);
+%!assert (csaps(x',y',1,x'), y', 10*eps);
+%!assert (csaps(x',y',1,x), y, 10*eps);
+%!assert (csaps(x,[y 2*y],1,x)', [y 2*y], 10*eps);
+
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <cde...@us...> - 2012-03-27 08:21:56
|
Revision: 10074
http://octave.svn.sourceforge.net/octave/?rev=10074&view=rev
Author: cdemills
Date: 2012-03-27 08:21:45 +0000 (Tue, 27 Mar 2012)
Log Message:
-----------
- added a comment about specific part of code
Modified Paths:
--------------
trunk/octave-forge/extra/dataframe/inst/@dataframe/dataframe.m
Modified: trunk/octave-forge/extra/dataframe/inst/@dataframe/dataframe.m
===================================================================
--- trunk/octave-forge/extra/dataframe/inst/@dataframe/dataframe.m 2012-03-27 08:17:47 UTC (rev 10073)
+++ trunk/octave-forge/extra/dataframe/inst/@dataframe/dataframe.m 2012-03-27 08:21:45 UTC (rev 10074)
@@ -272,6 +272,8 @@
the_line = cellfun (@(x) sscanf (x, "%f", locales), dummy, \
'UniformOutput', false);
else
+ %# this code require a patch to src/file-io.cc in main
+ %# Octave tree
the_line = sscanf (dummy, "%f", locales);
the_line = cellfun (@(x) x{1}, the_line, 'UniformOutput', false);
endif
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.
