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: <i7...@us...> - 2012-06-12 14:12:09
|
Revision: 10623
http://octave.svn.sourceforge.net/octave/?rev=10623&view=rev
Author: i7tiol
Date: 2012-06-12 14:11:58 +0000 (Tue, 12 Jun 2012)
Log Message:
-----------
Add missing NEWS entry for Makefile.
Modified Paths:
--------------
trunk/octave-forge/main/optim/NEWS
Modified: trunk/octave-forge/main/optim/NEWS
===================================================================
--- trunk/octave-forge/main/optim/NEWS 2012-06-12 14:08:57 UTC (rev 10622)
+++ trunk/octave-forge/main/optim/NEWS 2012-06-12 14:11:58 UTC (rev 10623)
@@ -20,7 +20,10 @@
honour an index of constraints has been removed, the respective
feature is still available by setting some options.
+ ** Makefile fixed to work with non-standard linker options e.g on
+ Apple.
+
Summary of important user-visible changes for optim 1.1.0:
-------------------------------------------------------------------
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <i7...@us...> - 2012-06-12 14:09:04
|
Revision: 10622
http://octave.svn.sourceforge.net/octave/?rev=10622&view=rev
Author: i7tiol
Date: 2012-06-12 14:08:57 +0000 (Tue, 12 Jun 2012)
Log Message:
-----------
Completed optimization suite with new backend lm_feasible.
Modified Paths:
--------------
trunk/octave-forge/main/optim/DESCRIPTION
trunk/octave-forge/main/optim/NEWS
trunk/octave-forge/main/optim/inst/nonlin_curvefit.m
trunk/octave-forge/main/optim/inst/nonlin_min.m
trunk/octave-forge/main/optim/inst/nonlin_residmin.m
trunk/octave-forge/main/optim/inst/private/__lm_svd__.m
trunk/octave-forge/main/optim/inst/private/__nonlin_residmin__.m
Added Paths:
-----------
trunk/octave-forge/main/optim/inst/private/__lm_feasible__.m
Modified: trunk/octave-forge/main/optim/DESCRIPTION
===================================================================
--- trunk/octave-forge/main/optim/DESCRIPTION 2012-06-12 13:42:55 UTC (rev 10621)
+++ trunk/octave-forge/main/optim/DESCRIPTION 2012-06-12 14:08:57 UTC (rev 10622)
@@ -1,6 +1,6 @@
Name: Optim
-Version: 1.1.1
-Date: 2012-06-07
+Version: 1.2.0
+Date: 2012-06-12
Author: various authors
Maintainer: Octave-Forge community <oct...@li...>
Title: Optimization.
Modified: trunk/octave-forge/main/optim/NEWS
===================================================================
--- trunk/octave-forge/main/optim/NEWS 2012-06-12 13:42:55 UTC (rev 10621)
+++ trunk/octave-forge/main/optim/NEWS 2012-06-12 14:08:57 UTC (rev 10622)
@@ -1,3 +1,26 @@
+Summary of important user-visible changes for optim 1.2.0:
+-------------------------------------------------------------------
+
+ ** Together with the new backend "lm_feasible" there is now a
+ complete suite of backends for optimization with linear and
+ general equality and inequality constraints, for scalar valued
+ objective functions and for array valued model function, which
+ features, a.o., honouring of constraints throughout optimization
+ and handling of structure-based parameters. The respective user
+ functions (frontends) are
+
+ nonlin_min nonlin_residmin nonlin_curvefit
+
+ together with a user function for statistics
+
+ residmin_stat
+
+ ** The requirement of nonlin_min, nonlin_residmin, and
+ nonlin_curvefit for the general constraint functions being able to
+ honour an index of constraints has been removed, the respective
+ feature is still available by setting some options.
+
+
Summary of important user-visible changes for optim 1.1.0:
-------------------------------------------------------------------
Modified: trunk/octave-forge/main/optim/inst/nonlin_curvefit.m
===================================================================
--- trunk/octave-forge/main/optim/inst/nonlin_curvefit.m 2012-06-12 13:42:55 UTC (rev 10621)
+++ trunk/octave-forge/main/optim/inst/nonlin_curvefit.m 2012-06-12 14:08:57 UTC (rev 10622)
@@ -81,3 +81,21 @@
endif
endfunction
+
+%!demo
+%! ## Example for linear inequality constraints
+%! ## (see also the same example in 'demo nonlin_residmin')
+%!
+%! ## independents and observations
+%! indep = 1:5;
+%! obs = [1, 2, 4, 7, 14];
+%! ## model function:
+%! f = @ (p, x) p(1) * exp (p(2) * x);
+%! ## initial values:
+%! init = [.25; .25];
+%! ## linear constraints, A.' * parametervector + B >= 0
+%! A = [1; -1]; B = 0; # p(1) >= p(2);
+%! settings = optimset ("inequc", {A, B});
+%!
+%! ## start optimization
+%! [p, model_values, cvg, outp] = nonlin_curvefit (f, init, indep, obs, settings)
Modified: trunk/octave-forge/main/optim/inst/nonlin_min.m
===================================================================
--- trunk/octave-forge/main/optim/inst/nonlin_min.m 2012-06-12 13:42:55 UTC (rev 10621)
+++ trunk/octave-forge/main/optim/inst/nonlin_min.m 2012-06-12 14:08:57 UTC (rev 10622)
@@ -22,10 +22,6 @@
## interface; any additionally needed constants can be supplied by
## wrapping the user functions into anonymous functions.
##
-## At the moment there is only a stochastic backend implemented, so that
-## some settings in the frontend have no use. A backend for
-## feasible-path sequential quadratic programming is being worked on.
-##
## The following description applies to usage with vector-based
## parameter handling. Differences in usage for structure-based
## parameter handling will be explained in a separate section below.
@@ -52,14 +48,14 @@
## currently can be @code{0} (maximum number of iterations exceeded),
## @code{1} (fixed number of iterations completed, e.g. in stochastic
## optimizers), @code{2} (parameter change less than specified precision
-## in two consecutive iterations), or @code{3} (improvement in objective
-## function less than specified).
+## in two consecutive iterations), @code{3} (improvement in objective
+## function less than specified), or @code{-4} (algorithm got stuck).
##
## @var{settings}:
##
-## @code{Algorithm}: String specifying the backend. Currently there is
-## no default, since only a stochastic backend is implemented
-## (@code{"siman"}). It is described in a separate section below.
+## @code{Algorithm}: String specifying the backend. Currently available
+## are @{"lm_feasible"} (default) and @code{"siman"}. They are described
+## in separate sections below.
##
## @code{objf_grad}: Function computing the gradient of the objective
## function with respect to the parameters, assuming residuals are
@@ -138,24 +134,29 @@
## @code{h (p[, idx]) >= 0}; @code{p} is the column vector of optimized
## paraters and the optional argument @code{idx} is a logical index.
## @code{h} has to return the values of all constraints if @code{idx} is
-## not given, and has to return only the indexed constraints if
+## not given. It may choose to return only the indexed constraints if
## @code{idx} is given (so computation of the other constraints can be
-## spared). In gradient determination, this function may be called with
-## an informational third argument, whose content depends on the
-## function for gradient determination. If a second entry for general
-## inequality constraints is given, it must be a function computing the
-## jacobian of the constraints with respect to the parameters. For this
-## function, the description of @code{dfdp} above applies, except that
-## it is called with 3 arguments since it has an additional argument
-## @code{idx} --- a logical index --- at second position, indicating
-## which rows of the jacobian must be returned, and except that the
-## default function calls @code{h} with 3 arguments, since the argument
-## @code{idx} is also supplied. Note that specifying linear constraints
-## as general constraints will generally waste performance, even if
-## further, non-linear, general constraints are also specified.
+## spared); in this case, the additional setting @code{inequc_f_idx} has
+## to be set to @code{true}. In gradient determination, this function
+## may be called with an informational third argument, whose content
+## depends on the function for gradient determination. If a second entry
+## for general inequality constraints is given, it must be a function
+## computing the jacobian of the constraints with respect to the
+## parameters. For this function, the description of @code{dfdp} above
+## applies, with 2 exceptions: 1) it is called with 3 arguments since it
+## has an additional argument @code{idx} --- a logical index --- at
+## second position, indicating which rows of the jacobian must be
+## returned (if the function chooses to return only indexed rows, the
+## additional setting @code{inequc_df_idx} has to be set to
+## @code{true}). 2) the default jacobian function calls @code{h} with 3
+## arguments, since the argument @code{idx} is also supplied. Note that
+## specifying linear constraints as general constraints will generally
+## waste performance, even if further, non-linear, general constraints
+## are also specified.
##
## @code{equc}: Equality constraints. Specified the same way as
-## inequality constraints (see @code{inequc}).
+## inequality constraints (see @code{inequc}). The respective additional
+## settings are named @code{equc_f_idx} and @code{equc_df_idx}.
##
## @code{cpiv}: Function for complementary pivoting, usable in
## algorithms for constraints. Default: @ cpiv_bard. Only the default
@@ -242,8 +243,30 @@
## parameters must contain arrays with dimensions reshapable to those of
## the respective parameters.
##
-## Description of backends (currently only one)
+## Description of backends
##
+## "lm_feasible"
+##
+## A Levenberg/Marquardt-like optimizer, attempting to honour
+## constraints throughout the course of optimization. This means that
+## the initial parameters must not violate constraints (to find an
+## initial feasible set of parameters, e.g. Octaves @code{sqp} can be
+## used, by specifying an objective function which is constant or which
+## returns the quadratic distance to the initial values). If the
+## constraints need only be honoured in the result of the optimization,
+## Octaves @code{sqp} may be preferable. The Hessian is either supplied
+## by the user or is approximated by the BFGS algorithm.
+##
+## Returned value @var{cvg} will be @code{2} or @code{3} for success and
+## @code{0} or @code{-4} for failure (see above for meaning).
+##
+## Backend-specific defaults are: @code{MaxIter}: 20, @code{fract_prec}:
+## @code{zeros (size (parameters))}, @code{max_fract_change}: @code{Inf}
+## for all parameters.
+##
+## Interpretation of @code{Display}: if set to @code{"iter"}, currently
+## only information on applying @code{max_fract_change} is printed.
+##
## "siman"
##
## A simulated annealing (stochastic) optimizer, changing all parameters
@@ -335,10 +358,14 @@
"fixed", [], \
"inequc", [], \
"equc", [], \
+ "inequc_f_idx", false, \
+ "inequc_df_idx", false, \
+ "equc_f_idx", false, \
+ "equc_df_idx", false, \
"TolFun", stol_default, \
"MaxIter", [], \
"Display", "off", \
- "Algorithm", [], \ # set reasonable default, once available
+ "Algorithm", "lm_feasible", \
"T_init", .01, \
"T_min", 1.0e-5, \
"mu_T", 1.005, \
@@ -351,10 +378,12 @@
return;
endif
- assign = @ assign; # Is this faster in repeated calls?
+ if (nargin < 2 || nargin > 3)
+ print_usage ();
+ endif
- if (nargin != 3)
- error ("incorrect number of arguments");
+ if (nargin == 2)
+ settings = struct ();
endif
if (ischar (f))
@@ -418,13 +447,36 @@
mc_struct = isstruct (mc);
emc_struct = isstruct (emc);
- ## correct "_pstruct" settings if functions are not supplied
+ ## correct "_pstruct" settings if functions are not supplied, handle
+ ## constraint functions not honoring indices
if (isempty (dfdp)) dfdp_pstruct = false; endif
if (isempty (hessian)) hessian_pstruct = false; endif
- if (isempty (f_genicstr)) f_inequc_pstruct = false; endif
- if (isempty (f_genecstr)) f_equc_pstruct = false; endif
- if (! user_df_gencstr) df_inequc_pstruct = false; endif
- if (! user_df_genecstr) df_equc_pstruct = false; endif
+ if (isempty (f_genicstr))
+ f_inequc_pstruct = false;
+ elseif (! optimget (settings, "inequc_f_idx", false))
+ f_genicstr = @ (p, varargin) apply_idx_if_given \
+ (f_genicstr (p, varargin{:}), varargin{:});
+ endif
+ if (isempty (f_genecstr))
+ f_equc_pstruct = false;
+ elseif (! optimget (settings, "equc_f_idx", false))
+ f_genecstr = @ (p, varargin) apply_idx_if_given \
+ (f_genecstr (p, varargin{:}), varargin{:});
+ endif
+ if (user_df_gencstr)
+ if (! optimget (settings, "inequc_df_idx", false))
+ df_gencstr = @ (varargin) df_gencstr (varargin{:})(varargin{2}, :);
+ endif
+ else
+ df_inequc_pstruct = false;
+ endif
+ if (user_df_genecstr)
+ if (! optimget (settings, "equc_df_idx", false))
+ df_genecstr = @ (varargin) df_genecstr (varargin{:})(varargin{2}, :);
+ endif
+ else
+ df_equc_pstruct = false;
+ endif
## some settings require a parameter order
if (pin_struct || ! isempty (pconf) || f_inequc_pstruct || \
@@ -988,7 +1040,6 @@
endif
hook.Display = optimget (settings, "Display", "off");
hook.testing = optimget (settings, "debug", false);
- hook.new_s = optimget (settings, "lm_svd_feasible_alt_s", false);
hook.siman.T_init = optimget (settings, "T_init", .01);
hook.siman.T_min = optimget (settings, "T_min", 1.0e-5);
hook.siman.mu_T = optimget (settings, "mu_T", 1.005);
@@ -999,10 +1050,7 @@
optimget (settings, "trace_steps", false);
hook.siman_log = \
optimget (settings, "siman_log", false);
- backend = optimget (settings, "Algorithm");
- if (isempty (backend))
- error ("At the moment, there is no default algorithm, it must be specified.");
- endif
+ backend = optimget (settings, "Algorithm", "lm_feasible");
backend = map_matlab_algorithm_names (backend);
backend = map_backend (backend);
@@ -1246,7 +1294,7 @@
## passed function for complementary pivoting
## hook.cpiv = cpiv; # set before
- ## passed value of residual function for initial parameters
+ ## passed value of objective function for initial parameters
hook.f_pin = f_pin;
## passed options
@@ -1283,12 +1331,12 @@
function backend = map_backend (backend)
switch (backend)
- case "sqp_infeasible"
- backend = "__sqp__";
- case "sqp"
- backend = "__sqp__";
- case "sqp_feasible"
- backend = "__sqp_feasible__";
+ ## case "sqp_infeasible"
+ ## backend = "__sqp__";
+ ## case "sqp"
+ ## backend = "__sqp__";
+ case "lm_feasible"
+ backend = "__lm_feasible__";
case "siman"
backend = "__siman__";
otherwise
@@ -1336,7 +1384,24 @@
endfunction
+function ret = apply_idx_if_given (ret, varargin)
+
+ if (nargin > 1)
+ ret = ret(varargin{1});
+ endif
+
+endfunction
+
%!demo
+%! ## Example for default optimization (Levenberg/Marquardt with
+%! ## BFGS), one non-linear equality constraint. Constrained optimum is
+%! ## at p = [0; 1].
+%! objective_function = @ (p) p(1)^2 + p(2)^2;
+%! pin = [-2; 5];
+%! constraint_function = @ (p) p(1)^2 + 1 - p(2);
+%! [p, objf, cvg, outp] = nonlin_min (objective_function, pin, optimset ("equc", {constraint_function}))
+
+%!demo
%! ## Example for simulated annealing, two parameters, "trace_steps"
%! ## is true;
%! t_init = .2;
Modified: trunk/octave-forge/main/optim/inst/nonlin_residmin.m
===================================================================
--- trunk/octave-forge/main/optim/inst/nonlin_residmin.m 2012-06-12 13:42:55 UTC (rev 10621)
+++ trunk/octave-forge/main/optim/inst/nonlin_residmin.m 2012-06-12 14:08:57 UTC (rev 10622)
@@ -129,21 +129,25 @@
## @code{h (p[, idx]) >= 0}; @code{p} is the column vector of optimized
## paraters and the optional argument @code{idx} is a logical index.
## @code{h} has to return the values of all constraints if @code{idx} is
-## not given, and has to return only the indexed constraints if
+## not given. It may choose to return only the indexed constraints if
## @code{idx} is given (so computation of the other constraints can be
-## spared). In gradient determination, this function may be called with
-## an informational third argument, whose content depends on the
-## function for gradient determination. If a second entry for general
-## inequality constraints is given, it must be a function computing the
-## jacobian of the constraints with respect to the parameters. For this
-## function, the description of @code{dfdp} above applies, except that
-## it is called with 3 arguments since it has an additional argument
-## @code{idx} --- a logical index --- at second position, indicating
-## which rows of the jacobian must be returned, and except that the
-## default function calls @code{h} with 3 arguments, since the argument
-## @code{idx} is also supplied. Note that specifying linear constraints
-## as general constraints will generally waste performance, even if
-## further, non-linear, general constraints are also specified.
+## spared); in this case, the additional setting @code{inequc_f_idx} has
+## to be set to @code{true}. In gradient determination, this function
+## may be called with an informational third argument, whose content
+## depends on the function for gradient determination. If a second entry
+## for general inequality constraints is given, it must be a function
+## computing the jacobian of the constraints with respect to the
+## parameters. For this function, the description of @code{dfdp} above
+## applies, with 2 exceptions: 1) it is called with 3 arguments since it
+## has an additional argument @code{idx} --- a logical index --- at
+## second position, indicating which rows of the jacobian must be
+## returned (if the function chooses to return only indexed rows, the
+## additional setting @code{inequc_df_idx} has to be set to
+## @code{true}). 2) the default jacobian function calls @code{h} with 3
+## arguments, since the argument @code{idx} is also supplied. Note that
+## specifying linear constraints as general constraints will generally
+## waste performance, even if further, non-linear, general constraints
+## are also specified.
##
## @code{equc}: Equality constraints. Specified the same way as
## inequality constraints (see @code{inequc}).
@@ -281,3 +285,20 @@
[p, resid, cvg, outp] = __nonlin_residmin__ (varargin{:});
endfunction
+
+%!demo
+%! ## Example for linear inequality constraints
+%! ## (see also the same example in 'demo nonlin_curvefit')
+%!
+%! ## independents
+%! indep = 1:5;
+%! ## residual function:
+%! f = @ (p) p(1) * exp (p(2) * indep) - [1, 2, 4, 7, 14];
+%! ## initial values:
+%! init = [.25; .25];
+%! ## linear constraints, A.' * parametervector + B >= 0
+%! A = [1; -1]; B = 0; # p(1) >= p(2);
+%! settings = optimset ("inequc", {A, B});
+%!
+%! ## start optimization
+%! [p, residuals, cvg, outp] = nonlin_residmin (f, init, settings)
Added: trunk/octave-forge/main/optim/inst/private/__lm_feasible__.m
===================================================================
--- trunk/octave-forge/main/optim/inst/private/__lm_feasible__.m (rev 0)
+++ trunk/octave-forge/main/optim/inst/private/__lm_feasible__.m 2012-06-12 14:08:57 UTC (rev 10622)
@@ -0,0 +1,505 @@
+## Copyright (C) 2012 Olaf Till <i7...@t-...>
+##
+## 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/>.
+
+function [p_res, objf, cvg, outp] = __lm_feasible__ (f, pin, hook)
+
+ ## some backend specific defaults
+ fract_prec_default = 0;
+ max_fract_step_default = Inf;
+
+ ## needed for some anonymous functions
+ if (exist ("ifelse") != 5)
+ ifelse = @ scalar_ifelse;
+ endif
+
+ n = length (pin);
+
+ ## passed constraints
+ mc = hook.mc; # matrix of linear constraints
+ vc = hook.vc; # vector of linear constraints
+ f_cstr = hook.f_cstr; # function of all constraints
+ df_cstr = hook.df_cstr; # function of derivatives of all constraints
+ n_gencstr = hook.n_gencstr; # number of non-linear constraints
+ eq_idx = hook.eq_idx; # logical index of equality constraints in all
+ # constraints
+ lbound = hook.lbound; # bounds, subset of linear inequality
+ ubound = hook.ubound; # constraints in mc and vc
+
+ ## passed values of constraints for initial parameters
+ pin_cstr = hook.pin_cstr;
+
+ ## passed return value of f for initial parameters
+ f_pin = hook.f_pin;
+
+ ## passed function for gradient of objective function
+ grad_f = hook.dfdp;
+
+ ## passed function for hessian of objective function
+ if (isempty (hessian = hook.hessian))
+ user_hessian = false;
+ A = eye (n);
+ else
+ user_hessian = true;
+ endif
+
+ ## passed function for complementary pivoting
+ cpiv = hook.cpiv;
+
+ ## passed options
+ ftol = hook.TolFun;
+ if (isempty (niter = hook.MaxIter)) niter = 20; endif
+ fixed = hook.fixed;
+ maxstep = hook.max_fract_change;
+ maxstep(isna (maxstep)) = max_fract_step_default;
+ pprec = hook.fract_prec;
+ pprec(isna (pprec)) = fract_prec_default;
+ verbose = strcmp (hook.Display, "iter");
+
+ ## some useful variables derived from passed variables
+ n_lcstr = size (vc, 1);
+ have_constraints_except_bounds = \
+ n_lcstr + n_gencstr > \
+ sum (lbound != -Inf) + sum (ubound != Inf);
+ ac_idx = true (n_lcstr + n_gencstr, 1); # index of all constraints
+ nc_idx = false (n_lcstr + n_gencstr, 1); # none of all constraints
+ gc_idx = cat (1, false (n_lcstr, 1), true (n_gencstr, 1)); # gen. constr.
+
+ nz = 20 * eps; # This is arbitrary. Accuracy of equality constraints.
+
+ ## backend-specific checking of options and constraints
+ ##
+ if (any (pin < lbound | pin > ubound) ||
+ any (pin_cstr.inequ.lin_except_bounds < 0) ||
+ any (pin_cstr.inequ.gen < 0) ||
+ any (abs (pin_cstr.equ.lin)) >= nz ||
+ any (abs (pin_cstr.equ.gen)) >= nz)
+ error ("Initial parameters violate constraints.");
+ endif
+ ##
+ idx = lbound == ubound;
+ if (any (idx))
+ warning ("lower and upper bounds identical for some parameters, fixing the respective parameters");
+ fixed(idx) = true;
+ endif
+ if (all (fixed))
+ error ("no free parameters");
+ endif
+ if (n_gencstr > 0 && any (! isinf (maxstep)))
+ warning ("setting both a maximum fractional step change of parameters and general constraints may result in inefficiency and failure");
+ endif
+
+ ## fill constant fields of hook for derivative-functions; some fields
+ ## may be backend-specific
+ dfdp_hook.fixed = fixed; # this may be handled by the frontend, but
+ # the backend still may add to it
+
+ ## set up for iterations
+ p = pbest = pin;
+ vf = fbest = f_pin;
+ iter = 0;
+ done = false;
+ ll = 1;
+ ltab = [.1, 1, 1e2, 1e4, 1e6];
+ chgprev = Inf (n, 1);
+ df = [];
+ c_act = false (n, 1);
+ dca = zeros (n, 0);
+
+ while (! done)
+
+ iter++;
+
+ ## gradient of objective function
+ old_df = df;
+ df = grad_f (p, setfield (dfdp_hook, "f", f))(:);
+
+ ## constraints, preparation of some constants
+ v_cstr = f_cstr (p, ac_idx);
+ old_c_act = c_act;
+ old_dca = dca;
+ c_act = v_cstr < nz | eq_idx; # index of active constraints
+ if (any (c_act))
+
+ if (n_gencstr)
+ ## full gradient is needed later
+ dct = df_cstr (p, ac_idx, setfield (dfdp_hook, "f", v_cstr));
+ dct(:, fixed) = 0; # for user supplied dfdp; necessary?
+ dcat = dct(c_act, :);
+ else
+ dcat = df_cstr (p, c_act, setfield (dfdp_hook, "f", v_cstr));
+ dcat(:, fixed) = 0; # for user supplied dfdp; necessary?
+ endif
+
+ dca = dcat.';
+
+ a_eq_idx = eq_idx(c_act);
+
+ else
+
+ dca = zeros (n, 0);
+
+ endif
+
+ ## hessian of objectiv function
+ if (user_hessian)
+
+ A = hessian (p);
+ idx = isnan (A);
+ A(idx) = A.'(idx);
+ if (any (isnan (A(:))))
+ error ("some second derivatives undefined by user function");
+ endif
+ if (! isreal (A))
+ error ("second derivatives given by user function not real");
+ endif
+ if (! issymmetric (A))
+ error ("Hessian returned by user function not symmetric");
+ endif
+
+ elseif (iter > 1)
+
+ if (any (chg))
+
+ ## approximate Hessian of Lagrangian
+
+ ## I wonder if this hassle here and above with accounting for
+ ## changing active sets is indeed better than just approximating
+ ## the Hessian only of the objective function.
+ ##
+ ## index, over all constraints, of constraints active both
+ ## previously and currently
+ s_c_act = old_c_act & c_act;
+ ## index, over currently active constraints, of constraints
+ ## active both previously and currently
+ id_new = s_c_act(c_act);
+ ## index, over previously active constraints, of constraints
+ ## active both previously and currently
+ id_old = s_c_act(old_c_act);
+ ## gradients of currently active constraints which were also
+ ## active previously
+ dca_new_id = dca(:, id_new);
+ ## gradients of previously active constraints which are also
+ ## active currently
+ dca_old_id = old_dca(:, id_old);
+ ## index, over constraints active both previously and currently,
+ ## of (old) non-zero multipliers (bidx set below previously)
+ bidx_old_id = bidx(id_old);
+ ## index, over (old) non-zero multipliers, of constraints active
+ ## both previously and currently (bidx set below previously)
+ old_l_idx = id_old(bidx);
+
+ ## difference of derivatives of new and old active constraints,
+ ## multiplied by multipliers, as used for BFGS update (lb set
+ ## below previously)
+ dch = (dca_new_id(:, bidx_old_id) - \
+ dca_old_id(:, bidx_old_id)) * \
+ lb(old_l_idx);
+
+ y = df - old_df - dch;
+
+ ## Damped BFGS according to Nocedal & Wright, 2nd edition,
+ ## procedure 18.2.
+ chgt = chg.';
+ sAs = chgt * A * chg;
+ cy = chgt * y;
+ if (cy >= .2 * sAs)
+ th = 1;
+ else
+ if ((den1 = sAs - cy) == 0)
+ cvg = -4;
+ break;
+ endif
+ th = .8 * sAs / den1;
+ endif
+ Ac = A * chg;
+ r = th * y + (1 - th) * Ac;
+
+ if ((den2 = chgt * r) == 0 || sAs == 0)
+ cvg = -4;
+ break;
+ endif
+ A += r * r.' / den2 - Ac * Ac.' / sAs;
+
+ endif
+
+ endif
+
+ ## Inverse scaled decomposition A = G * (1 ./ L) * G.'
+ ##
+ ## make matrix Binv for scaling
+ Binv = diag (A);
+ nidx = ! (idx = Binv == 0);
+ Binv(nidx) = 1 ./ sqrt (abs (Binv(nidx)));
+ Binv(idx) = 1;
+ Binv = diag (Binv);
+ ## eigendecomposition of scaled A
+ [V, L] = eig (Binv * A * Binv);
+ L = diag (L);
+ ## A is symmetric, so V and L are real, delete any imaginary parts,
+ ## which might occur due to inaccuracy
+ V = real (V);
+ L = real (L);
+ ##
+ nminL = - min (L) * 1.1 / ltab(1);
+ G = Binv * V;
+
+ ## Levenberg/Marquardt
+ fgoal = (1 - ftol) * vf;
+ for l = ltab
+
+ ll = max (ll, nminL);
+ l = max (1e-7, ll * l);
+
+ R = G * diag (1 ./ (L + l)) * G.';
+
+ ## step computation
+ if (any (c_act))
+
+ ## some constraints are active, quadratic programming
+
+ tp = dcat * R;
+ [lb, bidx, ridx, tbl] = cpiv (- tp * df, tp * dca, a_eq_idx);
+ chg = R * (dca(:, bidx) * lb - df); # step direction
+
+ ## indices for different types of constraints
+ c_inact = ! c_act; # inactive constraints
+ c_binding = c_unbinding = nc_idx;
+ c_binding(c_act) = bidx; # constraints selected binding
+ c_unbinding(c_act) = ridx; # constraints unselected binding
+ c_nonbinding = c_act & ! (c_binding | c_unbinding); #
+ #constraints selected non-binding
+
+ else
+
+ ## no constraints are active, chg is the Levenberg/Marquardt step
+
+ chg = - R * df; # step direction
+
+ lb = zeros (0, 1);
+ bidx = false (0, 1);
+
+ ## indices for different types of constraints (meaning see above)
+ c_inact = ac_idx;
+ c_binding = nc_idx;
+ c_unbinding = nc_idx;
+ c_nonbinding = nc_idx;
+
+ endif
+
+ ## apply inactive and non-binding constraints to step width
+ ##
+ ## linear constraints
+ k = 1;
+ c_tp = c_inact(1:n_lcstr);
+ mcit = mc(:, c_tp).';
+ vci = vc(c_tp);
+ hstep = mcit * chg;
+ idx = hstep < 0;
+ if (any (idx))
+ k = min (1, min (- (vci(idx) + mcit(idx, :) * p) ./ \
+ hstep(idx)));
+ endif
+ ##
+ ## general constraints
+ if (n_gencstr)
+ c_tp = gc_idx & (c_nonbinding | c_inact);
+ if (any (c_tp) && any (f_cstr (p + k * chg, c_tp) < 0))
+ [k, fval, info] = \
+ fzero (@ (x) min (cat (1, \
+ f_cstr (p + x * chg, c_tp), \
+ k - x, \
+ ifelse (x < 0, -Inf, Inf))), \
+ 0);
+ if (info != 1 || abs (fval) >= nz)
+ error ("could not find stepwidth to satisfy inactive and non-binding general inequality constraints");
+ endif
+ endif
+ endif
+ ##
+ chg = k * chg;
+
+ ## if necessary, regain binding constraints and one of the
+ ## possibly active previously inactive or non-binding constraints
+ if (any (gc_idx & c_binding)) # none selected binding => none
+ # unselected binding
+ ptp1 = p + chg;
+
+ tp = true;
+ nt_nosuc = true;
+ lim = 20;
+ while (nt_nosuc && lim >= 0)
+ ## we keep d_p.' * inv (R) * d_p minimal in each step of the
+ ## inner loop
+ c_tp0 = c_inact | c_nonbinding;
+ c_tp1 = c_inact | (gc_idx & c_nonbinding);
+ btbl = tbl(bidx, bidx);
+ c_tp2 = c_binding;
+ ## once (any(tp)==false), it would not get true again even
+ ## with the following assignment
+ if (any (tp) && \
+ any (tp = f_cstr (ptp1, c_tp1) < nz))
+ ## keep only the first true entry in tp
+ tp(tp) = logical (cat (1, 1, zeros (sum (tp) - 1, 1)));
+ ## supplement binding index with one (the first) getting
+ ## binding in c_tp1
+ c_tp2(c_tp1) = tp;
+ ## gradient of this added constraint
+ caddt = dct(c_tp2 & ! c_binding, :);
+ cadd = caddt.';
+ C = dct(c_binding, :) * R * cadd;
+ Ct = C.';
+ T = [btbl, btbl * C; \
+ -Ct * btbl, caddt * R * cadd - Ct * btbl * C];
+ btbl = gjp (T, size (T, 1));
+ endif
+ dcbt = dct(c_tp2, :);
+ mfc = - R * dcbt.' * btbl;
+
+ ptp2 = ptp1;
+ nt_niter = nt_niter_start = 100;
+ while (nt_nosuc && nt_niter >= 0)
+ hv = f_cstr (ptp2, c_tp2);
+ if (all (abs (hv) < nz))
+ nt_nosuc = false;
+ chg = ptp2 - p;
+ else
+ ptp2 = ptp2 + mfc * hv; # step should be zero for each
+ # component for which the parameter is
+ # "fixed"
+ endif
+ nt_niter--;
+ endwhile
+
+ if (nt_nosuc || \
+ any (abs (chg) > abs (p .* maxstep)) || \
+ any (f_cstr (ptp2, c_tp0) < -nz))
+ ## if (nt_nosuc), regaining did not converge, else,
+ ## regaining violated type 3 and 4.
+ nt_nosuc = true;
+ ptp1 = (p + ptp1) / 2;
+ endif
+ if (! nt_nosuc && \
+ any ((tp = f_cstr (ptp2, c_unbinding)) < 0))
+ [discarded, id] = min(tp);
+ tid = find (ridx);
+ id = tid(id); # index within active constraints
+ unsuccessful_exchange = false;
+ if (abs (tbl(id, id)) < nz) # Bard: not absolute value
+ ## exchange this unselected binding constraint against a
+ ## binding constraint, but not against an equality
+ ## constraint
+ tbidx = bidx & ! a_eq_idx;
+ if (! any (tbidx))
+ unsuccessful_exchange = true;
+ else
+ [discarded, idm] = max (abs (tbl(tbidx, id)));
+ tid = find (tbidx);
+ idm = tid(idm); # -> index within active constraints
+ tbl = gjp (tbl, idm);
+ bidx(idm) = false;
+ ridx(idm) = true;
+ endif
+ endif
+ if (unsuccessful_exchange)
+ ## It probably doesn't look good now; this desperate last
+ ## attempt is not in the original algortithm, since that
+ ## didn't account for equality constraints.
+ ptp1 = (p + ptp1) / 2;
+ else
+ tbl = gjp (tbl, id);
+ bidx(id) = true;
+ ridx(id) = false;
+ c_binding = nc_idx;
+ c_binding(c_act) = bidx;
+ c_unbinding = nc_idx;
+ c_unbinding(c_act) = ridx;
+ endif
+ ## regaining violated type 2 constraints
+ nt_nosuc = true;
+ endif
+ lim--;
+ endwhile
+ if (nt_nosuc)
+ error ("could not regain binding constraints");
+ endif
+ else
+ ## check the maximal stepwidth and apply as necessary
+ ochg = chg;
+ idx = ! isinf (maxstep);
+ limit = abs (maxstep(idx) .* p(idx));
+ chg(idx) = min (max (chg(idx), - limit), limit);
+ if (verbose && any (ochg != chg))
+ printf ("Change in parameter(s): %s:maximal fractional stepwidth enforced", \
+ sprintf ("%d ", find (ochg != chg)));
+ endif
+ endif # regaining
+
+ aprec = abs (pprec .* pbest);
+ if (any (abs (chg) > 0.1 * aprec)) # only worth evaluating
+ # function if there is some
+ # non-miniscule change
+ p_chg = p + chg;
+ ## since the projection method may have slightly violated
+ ## constraints due to inaccuracy, correct parameters to bounds
+ ## --- but only if no further constraints are given, otherwise
+ ## the inaccuracy in honoring them might increase by this
+ if (! have_constraints_except_bounds)
+ lidx = p_chg < lbound;
+ uidx = p_chg > ubound;
+ p_chg(lidx, 1) = lbound(lidx, 1);
+ p_chg(uidx, 1) = ubound(uidx, 1);
+ chg(lidx, 1) = p_chg(lidx, 1) - p(lidx, 1);
+ chg(uidx, 1) = p_chg(uidx, 1) - p(uidx, 1);
+ endif
+ ##
+ if (! isreal (vf_chg = f (p_chg)))
+ error ("objective function not real");
+ endif
+ if (vf_chg < fbest)
+ pbest = p_chg;
+ fbest = vf_chg;
+ endif
+ if (vf_chg <= fgoal)
+ p = p_chg;
+ vf = vf_chg;
+ break;
+ endif
+ endif
+ endfor
+
+ ll = l;
+
+ aprec = abs (pprec .* pbest);
+ if (vf_chg < eps || vf_chg > fgoal)
+ cvg = 3;
+ done = true;
+ elseif (all (abs (chg) <= aprec) && all (abs (chgprev) <= aprec))
+ cvg = 2;
+ done = true;
+ elseif (iter == niter)
+ cvg = 0;
+ done = true;
+ else
+ chgprev = chg;
+ endif
+
+ endwhile
+
+ ## return result
+ p_res = pbest;
+ objf = fbest;
+ outp.niter = iter;
+
+endfunction
Modified: trunk/octave-forge/main/optim/inst/private/__lm_svd__.m
===================================================================
--- trunk/octave-forge/main/optim/inst/private/__lm_svd__.m 2012-06-12 13:42:55 UTC (rev 10621)
+++ trunk/octave-forge/main/optim/inst/private/__lm_svd__.m 2012-06-12 14:08:57 UTC (rev 10622)
@@ -145,7 +145,7 @@
epsLlast=1;
epstab=[.1, 1, 1e2, 1e4, 1e6];
ac_idx = true (n_lcstr + n_gencstr, 1); % all constraints
- nc_idx = false (n_lcstr + n_gencstr, 1); % non of all constraints
+ nc_idx = false (n_lcstr + n_gencstr, 1); % none of all constraints
gc_idx = cat (1, false (n_lcstr, 1), true (n_gencstr, 1)); % gen. constr.
lc_idx = ~gc_idx;
@@ -159,10 +159,10 @@
# violated at start
if (any (c_act))
if (n_gencstr > 0)
+ %% full gradient is needed later
dct = df_cstr (p, ac_idx, ...
setfield (dfdp_hook, 'f', v_cstr));
dct(:, fixed) = 0; % for user supplied dfdp; necessary?
- dc = dct.';
dcat = dct(c_act, :);
else
dcat = df_cstr (p, c_act, ...
@@ -420,6 +420,12 @@
if (any(abs(chg) > 0.1*aprec))%--- % only worth evaluating
% function if there is some non-miniscule
% change
+ %% In the code of the outer loop before the inner loop pbest is
+ %% actually identical to p, since once they deviate, the outer
+ %% loop will not be repeated. Though the inner loop can still be
+ %% repeated in this case, pbest is not used in it. Since pprev
+ %% is set from pbest in the outer loop before the inner loop, it
+ %% is also identical to p up to here.
p=chg+pprev;
%% since the projection method may have slightly violated
%% constraints due to inaccuracy, correct parameters to bounds
Modified: trunk/octave-forge/main/optim/inst/private/__nonlin_residmin__.m
===================================================================
--- trunk/octave-forge/main/optim/inst/private/__nonlin_residmin__.m 2012-06-12 13:42:55 UTC (rev 10621)
+++ trunk/octave-forge/main/optim/inst/private/__nonlin_residmin__.m 2012-06-12 14:08:57 UTC (rev 10622)
@@ -69,6 +69,10 @@
"fixed", [], \
"inequc", [], \
"equc", [], \
+ "inequc_f_idx", false, \
+ "inequc_df_idx", false, \
+ "equc_f_idx", false, \
+ "equc_df_idx", false, \
"weights", [], \
"TolFun", stol_default, \
"MaxIter", [], \
@@ -142,12 +146,35 @@
mc_struct = isstruct (mc);
emc_struct = isstruct (emc);
- ## correct "_pstruct" settings if functions are not supplied
+ ## correct "_pstruct" settings if functions are not supplied, handle
+ ## constraint functions not honoring indices
if (isempty (dfdp)) dfdp_pstruct = false; endif
- if (isempty (f_genicstr)) f_inequc_pstruct = false; endif
- if (isempty (f_genecstr)) f_equc_pstruct = false; endif
- if (! user_df_gencstr) df_inequc_pstruct = false; endif
- if (! user_df_genecstr) df_equc_pstruct = false; endif
+ if (isempty (f_genicstr))
+ f_inequc_pstruct = false;
+ elseif (! optimget (settings, "inequc_f_idx", false))
+ f_genicstr = @ (p, varargin) apply_idx_if_given \
+ (f_genicstr (p, varargin{:}), varargin{:});
+ endif
+ if (isempty (f_genecstr))
+ f_equc_pstruct = false;
+ elseif (! optimget (settings, "equc_f_idx", false))
+ f_genecstr = @ (p, varargin) apply_idx_if_given \
+ (f_genecstr (p, varargin{:}), varargin{:});
+ endif
+ if (user_df_gencstr)
+ if (! optimget (settings, "inequc_df_idx", false))
+ df_gencstr = @ (varargin) df_gencstr (varargin{:})(varargin{2}, :);
+ endif
+ else
+ df_inequc_pstruct = false;
+ endif
+ if (user_df_genecstr)
+ if (! optimget (settings, "equc_df_idx", false))
+ df_genecstr = @ (varargin) df_genecstr (varargin{:})(varargin{2}, :);
+ endif
+ else
+ df_equc_pstruct = false;
+ endif
## some settings require a parameter order
if (pin_struct || ! isempty (pconf) || f_inequc_pstruct || \
@@ -995,3 +1022,11 @@
endif
endfunction
+
+function ret = apply_idx_if_given (ret, varargin)
+
+ if (nargin > 1)
+ ret = ret(varargin{1});
+ endif
+
+endfunction
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <i7...@us...> - 2012-06-12 13:43:06
|
Revision: 10621
http://octave.svn.sourceforge.net/octave/?rev=10621&view=rev
Author: i7tiol
Date: 2012-06-12 13:42:55 +0000 (Tue, 12 Jun 2012)
Log Message:
-----------
Avoid recursive variables in Makefile.
Modified Paths:
--------------
trunk/octave-forge/main/optim/src/Makefile
Modified: trunk/octave-forge/main/optim/src/Makefile
===================================================================
--- trunk/octave-forge/main/optim/src/Makefile 2012-06-12 10:39:39 UTC (rev 10620)
+++ trunk/octave-forge/main/optim/src/Makefile 2012-06-12 13:42:55 UTC (rev 10621)
@@ -1,10 +1,14 @@
MKOCTFILE ?= mkoctfile
-LAPACK_LIBS ?= $(shell $(MKOCTFILE) -p LAPACK_LIBS)
+ifndef LAPACK_LIBS
+LAPACK_LIBS := $(shell $(MKOCTFILE) -p LAPACK_LIBS)
+endif
OCTAVE_LAPACK_LIBS := $(shell $(MKOCTFILE) -p LAPACK_LIBS)
# reported necessary for Apple's VecLib framework by Carlo de Falco
# <car...@gm...>
-BLAS_LIBS ?= $(shell $(MKOCTFILE) -p BLAS_LIBS)
+ifndef BLAS_LIBS
+BLAS_LIBS := $(shell $(MKOCTFILE) -p BLAS_LIBS)
+endif
OCTAVE_BLAS_LIBS := $(shell $(MKOCTFILE) -p BLAS_LIBS)
# Passing LFLAGS, supplemented with LAPACK_LIBS and BLAS_LIBS, in the
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-06-12 10:39:49
|
Revision: 10620
http://octave.svn.sourceforge.net/octave/?rev=10620&view=rev
Author: paramaniac
Date: 2012-06-12 10:39:39 +0000 (Tue, 12 Jun 2012)
Log Message:
-----------
control: minor style fix in margin.m
Modified Paths:
--------------
trunk/octave-forge/main/control/inst/margin.m
Modified: trunk/octave-forge/main/control/inst/margin.m
===================================================================
--- trunk/octave-forge/main/control/inst/margin.m 2012-06-12 05:12:31 UTC (rev 10619)
+++ trunk/octave-forge/main/control/inst/margin.m 2012-06-12 10:39:39 UTC (rev 10620)
@@ -275,9 +275,9 @@
gamma_db = 20 * log10 (gamma);
wv = [min(w), max(w)];
- ax_vec_mag = __axis_limits__ ([reshape(w, [], 1), reshape(mag_db, [], 1)]);
+ ax_vec_mag = __axis_limits__ ([w(:), mag_db(:)]);
ax_vec_mag(1:2) = wv;
- ax_vec_pha = __axis_limits__ ([reshape(w, [], 1), reshape(pha, [], 1)]);
+ ax_vec_pha = __axis_limits__ ([w(:), pha(:)]);
ax_vec_pha(1:2) = wv;
wgm = [w_gamma, w_gamma];
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-06-12 05:12:37
|
Revision: 10619
http://octave.svn.sourceforge.net/octave/?rev=10619&view=rev
Author: paramaniac
Date: 2012-06-12 05:12:31 +0000 (Tue, 12 Jun 2012)
Log Message:
-----------
control: add the sage worksheet I used to compute mimo transfer function interconnections
Added Paths:
-----------
trunk/octave-forge/main/control/devel/mimo_transfer_function_sage_worksheet.zip
Property Changed:
----------------
trunk/octave-forge/main/control/devel/
Property changes on: trunk/octave-forge/main/control/devel
___________________________________________________________________
Deleted: svn:ignore
- *.zip
Added: trunk/octave-forge/main/control/devel/mimo_transfer_function_sage_worksheet.zip
===================================================================
(Binary files differ)
Property changes on: trunk/octave-forge/main/control/devel/mimo_transfer_function_sage_worksheet.zip
___________________________________________________________________
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: <et...@us...> - 2012-06-12 01:51:36
|
Revision: 10618
http://octave.svn.sourceforge.net/octave/?rev=10618&view=rev
Author: etienne
Date: 2012-06-12 01:51:29 +0000 (Tue, 12 Jun 2012)
Log Message:
-----------
Credits
Modified Paths:
--------------
trunk/octave-forge/main/vrml/NEWS
Modified: trunk/octave-forge/main/vrml/NEWS
===================================================================
--- trunk/octave-forge/main/vrml/NEWS 2012-06-12 01:46:12 UTC (rev 10617)
+++ trunk/octave-forge/main/vrml/NEWS 2012-06-12 01:51:29 UTC (rev 10618)
@@ -1,8 +1,11 @@
Summary of important user-visible changes for vrml 1.0.13:
-------------------------------------------------------------------
+ ** Credits: Thank you to Carnë Draug <car...@gm...> for explaining me the basics of the pkg system.
+
** The following functions are new:
+ data2vrml
checker_color
vrml_Box
vrml_Sphere
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <et...@us...> - 2012-06-12 01:46:20
|
Revision: 10617
http://octave.svn.sourceforge.net/octave/?rev=10617&view=rev
Author: etienne
Date: 2012-06-12 01:46:12 +0000 (Tue, 12 Jun 2012)
Log Message:
-----------
Be more quiet
Modified Paths:
--------------
trunk/octave-forge/main/vrml/inst/vrml_surf.m
Modified: trunk/octave-forge/main/vrml/inst/vrml_surf.m
===================================================================
--- trunk/octave-forge/main/vrml/inst/vrml_surf.m 2012-06-12 01:45:45 UTC (rev 10616)
+++ trunk/octave-forge/main/vrml/inst/vrml_surf.m 2012-06-12 01:46:12 UTC (rev 10617)
@@ -180,7 +180,7 @@
C = C4;
if numel (size (col)) == 2 && all (size (col) == size (defaultCol)) && all (col == defaultCol)
- col = [col, 0.8*col, 0.9*col]
+ col = [col, 0.8*col, 0.9*col];
end
if numel (col) == 3
col = col(:);
@@ -194,9 +194,9 @@
sideCol = col(4:6)*ones(1,4);
elseif numel (col) == 9
col = col(:);
- topCol = col(1:3)
- botCol = col(7:9)
- sideCol = col(4:6)*ones(1,4)
+ topCol = col(1:3);
+ botCol = col(7:9);
+ sideCol = col(4:6)*ones(1,4);
end
col = ones(3, R-1, C-1);
for i=1:3
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <et...@us...> - 2012-06-12 01:45:51
|
Revision: 10616
http://octave.svn.sourceforge.net/octave/?rev=10616&view=rev
Author: etienne
Date: 2012-06-12 01:45:45 +0000 (Tue, 12 Jun 2012)
Log Message:
-----------
Add key, else pkg complains
Modified Paths:
--------------
trunk/octave-forge/main/vrml/inst/test_moving_surf.m
Modified: trunk/octave-forge/main/vrml/inst/test_moving_surf.m
===================================================================
--- trunk/octave-forge/main/vrml/inst/test_moving_surf.m 2012-06-11 16:51:02 UTC (rev 10615)
+++ trunk/octave-forge/main/vrml/inst/test_moving_surf.m 2012-06-12 01:45:45 UTC (rev 10616)
@@ -13,6 +13,10 @@
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.
+###key test_moving_surf
+##
+## Test vmesh.m
+
R = 15;
C = 23;
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <et...@us...> - 2012-06-11 16:51:12
|
Revision: 10615
http://octave.svn.sourceforge.net/octave/?rev=10615&view=rev
Author: etienne
Date: 2012-06-11 16:51:02 +0000 (Mon, 11 Jun 2012)
Log Message:
-----------
Updates and new version number
Modified Paths:
--------------
trunk/octave-forge/main/vrml/DESCRIPTION
trunk/octave-forge/main/vrml/INDEX
trunk/octave-forge/main/vrml/Makefile
trunk/octave-forge/main/vrml/doc/Makefile
Modified: trunk/octave-forge/main/vrml/DESCRIPTION
===================================================================
--- trunk/octave-forge/main/vrml/DESCRIPTION 2012-06-11 13:34:31 UTC (rev 10614)
+++ trunk/octave-forge/main/vrml/DESCRIPTION 2012-06-11 16:51:02 UTC (rev 10615)
@@ -1,6 +1,6 @@
Name: Vrml
-Version: 1.0.12
-Date: 2010-04-11
+Version: 1.0.13
+Date: 2012-06-04
Author: Etienne Grossmann <et...@eg...>
Maintainer: Etienne Grossmann <et...@eg...>
Title: VRML.
Modified: trunk/octave-forge/main/vrml/INDEX
===================================================================
--- trunk/octave-forge/main/vrml/INDEX 2012-06-11 13:34:31 UTC (rev 10614)
+++ trunk/octave-forge/main/vrml/INDEX 2012-06-11 16:51:02 UTC (rev 10615)
@@ -21,6 +21,7 @@
vrml_lines
vrml_group
Support Functions
+ data2vrml
save_vrml
best_dir
best_dir_cov
Modified: trunk/octave-forge/main/vrml/Makefile
===================================================================
--- trunk/octave-forge/main/vrml/Makefile 2012-06-11 13:34:31 UTC (rev 10614)
+++ trunk/octave-forge/main/vrml/Makefile 2012-06-11 16:51:02 UTC (rev 10615)
@@ -15,5 +15,10 @@
clean:
@for _dir in $(SUBDIRS); do \
- $(MAKE) -C $$_dir $(MAKECMDGOALS); \
+ $(MAKE) -C $$_dir clean ; \
done
+ find . \( -name \*~ -o -name \*.bak -o -name \#\* -o -name \.\#\* -\) -exec rm -f {} \;
+
+documentation:
+ make -C ./doc all
+
Modified: trunk/octave-forge/main/vrml/doc/Makefile
===================================================================
--- trunk/octave-forge/main/vrml/doc/Makefile 2012-06-11 13:34:31 UTC (rev 10614)
+++ trunk/octave-forge/main/vrml/doc/Makefile 2012-06-11 16:51:02 UTC (rev 10615)
@@ -2,13 +2,18 @@
TEX = vrml-mini-howto-1.tex
PDF = $(patsubst %.tex,%.pdf,$(TEX))
+DVIPDF = dvipdf
-all : $(PDF) html/index.html
+all : $(PDF)
%.pdf : %.tex
- latex -interaction=nonstopmode $< > /dev/null 2>&1
- latex -interaction=nonstopmode $< > /dev/null 2>&1
- $(DVIPDF) $(@:.pdf=.dvi)
+ if [ `which latex` ] && [ `which $(DVIPDF)` ] ; then \
+ latex -interaction=nonstopmode $< > /dev/null 2>&1 ; \
+ latex -interaction=nonstopmode $< > /dev/null 2>&1 ; \
+ $(DVIPDF) $(@:.pdf=.dvi) ; \
+ else \
+ echo "Lacking programs to build doc. Skipping" ; \
+ fi
# Note verbosity=0 as well as making latex2html quieter, has the side-effect
# of not including a url to the raw text, which it'll get wrong
@@ -23,7 +28,9 @@
clean:
rm -fr $(patsubst %.tex,%,$(TEX)) html *.log
- rm -f $(PDF) *~
rm -f $(patsubst %.tex,%.aux,$(TEX))
rm -f $(patsubst %.tex,%.out,$(TEX))
rm -f $(patsubst %.tex,%.dvi,$(TEX))
+
+realclean: clean
+ rm -f $(PDF) *~
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <aba...@us...> - 2012-06-11 13:34:37
|
Revision: 10614
http://octave.svn.sourceforge.net/octave/?rev=10614&view=rev
Author: abarth93
Date: 2012-06-11 13:34:31 +0000 (Mon, 11 Jun 2012)
Log Message:
-----------
Modified Paths:
--------------
trunk/octave-forge/main/octcdf/inst/ncread.m
Added Paths:
-----------
trunk/octave-forge/main/octcdf/inst/ncinfo.m
Added: trunk/octave-forge/main/octcdf/inst/ncinfo.m
===================================================================
--- trunk/octave-forge/main/octcdf/inst/ncinfo.m (rev 0)
+++ trunk/octave-forge/main/octcdf/inst/ncinfo.m 2012-06-11 13:34:31 UTC (rev 10614)
@@ -0,0 +1,57 @@
+% finfo = ncinfo(filename)
+% vinfo = ncinfo(filename,varname)
+% return information about variable varname in netCDF
+% file filename
+
+function info = ncinfo(filename,varname)
+
+nc = netcdf(filename,'r');
+if nargin == 1
+ info.Filename = filename;
+
+ variables = var(nc);
+ for i=1:length(variables)
+ info.Variables(i) = ncinfo_var(nc,filename,name(variables{i}));
+ end
+elseif nargin == 2
+ info = ncinfo_var(nc,filename,varname);
+end
+
+close(nc);
+end
+
+function vinfo = ncinfo_var(nc,filename,varname)
+
+
+nv = nc{varname};
+
+vinfo.Size = fliplr(size(nv));
+vinfo.Filename = filename;
+vinfo.Dimensions = {};
+vinfo.Name = varname;
+
+dims = fliplr(dim(nv));
+
+for i=1:length(dims)
+ vinfo.Dimensions{i} = name(dims{i});
+end
+
+
+na = att(nv);
+
+%vinfo.Attributes = [];
+
+for j=1:length(na)
+ tmp = struct();
+ nm = name(na{j});
+
+ tmp.Name = nm;
+ tmp.Value = na{j}(:);
+ vinfo.Attributes(j) = tmp;
+end
+
+vinfo.FillValue = fillval(nv);
+
+
+end
+
Modified: trunk/octave-forge/main/octcdf/inst/ncread.m
===================================================================
--- trunk/octave-forge/main/octcdf/inst/ncread.m 2012-06-11 10:31:23 UTC (rev 10613)
+++ trunk/octave-forge/main/octcdf/inst/ncread.m 2012-06-11 13:34:31 UTC (rev 10614)
@@ -7,8 +7,11 @@
nc = netcdf(filename,'r');
nv = nc{varname};
sz = size(nv); sz = sz(end:-1:1);
-nd = ndims(nv);
+% number of dimenions
+%nd = length(sz);
+nd = length(dim(nv));
+
if nargin == 2
start = ones(1,nd);
count = inf*ones(1,nd);
@@ -36,7 +39,7 @@
for i=1:nd
subsr.subs{nd-i+1} = start(i):stride(i):endi(i);
end
-start,endi
+%start,endi
x = subsref(nv,subsr);
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-06-11 10:31:34
|
Revision: 10613
http://octave.svn.sourceforge.net/octave/?rev=10613&view=rev
Author: paramaniac
Date: 2012-06-11 10:31:23 +0000 (Mon, 11 Jun 2012)
Log Message:
-----------
control-devel: minor changes in comments
Modified Paths:
--------------
trunk/octave-forge/extra/control-devel/inst/__slicot_identification__.m
Modified: trunk/octave-forge/extra/control-devel/inst/__slicot_identification__.m
===================================================================
--- trunk/octave-forge/extra/control-devel/inst/__slicot_identification__.m 2012-06-11 07:29:39 UTC (rev 10612)
+++ trunk/octave-forge/extra/control-devel/inst/__slicot_identification__.m 2012-06-11 10:31:23 UTC (rev 10613)
@@ -135,7 +135,9 @@
## perform system identification
[a, b, c, d, q, ry, s, k, x0] = slident (dat.y, dat.u, nobr, n, meth, alg, conct, ctrl, rcond, tol);
- ## L L' = Ry, e = L v, v becomes white noise with identity covariance matrix
+ ## compute noise variance matrix factor L
+ ## L L' = Ry, e = L v
+ ## v becomes white noise with identity covariance matrix
l = chol (ry, "lower");
## assemble model
@@ -169,6 +171,7 @@
## state covariance matrix Q
## output covariance matrix Ry
## state-output cross-covariance matrix S
+ ## noise variance matrix factor L
info = struct ("K", k, "Q", q, "Ry", ry, "S", s, "L", l);
endfunction
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <cd...@us...> - 2012-06-11 07:29:46
|
Revision: 10612
http://octave.svn.sourceforge.net/octave/?rev=10612&view=rev
Author: cdf
Date: 2012-06-11 07:29:39 +0000 (Mon, 11 Jun 2012)
Log Message:
-----------
bug fix in newton solver
Modified Paths:
--------------
trunk/octave-forge/extra/secs1d/inst/secs1d_dd_newton.m
Modified: trunk/octave-forge/extra/secs1d/inst/secs1d_dd_newton.m
===================================================================
--- trunk/octave-forge/extra/secs1d/inst/secs1d_dd_newton.m 2012-06-09 15:38:32 UTC (rev 10611)
+++ trunk/octave-forge/extra/secs1d/inst/secs1d_dd_newton.m 2012-06-11 07:29:39 UTC (rev 10612)
@@ -128,12 +128,12 @@
A21 = - bim1a_laplacian (x, un .* nm, 1);
A22 = bim1a_advection_diffusion (x, un, 1, 1, V) + bim1a_reaction (x, 1, p .* fact);
A23 = bim1a_reaction (x, 1, n .* fact);
- R2 = A22 * n + bim1a_rhs (x, 1, (p .* n - theta .^ 2) .* fact);
+ R2 = A22 * n + bim1a_rhs (x, 1, (- theta .^ 2) .* fact);
A31 = bim1a_laplacian (x, up .* pm, 1);
A32 = bim1a_reaction (x, 1, p .* fact);
A33 = bim1a_advection_diffusion (x, up, 1, 1, -V) + bim1a_reaction (x, 1, n .* fact);
- R3 = A33 * p + bim1a_rhs (x, 1, (p .* n - theta .^ 2) .* fact);
+ R3 = A33 * p + bim1a_rhs (x, 1, (- theta .^ 2) .* fact);
res = [R1(2:end-1); R2(2:end-1); R3(2:end-1)];
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-06-09 15:38:38
|
Revision: 10611
http://octave.svn.sourceforge.net/octave/?rev=10611&view=rev
Author: paramaniac
Date: 2012-06-09 15:38:32 +0000 (Sat, 09 Jun 2012)
Log Message:
-----------
control-devel: add modified slident routine to compute singular values only
Modified Paths:
--------------
trunk/octave-forge/extra/control-devel/src/devel_slicot_functions.cc
Added Paths:
-----------
trunk/octave-forge/extra/control-devel/src/slib01ad.cc
Modified: trunk/octave-forge/extra/control-devel/src/devel_slicot_functions.cc
===================================================================
--- trunk/octave-forge/extra/control-devel/src/devel_slicot_functions.cc 2012-06-09 09:29:34 UTC (rev 10610)
+++ trunk/octave-forge/extra/control-devel/src/devel_slicot_functions.cc 2012-06-09 15:38:32 UTC (rev 10611)
@@ -1,6 +1,6 @@
-// #include "slib01ad.cc" // preprocess the input-output data
#include "slident.cc" // system identification
#include "slib01cd.cc" // compute initial state vector
+#include "slib01ad.cc" // compute singular values
#include "slident_a.cc"
#include "slident_b.cc"
Added: trunk/octave-forge/extra/control-devel/src/slib01ad.cc
===================================================================
--- trunk/octave-forge/extra/control-devel/src/slib01ad.cc (rev 0)
+++ trunk/octave-forge/extra/control-devel/src/slib01ad.cc 2012-06-09 15:38:32 UTC (rev 10611)
@@ -0,0 +1,383 @@
+/*
+
+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 <octave/f77-fcn.h>
+#include <octave/Cell.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);
+}
+
+// PKG_ADD: autoload ("slib01ad", "devel_slicot_functions.oct");
+DEFUN_DLD (slib01ad, 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 != 10)
+ {
+ print_usage ();
+ }
+ else
+ {
+////////////////////////////////////////////////////////////////////////////////////
+// SLICOT IB01AD - preprocess the input-output data //
+////////////////////////////////////////////////////////////////////////////////////
+
+ // arguments in
+ char meth_a;
+ char meth_b;
+ char alg;
+ char jobd;
+ char batch;
+ char conct;
+ char ctrl = 'N';
+
+ const Cell y_cell = args(0).cell_value ();
+ const Cell u_cell = args(1).cell_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 iconct = args(6).int_value ();
+ const int ictrl = args(7).int_value (); // ignored
+
+ double rcond = args(8).double_value ();
+ double tol_a = args(9).double_value ();
+
+ double tol_b = rcond;
+ double tol_c = rcond;
+
+
+ switch (imeth)
+ {
+ case 0:
+ meth_a = 'M';
+ meth_b = 'M';
+ break;
+ case 1:
+ meth_a = 'N';
+ meth_b = 'N';
+ break;
+ case 2:
+ meth_a = 'N'; // no typo here
+ meth_b = 'C';
+ 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_a == 'M')
+ jobd = 'M';
+ else // meth_a == 'N'
+ jobd = 'N'; // IB01AD.f says: This parameter is not relevant for METH = 'N'
+
+ if (iconct == 0)
+ conct = 'C';
+ else
+ conct = 'N';
+/*
+ if (ictrl == 0)
+ ctrl = 'C';
+ else
+ ctrl = 'N';
+*/
+ // m and l are equal for all experiments, checked by iddata class
+ int n_exp = y_cell.nelem (); // number of experiments
+ int m = u_cell.elem(0).columns (); // m: number of inputs
+ int l = y_cell.elem(0).columns (); // l: number of outputs
+ int nsmpl = 0; // total number of samples
+
+ // arguments out
+ int n;
+ int ldr;
+
+ if (meth_a == 'M' && jobd == 'M')
+ ldr = max (2*(m+l)*nobr, 3*m*nobr);
+ else if (meth_a == 'N' || (meth_a == '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);
+
+
+ // repeat for every experiment in the dataset
+ for (int i = 0; i < n_exp; i++)
+ {
+ if (n_exp == 1)
+ batch = 'O'; // one block only
+ else if (i == 0)
+ batch = 'F'; // first block
+ else if (i == n_exp-1)
+ batch = 'L'; // last block
+ else
+ batch = 'I'; // intermediate block
+
+ Matrix y = y_cell.elem(i).matrix_value ();
+ Matrix u = u_cell.elem(i).matrix_value ();
+
+ // y.rows == u.rows is checked by iddata class
+ // int m = u.columns (); // m: number of inputs
+ // int l = y.columns (); // l: number of outputs
+ int nsmp = y.rows (); // nsmp: number of samples in the current experiment
+ nsmpl += nsmp; // nsmpl: total number of samples of all experiments
+
+ // minimal nsmp size checked by __slicot_identification__.m
+ 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;
+
+ // workspace
+ int liwork_a;
+
+ if (meth_a == 'N') // if METH = 'N'
+ liwork_a = (m+l)*nobr;
+ else if (alg == 'F') // if METH = 'M' and ALG = 'F'
+ liwork_a = m+l;
+ else // if METH = 'M' and ALG = 'C' or 'Q'
+ liwork_a = 0;
+
+ // TODO: Handle 'k' for DWORK
+
+ int ldwork_a;
+ int ns = nsmp - 2*nobr + 1;
+
+ if (alg == 'C')
+ {
+ if (batch == 'F' || batch == 'I')
+ {
+ if (conct == 'C')
+ ldwork_a = (4*nobr-2)*(m+l);
+ else // (conct == 'N')
+ ldwork_a = 1;
+ }
+ else if (meth_a == 'M') // && (batch == 'L' || batch == 'O')
+ {
+ if (conct == 'C' && batch == 'L')
+ ldwork_a = max ((4*nobr-2)*(m+l), 5*l*nobr);
+ else if (jobd == 'M')
+ ldwork_a = max ((2*m-1)*nobr, (m+l)*nobr, 5*l*nobr);
+ else // (jobd == 'N')
+ ldwork_a = 5*l*nobr;
+ }
+ else // meth_b == 'N' && (batch == 'L' || batch == 'O')
+ {
+ ldwork_a = 5*(m+l)*nobr + 1;
+ }
+ }
+ else if (alg == 'F')
+ {
+ if (batch != 'O' && conct == 'C')
+ ldwork_a = (m+l)*2*nobr*(m+l+3);
+ else if (batch == 'F' || batch == 'I') // && conct == 'N'
+ ldwork_a = (m+l)*2*nobr*(m+l+1);
+ else // (batch == 'L' || '0' && conct == 'N')
+ ldwork_a = (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_a = 4*(m+l)*nobr;
+ }
+ else if (ldr >= ns && batch == 'O')
+ {
+ if (meth_a == 'M')
+ ldwork_a = max (4*(m+l)*nobr, 5*l*nobr);
+ else // (meth == 'N')
+ ldwork_a = 5*(m+l)*nobr + 1;
+ }
+ else if (conct == 'C' && (batch == 'I' || batch == 'L'))
+ {
+ ldwork_a = 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_a = 6*(m+l)*nobr;
+ }
+ }
+
+ /*
+ IB01AD.f Lines 438-445
+ C FURTHER COMMENTS
+ C
+ C For ALG = 'Q', BATCH = 'O' and LDR < NS, or BATCH <> 'O', the
+ C calculations could be rather inefficient if only minimal workspace
+ C (see argument LDWORK) is provided. It is advisable to provide as
+ C much workspace as possible. Almost optimal efficiency can be
+ C obtained for LDWORK = (NS+2)*(2*(M+L)*NOBR), assuming that the
+ C cache size is large enough to accommodate R, U, Y, and DWORK.
+ */
+
+ ldwork_a = max (ldwork_a, (ns+2)*(2*(m+l)*nobr));
+
+ /*
+ 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_a, liwork_a);
+ OCTAVE_LOCAL_BUFFER (double, dwork_a, ldwork_a);
+
+ // error indicators
+ int iwarn_a = 0;
+ int info_a = 0;
+
+
+ // SLICOT routine IB01AD
+ F77_XFCN (ib01ad, IB01AD,
+ (meth_a, 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_a,
+ iwork_a,
+ dwork_a, ldwork_a,
+ iwarn_a, info_a));
+
+
+ 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: IB01AD", info_a, 2, err_msg);
+ warning_msg ("ident: IB01AD", iwarn_a, 5, warn_msg);
+ }
+
+
+ // resize
+ int rs = 2*(m+l)*nobr;
+ r.resize (rs, rs);
+
+
+ // return values
+ retval(0) = sv;
+ retval(1) = octave_value (n);
+ }
+
+ return retval;
+}
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-06-09 09:29:40
|
Revision: 10610
http://octave.svn.sourceforge.net/octave/?rev=10610&view=rev
Author: paramaniac
Date: 2012-06-09 09:29:34 +0000 (Sat, 09 Jun 2012)
Log Message:
-----------
control: remove infinite poles, rename variable
Modified Paths:
--------------
trunk/octave-forge/main/control/inst/@lti/isstable.m
trunk/octave-forge/main/control/inst/__is_stable__.m
trunk/octave-forge/main/control/inst/isstabilizable.m
Modified: trunk/octave-forge/main/control/inst/@lti/isstable.m
===================================================================
--- trunk/octave-forge/main/control/inst/@lti/isstable.m 2012-06-09 09:02:40 UTC (rev 10609)
+++ trunk/octave-forge/main/control/inst/@lti/isstable.m 2012-06-09 09:29:34 UTC (rev 10610)
@@ -31,9 +31,9 @@
print_usage ();
endif
- eigw = pole (sys);
+ pol = pole (sys);
ct = isct (sys);
- bool = __is_stable__ (eigw, ct, tol);
+ bool = __is_stable__ (pol, ct, tol);
-endfunction
\ No newline at end of file
+endfunction
Modified: trunk/octave-forge/main/control/inst/__is_stable__.m
===================================================================
--- trunk/octave-forge/main/control/inst/__is_stable__.m 2012-06-09 09:02:40 UTC (rev 10609)
+++ trunk/octave-forge/main/control/inst/__is_stable__.m 2012-06-09 09:29:34 UTC (rev 10610)
@@ -22,12 +22,12 @@
## Created: December 2010
## Version: 0.1
-function bool = __is_stable__ (eigw, ct = true, tol = 0)
+function bool = __is_stable__ (pol, ct = true, tol = 0)
if (ct) # continuous-time
- bool = all (real (eigw) < -tol*(1 + abs (eigw)));
+ bool = all (real (pol) < -tol*(1 + abs (pol)));
else # discrete-time
- bool = all (abs (eigw) < 1 - tol);
+ bool = all (abs (pol) < 1 - tol);
endif
-endfunction
\ No newline at end of file
+endfunction
Modified: trunk/octave-forge/main/control/inst/isstabilizable.m
===================================================================
--- trunk/octave-forge/main/control/inst/isstabilizable.m 2012-06-09 09:02:40 UTC (rev 10609)
+++ trunk/octave-forge/main/control/inst/isstabilizable.m 2012-06-09 09:29:34 UTC (rev 10610)
@@ -1,4 +1,4 @@
-## Copyright (C) 2009, 2010 Lukas F. Reichlin
+## Copyright (C) 2009, 2010, 2012 Lukas F. Reichlin
##
## This file is part of LTI Syncope.
##
@@ -73,7 +73,7 @@
## Author: Lukas Reichlin <luk...@gm...>
## Created: October 2009
-## Version: 0.3.1
+## Version: 0.4
function bool = isstabilizable (a, b = [], e = [], tol = [], dflg = 0)
@@ -109,7 +109,7 @@
auncont = ac(uncont_idx, uncont_idx);
## calculate poles of uncontrollable part
- eigw = eig (auncont);
+ pol = eig (auncont);
else
## controllability staircase form - output matrix c has no influence
[ac, ec, ~, ~, ~, ~, ncont] = sltg01hd (a, e, b, zeros (1, columns (a)), tol);
@@ -120,10 +120,15 @@
euncont = ec(uncont_idx, uncont_idx);
## calculate poles of uncontrollable part
- eigw = eig (auncont, euncont);
+ pol = eig (auncont, euncont);
+
+ ## remove infinite poles
+ tolinf = norm ([auncont, euncont], 2);
+ idx = find (abs (pol) < tolinf/eps);
+ pol = pol(idx);
endif
## check whether uncontrollable poles are stable
- bool = __is_stable__ (eigw, ! dflg, tol);
+ bool = __is_stable__ (pol, ! dflg, tol);
endfunction
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-06-09 09:02:46
|
Revision: 10609
http://octave.svn.sourceforge.net/octave/?rev=10609&view=rev
Author: paramaniac
Date: 2012-06-09 09:02:40 +0000 (Sat, 09 Jun 2012)
Log Message:
-----------
control: relocate test
Modified Paths:
--------------
trunk/octave-forge/main/control/inst/@lti/subsref.m
trunk/octave-forge/main/control/inst/ltimodels.m
trunk/octave-forge/main/control/inst/test_control.m
Modified: trunk/octave-forge/main/control/inst/@lti/subsref.m
===================================================================
--- trunk/octave-forge/main/control/inst/@lti/subsref.m 2012-06-09 07:37:08 UTC (rev 10608)
+++ trunk/octave-forge/main/control/inst/@lti/subsref.m 2012-06-09 09:02:40 UTC (rev 10609)
@@ -21,7 +21,7 @@
## Author: Lukas Reichlin <luk...@gm...>
## Created: September 2009
-## Version: 0.3
+## Version: 0.4
function a = subsref (a, s)
@@ -50,3 +50,11 @@
a = subsref (a, s(2:end));
endfunction
+
+
+## lti: subsref
+%!shared a
+%! s = tf ("s");
+%! G = (s+1)*s*5/(s+1)/(s^2+s+1);
+%! a = G(1,1).num{1,1}(1);
+%!assert (a, 5, 1e-4);
Modified: trunk/octave-forge/main/control/inst/ltimodels.m
===================================================================
--- trunk/octave-forge/main/control/inst/ltimodels.m 2012-06-09 07:37:08 UTC (rev 10608)
+++ trunk/octave-forge/main/control/inst/ltimodels.m 2012-06-09 09:02:40 UTC (rev 10609)
@@ -97,15 +97,7 @@
%!assert (isdt (ltisys));
-## lti: subsref
-%!shared a
-%! s = tf ("s");
-%! G = (s+1)*s*5/(s+1)/(s^2+s+1);
-%! a = G(1,1).num{1,1}(1);
-%!assert (a, 5, 1e-4);
-
-
## ==============================================================================
## TF Tests
## ==============================================================================
Modified: trunk/octave-forge/main/control/inst/test_control.m
===================================================================
--- trunk/octave-forge/main/control/inst/test_control.m 2012-06-09 07:37:08 UTC (rev 10608)
+++ trunk/octave-forge/main/control/inst/test_control.m 2012-06-09 09:02:40 UTC (rev 10609)
@@ -54,6 +54,7 @@
test @lti/plus
test @lti/prescale
test @lti/sminreal
+test @lti/subsref
test @lti/zero
## robust control
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-06-09 07:37:15
|
Revision: 10608
http://octave.svn.sourceforge.net/octave/?rev=10608&view=rev
Author: paramaniac
Date: 2012-06-09 07:37:08 +0000 (Sat, 09 Jun 2012)
Log Message:
-----------
control-devel: add draft code
Added Paths:
-----------
trunk/octave-forge/extra/control-devel/devel/arx_siso.m
Added: trunk/octave-forge/extra/control-devel/devel/arx_siso.m
===================================================================
--- trunk/octave-forge/extra/control-devel/devel/arx_siso.m (rev 0)
+++ trunk/octave-forge/extra/control-devel/devel/arx_siso.m 2012-06-09 07:37:08 UTC (rev 10608)
@@ -0,0 +1,229 @@
+## Copyright (C) 2012 Lukas F. Reichlin
+##
+## This file is part of LTI Syncope.
+##
+## LTI Syncope is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## LTI Syncope is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{sys} =} arx (@var{dat}, @var{na}, @var{nb})
+## ARX
+## @end deftypefn
+
+## Author: Lukas Reichlin <luk...@gm...>
+## Created: April 2012
+## Version: 0.1
+
+function [sys, varargout] = arx (dat, varargin)
+
+ ## TODO: delays
+
+ if (nargin < 2)
+ print_usage ();
+ endif
+
+ if (! isa (dat, "iddata"))
+ error ("arx: first argument must be an iddata dataset");
+ endif
+
+% if (nargin > 2) # arx (dat, ...)
+ if (is_real_scalar (varargin{1})) # arx (dat, n, ...)
+ varargin = horzcat (varargin(2:end), {"na"}, varargin(1), {"nb"}, varargin(1));
+ endif
+ if (isstruct (varargin{1})) # arx (dat, opt, ...), arx (dat, n, opt, ...)
+ varargin = horzcat (__opt2cell__ (varargin{1}), varargin(2:end));
+ endif
+% endif
+
+ nkv = numel (varargin); # number of keys and values
+
+ if (rem (nkv, 2))
+ error ("arx: keys and values must come in pairs");
+ endif
+
+
+ ## p: outputs, m: inputs, ex: experiments
+ [~, p, m, ex] = size (dat); # dataset dimensions
+
+ ## extract data
+ Y = dat.y;
+ U = dat.u;
+ tsam = dat.tsam;
+
+ ## multi-experiment data requires equal sampling times
+ if (ex > 1 && ! isequal (tsam{:}))
+ error ("arx: require equally sampled experiments");
+ else
+ tsam = tsam{1};
+ endif
+
+
+ ## default arguments
+ na = [];
+ nb = []; % ???
+ nk = [];
+
+ ## handle keys and values
+ for k = 1 : 2 : nkv
+ key = lower (varargin{k});
+ val = varargin{k+1};
+ switch (key)
+ ## TODO: proper argument checking
+ case "na"
+ na = val;
+ case "nb"
+ nb = val;
+ case "nk"
+ error ("nk");
+ otherwise
+ warning ("arx: invalid property name '%s' ignored", key);
+ endswitch
+ endfor
+
+
+ if (is_real_scalar (na, nb))
+ na = repmat (na, p, 1); # na(p-by-1)
+ nb = repmat (nb, p, m); # nb(p-by-m)
+ elseif (! (is_real_vector (na) && is_real_matrix (nb) \
+ && rows (na) == p && rows (nb) == p && columns (nb) == m))
+ error ("arx: require na(%dx1) instead of (%dx%d) and nb(%dx%d) instead of (%dx%d)", \
+ p, rows (na), columns (na), p, m, rows (nb), columns (nb));
+ endif
+
+ max_nb = max (nb, [], 2); # one maximum for each row/output, max_nb(p-by-1)
+ n = max (na, max_nb); # n(p-by-1)
+
+ ## create empty cells for numerator and denominator polynomials
+ % num = cell (p, m+p);
+ % den = cell (p, m+p);
+ num = cell (p, m);
+ den = cell (p, m);
+
+ ## MIMO (p-by-m) models are identified as pm SISO models
+ ## For multi-experiment data, minimize the trace of the error
+ for i = 1 : p # for every output
+ for j = 1 : m # for every input
+ Phi = cell (ex, 1); # one regression matrix per experiment
+ for e = 1 : ex # for every experiment
+ ## avoid warning: toeplitz: column wins anti-diagonal conflict
+ ## therefore set first row element equal to y(1)
+ PhiY = toeplitz (Y{e}(1:end-1, i), [Y{e}(1, i); zeros(na(i,j)-1, 1)]);
+ PhiU = toeplitz (U{e}(1:end-1, j), [U{e}(1, j); zeros(nb(i,j)-1, 1)]);
+ Phi{e} = [-PhiY, PhiU](n(i):end, :);
+ endfor
+
+ ## compute parameter vector Theta
+ Theta = __theta__ (Phi, Y, i, n);
+
+ ## extract polynomials A and B from Theta
+ A = [1; Theta(1:na(i,j))]; # a0 = 1, a1 = Theta(1), an = Theta(n)
+ B = [0; Theta(na(i,j)+1:end)]; # b0 = 0 (leading zero required by filt)
+
+ num(i,j) = A;
+ den(i,j) = B;
+
+ %{
+ ## add error inputs
+ Be = repmat ({0}, 1, p); # there are as many error inputs as system outputs (p)
+ Be(i) = 1; # inputs m+1:m+p are zero, except m+i which is one
+ num(i, :) = [B, Be]; # numerator polynomials for output i, individual for each input
+ den(i, :) = repmat ({A}, 1, m+p); # in a row (output i), all inputs have the same denominator polynomial
+ %}
+ endfor
+ endfor
+
+ %{
+ ## A(q) y(t) = B(q) u(t) + e(t)
+ ## there is only one A per row
+ ## B(z) and A(z) are a Matrix Fraction Description (MFD)
+ ## y = A^-1(q) B(q) u(t) + A^-1(q) e(t)
+ ## since A(q) is a diagonal polynomial matrix, its inverse is trivial:
+ ## the corresponding transfer function has common row denominators.
+ %}
+
+ sys = filt (num, den, tsam); # filt creates a transfer function in z^-1
+
+ ## compute initial state vector x0 if requested
+ ## this makes only sense for state-space models, therefore convert TF to SS
+ if (nargout > 1)
+ sys = prescale (ss (sys(:,1:m)));
+ x0 = slib01cd (Y, U, sys.a, sys.b, sys.c, sys.d, 0.0);
+ ## return x0 as vector for single-experiment data
+ ## instead of a cell containing one vector
+ if (numel (x0) == 1)
+ x0 = x0{1};
+ endif
+ varargout{1} = x0;
+ endif
+
+endfunction
+
+
+function theta = __theta__ (phi, y, i, n)
+
+ if (numel (phi) == 1) # single-experiment dataset
+ ## use "square-root algorithm"
+ A = horzcat (phi{1}, y{1}(n(i)+1:end, i)); # [Phi, Y]
+ R0 = triu (qr (A, 0)); # 0 for economy-size R (without zero rows)
+ R1 = R0(1:end-1, 1:end-1); # R1 is triangular - can we exploit this in R1\R2?
+ R2 = R0(1:end-1, end);
+ theta = __ls_svd__ (R1, R2); # R1 \ R2
+
+ ## Theta = Phi \ Y(n+1:end, :); # naive formula
+ ## theta = __ls_svd__ (phi{1}, y{1}(n(i)+1:end, i));
+ else # multi-experiment dataset
+ ## TODO: find more sophisticated formula than
+ ## Theta = (Phi1' Phi + Phi2' Phi2 + ...) \ (Phi1' Y1 + Phi2' Y2 + ...)
+
+ ## covariance matrix C = (Phi1' Phi + Phi2' Phi2 + ...)
+ tmp = cellfun (@(Phi) Phi.' * Phi, phi, "uniformoutput", false);
+ % rc = cellfun (@rcond, tmp); # C auch noch testen? QR oder SVD?
+ C = plus (tmp{:});
+
+ ## PhiTY = (Phi1' Y1 + Phi2' Y2 + ...)
+ tmp = cellfun (@(Phi, Y) Phi.' * Y(n(i)+1:end, i), phi, y, "uniformoutput", false);
+ PhiTY = plus (tmp{:});
+
+ ## pseudoinverse Theta = C \ Phi'Y
+ theta = __ls_svd__ (C, PhiTY);
+ endif
+
+endfunction
+
+
+function x = __ls_svd__ (A, b)
+
+ ## solve the problem Ax=b
+ ## x = A\b would also work,
+ ## but this way we have better control and warnings
+
+ ## solve linear least squares problem by pseudoinverse
+ ## the pseudoinverse is computed by singular value decomposition
+ ## M = U S V* ---> M+ = V S+ U*
+ ## Th = Ph \ Y = Ph+ Y
+ ## Th = V S+ U* Y, S+ = 1 ./ diag (S)
+
+ [U, S, V] = svd (A, 0); # 0 for "economy size" decomposition
+ S = diag (S); # extract main diagonal
+ r = sum (S > eps*S(1));
+ if (r < length (S))
+ warning ("arx: rank-deficient coefficient matrix");
+ warning ("sampling time too small");
+ warning ("persistence of excitation");
+ endif
+ V = V(:, 1:r);
+ S = S(1:r);
+ U = U(:, 1:r);
+ x = V * (S .\ (U' * b)); # U' is the conjugate transpose
+
+endfunction
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <jpi...@us...> - 2012-06-09 00:53:37
|
Revision: 10607
http://octave.svn.sourceforge.net/octave/?rev=10607&view=rev
Author: jpicarbajal
Date: 2012-06-09 00:53:32 +0000 (Sat, 09 Jun 2012)
Log Message:
-----------
signal: Adding new functions. clustersegment is old and can be improved a lot.
Modified Paths:
--------------
trunk/octave-forge/main/signal/inst/clustersegment.m
Modified: trunk/octave-forge/main/signal/inst/clustersegment.m
===================================================================
--- trunk/octave-forge/main/signal/inst/clustersegment.m 2012-06-09 00:52:28 UTC (rev 10606)
+++ trunk/octave-forge/main/signal/inst/clustersegment.m 2012-06-09 00:53:32 UTC (rev 10607)
@@ -59,11 +59,11 @@
%!demo
%! xhi = [0 0 1 1 1 0 0 1 0 0 0 1 1];
-%! ranges = intervalSegment(xhi)
+%! ranges = clustersegment(xhi)
%!
%! % The first sequence of 1's in xhi is
%! xhi(ranges{1}(1,:))
%!demo
%! xhi = rand(3,10)>0.4
-%! ranges = intervalSegment(xhi)
+%! ranges = clustersegment(xhi)
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <jpi...@us...> - 2012-06-09 00:52:34
|
Revision: 10606
http://octave.svn.sourceforge.net/octave/?rev=10606&view=rev
Author: jpicarbajal
Date: 2012-06-09 00:52:28 +0000 (Sat, 09 Jun 2012)
Log Message:
-----------
signal: Adding new functions. clustersegment is old and can be improved a lot.
Modified Paths:
--------------
trunk/octave-forge/main/signal/INDEX
trunk/octave-forge/main/signal/NEWS
Added Paths:
-----------
trunk/octave-forge/main/signal/inst/clustersegment.m
trunk/octave-forge/main/signal/inst/movingrms.m
trunk/octave-forge/main/signal/inst/schtrig.m
Modified: trunk/octave-forge/main/signal/INDEX
===================================================================
--- trunk/octave-forge/main/signal/INDEX 2012-06-08 19:20:19 UTC (rev 10605)
+++ trunk/octave-forge/main/signal/INDEX 2012-06-09 00:52:28 UTC (rev 10606)
@@ -27,6 +27,7 @@
sgolayfilt
sosfilt
medfilt1
+ movingrms
Filter analysis
freqs freqs_plot
grpdelay
@@ -119,7 +120,8 @@
kaiser
nuttallwin
triang
- gaussian gausswin
+ gaussian
+ gausswin
tukeywin
rectwin
welchwin
@@ -160,3 +162,5 @@
wrev
zerocrossing
sampled2continuous
+ schtrig
+ clustersegment
Modified: trunk/octave-forge/main/signal/NEWS
===================================================================
--- trunk/octave-forge/main/signal/NEWS 2012-06-08 19:20:19 UTC (rev 10605)
+++ trunk/octave-forge/main/signal/NEWS 2012-06-09 00:52:28 UTC (rev 10606)
@@ -4,6 +4,9 @@
signal-1.1.4 Release Date: 2012-XX-XX Release Manager:
===============================================================================
+ ** The following functions are new:
+ movingrms schtrig clustersegment
+
** signal is no longer dependent on the image package.
** signal is now dependent on the general package.
Added: trunk/octave-forge/main/signal/inst/clustersegment.m
===================================================================
--- trunk/octave-forge/main/signal/inst/clustersegment.m (rev 0)
+++ trunk/octave-forge/main/signal/inst/clustersegment.m 2012-06-09 00:52:28 UTC (rev 10606)
@@ -0,0 +1,69 @@
+## Copyright (c) 2010 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
+## (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/>.
+
+%% contRange = clustersegment(xhi)
+% The function calculates the initial and end index values of the sequences of
+% 1's in the rows of xhi. The result is returned in a cell of size 1xNp, with Np
+% the numer of rows in xhi. Each element of the cell has two rows; the first row
+% is the inital index of a sequence of 1's and the second row is the end index
+% of that sequence.
+
+function contRange = clustersegment(xhi)
+
+ % Find discontinuities
+ bool_discon = diff(xhi,1,2);
+ [Np Na] = size(xhi);
+ contRange = cell(1,Np);
+
+ for i = 1:Np
+ idxUp = find(bool_discon(i,:)>0)+1;
+ idxDwn = find(bool_discon(i,:)<0);
+ tLen = length(idxUp) + length(idxDwn);
+
+ if xhi(i,1)==1
+ % first event was down
+ contRange{i}(1) = 1;
+ contRange{i}(2:2:tLen+1) = idxDwn;
+ contRange{i}(3:2:tLen+1) = idxUp;
+ else
+ % first event was up
+ contRange{i}(1:2:tLen) = idxUp;
+ contRange{i}(2:2:tLen) = idxDwn;
+ end
+
+ if xhi(i,end)==1
+ % last event was up
+ contRange{i}(end+1) = Na;
+ end
+
+ tLen = length(contRange{i});
+ if tLen ~=0
+ contRange{i}=reshape(contRange{i},2,tLen/2);
+ end
+
+ end
+
+endfunction
+
+%!demo
+%! xhi = [0 0 1 1 1 0 0 1 0 0 0 1 1];
+%! ranges = intervalSegment(xhi)
+%!
+%! % The first sequence of 1's in xhi is
+%! xhi(ranges{1}(1,:))
+
+%!demo
+%! xhi = rand(3,10)>0.4
+%! ranges = intervalSegment(xhi)
Property changes on: trunk/octave-forge/main/signal/inst/clustersegment.m
___________________________________________________________________
Added: svn:executable
+ *
Added: trunk/octave-forge/main/signal/inst/movingrms.m
===================================================================
--- trunk/octave-forge/main/signal/inst/movingrms.m (rev 0)
+++ trunk/octave-forge/main/signal/inst/movingrms.m 2012-06-09 00:52:28 UTC (rev 10606)
@@ -0,0 +1,96 @@
+## 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
+## (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{rmsx},@var{w}] =} movingrms (@var{x},@var{w},@var{rc},@var{Fs}=1)
+## Calculates moving RMS value of the signal in @var{x}.
+##
+## The signla is convoluted against a sigmoid window of width @var{w} and risetime
+## @var{rc}. The units of these to parameters are relative ot the vlaue of the sampling
+## frequency given in @var{Fs} (Default value = 1).
+##
+## Run @code{demo movingrms} to see an example.
+##
+## @seealso{sigmoid_train}
+## @end deftypefn
+
+function [rmsx w]= movingrms (x,width, risetime, Fs=1)
+
+ [N nc] = size (x);
+ if width*Fs > N/2
+ idx = [1 N];
+ w = ones(N,1);
+ else
+ idx = round ((N + width*Fs*[-1 1])/2);
+ w = sigmoid_train ((1:N)', idx, risetime*Fs);
+ end
+ fw = fft (w.^2);
+ fx = fft (x.^2);
+ rmsx = real(ifft (fx.*fw)/(N-1));
+
+ rmsx (rmsx < eps*max(rmsx(:))) = 0;
+
+ rmsx = circshift (sqrt (rmsx), round(mean(idx)));
+% w = circshift (w, -idx(1));
+
+endfunction
+
+%!demo
+%! N = 128;
+%! t = linspace(0,1,N)';
+%! x = sigmoid_train (t,[0.4 inf],1e-2).*(2*rand(size(t))-1);
+%!
+%! Fs = 1/diff(t(1:2));
+%! width = 0.05;
+%! rc = 5e-3;
+%! [wx w] = movingrms (zscore (x),width,rc,Fs);
+%!
+%! close all
+%! figure ()
+%!
+%! area (t,wx,'facecolor',[0.85 0.85 1],'edgecolor','b','linewidth',2);
+%! hold on;
+%! h = plot (t,x,'r-;Data;',t,w,'g-;Window;');
+%! set (h, 'linewidth', 2);
+%! hold off
+%!
+%! % ---------------------------------------------------------------------------
+%! % The shaded plot shows the local RMS of the Data: white noise with onset at
+%! % aprox. t== 0.4.
+%! % The observation window is also shown.
+
+%!demo
+%! N = 128;
+%! t = linspace(0,1,N)';
+%! x = exp(-((t-0.5)/0.1).^2) + 0.1*rand(N,1);
+%!
+%! Fs = 1/diff(t(1:2));
+%! width = 0.1;
+%! rc = 2e-3;
+%! [wx w] = movingrms (zscore (x),width,rc,Fs);
+%!
+%! close all
+%! figure ()
+%!
+%! area (t,wx,'facecolor',[0.85 0.85 1],'edgecolor','b','linewidth',2);
+%! hold on;
+%! h = plot (t,x,'r-;Data;',t,w,'g-;Window;');
+%! set (h, 'linewidth', 2);
+%! hold off
+%!
+%! % ---------------------------------------------------------------------------
+%! % The shaded plot shows the local RMS of the Data: Gausian with centered at
+%! % aprox. t== 0.5.
+%! % The observation window is also shown.
Added: trunk/octave-forge/main/signal/inst/schtrig.m
===================================================================
--- trunk/octave-forge/main/signal/inst/schtrig.m (rev 0)
+++ trunk/octave-forge/main/signal/inst/schtrig.m 2012-06-09 00:52:28 UTC (rev 10606)
@@ -0,0 +1,95 @@
+## 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
+## (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{rmsx},@var{w}] =} schtrig (@var{x},@var{lvl},@var{rst}=1)
+## Implements a multisignal Schmitt trigger with levels @var{lvl}.
+##
+## The triger works along the first dimension of the array @var{x}. When @code{@var{rst}==1}
+## the state of the trigger for all signals is set to the low state (i.e. 0).
+##
+## Run @code{demo schtrig} to see an example.
+##
+## @seealso{clustersegment}
+## @end deftypefn
+
+function v = schtrig (x, lvl, rst = 1)
+
+ persistent st0;
+
+ if length(lvl) == 1
+ lvl = abs (lvl)*[1 -1];
+ else
+ lvl = sort(lvl,'descend');
+ end
+
+ [nT nc] = size(x);
+
+ v = NA (nT, nc);
+
+ if rst || isempty(st0)
+ st0 = zeros(1,nc);
+ printf ("Trigger initialized!\n");
+ flush (stdout);
+ end
+
+ v(1,:) = st0;
+
+ % Signal is above up level
+ up = x > lvl(1);
+ v(up) = 1;
+
+ % Signal is below down level
+ dw = x < lvl(2);
+ v(dw) = 0;
+
+ % Resolve intermediate states
+ % Find data between the levels
+ idx = isnan (v);
+ ranges = clustersegment (idx');
+
+ for i=1:nc
+ % Record the state at the begining of the interval between levels
+ if !isempty (ranges{i})
+ prev = ranges{i}(1,:)-1;
+ prev(prev<1) = 1;
+ st0 = v(prev,i);
+
+ % Copy the initial state to the interval
+ ini_idx = ranges{i}(1,:);
+ end_idx = ranges{i}(2,:);
+ for j =1:length(ini_idx)
+ v(ini_idx(j):end_idx(j),i) = st0(j);
+ end
+ end
+ end
+
+ st0 = v(end,:);
+
+endfunction
+
+%!demo
+%! t = linspace(0,1,100)';
+%! x = sin (2*pi*2*t) + sin (2*pi*5*t).*[0.8 0.3];
+%!
+%! lvl = [0.8 0.25]';
+%! v = schtrig (x,lvl);
+%!
+%! h = plot(t,x,t,v);
+%! set (h([1 3]),'color','b');
+%! set (h([2 4]),'color',[0 1 0.5]);
+%! set (h,'linewidth',2);
+%! line([0; 1],lvl([1; 1]),'color','r');
+%! line([0;1],lvl([2;2]),'color','b')
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <be...@us...> - 2012-06-08 19:20:25
|
Revision: 10605
http://octave.svn.sourceforge.net/octave/?rev=10605&view=rev
Author: benjf5
Date: 2012-06-08 19:20:19 +0000 (Fri, 08 Jun 2012)
Log Message:
-----------
Sanity checks on fastlscomplex added, usage string improved.
Modified Paths:
--------------
trunk/octave-forge/extra/lssa/fastlscomplex.cc
Modified: trunk/octave-forge/extra/lssa/fastlscomplex.cc
===================================================================
--- trunk/octave-forge/extra/lssa/fastlscomplex.cc 2012-06-08 18:43:09 UTC (rev 10604)
+++ trunk/octave-forge/extra/lssa/fastlscomplex.cc 2012-06-08 19:20:19 UTC (rev 10605)
@@ -12,26 +12,42 @@
#include <exception>
ComplexRowVector flscomplex( RowVector tvec , ComplexRowVector xvec ,
- int coefficients, int octaves, double maxfreq );
+ double maxfreq , int octaves , int coefficients);
-DEFUN_DLD(fastlscomplex,args,nargout, "Takes the complex fast least-squares transform of a set of data.") {
+DEFUN_DLD(fastlscomplex,args,nargout, "fastlscomplex(time,magnitude,maximum_frequency,octaves,coefficients)") {
+ if ( args.length() != 5 ) {
+ print_usage();
+ return octave_value_list ();
+ }
RowVector tvals = args(0).row_vector_value();
ComplexRowVector xvals = args(1).complex_row_vector_value();
- int ncoeff = args(2).int_value();
+ double omegamax = args(2).double_value();
int noctaves = args(3).int_value();
- double omegamax = args(4).double_value();
+ int ncoeff = args(4).int_value();
+ if ( tvals.numel() != xvals.numel() ){
+ if ( tvals.numel() > xvals.numel() ) {
+ error("More time values than magnitude values.");
+ } else {
+ error("More magnitude values than time values.");
+ }
+ }
+ if ( ncoeff == 0 ) error("No coefficients to compute.");
+ if ( noctaves == 0 ) error("No octaves to compute over.");
+ if ( omegamax == 0 ) error("No difference between minimal and maximal frequency.");
octave_value_list retval;
if ( !error_state) {
- ComplexRowVector results = flscomplex(tvals,xvals,ncoeff,noctaves,omegamax);
+ ComplexRowVector results = flscomplex(tvals,xvals,omegamax,noctaves,ncoeff);
retval(0) = octave_value(results);
+ } else {
+ return octave_value_list ();
}
return retval;
}
ComplexRowVector flscomplex( RowVector tvec , ComplexRowVector xvec ,
- int coefficients, int octaves, double maxfreq ) {
+ double maxfreq, int octaves, int coefficients ) {
struct Precomputation_Record {
Precomputation_Record *next;
std::complex<double> power_series[12]; // I'm using 12 as a matter of compatibility, only.
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <jpi...@us...> - 2012-06-08 18:43:15
|
Revision: 10604
http://octave.svn.sourceforge.net/octave/?rev=10604&view=rev
Author: jpicarbajal
Date: 2012-06-08 18:43:09 +0000 (Fri, 08 Jun 2012)
Log Message:
-----------
statistics: Adding explanation ofr the difference between the formulas in the book GPML and the one sin the function.
Modified Paths:
--------------
trunk/octave-forge/main/statistics/inst/regress_gp.m
Modified: trunk/octave-forge/main/statistics/inst/regress_gp.m
===================================================================
--- trunk/octave-forge/main/statistics/inst/regress_gp.m 2012-06-08 18:39:56 UTC (rev 10603)
+++ trunk/octave-forge/main/statistics/inst/regress_gp.m 2012-06-08 18:43:09 UTC (rev 10604)
@@ -46,6 +46,19 @@
x = [ones(1,size(x,1)); x'];
+ ## Juan Pablo Carbajal <car...@if...>
+ ## Note that in the book the equation (below 2.11) for the A reads
+ ## A = (1/sy^2)*x*x' + inv (Vp);
+ ## where sy is the scalar variance of the of the residuals (i.e y = x' * w + epsilon)
+ ## and epsilon is drawn from N(0,sy^2). Vp is the variance of the parameters w.
+ ## Note that
+ ## (sy^2 * A)^{-1} = (1/sy^2)*A^{-1} = (x*x' + sy^2 * inv(Vp))^{-1};
+ ## and that the formula for the w mean is
+ ## (1/sy^2)*A^{-1}*x*y
+ ## Then one obtains
+ ## inv(x*x' + sy^2 * inv(Vp))*x*y
+ ## Looking at the formula bloew we see that Sp = (1/sy^2)*Vp
+ ## making the regression depend on only one parameter, Sp, and not two.
A = x*x' + inv (Sp);
K = inv (A);
wm = K*x*y;
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <be...@us...> - 2012-06-08 18:40:02
|
Revision: 10603
http://octave.svn.sourceforge.net/octave/?rev=10603&view=rev
Author: benjf5
Date: 2012-06-08 18:39:56 +0000 (Fri, 08 Jun 2012)
Log Message:
-----------
A few changes made to fastlscomplex, improved speed just slightly.
Modified Paths:
--------------
trunk/octave-forge/extra/lssa/fastlscomplex.cc
Modified: trunk/octave-forge/extra/lssa/fastlscomplex.cc
===================================================================
--- trunk/octave-forge/extra/lssa/fastlscomplex.cc 2012-06-08 18:12:34 UTC (rev 10602)
+++ trunk/octave-forge/extra/lssa/fastlscomplex.cc 2012-06-08 18:39:56 UTC (rev 10603)
@@ -51,6 +51,11 @@
for ( int array_iter = 0 ; array_iter < 12 ; array_iter++ ) {
temp_array[array_iter] = std::complex<double> ( 0 , 0 );
}
+ int factorial_array[12];
+ factorial_array[0] = 1;
+ for ( int i = 1 ; i < 12 ; i++ ) {
+ factorial_array[i] = factorial_array[i-1] * i;
+ }
n_inv = 1.0 / n;
mu = (0.5 * M_PI)/length; // Per the article; this is in place to improve numerical accuracy if desired.
/* Viz. the paper, in which Dtau = c / omega_max, and c is stated as pi/2 for floating point processors,
@@ -88,33 +93,29 @@
alpha.real() = xvec(k).real();
alpha.imag() = xvec(k).imag();
// alpha = std::complex<double> ( xvec(k).real() , xvec(k).imag() );
- if ( j == 0 ) {
- p = 1;
+ // alpha *= ( -1 * im * mu * ( tvec(k) - tau_h ) ) * p;
+ if ( !( j % 2 ) ) {
+ if ( ! ( j % 4 ) ) {
+ alpha.real() = xvec(k).real() * pow(mu,j) * pow(tvec(k)-tau_h,j) / factorial_array[j];
+ alpha.imag() = xvec(k).imag() * pow(mu,j) * pow(tvec(k)-tau_h,j) / factorial_array[j];
+ } else {
+ alpha.real() = -1 * xvec(k).real() * pow(mu,j) * pow(tvec(k)-tau_h,j) / factorial_array[j];
+ alpha.imag() = -1 * xvec(k).imag() * pow(mu,j) * pow(tvec(k)-tau_h,j) / factorial_array[j];
+ }
} else {
- p *= 1/(double)j;
- // alpha *= ( -1 * im * mu * ( tvec(k) - tau_h ) ) * p;
- if ( !( j % 2 ) ) {
- if ( ! ( j % 4 ) ) {
- alpha.real() = xvec(k).real() * pow(mu,j) * pow(tvec(k)-tau_h,j);
- alpha.imag() = xvec(k).imag() * pow(mu,j) * pow(tvec(k)-tau_h,j);
- } else {
- alpha.real() = -1 * xvec(k).real() * pow(mu,j) * pow(tvec(k)-tau_h,j);
- alpha.imag() = -1 * xvec(k).imag() * pow(mu,j) * pow(tvec(k)-tau_h,j);
- } } else {
- if ( ! ( j % 3 ) ) {
- alpha.real() = -1 * xvec(k).imag() * pow(mu,j) * pow(tvec(k)-tau_h,j);
- alpha.imag() = -1 * xvec(k).real() * pow(mu,j) * pow(tvec(k)-tau_h,j);
- } else {
- alpha.real() = xvec(k).imag() * pow(mu,j) * pow(tvec(k)-tau_h,j);
- alpha.imag() = xvec(k).real() * pow(mu,j) * pow(tvec(k)-tau_h,j);
- }
+ if ( ! ( j % 3 ) ) {
+ alpha.real() = -1 * xvec(k).imag() * pow(mu,j) * pow(tvec(k)-tau_h,j) / factorial_array[j];
+ alpha.imag() = -1 * xvec(k).real() * pow(mu,j) * pow(tvec(k)-tau_h,j) / factorial_array[j];
+ } else {
+ alpha.real() = xvec(k).imag() * pow(mu,j) * pow(tvec(k)-tau_h,j) / factorial_array[j];
+ alpha.imag() = xvec(k).real() * pow(mu,j) * pow(tvec(k)-tau_h,j) / factorial_array[j];
}
}
record_current->power_series[j].real() += alpha.real();
record_current->power_series[j].imag() += alpha.imag();
- // Computes each next step of the power series for the given power series element.
- // j was reused since it was a handy inner-loop variable, even though I used it twice here.
}
+ // Computes each next step of the power series for the given power series element.
+ // j was reused since it was a handy inner-loop variable, even though I used it twice here.
record_current->stored_data = true;
k++;
}
@@ -167,8 +168,9 @@
* Otherwise, merge with the next computation range.
*/
double exp_power_series_elements[12];
- for ( int r_iter = 0 ; r_iter < 12 ; r_iter++ ) {
- exp_power_series_elements[r_iter] = ( r_iter == 0 ) ? 1 : exp_power_series_elements[r_iter-1]
+ exp_power_series_elements[0] = 1;
+ for ( int r_iter = 1 ; r_iter < 12 ; r_iter++ ) {
+ exp_power_series_elements[r_iter] = exp_power_series_elements[r_iter-1]
* ( mu * loop_delta_tau) * ( 1.0 / ( (double) r_iter ) );
}
try{
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <et...@us...> - 2012-06-08 18:12:40
|
Revision: 10602
http://octave.svn.sourceforge.net/octave/?rev=10602&view=rev
Author: etienne
Date: 2012-06-08 18:12:34 +0000 (Fri, 08 Jun 2012)
Log Message:
-----------
Fix bugs
Modified Paths:
--------------
trunk/octave-forge/main/vrml/inst/test_moving_surf.m
trunk/octave-forge/main/vrml/inst/vrml_TimeSensor.m
Modified: trunk/octave-forge/main/vrml/inst/test_moving_surf.m
===================================================================
--- trunk/octave-forge/main/vrml/inst/test_moving_surf.m 2012-06-08 18:09:31 UTC (rev 10601)
+++ trunk/octave-forge/main/vrml/inst/test_moving_surf.m 2012-06-08 18:12:34 UTC (rev 10602)
@@ -33,8 +33,7 @@
s2 = vrml_anim ("Coordinate",[a,b,a],"foo.set_point",[0 0.5 1],5);
- s3 = vrml_faces ([-1 -1 1 1;-1 1 1 -1;0.1 0.1 0.1 0.1],\
- list ([1 2 3 4]),"tran",0.4,"col",[0.3 0.9 0.4]);
+ s3 = vrml_faces ([-1 -1 1 1;-1 1 1 -1;0.1 0.1 0.1 0.1], {[1 2 3 4]}, "tran",0.4,"col",[0.3 0.9 0.4]);
vrml_browse ([s1,s2,s3])
@@ -57,7 +56,7 @@
s4 = vrml_anim ("Coordinate",[a,b,a],"foo.set_point",[0 0.5 1],tn);
s3 = vrml_faces ([-1 -1 1 1;-1 1 1 -1;0.1 0.1 0.1 0.1],\
- list ([1 2 3 4]),"tran",0.4,"col",[0.9 0.4 0.4]);
+ {[1 2 3 4]},"tran",0.4,"col",[0.9 0.4 0.4]);
vrml_browse ([s1,s2,s3,s4])
printf ("Press a key. \n"); pause ();
Modified: trunk/octave-forge/main/vrml/inst/vrml_TimeSensor.m
===================================================================
--- trunk/octave-forge/main/vrml/inst/vrml_TimeSensor.m 2012-06-08 18:09:31 UTC (rev 10601)
+++ trunk/octave-forge/main/vrml/inst/vrml_TimeSensor.m 2012-06-08 18:12:34 UTC (rev 10602)
@@ -75,11 +75,10 @@
end
end
+DEF = 0;
if rem (length (varargin), 2), error ("vrml_TimeSensor : Odd n. of arguments"); end
-DEF = 0;
-
l = {"TimeSensor {\n"};
i = 1;
while i < length (varargin)
@@ -102,8 +101,7 @@
if verbose, printf ("vrml_TimeSensor : Defining node '%s'\n",v); end
if DEF, error ("vrml_TimeSensor : Multiple DEFs found"); end
- ##l = splice (l,1,0,{"DEF ",v," "});
- varargin = {"DEF ",v," ",l};
+ l = {sprintf("DEF %s ", v), l{:}};
DEF = 1;
else # Add data field
@@ -151,7 +149,14 @@
end
end
+
l{end+1} = "}\n";
-s = feval ("strcat", l{:});
+
+s = "";
+for i=1:numel(l)
+ s = [s, sprintf(l{i})];
+endfor
+### Stupid strcat removes trailing spaces in l's elements
+### s = strcat (l{:});
endfunction
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <et...@us...> - 2012-06-08 18:09:39
|
Revision: 10601
http://octave.svn.sourceforge.net/octave/?rev=10601&view=rev
Author: etienne
Date: 2012-06-08 18:09:31 +0000 (Fri, 08 Jun 2012)
Log Message:
-----------
Put back minimal data2vrml()
Added Paths:
-----------
trunk/octave-forge/main/vrml/inst/data2vrml.m
Added: trunk/octave-forge/main/vrml/inst/data2vrml.m
===================================================================
--- trunk/octave-forge/main/vrml/inst/data2vrml.m (rev 0)
+++ trunk/octave-forge/main/vrml/inst/data2vrml.m 2012-06-08 18:09:31 UTC (rev 10601)
@@ -0,0 +1,40 @@
+## Copyright (C) 2012 Etienne Grossmann <et...@eg...>
+##
+## 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/>.
+
+## s = data2vrml (typeStr, value) - Convert 'value' to VRML code of type typeStr
+##
+## TODO: Improve this function
+##
+## If typeStr is "SFBool", then s is "TRUE" or "FALSE"
+## If typeStr is "MFString", then s is sprintf ("%s", value)
+## otherwise s is sprintf ("%f", value)
+##
+function s = data2vrml (typeStr, value)
+
+if strcmp (typeStr, "SFBool")
+
+ if value, s = "TRUE";
+ else s = "FALSE";
+ endif
+
+elseif strcmp (typeStr, "MFSTring")
+
+ s = sprintf ("%s", value);
+
+else
+
+ s = sprintf ("%f", value);
+
+end
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-06-08 17:47:50
|
Revision: 10600
http://octave.svn.sourceforge.net/octave/?rev=10600&view=rev
Author: paramaniac
Date: 2012-06-08 17:47:43 +0000 (Fri, 08 Jun 2012)
Log Message:
-----------
control: remove cruft
Removed Paths:
-------------
trunk/octave-forge/main/control/devel/makefile_helpers.m
trunk/octave-forge/main/control/devel/minreal.m
trunk/octave-forge/main/control/devel/test_ctrbf.m
trunk/octave-forge/main/control/devel/tf2ss_draft.m
Deleted: trunk/octave-forge/main/control/devel/makefile_helpers.m
===================================================================
--- trunk/octave-forge/main/control/devel/makefile_helpers.m 2012-06-08 16:45:39 UTC (rev 10599)
+++ trunk/octave-forge/main/control/devel/makefile_helpers.m 2012-06-08 17:47:43 UTC (rev 10600)
@@ -1,23 +0,0 @@
-## ==============================================================================
-## Developer Makefile for OCT-files
-## ==============================================================================
-## USAGE: * fetch control from Octave-Forge by svn
-## * add control/inst, control/src and control/devel to your Octave path
-## * run makefile_*
-## ==============================================================================
-
-homedir = pwd ();
-develdir = fileparts (which ("makefile_helpers"));
-srcdir = [develdir, "/../src"];
-cd (srcdir);
-
-mkoctfile is_real_scalar.cc
-
-mkoctfile is_real_vector.cc
-
-mkoctfile is_real_matrix.cc
-
-mkoctfile is_real_square_matrix.cc
-
-system ("rm *.o");
-cd (homedir);
\ No newline at end of file
Deleted: trunk/octave-forge/main/control/devel/minreal.m
===================================================================
--- trunk/octave-forge/main/control/devel/minreal.m 2012-06-08 16:45:39 UTC (rev 10599)
+++ trunk/octave-forge/main/control/devel/minreal.m 2012-06-08 17:47:43 UTC (rev 10600)
@@ -1,45 +0,0 @@
-## Copyright (C) 2010 Benjamin Fernandez <ma...@be...>
-##
-## 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 Octave; see the file COPYING. If not, see
-## <http://www.gnu.org/licenses/>.
-
-## -*- texinfo -*-
-## @deftypefn{Function File} {[@var{Amin}, @var{Bmin}, @var{Cmin}] =} ctrbf (@var{A}, @var{B}, @var{C})
-## @deftypefnx{Function File} {[@var{Amin}, @var{Bmin}, @var{Cmin}] =} ctrbf (@var{A}, @var{B}, @var{C}, @var{TOL})
-## Minimal realization of the system (A,B,C).
-## If the system is controlable and observable, the syste es minimal.
-## If the system is not controlable or/and observable, a new system created
-## with the controlabe and observable subspace is equivalent
-## which means that Cco(sI-Aco)^(-1)Bco = C(sI-A)^(-1)B.
-## @end deftypefn
-
-## Author: Benjamin Fernandez <ma...@be...>
-## Created: 2010-07-10
-
-function [Amin,Bmin,Cmin] = minreal(A,B,C,TOL)
- if(nargin<3 || nargin>4)
- print_usage();
- endif
- n = length(A);
- if(nargin == 3)
- TOL = n*norm(A,1)*eps;
- endif
- [Abar,Bbar,Cbar,T,K,Ac,Bc,Cc,Cnc,lc] = ctrbf(A,B,C);
- [Abar,Bbar,Cbar,T,K,Amin,Bmin,Cmin,Cno,lo] = obsvf(Ac,Bc,Cc);
- l = lc+lo;
- disp('States reduced:');
- disp(l);
-
-endfunction
-
Deleted: trunk/octave-forge/main/control/devel/test_ctrbf.m
===================================================================
--- trunk/octave-forge/main/control/devel/test_ctrbf.m 2012-06-08 16:45:39 UTC (rev 10599)
+++ trunk/octave-forge/main/control/devel/test_ctrbf.m 2012-06-08 17:47:43 UTC (rev 10600)
@@ -1,84 +0,0 @@
-a = [ -1.0 0.0 0.0
- -2.0 -2.0 -2.0
- -1.0 0.0 -3.0 ];
-
-b = [ 1.0 0.0 0.0
- 0.0 2.0 1.0 ].';
-
-
-c = [ 0.0 2.0 1.0
- 1.0 0.0 0.0 ];
-
-
-[ac, bc, cc, z, ncont] = sltb01ud (a, b, c, 0.0)
-
-z\a*z
-
-a = [ 1 1
- 4 -2 ];
-
-b = [ 1 -1
- 1 -1 ];
-
-c = [ 1 0
- 0 1 ];
-
-[ac, bc, cc, z, ncont] = sltb01ud (a, b, c, 0.0)
-
-%{
- The transformed state dynamics matrix of a controllable realization is
- -3.0000 2.2361
- 0.0000 -1.0000
-
- and the dimensions of its diagonal blocks are
- 2
-
- The transformed input/state matrix B of a controllable realization is
- 0.0000 -2.2361
- 1.0000 0.0000
-
- The transformed output/state matrix C of a controllable realization is
- -2.2361 0.0000
- 0.0000 1.0000
-
- The controllability index of the transformed system representation = 1
-
- The similarity transformation matrix Z is
- 0.0000 1.0000 0.0000
- -0.8944 0.0000 -0.4472
- -0.4472 0.0000 0.8944
-%}
-%{
-A =
- 1 1
- 4 -2
-
-B =
- 1 -1
- 1 -1
-
-C =
- 1 0
- 0 1
-and locate the uncontrollable mode.
-
-[Abar,Bbar,Cbar,T,k]=ctrbf(A,B,C)
-
-Abar =
- -3.0000 0
- -3.0000 2.0000
-
-Bbar =
- 0.0000 0.0000
- 1.4142 -1.4142
-
-Cbar =
- -0.7071 0.7071
- 0.7071 0.7071
-
-T =
- -0.7071 0.7071
- 0.7071 0.7071
-k =
- 1 0
-%}
\ No newline at end of file
Deleted: trunk/octave-forge/main/control/devel/tf2ss_draft.m
===================================================================
--- trunk/octave-forge/main/control/devel/tf2ss_draft.m 2012-06-08 16:45:39 UTC (rev 10599)
+++ trunk/octave-forge/main/control/devel/tf2ss_draft.m 2012-06-08 17:47:43 UTC (rev 10600)
@@ -1,76 +0,0 @@
-function retsys = tf2ss_draft ()
-
-num = {[1, 5, 7], [1]; [1, 7], [1, 5, 5]};
-den = {[1, 5, 6], [1, 2]; [1, 8, 6], [1, 3, 2]};
-sys = tf (num, den);
-
-%{
-sys = tf (1, [1, 0])
-sys = tf (1, [1, 1])
-
-sys = tf (1, conv ([1, 1, 1], [1, 4, 6, 4, 1]))
-
-sys = tf ()
-sys = tf ("s")
-%}
-
-sys = tf (WestlandLynx);
-tol = sqrt (eps)
-
- [p, m] = size (sys);
- [num, den] = tfdata (sys);
-
- numc = cell (p, m);
- denc = cell (p, 1);
-
- ## multiply all denominators in a row and
- ## update each numerator accordingly
- for i = 1 : p
- denc(i) = __conv__ (den{i,:});
- for j = 1 : m
- idx = setdiff (1:m, j);
- numc(i,j) = __conv__ (num{i,j}, den{i,idx});
- endfor
- endfor
-
- len_numc = cellfun (@length, numc);
- len_denc = cellfun (@length, denc);
-
- ## tfpoly ensures that there are no leading zeros
- tmp = len_numc > repmat (len_denc, 1, m);
- if (any (tmp(:)))
- error ("tf: tf2ss: system must be proper");
- endif
-
- max_len_denc = max (len_denc(:));
- ucoeff = zeros (p, m, max_len_denc);
- dcoeff = zeros (p, max_len_denc);
- index = len_denc-1;
-
- for i = 1 : p
- len = len_denc(i);
- dcoeff(i, 1:len) = denc{i};
- for j = 1 : m
- ucoeff(i, j, len-len_numc(i,j)+1 : len) = numc{i,j};
- endfor
- endfor
-index, prod (index), eps*prod (index)
-min (sqrt (eps), eps*prod (index))
- [a, b, c, d] = sltd04ad (ucoeff, dcoeff, index, tol);
-
- retsys = ss (a, b, c, d);
-
-endfunction
-
-
-function vec = __conv__ (vec, varargin)
-
- if (nargin == 1)
- return;
- else
- for k = 1 : nargin-1
- vec = conv (vec, varargin{k});
- endfor
- endif
-
-endfunction
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <prn...@us...> - 2012-06-08 16:45:45
|
Revision: 10599
http://octave.svn.sourceforge.net/octave/?rev=10599&view=rev
Author: prnienhuis
Date: 2012-06-08 16:45:39 +0000 (Fri, 08 Jun 2012)
Log Message:
-----------
Updated to release 1.0.19
Modified Paths:
--------------
trunk/octave-forge/main/io/DESCRIPTION
trunk/octave-forge/main/io/NEWS
Modified: trunk/octave-forge/main/io/DESCRIPTION
===================================================================
--- trunk/octave-forge/main/io/DESCRIPTION 2012-06-08 16:39:09 UTC (rev 10598)
+++ trunk/octave-forge/main/io/DESCRIPTION 2012-06-08 16:45:39 UTC (rev 10599)
@@ -1,6 +1,6 @@
Name: io
-Version: 1.0.18
-Date: 2012-03-22
+Version: 1.0.19
+Date: 2012-06-08
Author: various authors
Maintainer: Philip Nienhuis <prn...@us...>
Title: Input/Output
Modified: trunk/octave-forge/main/io/NEWS
===================================================================
--- trunk/octave-forge/main/io/NEWS 2012-06-08 16:39:09 UTC (rev 10598)
+++ trunk/octave-forge/main/io/NEWS 2012-06-08 16:45:39 UTC (rev 10599)
@@ -1,6 +1,8 @@
Summary of important user-visible changes for releases of the io package
-post-1.0.18 (in SVN only)
+===============================================================================
+io-1.0.19 Release Date: 2012-06-08 Release Manager: Philip Nienhuis
+===============================================================================
** Bug fixes:
--- getusedrange subfunc getusedrange_jod: str2num applied to indices rather
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.
