Re: dynamically linked functions

[Phillip Smith]

>>>>> "Kaj" == Kaj Wiik <address@hidden> writes:

Kaj> On Thu, 30 Mar 1995, Phillip Smith wrote:
>> I just wanted to send kudos for the dynamically linked function
>> support in Octave 1.1.1 (I'm using a Sparc 10 running SunOS 4.1.3).
>> I recentedly finished writing several functions to support i/o to
>> an in-house datafile format.  Everything went smoothly, even though
>> this

Kaj> Very interesting! Could you please send function(/s) as an
Kaj> example for the list? It would be most helpful for ones like me
Kaj> with very little programming experience.

Kaj> Kaj Wiik HUT/Metsahovi

It's not particularly pretty, but here is an example routine to
perform linear interpolation.  I'm using gcc 2.6.3 and libg++ 2.6.2.
The following lines compile the source on my machine and copy it to
one of the default directories for oct files:

c++ -DHAVE_CONFIG_H -fno-implicit-templates \
 -I/home/windchime/u/spock/packages/math/octave-1.1.1 \
 -I/home/windchime/u/spock/packages/math/octave-1.1.1/src \
 -I/usr/local/include/octave -g -target sun4 -c xinterp.cc
cp xinterp2.o /usr/local/lib/octave/site/oct/sparc-sun-sunos4.1.3/xinterp2.oct

I hope this helps ... 

Regards,
Phillip
address@hidden

---------------------------- begin enclosure --------------------------
#include <iostream.h>

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include "defun-dld.h"
#include "tree-const.h"
#include "help.h"
#include "utils.h"
#include "error.h"
#include "gripes.h"

DEFUN_DLD ("xinterp", Fxinterp, Sxinterp, FSxinterp, -1, -1,
"[z1,z2,...,zn] = xinterp (x,xi,y1,y2,...,yn)\n\
\n\
given y=f(x), linearly interpolate to find z=f(xi).\n\
Both x and xi must be monotonic increasing.  If some \n\
values of xi are outside the bounds of x, a warning message \n\
is printed and the corresponding z values are computed by \n\
linear extrapolation.")
{
  Octave_object retval;
  int i,j,k,l;
  int nr, nc;
  int lx, lxi;
  double a;
  double dx;

  // check for correct number of arguments
  
  int nargin = args.length ();
  
  if ( nargin < 3 )
    {
      print_usage ("xinterp");
      return retval;
    }
  
  if ( nargin==3 && nargout == 0 )
    nargout = 1;
  
  if ( nargout != nargin-2 )
    {
      print_usage ("xinterp");
      return retval;
    }

  // check that x is a monotonic increasing vector
  
  int nr_x = args(0).rows ();
  int nc_x = args(0).columns ();

  if ( ! args(0).is_real_type () )
    {
      error("xinterp: x isn't real",i);
      return retval;
    }
  
  if ( nr_x == 0 || nc_x == 0 )
    {
      error("xinterp: x has a zero dimension");
      return retval;
    }
  
  if ( nr_x == 1 && nc_x == 1 )
    {
      error("xinterp: x cannot be a scalar");
      return retval;
    }

  if ( nr_x != 1 && nc_x != 1 )
    {
      error("xinterp: x must be a vector");
      return retval;
    }

  ColumnVector x = args(0).vector_value ();
  lx = x.length ();

  double min_x = x (0);
  double max_x = x (0);

  for ( i=1; i<lx; i++ ) 
    {
      a = x (i);
      
      if ( a < min_x )
        min_x = a;
      else if ( a > max_x )
        max_x = a;
      
      if ( a <= x (i-1) ) 
        {
          error("xinterp: x is not monotonic increasing");
          return retval;
        }
    }

  // check that xi is a monotonic increasing vector
  
  int nr_xi = args(1).rows ();
  int nc_xi = args(1).columns ();

  if ( ! args(1).is_real_type () )
    {
      error("xinterp: xi isn't real",i);
      return retval;
    }
  
  if ( nr_xi == 0 || nc_xi == 0 )
    {
      error("xinterp: xi has a zero dimension");
      return retval;
    }
  
  if ( nr_xi != 1 && nc_xi != 1 )
    {
      error("xinterp: xi must be a vector");
      return retval;
    }

  ColumnVector xi = args(1).vector_value ();
  lxi = xi.length ();
  
  double min_xi = xi (0);
  double max_xi = xi (0);

  for ( i=1; i<lxi; i++ ) 
    {
      a = xi (i);
      
      if ( a <= xi (i-1) ) 
        {
          error("xinterp: xi is not monotonic increasing");
          return retval;
        }

      if ( a < min_xi )
        min_xi = a;
      else if ( a > max_xi )
        max_xi = a;
    }

  if ( min_xi < min_x || max_xi > max_x )
    warning("xinterp: some values of xi are outside the bounds of x");

  // check that all y are the same length as x

  for ( i=2; i<nargin; i++ ) 
    {

      if ( ! args(i).is_real_type () )
        {
          error("xinterp: y%d isn't real",i-1);
          return retval;
        }

      int nr = args(i).rows ();
      int nc = args(i).columns ();
      
      if ( nr == 1 || nc == 1 ) 
        {
          if ( nr*nc != lx )
            {
              error("xinterp: y%d(%d) isn't the same length as 
x(%d)",i-1,nr*nc,lx);
              return retval;
            }
        } 
      else
        { 
          if ( nr != lx )
            {
              error("xinterp: y%d(%d,%d) doesn't conform to 
x(%d)",i-1,nr,nc,lx);
              return retval;
            }
        }
          
    }

  // compute cross-reference vector which maps from xi into x, i.e., 
  // find k=xref(j) such that x(k) <= xi(j) < x(k+1). If xi(j)>max(x)
  // then xref(j)=length(x)-2, i.e., extrapolate based on final two 
  // samples. If xi(j)<min(x) then xref(j)=0, i.e., extrapolate based 
  // on initial two samples.

  ColumnVector xref (lxi);

  int ndx = lx-2;   // ensure that ndx+1 is always a valid index value for x

  for ( i=lxi-1; i>=0; i-- )
    {
      while ( x (ndx) > xi (i)  && ndx>0 )
        ndx--;
      xref (i) = ndx;
    }

  // perform the interpolation

  for ( i=2; i<nargin; i++ ) 
    {
      
      int nr = args(i).rows ();
      int nc = args(i).columns ();
      if ( nr == 1 )
        nc = lxi;
      else if ( nc == 1 )
        nr = lxi;
      else 
        nr = lxi;
      
      Matrix tmp (nr,nc,0.0);
      Matrix y = args(i).matrix_value();

      if ( nr == 1 )
        {
          for ( j=0; j<lxi; j++ )
            {
              k = NINT( xref (j) );
              tmp.elem (0,j) = y.elem(0,k) + ( y.elem(0,k+1)-y.elem(0,k) ) /
                ( x(k+1) - x(k) ) * ( xi(j) - x(k) );
            }
        }
      else if ( nc == 1 )
        {
          for ( j=0; j<lxi; j++ )
            {
              k = NINT( xref (j) );
              tmp.elem (j,0) = y.elem(k,0) + ( y.elem(k+1,0)-y.elem(k,0) ) /
                ( x(k+1) - x(k) ) * ( xi(j) - x(k) );
            }
        }
      else
        {
          for ( j=0; j<lxi; j++ )
            {
              k = NINT( xref(j) );
              dx = ( xi(j) - x(k) ) / ( x(k+1) - x(k) );
              for ( l=0; l<nc; l++ ) 
                tmp.elem(j,l) = y.elem(k,l) + (y.elem(k+1,l)-y.elem(k,l))*dx;
            }
        }
      retval (i-2) = tmp;
    }

  return retval;
}
----------------------------- end enclosure ---------------------------