changeset 18:b1399f66b364

Moved some functions to general/numerical/scalar, some to stats

--- a/general/numerical/frint.m	Thu Jan 17 13:21:46 2013 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,22 +0,0 @@
-function y=frint(r,x,dc)
-% frint - Fractional integration
-%
-% y=frint(r,x[,dc])
-% Returns the rth order integral of x, where r can be fractional
-% if dc=0 (default is 1) dc term is not added back in
-
-x=x(:);
-n=length(x);
-z=fft(x); z0=z(1); z(1)=0;
-K=(0:n-1)';
-D=exp(2*pi*i*K/n)-1;
-D(1)=1; % to avoid divide by zero
-z2=z./D.^r;
-y=real(ifft(z2)); % integral minus DC term
-
-% put DC term  back as a ramp
-if nargin<3 || dc==1,
-	y=y+sign(z0)*(K*abs(z0)/n).^r;
-end
-
-

--- a/general/numerical/hzeta.c	Thu Jan 17 13:21:46 2013 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,96 +0,0 @@
-/*=================================================================
- * hzeta.c 
- *
- * Wrapper function for gsl_sf_hzeta_e
- *
- * This MEX function computes the Hurwitz-zeta-function by calling 
- * the function gsl_sf_hzeta_e from the GNU Scientific Library. For 
- * x>1 and a>=0 the Hurwitz-zeta-function is defined as
- * 
- *                        oo 
- *                       -----
- *                        \         1
- *       hzeta(x, a) =     )     -------- .
- *                        /             x
- *                       -----   (i + a)
- *                       i = 0
- *
- * The MEX function takes two arguments, x and a. If both arguments 
- * are matrices of the same size, then the MEX function will compute
- * the Hurwitz-zeta-function elementwise and the result of the MEX 
- * function will be a matrix of the same size as the inputs. If one 
- * input argument is a scalar s, then the function behaves as if this
- * argument would be a matrix filled with the value s of the same 
- * size as the other argument.
- *
- * Compile this MEX function with "mex hzeta.c -lgsl". See Matlab
- * documentation for additional information on compiling MEX functions.
- *
- * (C) by Heiko Bauke, 2007, [e-mail]
- *
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU General Public License in
- * version 2 as published by the Free Software Foundation.
- *
- * 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, write to the Free Software
- * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
- * 02111-1307, USA.
- *
- *=================================================================*/
-
-#include "mex.h"
-#include <gsl/gsl_errno.h>
-#include <gsl/gsl_math.h>
-#include <gsl/gsl_sf_zeta.h>
-
-void mexFunction(int nlhs, mxArray *plhs[], 
-                 int nrhs, const mxArray *prhs[]) {
-  double *x, *y, *a;
-  int nrows_x, ncols_x, i_x, 
-      nrows_a, ncols_a, i_a, 
-      nrows_y, ncols_y, i_y, 
-      status;
-  gsl_sf_result result;
-
-  /* check input arguments */
-  if (nrhs!=2)
-    mexErrMsgTxt("Two input arguments required.\n");
-  if (!mxIsDouble(prhs[0]))
-    mexErrMsgTxt("1st argument not a matrix of double."); 
-  nrows_x=mxGetM(prhs[0]);
-  ncols_x=mxGetN(prhs[0]);
-  if (!mxIsDouble(prhs[1]))
-    mexErrMsgTxt("2nd argument not a matrix of double."); 
-  nrows_a=mxGetM(prhs[1]);
-  ncols_a=mxGetN(prhs[1]);
-  if (!( (nrows_a==nrows_x && ncols_a==ncols_x) ||
-	 (nrows_x==1 && ncols_x==1) ||
-	 (nrows_a==1 && ncols_a==1) ) ) 
-    mexErrMsgTxt("arguments must be matrices of the same size or scalars."); 
-  /* create matrix for the return argument */ 
-  nrows_y=nrows_x>nrows_a ? nrows_x : nrows_a;
-  ncols_y=ncols_x>ncols_a ? ncols_x : ncols_a;
-  plhs[0]=mxCreateDoubleMatrix(nrows_y, ncols_y, mxREAL); 
-  /* assign pointers to each input and output */ 
-  x=mxGetPr(prhs[0]); 
-  a=mxGetPr(prhs[1]);
-  y=mxGetPr(plhs[0]); 
-  /* compute zeta function */
-  gsl_set_error_handler_off();
-  for (i_x=i_a=i_y=0; i_y<nrows_y*ncols_y; ++i_y) {
-    if (gsl_sf_hzeta_e(x[i_x], a[i_a], &result)==GSL_SUCCESS)
-      y[i_y]=result.val;
-    else
-      y[i_y]=GSL_NAN;
-    if (!(nrows_x==1 && ncols_x==1))
-      ++i_x;
-    if (!(nrows_a==1 && ncols_a==1))
-      ++i_a;
-  }
-}

Binary file general/numerical/hzeta.mexmac has changed

--- a/general/numerical/measure.m	Thu Jan 17 13:21:46 2013 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,8 +0,0 @@
-function b=measure(x,l)
-% measure - returns vector in given direction but with a different length
-%
-% measure :: [[N]] -> [[N]] ~'unit vector of length 1'.
-% measure :: [[N]], L:nonneg -> [[N]] ~'vector of length L'.
-
-if nargin<2, l=1; end
-b=(l/norm(x))*x;

--- a/general/numerical/safe_exp.m	Thu Jan 17 13:21:46 2013 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,11 +0,0 @@
-function [y,k]=safe_exp(x,l)
-% safe_exp - high dynamic range exponential
-%
-% safe_exp :: X:[[N,M]] -> [[N,M]], [[1,M]].
-%
-% returns y and k such that exp(X) = y * exp(k) and
-% maximum value in y is 1.
-
-k=max(x,[],1);
-y=exp(x-repmat(k,size(x,1),1));
-	

--- a/general/numerical/safe_exp2.m	Thu Jan 17 13:21:46 2013 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,11 +0,0 @@
-function [y,k]=safe_exp2(x,l)
-% safe_exp2 - high dynamic range exponential
-%
-% safe_exp2 :: X:[[N,M]], L:real -> [[N,M]], [[1,M]].
-%
-% returns y and k such that exp(X) = y * exp(k) and
-% maximum value in y is exp(L).
-% Maximum sensible value of L is log(realmax)=709.78.
-
-k=max(x,[],1)-l; 
-y=exp(x-repmat(k,size(x,1),1));