octave-maintainers
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Moving code from octave-forge to octave [Was: polyderiv problem?]


From: David Bateman
Subject: Moving code from octave-forge to octave [Was: polyderiv problem?]
Date: Wed, 09 Feb 2005 15:52:38 +0100
User-agent: Mozilla Thunderbird 0.8 (X11/20040923)

I'm move this discussion to maintainers, as that seems appropriate at this point.

Paul Kienzle wrote:


On Feb 9, 2005, at 6:10 AM, David Bateman wrote:

John W. Eaton wrote:

It looks like you are using the polyderiv.m from octave-forge, but in
any case, Octave's version of this function has the same bug.  Please
try the following patch.  I haven't looked at the version from
octave-forge to see if it is significantly different from the one in
Octave, but if it is, it would be nice if someone would reconcile the
differences and submit a patch.  If they are the same, then it might
be good to remove the version from octave-forge.

Thanks,


If you would accept some of the other stuff in FIXES and patches, then I'd be willing to feed them to you in small digestable pieces...


The README says why everything is there.

We should probably drop hankel, toeplitz, triu, tril, but make sure the test cases get into octave first.

The need for these function points to a problem. I suspect the proper fix is to code some of them as oct-files. In fact I include here an implementation of triu/tril as an oct-file, that is against the sparse-merge branch. This version is implemented in a manner similar to the octave version but is twice as fast as the octave-forge version and uses 5 times less memory. Perhaps this should replace the octave version and we should dump the existing octave-forge and octave versions. toeplitz and hankel should have the same treatment.

As for the test cases, triu/tril don't have any. But toeplitz and hankel do and should be ported.

tf2zp and zp2tf use a very simple algorithm. Confirm with control systems toolbox folks that it is sufficient for their needs. As far as I can tell it reduces to the same thing, but there is a lot to trace through in the control systems toolbox to be sure of this.

Ok, so someone needs to answer that one.

grid is fine.  I see no reason not to include it.

Rand will require some work if John wants to preserve the existing sequences for the existing seeds. Using the 'seed'/'state' distinction will allow that.

I'd vote to get rid of randpak completely and just dump the existing sequences with existing seeds. If we were to go this path, it would be better to try and get the same behavior with the same keys as matlab. However that means we have to figure out what the uniform generator that is used as a base to matlab is... I see this as hard to impossible, so again I'd prefer to not generate the same sequences as matlab.

lin2mu and mu2lin change the interface.

They don't seem to be basically incompatiable though. So again perhaps they should be included without change.

poly* you've already done.

Some shadow functions exist in other parts of octave-forge, usually because Octave forge hacked something together which worked before John had a chance to add it properly to Octave.
admin/make_index should tell you which.

These should also be treated in a general de-crufting of octave-forge...

D.

--
David Bateman                                address@hidden
Motorola Labs - Paris +33 1 69 35 48 04 (Ph) Parc Les Algorithmes, Commune de St Aubin +33 1 69 35 77 01 (Fax) 91193 Gif-Sur-Yvette FRANCE

The information contained in this communication has been classified as: [x] General Business Information [ ] Motorola Internal Use Only [ ] Motorola Confidential Proprietary

/*

Copyright (C) 2004 David Bateman

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, 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, write to the Free
Software Foundation, 59 Temple Place - Suite 330, Boston, MA
02111-1307, USA.

*/

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include "dNDArray.h"
#include "CNDArray.h"
#include "lo-mappers.h"

#include "defun-dld.h"
#include "error.h"
#include "oct-obj.h"

DEFUN_DLD (tril, args, nargout,
  "-*- texinfo -*-\n\
@deftypefn {Function File} {} tril (@var{a}, @var{k})\n\
@deftypefnx {Function File} {} triu (@var{a}, @var{k})\n\
Return a new matrix formed by extracting extract the lower (@code{tril})\n\
or upper (@code{triu}) triangular part of the matrix @var{a}, and\n\
setting all other elements to zero.  The second argument is optional,\n\
and specifies how many diagonals above or below the main diagonal should\n\
also be set to zero.\n\
\n\
The default value of @var{k} is zero, so that @code{triu} and\n\
@code{tril} normally include the main diagonal as part of the result\n\
matrix.\n\
\n\
If the value of @var{k} is negative, additional elements above (for\n\
@code{tril}) or below (for @code{triu}) the main diagonal are also\n\
selected.\n\
\n\
The absolute value of @var{k} must not be greater than the number of\n\
sub- or super-diagonals.\n\
\n\
For example,\n\
\n\
@example\n\
@group\n\
tril (ones (3), -1)\n\
     @result{}  0  0  0\n\
         1  0  0\n\
         1  1  0\n\
@end group\n\
@end example\n\
\n\
@noindent\n\
and\n\
\n\
@example\n\
@group\n\
tril (ones (3), 1)\n\
     @result{}  1  1  0\n\
         1  1  1\n\
         1  1  1\n\
@end group\n\
@end example\n\
@end deftypefn\n\
\n\
@seealso{triu, diag}")
{

  octave_value retval;
  int nargin = args.length ();
  int k = 0;

  if (nargin == 2)
    {
      k = args(1).int_value();
      
      if (error_state)
        return retval;
    }

  if (nargin < 1 || nargin > 2)
    usage ("tril");
  else if (args(0).class_name () == "sparse")
    {
      if (args(0).is_complex_type())
        {
          SparseComplexMatrix m = args(0).sparse_complex_matrix_value();

          if (!error_state)
            {
              int nr = m.rows();
              int nc = m.cols();
              if ((k > 0 && k > nc) || (k < 0 && k < -nr))
                {
                  error ("tril: requested diagonal out of range");
                  return retval;
                }

              for (int j = 0; j < nc; j++)
                for (int i = m.cidx(j); i < m.cidx(j+1); i++)
                  if (m.ridx(i) < j-k)
                    m.data(i) = 0.;

              m.maybe_compress (true);
              retval = m;
            }
        }
      else
        {
          SparseMatrix m = args(0).sparse_matrix_value();

          if (!error_state)
            {
              int nr = m.rows();
              int nc = m.cols();
              if ((k > 0 && k > nc) || (k < 0 && k < -nr))
                {
                  error ("tril: requested diagonal out of range");
                  return retval;
                }

              for (int j = 0; j < nc; j++)
                for (int i = m.cidx(j); i < m.cidx(j+1); i++)
                  if (m.ridx(i) < j-k)
                    m.data(i) = 0.;

              m.maybe_compress (true);
              retval = m;
            }
        }
    }
  else
    {
      if (args(0).is_complex_type())
        {
          ComplexMatrix m = args(0).complex_matrix_value();

          if (!error_state)
            {
              Complex *m_vec = m.fortran_vec();
              int nr = m.rows();
              int nc = m.cols();
              if ((k > 0 && k > nc) || (k < 0 && k < -nr))
                {
                  error ("tril: requested diagonal out of range");
                  return retval;
                }

              for (int j = 0; j < nc; j++)
                for (int i = 0; i < (j-k < nr ? j-k : nr); i++)
                  m_vec[i+j*nr] = 0.;

              retval = m;
            }
        }
      else
        {
          Matrix m = args(0).matrix_value();

          if (!error_state)
            {
              double *m_vec = m.fortran_vec();
              int nr = m.rows();
              int nc = m.cols();
              if ((k > 0 && k > nc) || (k < 0 && k < -nr))
                {
                  error ("tril: requested diagonal out of range");
                  return retval;
                }

              for (int j = 0; j < nc; j++)
                for (int i = 0; i < (j-k < nr ? j-k : nr); i++)
                  m_vec[i+j*nr] = 0.;

              retval = m;
            }
        }
    }

  return retval;
}

DEFUN_DLD (triu, args, nargout,
  "-*- texinfo -*-\n\
@deftypefn {Function File} {} triu (@var{a}, @var{k})\n\
See tril.\n\
@deftypefn")
{
  octave_value retval;
  int nargin = args.length ();
  int k = 0;

  if (nargin == 2)
    {
      k = args(1).int_value();
      
      if (error_state)
        return retval;
    }

  if (nargin < 1 || nargin > 2)
    usage ("tril");
  else if (args(0).class_name () == "sparse")
    {
      if (args(0).is_complex_type())
        {
          SparseComplexMatrix m = args(0).sparse_complex_matrix_value();

          if (!error_state)
            {
              int nr = m.rows();
              int nc = m.cols();
              if ((k > 0 && k > nc) || (k < 0 && k < -nr))
                {
                  error ("triu: requested diagonal out of range");
                  return retval;
                }

              for (int j = 0; j < nc; j++)
                for (int i = m.cidx(j); i < m.cidx(j+1); i++)
                  if (m.ridx(i) > j-k)
                    m.data(i) = 0.;

              m.maybe_compress (true);
              retval = m;
            }
        }
      else
        {
          SparseMatrix m = args(0).sparse_matrix_value();

          if (!error_state)
            {
              int nr = m.rows();
              int nc = m.cols();
              if ((k > 0 && k > nc) || (k < 0 && k < -nr))
                {
                  error ("triu: requested diagonal out of range");
                  return retval;
                }

              for (int j = 0; j < nc; j++)
                for (int i = m.cidx(j); i < m.cidx(j+1); i++)
                  if (m.ridx(i) > j-k)
                    m.data(i) = 0.;

              m.maybe_compress (true);
              retval = m;
            }
        }
    }
  else
    {
      if (args(0).is_complex_type())
        {
          ComplexMatrix m = args(0).complex_matrix_value();

          if (!error_state)
            {
              Complex *m_vec = m.fortran_vec();
              int nr = m.rows();
              int nc = m.cols();
              if ((k > 0 && k > nc) || (k < 0 && k < -nr))
                {
                  error ("triu: requested diagonal out of range");
                  return retval;
                }

              for (int j = 0; j < nc; j++)
                for (int i = (j-k+1 > 0 ? j-k+1 : 0); i < nr; i++)
                  m_vec[i+j*nr] = 0.;
              
              retval = m;
            }
        }
      else
        {
          Matrix m = args(0).matrix_value();

          if (!error_state)
            {
              double *m_vec = m.fortran_vec();
              int nr = m.rows();
              int nc = m.cols();
              if ((k > 0 && k > nc) || (k < 0 && k < -nr))
                {
                  error ("triu: requested diagonal out of range");
                  return retval;
                }

              for (int j = 0; j < nc; j++)
                for (int i = (j-k+1 > 0 ? j-k+1 : 0); i < nr; i++)
                  m_vec[i+j*nr] = 0.;
              
              retval = m;
            }
        }
    }

  return retval;
}

/*
;;; Local Variables: ***
;;; mode: C++ ***
;;; End: ***
*/

reply via email to

[Prev in Thread] Current Thread [Next in Thread]