From: <hi...@us...> - 2008-02-11 07:32:35
|
Revision: 4617 http://octave.svn.sourceforge.net/octave/?rev=4617&view=rev Author: highegg Date: 2008-02-10 23:32:41 -0800 (Sun, 10 Feb 2008) Log Message: ----------- OctGPR: a dozen of initial corrections Modified Paths: -------------- trunk/octave-forge/main/octgpr/TODO trunk/octave-forge/main/octgpr/inst/demo_octgpr.m trunk/octave-forge/main/octgpr/src/gpr_train.cc trunk/octave-forge/main/octgpr/src/optdrv.f trunk/octave-forge/main/octgpr/src/train.c Added Paths: ----------- trunk/octave-forge/main/octgpr/src/ChangeLog Removed Paths: ------------- trunk/octave-forge/main/octgpr/src/trtest.f trunk/octave-forge/main/octgpr/src/vmcmp.f trunk/octave-forge/main/octgpr/src/vmfac.f Modified: trunk/octave-forge/main/octgpr/TODO =================================================================== --- trunk/octave-forge/main/octgpr/TODO 2008-02-11 05:24:53 UTC (rev 4616) +++ trunk/octave-forge/main/octgpr/TODO 2008-02-11 07:32:41 UTC (rev 4617) @@ -3,5 +3,6 @@ o add documentation where missing o add more demos o add RBF models with GCV training (in progress) -o improve trust-region training optimizer +o improve trust-region training optimizer - stopping criterion + basen on objective decrease Modified: trunk/octave-forge/main/octgpr/inst/demo_octgpr.m =================================================================== --- trunk/octave-forge/main/octgpr/inst/demo_octgpr.m 2008-02-11 05:24:53 UTC (rev 4616) +++ trunk/octave-forge/main/octgpr/inst/demo_octgpr.m 2008-02-11 07:32:41 UTC (rev 4617) @@ -1,3 +1,23 @@ +% Copyright (C) 2008 VZLU Prague, a.s., Czech Republic +% +% Author: Jaroslav Hajek <hi...@gm...> +% +% This file is part of OctGPR. +% +% OctGPR is free software; you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation; either version 2 of the License, or +% (at your option) any later version. +% +% This program is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this software; see the file COPYING. If not, see +% <http://www.gnu.org/licenses/>. +% % demo of Gaussian Process Regression package 1; Added: trunk/octave-forge/main/octgpr/src/ChangeLog =================================================================== --- trunk/octave-forge/main/octgpr/src/ChangeLog (rev 0) +++ trunk/octave-forge/main/octgpr/src/ChangeLog 2008-02-11 07:32:41 UTC (rev 4617) @@ -0,0 +1,13 @@ +2008-02-07 Jaroslav Hajek <hi...@gm...> + + * gpr_train.cc: corrected inline Texinfo documentation + * ChangeLog: added this file + +2008-02-08 Jaroslav Hajek <hi...@gm...> + + * train.c: corrected C->Fortran calls in to use proper #defines from forsubs.h + * vmfac.f vmcmp.f trstep.f: deleted (no use) + * optdrv.f: cosmetic changes + * ../inst/demo_octgpr.m: added copyright + * ../TODO: modified + Modified: trunk/octave-forge/main/octgpr/src/gpr_train.cc =================================================================== --- trunk/octave-forge/main/octgpr/src/gpr_train.cc 2008-02-11 05:24:53 UTC (rev 4616) +++ trunk/octave-forge/main/octgpr/src/gpr_train.cc 2008-02-11 07:32:41 UTC (rev 4617) @@ -66,7 +66,7 @@ @item imq: @code{f(t) = 1/sqrt(1+t^2)} \n\ @end enumerate \n\ \n\ -@var{opts} is a cell array in the form {\"option name\",option value,...}. \n\ +@var{opts} is a cell array in the form @{\"option name\",option value,...@}. \n\ Possible options: \n\ maxev: maximum number of factorizations to be used during training. default 500.\n\ tol: stopping tolerance (minimum trust-region radius). default 1e-6. \n\ Modified: trunk/octave-forge/main/octgpr/src/optdrv.f =================================================================== --- trunk/octave-forge/main/octgpr/src/optdrv.f 2008-02-11 05:24:53 UTC (rev 4616) +++ trunk/octave-forge/main/octgpr/src/optdrv.f 2008-02-11 07:32:41 UTC (rev 4617) @@ -196,8 +196,8 @@ do i = 1,ndim theta(i) = theta0(i) + stp(i) end do + nu = nu0 + stp(ndim+1) c normal return - ready for next evaluation - nu = nu0 + stp(ndim+1) info = 0 end subroutine Modified: trunk/octave-forge/main/octgpr/src/train.c =================================================================== --- trunk/octave-forge/main/octgpr/src/train.c 2008-02-11 05:24:53 UTC (rev 4616) +++ trunk/octave-forge/main/octgpr/src/train.c 2008-02-11 07:32:41 UTC (rev 4617) @@ -50,7 +50,7 @@ int i,code,info,l2nu = 1; /* setup scale factors */ - stheta_(&ndim,&nx,X,scal); + F77_stheta(&ndim,&nx,X,scal); info = 0; @@ -63,7 +63,7 @@ while(1) { /* call the optimization driver */ - optdrv_(&ndim,theta,nu,nll,dtheta,dnu, + F77_optdrv(&ndim,theta,nu,nll,dtheta,dnu, theta0,nu0,nll0,dtheta0,dnu0, &info,scal,&l2nu,VM,CP,IC); if (info == 1) { @@ -75,7 +75,7 @@ } /* evaluate objective */ - nllgpr_(&ndim,&nx,X,y,theta,nu,var,&nlin, + F77_nllgpr(&ndim,&nx,X,y,theta,nu,var,&nlin, mu,R,nll,corf,&info); if (IC[0] > opts->maxev || (opts->monitor && @@ -86,7 +86,7 @@ if (!info) /* evaluate gradient */ - nldgpr_(&ndim,&nx,X,theta,nu,var,R,dtheta,dnu,&info); + F77_nldgpr(&ndim,&nx,X,theta,nu,var,R,dtheta,dnu,&info); /* mark whether the step was successful */ if (info) info = 2; Deleted: trunk/octave-forge/main/octgpr/src/trtest.f =================================================================== --- trunk/octave-forge/main/octgpr/src/trtest.f 2008-02-11 05:24:53 UTC (rev 4616) +++ trunk/octave-forge/main/octgpr/src/trtest.f 2008-02-11 07:32:41 UTC (rev 4617) @@ -1,44 +0,0 @@ -c Copyright (C) 2008 VZLU Prague, a.s., Czech Republic -c -c Author: Jaroslav Hajek <hi...@gm...> -c -c This file is part of OctGPR. -c -c OctGPR is free software; you can redistribute it and/or modify -c it under the terms of the GNU General Public License as published by -c the Free Software Foundation; either version 2 of the License, or -c (at your option) any later version. -c -c This program is distributed in the hope that it will be useful, -c but WITHOUT ANY WARRANTY; without even the implied warranty of -c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -c GNU General Public License for more details. -c -c You should have received a copy of the GNU General Public License -c along with this software; see the file COPYING. If not, see -c <http://www.gnu.org/licenses/>. -c - program test - integer n - parameter(n = 2) - double precision w(n),ba(n),a,g(n),ga,s(n),sa,del,sal,sau - integer i,j,info - open(10,file='testcm.dat') - read(10,*) i - if (i /= n) stop 'dim mismatch' - read(10,*) w - read(10,*) ba,a - read(10,*) g,ga - read(10,*) del,sal,sau - read(10,*) s,sa - write(*,*) w - write(*,*) ba,a - write(*,*) g,ga - write(*,*) del,sal,sau - write(*,*) s,sa - print *,'mval = ',sum(w*s**2)/2+sum((g+sa*ba)*s)+ga*sa+a*sa**2/2 - close(10) - call trstp(n,w,ba,a,g,ga,s,sa,del,sal,sau,info) - print *,'s = ',s,' sa = ',sa,' info = ',info - print *,'mval = ',sum(w*s**2)/2+sum((g+sa*ba)*s)+ga*sa+a*sa**2/2 - end program Deleted: trunk/octave-forge/main/octgpr/src/vmcmp.f =================================================================== --- trunk/octave-forge/main/octgpr/src/vmcmp.f 2008-02-11 05:24:53 UTC (rev 4616) +++ trunk/octave-forge/main/octgpr/src/vmcmp.f 2008-02-11 07:32:41 UTC (rev 4617) @@ -1,54 +0,0 @@ -c Copyright (C) 2008 VZLU Prague, a.s., Czech Republic -c -c Author: Jaroslav Hajek <hi...@gm...> -c -c This file is part of OctGPR. -c -c OctGPR is free software; you can redistribute it and/or modify -c it under the terms of the GNU General Public License as published by -c the Free Software Foundation; either version 2 of the License, or -c (at your option) any later version. -c -c This program is distributed in the hope that it will be useful, -c but WITHOUT ANY WARRANTY; without even the implied warranty of -c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -c GNU General Public License for more details. -c -c You should have received a copy of the GNU General Public License -c along with this software; see the file COPYING. If not, see -c <http://www.gnu.org/licenses/>. -c - subroutine vmcmp(n,B,D,Z,ba,a,g,ga,s,sa,ssiz,redc) -c purpose: completes the solution of the mixed-norm -c trust-region subproblem initiated by vmfac and -c trstp, computes the actual step size and -c expected reduction according to the model. -c arguments: -c n (in) problem dimension. -c B (in) variable metric matrix in packed storage -c D (in) diagonal scaling matrix -c Z (in) variable metric polar base -c ba,a (in) quadratic coefficients -c g,ga (in) linear coefficients -c s,sa (io) solution (sa is not modified) -c ssiz (out) norm(s) -c redc (out) model reduction -c - integer n - double precision B(*),D(n),Z(n,n),ba(n),a,g(n),ga, - +s(n),sa,ssiz,redc - double precision tmp(n) - integer i,j -c transform from polar basis - call dgemv('N',n,n,1.d0,Z,n,s,1,0.d0,tmp,1) - do i = 1,n - s(i) = tmp(i)/D(i) - end do -c compute the actual step size - ssiz = dnrm2(n,tmp,1) -c compute the expected reduction - call dspmv('U',n,1.d0,B,tmp,1,0.d0,s,1) - redc = 0.5d0*ddot(n,tmp,1,s,1) + sa*ddot(n,ba,1,s,1) + - +0.5d0*a*sa**2 + ddot(n,g,1,s,1) + ga*sa - end subroutine - Deleted: trunk/octave-forge/main/octgpr/src/vmfac.f =================================================================== --- trunk/octave-forge/main/octgpr/src/vmfac.f 2008-02-11 05:24:53 UTC (rev 4616) +++ trunk/octave-forge/main/octgpr/src/vmfac.f 2008-02-11 07:32:41 UTC (rev 4617) @@ -1,58 +0,0 @@ -c Copyright (C) 2008 VZLU Prague, a.s., Czech Republic -c -c Author: Jaroslav Hajek <hi...@gm...> -c -c This file is part of OctGPR. -c -c OctGPR is free software; you can redistribute it and/or modify -c it under the terms of the GNU General Public License as published by -c the Free Software Foundation; either version 2 of the License, or -c (at your option) any later version. -c -c This program is distributed in the hope that it will be useful, -c but WITHOUT ANY WARRANTY; without even the implied warranty of -c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -c GNU General Public License for more details. -c -c You should have received a copy of the GNU General Public License -c along with this software; see the file COPYING. If not, see -c <http://www.gnu.org/licenses/>. -c - subroutine vmfac(n,B,D,WB,ZB,g,g1,ba,ba1) -c purpose: store in W and Z the eigendecomposition of -c B -c arguments: -c n (in) problem dimension -c B (in) variable metric matrix in packed storage -c D (in) diagonal scaling matrix -c WB (out) variable metric eigenvalues -c ZB (out) variable metric polar base -c g (in) the gradient -c g1 (out) transformed gradient -c ba (in) the theta-nu part of hessian model -c ba1 (out) transformed ba -c - integer n - double precision B(*),D(n),WB(n),ZB(n,n),g(n),g1(n),ba(n),ba1(n) - double precision B2(n*(n+1)/2),work(n,3) - integer info,i,j,k -c form the scaled matrix diag(D)*B*diag(D) - k = 0 - do j = 1,n - do i = 1,j - k = k + 1 - B2(k) = B(k) * D(i)*D(j) - end do - end do -c get the polar decomposition - call dspev('V','U',n,B2,WB,ZB,n,work,info) -c scale vectors - do i = 1,n - work(i,1) = g(i)/D(i) - work(i,2) = ba(i)/D(i) - end do -c transform to polar basis - call dgemv('T',n,n,1.d0,ZB,n,work(1,1),1,0.d0,g1,1) - call dgemv('T',n,n,1.d0,ZB,n,work(2,1),1,0.d0,ba1,1) - end subroutine - This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <hi...@us...> - 2008-02-14 09:19:33
|
Revision: 4626 http://octave.svn.sourceforge.net/octave/?rev=4626&view=rev Author: highegg Date: 2008-02-14 01:19:32 -0800 (Thu, 14 Feb 2008) Log Message: ----------- more doc fixes Modified Paths: -------------- trunk/octave-forge/main/octgpr/DESCRIPTION trunk/octave-forge/main/octgpr/src/gpr_predict.cc trunk/octave-forge/main/octgpr/src/gpr_train.cc Modified: trunk/octave-forge/main/octgpr/DESCRIPTION =================================================================== --- trunk/octave-forge/main/octgpr/DESCRIPTION 2008-02-14 04:09:10 UTC (rev 4625) +++ trunk/octave-forge/main/octgpr/DESCRIPTION 2008-02-14 09:19:32 UTC (rev 4626) @@ -1,5 +1,5 @@ Name: OctGPR -Version: 1.1.0 +Version: 1.1.1 Date: 2008-08-02 Author: Jaroslav Hajek (hi...@gm...) Title: Package for full dense Gaussian Process Regression Modified: trunk/octave-forge/main/octgpr/src/gpr_predict.cc =================================================================== --- trunk/octave-forge/main/octgpr/src/gpr_predict.cc 2008-02-14 04:09:10 UTC (rev 4625) +++ trunk/octave-forge/main/octgpr/src/gpr_predict.cc 2008-02-14 09:19:32 UTC (rev 4626) @@ -44,9 +44,7 @@ (The organization is determined by GPM.theta, as in @code{gpr_train}). \n\ \n\ @var{y} is set to the predicted dependent variable values. \n\ -\n\ If @var{sig} is requested, it is set to the estimated prediction deviations. \n\ -\n\ If @var{dy} is requested, it is populated with the prediction gradients. \n\ \n\ @seealso{gpr_train, gpr_setup}\n\ Modified: trunk/octave-forge/main/octgpr/src/gpr_train.cc =================================================================== --- trunk/octave-forge/main/octgpr/src/gpr_train.cc 2008-02-14 04:09:10 UTC (rev 4625) +++ trunk/octave-forge/main/octgpr/src/gpr_train.cc 2008-02-14 09:19:32 UTC (rev 4626) @@ -45,7 +45,7 @@ If requested, estimates the hyperparameters for Gaussian Process Regression (inverse\n\ length scales and relative noise) via reduced maximum likelihood, and then \n\ sets up the model for inference (prediction), storing necessary information in \n\ -the structure @var{GPM}, intended for use with gpr_train. \n\ +the structure @var{GPM}, intended for use with @code{gpr_predict}. \n\ \n\ @var{X} is the matrix of independent variables of the observations, \n\ @var{y} is a vector containing the dependent variables, \n\ @@ -60,21 +60,28 @@ \n\ @var{corf} specifies the decreasing function type for correlation function:\n\ @code{corr(x,y) = f((x-y).^2 * theta\')} (x,y,theta row vectors). Possible values:\n\ -@enumerate \n\ -@item gau: @code{f(t) = exp(-t)} \n\ -@item exp: @code{f(t) = exp(-sqrt(t))} \n\ -@item imq: @code{f(t) = 1/sqrt(1+t^2)} \n\ -@end enumerate \n\ \n\ -@var{opts} is a cell array in the form @{\"option name\",option value,...@}. \n\ -Possible options: \n\ +@table @option\n\ +@item gau\n\ +@code{f(t) = exp(-t)} \n\ +@item exp\n\ +@code{f(t) = exp(-sqrt(t))} \n\ +@item imq\n\ +@code{f(t) = 1/sqrt(1+t^2)} \n\ +@end table \n\ \n\ - maxev: maximum number of factorizations to be used during training. default 500.\n\ +@var{opts} is a cell array in the form @{\"option name\",option value,...@}.\n\ +Possible options:\n\ \n\ - tol: stopping tolerance (minimum trust-region radius). default 1e-6. \n\ +@table @option\n\ +@item maxev\n\ +maximum number of factorizations to be used during training. default 500.\n\ +@item tol\n\ +stopping tolerance (minimum trust-region radius). default 1e-6.\n\ +@item numin\n\ +minimum allowable noise. Default is @code{sqrt(1e1*eps)}.\n\ +@end table\n\ \n\ - numin: minimum allowable noise. Default is @code{sqrt(1e1*eps)}. \n\ -\n\ Training cell array @var{opts} is recognized even if other arguments are omitted. \n\ If it is not supplied (the last argument is not a cell array), training is skipped. \n\ \n\ This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <hi...@us...> - 2008-02-15 07:02:02
|
Revision: 4628 http://octave.svn.sourceforge.net/octave/?rev=4628&view=rev Author: highegg Date: 2008-02-14 23:02:05 -0800 (Thu, 14 Feb 2008) Log Message: ----------- improved demo Modified Paths: -------------- trunk/octave-forge/main/octgpr/inst/demo_octgpr.m trunk/octave-forge/main/octgpr/src/ChangeLog Modified: trunk/octave-forge/main/octgpr/inst/demo_octgpr.m =================================================================== --- trunk/octave-forge/main/octgpr/inst/demo_octgpr.m 2008-02-14 11:54:47 UTC (rev 4627) +++ trunk/octave-forge/main/octgpr/inst/demo_octgpr.m 2008-02-15 07:02:05 UTC (rev 4628) @@ -19,9 +19,9 @@ % <http://www.gnu.org/licenses/>. % % demo of Gaussian Process Regression package -1; +disp("2-dimensional GPR demo"); +disp("(set global variable nsamp to the number of random samples)"); -disp("2-dimensional GPR demo"); % define the test function (the well-known matlab "peaks") function z = testfun(x,y) z = 4 + 3 * (1-x).^2 .* exp(-(x.^2) - (y+1).^2) - ... @@ -30,7 +30,8 @@ end -disp("matlab peaks surface..."); +tit = "matlab ""peaks"" surface"; +disp(tit); % create the mesh onto which to interpolate t = linspace(-3,3,50); [xi,yi] = meshgrid(t,t); @@ -40,11 +41,16 @@ zimax = max(vec(zi)); zimin = min(vec(zi)); subplot(1,2,1); mesh(xi,yi,zi); +title(tit); pause; -disp("sampled at 400 random points"); -% create 400 random samples -xs = rand(400,1); ys = rand(400,1); +if (!exist("nsamp","var") || !isnumeric(nsamp)) + nsamp = 150 +end +tit = sprintf("sampled at %d random points",nsamp); +disp(tit); +% create random samples +xs = rand(nsamp,1); ys = rand(nsamp,1); xs = 6*xs-3; ys = 6*ys - 3; % evaluate at random samples zs = testfun(xs,ys); @@ -52,19 +58,23 @@ subplot(1,2,2); plot3(xs,ys,zs,".+"); +title(tit); pause; -disp("GPR model with heuristic hypers"); -ths = 0.3 ./ std(xys); +tit = "GPR model with heuristic hypers"; +disp(tit); +ths = 1 ./ std(xys); GPM = gpr_train(xys,zs,ths,1e-5); zm = gpr_predict(GPM,[vec(xi) vec(yi)]); zm = reshape(zm,size(zi)); zm = min(zm,zimax); zm = max(zm,zimin); subplot(1,2,2); mesh(xi,yi,zm); +title(tit); pause; -disp("GPR model with MLE training"); +tit = "GPR model with MLE training"; +disp(tit); fflush(stdout); GPM = gpr_train(xys,zs,ths,1e-5,{"tol",1e-5,"maxev",400,"numin",1e-8}); zm = gpr_predict(GPM,[vec(xi) vec(yi)]); @@ -72,6 +82,7 @@ zm = min(zm,zimax); zm = max(zm,zimin); subplot(1,2,2); mesh(xi,yi,zm); +title(tit); pause; close Modified: trunk/octave-forge/main/octgpr/src/ChangeLog =================================================================== --- trunk/octave-forge/main/octgpr/src/ChangeLog 2008-02-14 11:54:47 UTC (rev 4627) +++ trunk/octave-forge/main/octgpr/src/ChangeLog 2008-02-15 07:02:05 UTC (rev 4628) @@ -1,3 +1,7 @@ +2008-02-15 Jaroslav Hajek <hi...@gm...> + + * octgpr.m: improved demo. + 2008-02-13 Jaroslav Hajek <hi...@gm...> * train.c: correct allocation and copying result values. This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <hi...@us...> - 2008-02-24 17:58:17
|
Revision: 4673 http://octave.svn.sourceforge.net/octave/?rev=4673&view=rev Author: highegg Date: 2008-02-24 09:58:06 -0800 (Sun, 24 Feb 2008) Log Message: ----------- new stopping crit for training, ver 1.1.2 Modified Paths: -------------- trunk/octave-forge/main/octgpr/DESCRIPTION trunk/octave-forge/main/octgpr/src/ChangeLog trunk/octave-forge/main/octgpr/src/forsubs.h trunk/octave-forge/main/octgpr/src/gpr_train.cc trunk/octave-forge/main/octgpr/src/gprmod.h trunk/octave-forge/main/octgpr/src/nl0gpr.f trunk/octave-forge/main/octgpr/src/optdrv.f trunk/octave-forge/main/octgpr/src/train.c Modified: trunk/octave-forge/main/octgpr/DESCRIPTION =================================================================== --- trunk/octave-forge/main/octgpr/DESCRIPTION 2008-02-24 16:24:19 UTC (rev 4672) +++ trunk/octave-forge/main/octgpr/DESCRIPTION 2008-02-24 17:58:06 UTC (rev 4673) @@ -1,6 +1,6 @@ Name: OctGPR -Version: 1.1.1 -Date: 2008-08-02 +Version: 1.1.2 +Date: 2008-24-02 Author: Jaroslav Hajek (hi...@gm...) Title: Package for full dense Gaussian Process Regression Maintainer: Jaroslav Hajek (hi...@gm...) Modified: trunk/octave-forge/main/octgpr/src/ChangeLog =================================================================== --- trunk/octave-forge/main/octgpr/src/ChangeLog 2008-02-24 16:24:19 UTC (rev 4672) +++ trunk/octave-forge/main/octgpr/src/ChangeLog 2008-02-24 17:58:06 UTC (rev 4673) @@ -1,3 +1,9 @@ +2008-02-24 Jaroslav Hajek <hi...@gm...> + + * optdrv.f: implement objective reduction stopping criterion + * gprmod.h, train.c, gpr_train.cc: add support for the new feature + * forsubs.h nl0gpr.f: correct subroutine name + 2008-02-19 Jaroslav Hajek <hi...@gm...> * *.h *c *.cc: adjusted C/C++ sources to meet GNU coding standards @@ -6,7 +12,7 @@ * forsubs.h gprmod.h predict.c train.c setup.c: added const modifiers where appropriate. - * grp_train,gpr_predict: replaced Array::fortran_vec with Array::data + * grp_train, gpr_predict: replaced Array::fortran_vec with Array::data where appropriate, to prevent unnecessary copying. 2008-02-15 Jaroslav Hajek <hi...@gm...> Modified: trunk/octave-forge/main/octgpr/src/forsubs.h =================================================================== --- trunk/octave-forge/main/octgpr/src/forsubs.h 2008-02-24 16:24:19 UTC (rev 4672) +++ trunk/octave-forge/main/octgpr/src/forsubs.h 2008-02-24 17:58:06 UTC (rev 4673) @@ -33,6 +33,7 @@ #define F77_corimq F77_FUNC(corimq,CORIMQ) #define F77_nllgpr F77_FUNC(nllgpr,NLLGPR) #define F77_nldgpr F77_FUNC(nldgpr,NLDGPR) +#define F77_nl0gpr F77_FUNC(nl0gpr,NL0GPR) #define F77_nllbnd F77_FUNC(nllbnd,NLLBND) #define F77_infgpr F77_FUNC(infgpr,INFGPR) #define F77_stheta F77_FUNC(stheta,STHETA) Modified: trunk/octave-forge/main/octgpr/src/gpr_train.cc =================================================================== --- trunk/octave-forge/main/octgpr/src/gpr_train.cc 2008-02-24 16:24:19 UTC (rev 4672) +++ trunk/octave-forge/main/octgpr/src/gpr_train.cc 2008-02-24 17:58:06 UTC (rev 4673) @@ -43,32 +43,32 @@ @var{nu}, @var{nlin}, @var{corf}, @var{opts})\n\ @cindex Gaussian Process Regression model training.\n\ If requested, estimates the hyperparameters for Gaussian Process Regression (inverse\n\ -length scales and relative noise) via reduced maximum likelihood, and then \n\ -sets up the model for inference (prediction), storing necessary information in \n\ -the structure @var{GPM}, intended for use with @code{gpr_predict}. \n\ +length scales and relative noise) via reduced maximum likelihood, and then\n\ +sets up the model for inference (prediction), storing necessary information in\n\ +the structure @var{GPM}, intended for use with @code{gpr_predict}.\n\ \n\ -@var{X} is the matrix of independent variables of the observations, \n\ -@var{y} is a vector containing the dependent variables, \n\ -@var{theta} contains the (initial) inverse length scales for the regression model. \n\ -If @var{theta} is a row vector, rows of @var{X} correspond to observations, columns to \n\ -variables. Otherwise, it is the other way around. \n\ +@var{X} is the matrix of independent variables of the observations,\n\ +@var{y} is a vector containing the dependent variables,\n\ +@var{theta} contains the (initial) inverse length scales for the regression model.\n\ +If @var{theta} is a row vector, rows of @var{X} correspond to observations, columns to\n\ +variables. Otherwise, it is the other way around.\n\ \n\ -@var{nu} specifies the (initial) relative noise level. If not supplied, it defaults \n\ -to 1e-5. \n\ -@var{nlin} specifies the number of leading variables to include in linear \n\ -underlying trend. If not supplied, it defaults to 0 (constant trend). \n\ +@var{nu} specifies the (initial) relative noise level. If not supplied, it defaults\n\ +to 1e-5.\n\ +@var{nlin} specifies the number of leading variables to include in linear\n\ +underlying trend. If not supplied, it defaults to 0 (constant trend).\n\ \n\ @var{corf} specifies the decreasing function type for correlation function:\n\ @code{corr(x,y) = f((x-y).^2 * theta\')} (x,y,theta row vectors). Possible values:\n\ \n\ @table @option\n\ @item gau\n\ -@code{f(t) = exp(-t)} \n\ +@code{f(t) = exp(-t)}\n\ @item exp\n\ -@code{f(t) = exp(-sqrt(t))} \n\ +@code{f(t) = exp(-sqrt(t))}\n\ @item imq\n\ -@code{f(t) = 1/sqrt(1+t^2)} \n\ -@end table \n\ +@code{f(t) = 1/sqrt(1+t^2)}\n\ +@end table\n\ \n\ @var{opts} is a cell array in the form @{\"option name\",option value,...@}.\n\ Possible options:\n\ @@ -78,15 +78,20 @@ maximum number of factorizations to be used during training. default 500.\n\ @item tol\n\ stopping tolerance (minimum trust-region radius). default 1e-6.\n\ +the iteration terminates if the trust region gets below tol.\n\ +@item ftol\n\ +stopping tolerance (minimum objective reduction). default 1e-4.\n\ +the iteration terminates if the relative reduction of two successive\n\ +downhill steps gets below ftol and the second one is smaller.\n\ @item numin\n\ minimum allowable noise. Default is @code{sqrt(1e1*eps)}.\n\ @end table\n\ \n\ -Training cell array @var{opts} is recognized even if other arguments are omitted. \n\ -If it is not supplied (the last argument is not a cell array), training is skipped. \n\ +Training cell array @var{opts} is recognized even if other arguments are omitted.\n\ +If it is not supplied (the last argument is not a cell array), training is skipped.\n\ \n\ -On return the function creates the @var{GPM} structure, \n\ -which can subsequently be used for predictions with @code{gpr_predict}. \n\ +On return the function creates the @var{GPM} structure,\n\ +which can subsequently be used for predictions with @code{gpr_predict}.\n\ If @var{nll} is present, it is set to the resulting negative log likelihood.\n\ @seealso{gpr_predict}\n\ @end deftypefn") @@ -213,7 +218,8 @@ struct GPR_train_opts topts; // setup initial values topts.maxev = 500; - topts.tol= 1e-5; + topts.tol= 1e-6; + topts.ftol= 1e-4; topts.numin = 5e-8; topts.monitor = &progress_monitor; topts.instance = 0; @@ -239,6 +245,10 @@ { topts.tol = opts (++iopt).scalar_value (); } + else if (oname == "ftol") + { + topts.ftol = opts (++iopt).scalar_value (); + } else if (oname == "numin") { topts.numin = opts (++iopt).scalar_value (); Modified: trunk/octave-forge/main/octgpr/src/gprmod.h =================================================================== --- trunk/octave-forge/main/octgpr/src/gprmod.h 2008-02-24 16:24:19 UTC (rev 4672) +++ trunk/octave-forge/main/octgpr/src/gprmod.h 2008-02-24 17:58:06 UTC (rev 4673) @@ -33,7 +33,7 @@ /* training options for GPR */ struct GPR_train_opts { - double numin, tol; + double numin, tol, ftol; int maxev; int (*monitor)(void *instance, int num, double *nll); void *instance; Modified: trunk/octave-forge/main/octgpr/src/nl0gpr.f =================================================================== --- trunk/octave-forge/main/octgpr/src/nl0gpr.f 2008-02-24 16:24:19 UTC (rev 4672) +++ trunk/octave-forge/main/octgpr/src/nl0gpr.f 2008-02-24 17:58:06 UTC (rev 4673) @@ -18,7 +18,7 @@ c along with this software; see the file COPYING. If not, see c <http://www.gnu.org/licenses/>. c - subroutine nllbnd(nx,y,nu,nll0,nllinf) + subroutine nl0gpr(nx,y,nu,nll0,nllinf) c purpose: this function evaluates the negative log likelihood c of a GPR process regressor at boundaries: c nll0 is a common value for theta = 0 and Modified: trunk/octave-forge/main/octgpr/src/optdrv.f =================================================================== --- trunk/octave-forge/main/octgpr/src/optdrv.f 2008-02-24 16:24:19 UTC (rev 4672) +++ trunk/octave-forge/main/octgpr/src/optdrv.f 2008-02-24 17:58:06 UTC (rev 4673) @@ -56,6 +56,9 @@ c CP(8) last spatial step size c CP(9) last noise step size c CP(10) model-predicted reduction in last step +c CP(11) the objective upper limit +c CP(12) the objective reduction tolerance +c CP(13) last relative reduction of objective c c IC (in) integer counters c IC(1) number of evaluations @@ -69,7 +72,7 @@ logical l2nu double precision VM(*),CP(13) double precision stp(ndim+1),grd(ndim+1),Zg(ndim),Zba(ndim), - +snlo,snup,eps,ss + +snlo,snup,eps,ss,relr double precision dlamch,dnrm2,ddot external dlamch,dspmid,dtrstp,dnrm2,ddot integer iB,iba,iW,iZ,i @@ -81,7 +84,12 @@ if (info == 1) then c go directly to downhill if there is no last step - if (IC(1) == 0) goto 150 + if (IC(1) == 0) then +c if the first objective is better than the limit value, +c replace the limit value. + CP(11) = min(CP(11),nll) + goto 150 + end if c compute the reduction success ratio CP(10) = (nll - nll0) / CP(10) if (CP(10) < 0.1d0) then @@ -117,11 +125,13 @@ eps = sqrt(dlamch('E')) if (CP(1) < eps) CP(1) = eps if (CP(2) == 0d0) CP(2) = 1d-4 - if (CP(3) == 0d0) CP(3) = 1d-5 + if (CP(3) < 0d0) CP(3) = 1d-6 if (CP(4) == 0d0) CP(4) = 1d+2 if (CP(5) == 0d0) CP(5) = 0.6d0 if (CP(6) == 0d0) CP(6) = 1.5d0 if (CP(7) == 0d0) CP(7) = max(dxnrm2(ndim,scal,theta),nu/CP(2)) + if (CP(12) < 0d0) CP(12) = 1d-6 + CP(13) = CP(12) c set VM matrix to a multiple of identity call dspmid('U',ndim+1,1d-2,VM(iB)) nll0 = dlamch('O') @@ -149,8 +159,15 @@ 150 continue c HANDLE SUCCESSFUL STEP IC(1) = IC(1) + 1 -c if the step was downhill, replace stored values +c check downhill step if (nll < nll0) then +c compute the relative progress + if (nll0 < CP(11)) then + relr = (nll - nll0) / (nll0 - CP(11)) + else + relr = 0d0 + end if +c process downhill step IC(2) = IC(2) + 1 do i = 1,ndim theta0(i) = theta(i) @@ -160,6 +177,13 @@ dnu0(1) = dnu(1) if (l2nu) dnu0(2) = dnu(2) nll0 = nll +c check second termination criterion + if (IC(2) > 1 .and. relr < CP(13) .and. CP(13) < CP(12)) then +c indicate termination + info = 2 + return + end if + CP(13) = relr end if c overwrite the second derivative with the exact value VM(iba+ndim) = dnu0(2) Modified: trunk/octave-forge/main/octgpr/src/train.c =================================================================== --- trunk/octave-forge/main/octgpr/src/train.c 2008-02-24 16:24:19 UTC (rev 4672) +++ trunk/octave-forge/main/octgpr/src/train.c 2008-02-24 17:58:06 UTC (rev 4673) @@ -45,6 +45,7 @@ double *R = malloc (nx*(nx+2+nlin)*sizeof (double)); double *mu = malloc ((nlin+1)*(nlin+2)*sizeof (double)); double CP[10]; + double dummy; int IC[3]; int i, code, info, l2nu = 1; @@ -60,14 +61,18 @@ CP[1] = 1e-4; /* noise scale factor */ CP[2] = opts->tol; CP[3] = 1e+2; /* initial step size */ + CP[11] = opts->ftol; + /* setup the reference objective value */ + F77_nl0gpr(&nx,y,nu,&dummy,&CP[10]); + while (1) { /* call the optimization driver */ F77_optdrv (&ndim, theta, nu, nll, dtheta, dnu, theta0, nu0, nll0, dtheta0, dnu0, &info, scal, &l2nu, VM, CP, IC); - if (info == 1) + if (info == 1 || info == 2) { code = TRAIN_CONV; break; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <hi...@us...> - 2008-04-04 08:39:11
|
Revision: 4853 http://octave.svn.sourceforge.net/octave/?rev=4853&view=rev Author: highegg Date: 2008-04-04 01:39:10 -0700 (Fri, 04 Apr 2008) Log Message: ----------- replace normal equations by QR decomposition Modified Paths: -------------- trunk/octave-forge/main/octgpr/TODO trunk/octave-forge/main/octgpr/src/ChangeLog trunk/octave-forge/main/octgpr/src/configure trunk/octave-forge/main/octgpr/src/configure.in trunk/octave-forge/main/octgpr/src/nllgpr.f trunk/octave-forge/main/octgpr/src/setup.c trunk/octave-forge/main/octgpr/src/train.c Modified: trunk/octave-forge/main/octgpr/TODO =================================================================== --- trunk/octave-forge/main/octgpr/TODO 2008-04-03 13:40:48 UTC (rev 4852) +++ trunk/octave-forge/main/octgpr/TODO 2008-04-04 08:39:10 UTC (rev 4853) @@ -3,6 +3,4 @@ o add documentation where missing o add more demos o add RBF models with GCV training (in progress) -o improve trust-region training optimizer - stopping criterion - basen on objective decrease Modified: trunk/octave-forge/main/octgpr/src/ChangeLog =================================================================== --- trunk/octave-forge/main/octgpr/src/ChangeLog 2008-04-03 13:40:48 UTC (rev 4852) +++ trunk/octave-forge/main/octgpr/src/ChangeLog 2008-04-04 08:39:10 UTC (rev 4853) @@ -1,3 +1,11 @@ +2008-04-04 Jaroslav Hajek <hi...@gm...> + + * nllgpr.f: use QR decomposition instead of normal equations + for mean fitting. + * train.c, setup.c: change workspace size for nllgpr + * configure.in: update the LAPACK routines checked + * configure: regenerate + 2008-03-26 Jaroslav Hajek <hi...@gm...> * config.log, config.status, config.h, Makefile: remove. Modified: trunk/octave-forge/main/octgpr/src/configure =================================================================== --- trunk/octave-forge/main/octgpr/src/configure 2008-04-03 13:40:48 UTC (rev 4852) +++ trunk/octave-forge/main/octgpr/src/configure 2008-04-04 08:39:10 UTC (rev 4853) @@ -5229,15 +5229,15 @@ fi -{ echo "$as_me:$LINENO: checking for library containing dpotrs" >&5 -echo $ECHO_N "checking for library containing dpotrs... $ECHO_C" >&6; } -if test "${ac_cv_search_dpotrs+set}" = set; then +{ echo "$as_me:$LINENO: checking for library containing dgeqr2" >&5 +echo $ECHO_N "checking for library containing dgeqr2... $ECHO_C" >&6; } +if test "${ac_cv_search_dgeqr2+set}" = set; then echo $ECHO_N "(cached) $ECHO_C" >&6 else ac_func_search_save_LIBS=$LIBS cat >conftest.$ac_ext <<_ACEOF program main - call dpotrs + call dgeqr2 end _ACEOF for ac_lib in '' goto atlas lapack; do @@ -5265,7 +5265,7 @@ test ! -s conftest.err } && test -s conftest$ac_exeext && $as_test_x conftest$ac_exeext; then - ac_cv_search_dpotrs=$ac_res + ac_cv_search_dgeqr2=$ac_res else echo "$as_me: failed program was:" >&5 sed 's/^/| /' conftest.$ac_ext >&5 @@ -5275,26 +5275,158 @@ rm -f core conftest.err conftest.$ac_objext conftest_ipa8_conftest.oo \ conftest$ac_exeext - if test "${ac_cv_search_dpotrs+set}" = set; then + if test "${ac_cv_search_dgeqr2+set}" = set; then break fi done -if test "${ac_cv_search_dpotrs+set}" = set; then +if test "${ac_cv_search_dgeqr2+set}" = set; then : else - ac_cv_search_dpotrs=no + ac_cv_search_dgeqr2=no fi rm conftest.$ac_ext LIBS=$ac_func_search_save_LIBS fi -{ echo "$as_me:$LINENO: result: $ac_cv_search_dpotrs" >&5 -echo "${ECHO_T}$ac_cv_search_dpotrs" >&6; } -ac_res=$ac_cv_search_dpotrs +{ echo "$as_me:$LINENO: result: $ac_cv_search_dgeqr2" >&5 +echo "${ECHO_T}$ac_cv_search_dgeqr2" >&6; } +ac_res=$ac_cv_search_dgeqr2 if test "$ac_res" != no; then test "$ac_res" = "none required" || LIBS="$ac_res $LIBS" fi +{ echo "$as_me:$LINENO: checking for library containing dorm2r" >&5 +echo $ECHO_N "checking for library containing dorm2r... $ECHO_C" >&6; } +if test "${ac_cv_search_dorm2r+set}" = set; then + echo $ECHO_N "(cached) $ECHO_C" >&6 +else + ac_func_search_save_LIBS=$LIBS +cat >conftest.$ac_ext <<_ACEOF + program main + call dorm2r + end +_ACEOF +for ac_lib in '' goto atlas lapack; do + if test -z "$ac_lib"; then + ac_res="none required" + else + ac_res=-l$ac_lib + LIBS="-l$ac_lib $ac_func_search_save_LIBS" + fi + rm -f conftest.$ac_objext conftest$ac_exeext +if { (ac_try="$ac_link" +case "(($ac_try" in + *\"* | *\`* | *\\*) ac_try_echo=\$ac_try;; + *) ac_try_echo=$ac_try;; +esac +eval "echo \"\$as_me:$LINENO: $ac_try_echo\"") >&5 + (eval "$ac_link") 2>conftest.er1 + ac_status=$? + grep -v '^ *+' conftest.er1 >conftest.err + rm -f conftest.er1 + cat conftest.err >&5 + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); } && { + test -z "$ac_f77_werror_flag" || + test ! -s conftest.err + } && test -s conftest$ac_exeext && + $as_test_x conftest$ac_exeext; then + ac_cv_search_dorm2r=$ac_res +else + echo "$as_me: failed program was:" >&5 +sed 's/^/| /' conftest.$ac_ext >&5 + + +fi + +rm -f core conftest.err conftest.$ac_objext conftest_ipa8_conftest.oo \ + conftest$ac_exeext + if test "${ac_cv_search_dorm2r+set}" = set; then + break +fi +done +if test "${ac_cv_search_dorm2r+set}" = set; then + : +else + ac_cv_search_dorm2r=no +fi +rm conftest.$ac_ext +LIBS=$ac_func_search_save_LIBS +fi +{ echo "$as_me:$LINENO: result: $ac_cv_search_dorm2r" >&5 +echo "${ECHO_T}$ac_cv_search_dorm2r" >&6; } +ac_res=$ac_cv_search_dorm2r +if test "$ac_res" != no; then + test "$ac_res" = "none required" || LIBS="$ac_res $LIBS" + +fi + +{ echo "$as_me:$LINENO: checking for library containing dtrtrs" >&5 +echo $ECHO_N "checking for library containing dtrtrs... $ECHO_C" >&6; } +if test "${ac_cv_search_dtrtrs+set}" = set; then + echo $ECHO_N "(cached) $ECHO_C" >&6 +else + ac_func_search_save_LIBS=$LIBS +cat >conftest.$ac_ext <<_ACEOF + program main + call dtrtrs + end +_ACEOF +for ac_lib in '' goto atlas lapack; do + if test -z "$ac_lib"; then + ac_res="none required" + else + ac_res=-l$ac_lib + LIBS="-l$ac_lib $ac_func_search_save_LIBS" + fi + rm -f conftest.$ac_objext conftest$ac_exeext +if { (ac_try="$ac_link" +case "(($ac_try" in + *\"* | *\`* | *\\*) ac_try_echo=\$ac_try;; + *) ac_try_echo=$ac_try;; +esac +eval "echo \"\$as_me:$LINENO: $ac_try_echo\"") >&5 + (eval "$ac_link") 2>conftest.er1 + ac_status=$? + grep -v '^ *+' conftest.er1 >conftest.err + rm -f conftest.er1 + cat conftest.err >&5 + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); } && { + test -z "$ac_f77_werror_flag" || + test ! -s conftest.err + } && test -s conftest$ac_exeext && + $as_test_x conftest$ac_exeext; then + ac_cv_search_dtrtrs=$ac_res +else + echo "$as_me: failed program was:" >&5 +sed 's/^/| /' conftest.$ac_ext >&5 + + +fi + +rm -f core conftest.err conftest.$ac_objext conftest_ipa8_conftest.oo \ + conftest$ac_exeext + if test "${ac_cv_search_dtrtrs+set}" = set; then + break +fi +done +if test "${ac_cv_search_dtrtrs+set}" = set; then + : +else + ac_cv_search_dtrtrs=no +fi +rm conftest.$ac_ext +LIBS=$ac_func_search_save_LIBS +fi +{ echo "$as_me:$LINENO: result: $ac_cv_search_dtrtrs" >&5 +echo "${ECHO_T}$ac_cv_search_dtrtrs" >&6; } +ac_res=$ac_cv_search_dtrtrs +if test "$ac_res" != no; then + test "$ac_res" = "none required" || LIBS="$ac_res $LIBS" + +fi + { echo "$as_me:$LINENO: checking for library containing dpotri" >&5 echo $ECHO_N "checking for library containing dpotri... $ECHO_C" >&6; } if test "${ac_cv_search_dpotri+set}" = set; then Modified: trunk/octave-forge/main/octgpr/src/configure.in =================================================================== --- trunk/octave-forge/main/octgpr/src/configure.in 2008-04-03 13:40:48 UTC (rev 4852) +++ trunk/octave-forge/main/octgpr/src/configure.in 2008-04-04 08:39:10 UTC (rev 4853) @@ -86,7 +86,9 @@ # NOTE: BLAS should go first, so that it gets listed earlier. # many optimized blas libraries offer optimized xPOTRF AC_SEARCH_LIBS(dpotrf,goto atlas lapack) -AC_SEARCH_LIBS(dpotrs,goto atlas lapack) +AC_SEARCH_LIBS(dgeqr2,goto atlas lapack) +AC_SEARCH_LIBS(dorm2r,goto atlas lapack) +AC_SEARCH_LIBS(dtrtrs,goto atlas lapack) AC_SEARCH_LIBS(dpotri,goto atlas lapack) AC_SEARCH_LIBS(dsyev,goto atlas lapack) Modified: trunk/octave-forge/main/octgpr/src/nllgpr.f =================================================================== --- trunk/octave-forge/main/octgpr/src/nllgpr.f 2008-04-03 13:40:48 UTC (rev 4852) +++ trunk/octave-forge/main/octgpr/src/nllgpr.f 2008-04-04 08:39:10 UTC (rev 4853) @@ -36,9 +36,9 @@ c var (out) MLE estimated global variance c nlin (in) number of leading input variables to include in linear c mean trend. Set nlin = 0 to use a constant mean. -c mu (out) at least (nlin+1)*(nlin+2). The first nlin+1 components +c mu (out) at least (nlin+1)*3. The first nlin+1 components c contain the components of the linear trend (constant -c first). Rest is factorized matrix - see code. +c first). c R (out) at least (nx)*(nx+2+nlin). The first nx*(nx+1) components c contain useful factorization details and are reused in c nldgpr @@ -61,10 +61,10 @@ c increase nx. integer ndim,nx,nlin,info real*8 X(ndim,nx),y(nx),theta(ndim),nu,var - real*8 mu(0:nlin,0:nlin+1),R(nx,0:nx+1+nlin),nll + real*8 mu(0:nlin,0:2),R(nx,0:nx+1+nlin),nll external corr - external dwdis2,dpotrf,dcopy,dtrsv,dtrsm,dsyrk,dgemv,dpotrs, - +daxpy,xerbla + external dwdis2,dcopy,daxpy,dscal,dtrsv,dtrsm,xerbla, + +dpotrf,dtrtrs,dgeqr2,dorm2r integer i,j,k,nl1 real*8 sums,sum,dwdis2 @@ -117,28 +117,27 @@ if (nlin > 0) then nl1 = nlin + 1 -c fit a linear model using normal equations. Using QR decomposition can -c also be considered, but it's slower and uses more workspace. I believe -c that precision is not an issue here, as typically nx >> ndim +c fit a linear model using QR decomposition c form L \ M call dtrsm('L','L','N','N',nx,nl1, + 1.d0,R(1,1),nx,R(1,nx+1),nx) -c form MM = (M'*inv(R)*M) - call dsyrk('U','T',nl1,nx, - + 1.d0,R(1,nx+1),nx,0.d0,mu(0,1),nl1) -c factorize MM - call dpotrf('U',nl1,mu(0,1),nl1,info) +c factorize (unblocked code) Q*R = L \ M + call dgeqr2(nx,nl1,R(1,nx+1),nx,mu(0,1),mu(0,2),info) +c form yy = Q' * (L \ y) + call dorm2r('L','T',nx,1,nl1,R(1,nx+1),nx,mu(0,1), + + R(1,0),nx,mu(0,2),info) +c copy yy(1:nlin+1) + call dcopy(nl1,R(1,0),1,mu(0,0),1) +c set to zero to get residual yy(1:nlin+1) = 0 + call dscal(nl1,0d0,R(1,0),1) +c project to residual space: Q*Q' * yy + call dorm2r('L','N',nx,1,nl1,R(1,nx+1),nx,mu(0,1), + + R(1,0),nx,mu(0,2),info) +c solve R \ yy(1:nlin+1) for linear trend + call dtrtrs('U','N','N',nl1,1,R(1,nx+1),nx,mu(0,0), + + nl1,info) if (info /= 0) goto 502 -c form rhs = (L \ M)' * (L \ y) - call dgemv('T',nx,nl1, - + 1.d0,R(1,nx+1),nx,R(1,0),1,0.d0,mu(0,0),1) -c solve MM \ rhs - call dpotrs('U',nl1,1,mu(0,1),nl1,mu(0,0),nl1,info) - if (info /= 0) goto 502 -c remove the linear trend from L \ y - call dgemv('N',nx,nl1,-1.d0,R(1,nx+1),nx, - + mu(0,0),1,1.d0,R(1,0),1) else c specialized for constant fitting - for speed (and clarity, as this was c the initial version @@ -146,6 +145,7 @@ call dtrsv('L','N','N',nx,R(1,1),nx,R(1,nx+1),1) c accumulate (m'*inv(R)*m) and (y'*inv(R)*m) call dsdacc(nx,R(1,nx+1),R(1,0),sums,sum) + if (sums == 0d0) goto 502 c store the mean and matrix (1x1) mu(0,1) = sums mu(0,0) = sum / sums Modified: trunk/octave-forge/main/octgpr/src/setup.c =================================================================== --- trunk/octave-forge/main/octgpr/src/setup.c 2008-04-03 13:40:48 UTC (rev 4852) +++ trunk/octave-forge/main/octgpr/src/setup.c 2008-04-04 08:39:10 UTC (rev 4853) @@ -31,7 +31,7 @@ { /* allocate workspace */ double *R = malloc (nx*(nx+2+nlin)*DSIZE); - double *mmu = malloc ((nlin+1)*(nlin+2)*DSIZE); + double *mmu = malloc ((nlin+1)*3*DSIZE); int ierr; Modified: trunk/octave-forge/main/octgpr/src/train.c =================================================================== --- trunk/octave-forge/main/octgpr/src/train.c 2008-04-03 13:40:48 UTC (rev 4852) +++ trunk/octave-forge/main/octgpr/src/train.c 2008-04-04 08:39:10 UTC (rev 4853) @@ -24,12 +24,13 @@ #include "gprmod.h" +#define DSIZE sizeof (double) int GPR_train (int ndim, int nx, const double *X, const double *y, double *theta, double *nu, double *nll, int nlin, corfptr corf, struct GPR_train_opts *opts) { - double *wrk = malloc ((4*ndim+7)*sizeof (double)); + double *wrk = malloc ((4*ndim+7)*DSIZE); double *scal = wrk; double *theta0 = scal + ndim; double *nu0 = theta0 + ndim; @@ -39,11 +40,11 @@ double *dnu = dtheta + ndim; double *dtheta0 = dnu + 2; double *dnu0 = dtheta0 + ndim; - double *VM = malloc ((ndim+1)*(3*ndim+2)/2*sizeof (double)); + double *VM = malloc ((ndim+1)*(3*ndim+2)/2*DSIZE); /* workspace for nllgpr, nldgpr */ - double *R = malloc (nx*(nx+2+nlin)*sizeof (double)); - double *mu = malloc ((nlin+1)*(nlin+2)*sizeof (double)); + double *R = malloc (nx*(nx+2+nlin)*DSIZE); + double *mu = malloc ((nlin+1)*3*DSIZE); double CP[20]; double dummy; int IC[3]; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <hi...@us...> - 2008-04-17 13:53:48
|
Revision: 4949 http://octave.svn.sourceforge.net/octave/?rev=4949&view=rev Author: highegg Date: 2008-04-17 06:53:44 -0700 (Thu, 17 Apr 2008) Log Message: ----------- selecting RBF centers via k-means clustering Modified Paths: -------------- trunk/octave-forge/main/octgpr/DESCRIPTION trunk/octave-forge/main/octgpr/inst/demo_octgpr.m trunk/octave-forge/main/octgpr/src/Makefile.in trunk/octave-forge/main/octgpr/src/configure.in Added Paths: ----------- trunk/octave-forge/main/octgpr/ChangeLog trunk/octave-forge/main/octgpr/inst/rbf_centers.m trunk/octave-forge/main/octgpr/src/pdist2_mw.cc Removed Paths: ------------- trunk/octave-forge/main/octgpr/src/ChangeLog Added: trunk/octave-forge/main/octgpr/ChangeLog =================================================================== --- trunk/octave-forge/main/octgpr/ChangeLog (rev 0) +++ trunk/octave-forge/main/octgpr/ChangeLog 2008-04-17 13:53:44 UTC (rev 4949) @@ -0,0 +1,81 @@ +2008-04-17 Jaroslav Hajek <hi...@gm...> + + * inst/demo_octgpr.m: refactor demo + * inst/rbf_centers.m: use kmeans++ initialization due to Arthur & + Vassilvitskii + * src/Changelog -> ChangeLog: transform to package root dir + +2008-04-15 Jaroslav Hajek <hi...@gm...> + + * src/pdist2_mw.cc: new source. + * inst/rbf_centers.m: new function. + +2008-04-04 Jaroslav Hajek <hi...@gm...> + + * src/nllgpr.f: use QR decomposition instead of normal equations + for mean fitting. + * src/train.c, src/setup.c: change workspace size for nllgpr + * src/configure.in: update the LAPACK routines checked + * src/configure: regenerate + +2008-03-26 Jaroslav Hajek <hi...@gm...> + + * src/config.log, src/config.status, src/config.h, src/Makefile: remove. + * src/configure.in: add support for -fPIC flags. + * src/configure, src/config.h.in: regenerate. + * src/Makefile.in: add distclean target + +2008-02-28 Jaroslav Hajek <hi...@gm...> + + * src/corrf.f: added Matern-3/2 and Matern-5/2 correlation funcs + * src/get_corrf.c, src/forsubs.h: included the new correlation funcs + * src/gpr_train.cc: documentation for the new correlation funcs + +2008-02-24 Jaroslav Hajek <hi...@gm...> + + * src/optdrv.f: implement objective reduction stopping criterion + * src/gprmod.h, src/train.c, src/gpr_train.cc: add support for the new + feature + * src/forsubs.h src/nl0gpr.f: correct subroutine name + * src/train.c: correct CP array size + +2008-02-19 Jaroslav Hajek <hi...@gm...> + + * src/*.{h,c,cc}: adjusted C/C++ sources to meet GNU coding standards + +2008-02-18 Jaroslav Hajek <hi...@gm...> + + * src/forsubs.h src/gprmod.h src/predict.c src/train.c src/setup.c: + added const modifiers where appropriate. + * src/grp_train, src/gpr_predict: replaced Array::fortran_vec with + Array::data where appropriate, to prevent unnecessary copying. + +2008-02-15 Jaroslav Hajek <hi...@gm...> + + * inst/demo_octgpr.m: improved demo. + +2008-02-13 Jaroslav Hajek <hi...@gm...> + + * src/train.c: correct allocation and copying result values. + * src/nldgpr.f: compute 2nd derivative by a better formula. + * src/gpr_train.cc: improved documentation + * src/gpr_predict.cc: improved documentation + +2008-02-12 Jaroslav Hajek <hi...@gm...> + + * src/dscrot.f: correct the W*D*x and W'*D*x sequence. + * src/optdrv.f: correct calls to dscrot.f + +2008-02-08 Jaroslav Hajek <hi...@gm...> + + * src/train.c: corrected C->Fortran calls in to use proper #defines from forsubs.h + * src/vmfac.f vmcmp.f trstep.f: deleted (no use) + * src/optdrv.f: cosmetic changes + * inst/demo_octgpr.m: added copyright + * TODO: modified + +2008-02-07 Jaroslav Hajek <hi...@gm...> + + * src/gpr_train.cc: corrected inline Texinfo documentation + * src/ChangeLog: added this file + Modified: trunk/octave-forge/main/octgpr/DESCRIPTION =================================================================== --- trunk/octave-forge/main/octgpr/DESCRIPTION 2008-04-17 08:09:27 UTC (rev 4948) +++ trunk/octave-forge/main/octgpr/DESCRIPTION 2008-04-17 13:53:44 UTC (rev 4949) @@ -1,5 +1,5 @@ Name: OctGPR -Version: 1.1.2 +Version: 1.1.3 Date: 2008-24-02 Author: Jaroslav Hajek (hi...@gm...) Title: Package for full dense Gaussian Process Regression Modified: trunk/octave-forge/main/octgpr/inst/demo_octgpr.m =================================================================== --- trunk/octave-forge/main/octgpr/inst/demo_octgpr.m 2008-04-17 08:09:27 UTC (rev 4948) +++ trunk/octave-forge/main/octgpr/inst/demo_octgpr.m 2008-04-17 13:53:44 UTC (rev 4949) @@ -18,71 +18,140 @@ % along with this software; see the file COPYING. If not, see % <http://www.gnu.org/licenses/>. % -% demo of Gaussian Process Regression package -disp("2-dimensional GPR demo"); -disp("(set global variable nsamp to the number of random samples)"); -% define the test function (the well-known matlab "peaks") -function z = testfun(x,y) - z = 4 + 3 * (1-x).^2 .* exp(-(x.^2) - (y+1).^2) - ... - 10 * (x/5 - x.^3 - y.^5) .* exp(-x.^2 - y.^2)- ... - 1/3 * exp(-(x+1).^2 - y.^2); +% -*- texinfo -*- +% @deftypefn {Function File} demo_octgpr (1, nsamp = 150) +% @deftypefnx {Function File} demo_octgpr (2, ncnt = 20, npt = 500) +% OctGPR package demo function. +% First argument selects available demos: +% +% @itemize +% @item 1. GPR regression demo @* +% A function is sampled (with small noise), then reconstructed using GPR +% regression. @var{nsamp} specifies the number of samples. +% @seealso{gpr_train, gpr_predict} +% @item 2. RBF centers selection demo @* +% Radial basis centers are selected amongst random points. +% @var{ncnt} specifies number of centers, @var{npt} number of points. +% @seealso{rbf_centers} +% @end itemize +% @end deftypefn +function demo_octgpr (number, varargin) + switch (number) + case 1 + demo_octgpr1 (varargin{:}) + case 2 + demo_octgpr2 (varargin{:}) + otherwise + error ("demo_octgpr: invalid demo number") + endswitch +endfunction -end +% define the test function (the well-known matlab "peaks" plus some sines) +function z = testfun1 (x, y) + z = 4 + 3 * (1-x).^2 .* exp(-(x.^2) - (y+1).^2) ... + + 10 * (x/5 - x.^3 - y.^5) .* exp(-x.^2 - y.^2) ... + - 1/3 * exp(-(x+1).^2 - y.^2) ... + + 2*sin (x + y + 1e-1*x.*y); +endfunction -tit = "matlab ""peaks"" surface"; -disp(tit); -% create the mesh onto which to interpolate -t = linspace(-3,3,50); -[xi,yi] = meshgrid(t,t); +function demo_octgpr1 (nsamp = 150) + tit = "a peaked surface"; + disp (tit); -% evaluate -zi = testfun(xi,yi); -zimax = max(vec(zi)); zimin = min(vec(zi)); -subplot(1,2,1); -mesh(xi,yi,zi); -title(tit); -pause; + % create the mesh onto which to interpolate + t = linspace (-3, 3, 50); + [xi,yi] = meshgrid (t, t); -if (!exist("nsamp","var") || !isnumeric(nsamp)) - nsamp = 150 -end -tit = sprintf("sampled at %d random points",nsamp); -disp(tit); -% create random samples -xs = rand(nsamp,1); ys = rand(nsamp,1); -xs = 6*xs-3; ys = 6*ys - 3; -% evaluate at random samples -zs = testfun(xs,ys); -xys = [xs ys]; + % evaluate + zi = testfun1 (xi, yi); + zimax = max (vec (zi)); zimin = min (vec (zi)); + subplot (2, 2, 1); + mesh (xi, yi, zi); + title (tit); + subplot (2, 2, 3); + contourf (xi, yi, zi, 20); + pause; -subplot(1,2,2); -plot3(xs,ys,zs,".+"); -title(tit); -pause; + if (!exist ("nsamp", "var") || !isnumeric (nsamp)) + nsamp = 150; + endif -tit = "GPR model with heuristic hypers"; -disp(tit); -ths = 1 ./ std(xys); -GPM = gpr_train(xys,zs,ths,1e-5); -zm = gpr_predict(GPM,[vec(xi) vec(yi)]); -zm = reshape(zm,size(zi)); -zm = min(zm,zimax); zm = max(zm,zimin); -subplot(1,2,2); -mesh(xi,yi,zm); -title(tit); -pause; + tit = sprintf ("sampled at %d random points", nsamp); + disp (tit); + % create random samples + xs = rand (nsamp,1); ys = rand (nsamp,1); + xs = 6*xs-3; ys = 6*ys - 3; + % evaluate at random samples + zs = testfun1 (xs, ys); + xys = [xs ys]; -tit = "GPR model with MLE training"; -disp(tit); -fflush(stdout); -GPM = gpr_train(xys,zs,ths,1e-5,{"tol",1e-5,"maxev",400,"numin",1e-8}); -zm = gpr_predict(GPM,[vec(xi) vec(yi)]); -zm = reshape(zm,size(zi)); -zm = min(zm,zimax); zm = max(zm,zimin); -subplot(1,2,2); -mesh(xi,yi,zm); -title(tit); -pause; + subplot (2, 2, 2); + plot3 (xs, ys, zs, ".+"); + title (tit); + subplot (2, 2, 4); + plot (xs, ys, ".+"); + pause -close + tit = "GPR model with heuristic hypers"; + disp (tit); + ths = 1 ./ std (xys); + GPM = gpr_train (xys, zs, ths, 1e-5); + zm = gpr_predict (GPM, [vec(xi) vec(yi)]); + zm = reshape (zm, size(zi)); + zm = min (zm, zimax); zm = max (zm, zimin); + subplot (2, 2, 2); + mesh (xi, yi, zm); + title (tit); + subplot(2, 2, 4) + hold on + contourf (xi, yi, zm, 20); + plot (xs, ys, "+6"); + hold off + pause + + tit = "GPR model with MLE training"; + disp (tit); + fflush (stdout); + GPM = gpr_train (xys, zs, ths, 1e-5, {"tol", 1e-5, "maxev", 400, "numin", 1e-8}); + zm = gpr_predict (GPM, [vec(xi) vec(yi)]); + zm = reshape (zm, size (zi)); + zm = min (zm, zimax); zm = max (zm, zimin); + subplot (2, 2, 2); + mesh (xi, yi, zm); + title (tit); + subplot(2, 2, 4) + hold on + contourf (xi, yi, zm, 20); + plot (xs, ys, "+6"); + hold off + pause + + close +endfunction + +function demo_octgpr2 (ncnt = 50, npt = 500) + + npt = ncnt*ceil (npt/ncnt); + U = rand (ncnt, 2); + cs = min (pdist2_mw (U, 2) + diag (Inf (ncnt, 1))); + X = repmat (U, npt/ncnt, 1) + repmat (cs', npt/ncnt, 2) .* randn (npt, 2); + disp ("slightly clustered random points") + plot (X(:,1), X(:,2), "+"); + pause + + [U, ur] = rbf_centers(X, ncnt); + + fi = linspace (0, 2*pi, 20); + ncolors = rows (colormap); + hold on + for i = 1:rows (U) + xc = U(i,1) + ur(i) * cos (fi); + yc = U(i,2) + ur(i) * sin (fi); + line (xc, yc); + endfor + hold off + pause + close + +endfunction Added: trunk/octave-forge/main/octgpr/inst/rbf_centers.m =================================================================== --- trunk/octave-forge/main/octgpr/inst/rbf_centers.m (rev 0) +++ trunk/octave-forge/main/octgpr/inst/rbf_centers.m 2008-04-17 13:53:44 UTC (rev 4949) @@ -0,0 +1,97 @@ +% Copyright (C) 2008 VZLU Prague, a.s., Czech Republic +% +% Author: Jaroslav Hajek <hi...@gm...> +% +% This file is part of OctGPR. +% +% OctGPR is free software; you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation; either version 2 of the License, or +% (at your option) any later version. +% +% This program is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this software; see the file COPYING. If not, see +% <http://www.gnu.org/licenses/>. +% + +% -*- texinfo -*- +% @deftypefn {Function File} {[U, ur, iu]} = rbf_centers (@var{X}, @var{nu}, @var{theta}) +% Selects a given number of RBF centers based on Lloyd's clustering algorithm. +% +% @end deftypefn +function [U, ur, iu] = rbf_centers (X, nu, theta) + + pso_old = page_screen_output (0); + + if (nargin == 3) + X = dmult (X, theta); + elseif (nargin != 2) + print_usage (); + endif + + disp ("initializing ..."); + % the D^2 weighting initialization + + D = Inf; + kk = 1:rows (X); + cp = kk; + + for i = 1:nu + jj = sum (rand() * cp(end) < cp); + k(i) = kk(jj); + kk(jj) = []; + U = X(k(i),:); + D = min (D, pdist2_mw(X, U, 'ssq')'); + cp = cumsum (D(kk)); + endfor + + + % now perform the k-means algorithm + + U = X(k,:); + D = pdist2_mw(U, X, 'ssq'); + [xx, j] = min (D); + + it = 0; + do + for i = 1:nu + ij = find(j == i); + if (!isempty (ij)) + U(i,:) = mean (X(ij,:)); + else + U(i,:) = X(ceil (rand () * rows (X)), :); + endif + endfor + j1 = j; + D = pdist2_mw (U, X, 'ssq'); + [xx, j] = min (D); + printf ("k-means iteration %d\r", ++it); + fflush (stdout); + until (all (j == j1)) + printf ("\n"); + + if (nargout > 2) + iu = j; + endif + + if (nargout > 1) + ur = zeros (nu, 1); + for i = 1:nu + ij = (j == i); + ur(i) = sqrt (max (D(i,ij))); + endfor + endif + + if (nargin == 3) + U = dmult (U, 1./theta); + if (any(theta == 0)) + U(:,theta == 0) = 0; + endif + endif + +endfunction Deleted: trunk/octave-forge/main/octgpr/src/ChangeLog =================================================================== --- trunk/octave-forge/main/octgpr/src/ChangeLog 2008-04-17 08:09:27 UTC (rev 4948) +++ trunk/octave-forge/main/octgpr/src/ChangeLog 2008-04-17 13:53:44 UTC (rev 4949) @@ -1,68 +0,0 @@ -2008-04-04 Jaroslav Hajek <hi...@gm...> - - * nllgpr.f: use QR decomposition instead of normal equations - for mean fitting. - * train.c, setup.c: change workspace size for nllgpr - * configure.in: update the LAPACK routines checked - * configure: regenerate - -2008-03-26 Jaroslav Hajek <hi...@gm...> - - * config.log, config.status, config.h, Makefile: remove. - * configure.in: add support for -fPIC flags. - * configure, config.h.in: regenerate. - * Makefile.in: add distclean target - -2008-02-28 Jaroslav Hajek <hi...@gm...> - - * corrf.f: added Matern-3/2 and Matern-5/2 correlation funcs - * get_corrf.c, forsubs.h: included the new correlation funcs - * gpr_train.cc: documentation for the new correlation funcs - -2008-02-24 Jaroslav Hajek <hi...@gm...> - - * optdrv.f: implement objective reduction stopping criterion - * gprmod.h, train.c, gpr_train.cc: add support for the new feature - * forsubs.h nl0gpr.f: correct subroutine name - * train.c: correct CP array size - -2008-02-19 Jaroslav Hajek <hi...@gm...> - - * *.h *c *.cc: adjusted C/C++ sources to meet GNU coding standards - -2008-02-18 Jaroslav Hajek <hi...@gm...> - - * forsubs.h gprmod.h predict.c train.c setup.c: added const modifiers - where appropriate. - * grp_train, gpr_predict: replaced Array::fortran_vec with Array::data - where appropriate, to prevent unnecessary copying. - -2008-02-15 Jaroslav Hajek <hi...@gm...> - - * octgpr.m: improved demo. - -2008-02-13 Jaroslav Hajek <hi...@gm...> - - * train.c: correct allocation and copying result values. - * nldgpr.f: compute 2nd derivative by a better formula. - * gpr_train.cc: improved documentation - * gpr_predict.cc: improved documentation - -2008-02-12 Jaroslav Hajek <hi...@gm...> - - * dscrot.f: correct the W*D*x and W'*D*x sequence. - * optdrv.f: correct calls to dscrot.f - -2008-02-08 Jaroslav Hajek <hi...@gm...> - - * train.c: corrected C->Fortran calls in to use proper #defines from forsubs.h - * vmfac.f vmcmp.f trstep.f: deleted (no use) - * optdrv.f: cosmetic changes - * ../inst/demo_octgpr.m: added copyright - * ../TODO: modified - -2008-02-07 Jaroslav Hajek <hi...@gm...> - - * gpr_train.cc: corrected inline Texinfo documentation - * ChangeLog: added this file - Modified: trunk/octave-forge/main/octgpr/src/Makefile.in =================================================================== --- trunk/octave-forge/main/octgpr/src/Makefile.in 2008-04-17 08:09:27 UTC (rev 4948) +++ trunk/octave-forge/main/octgpr/src/Makefile.in 2008-04-17 13:53:44 UTC (rev 4949) @@ -35,7 +35,7 @@ OBJS_GPR_PRED=dsdacc.o dwdis2.o infgpr.o corrf.o \ get_corrf.o predict.o -all: gpr_train.oct gpr_predict.oct +all: gpr_train.oct gpr_predict.oct pdist2_mw.oct %.o: %.f $(F77) $(FFLAGS) -c $< @@ -51,6 +51,10 @@ $(MKOCTFILE) -o $@ gpr_train.o $(OBJS_GPR_TRAIN) $(LIBS) gpr_predict.oct: gpr_predict.o $(OBJS_GPR_PRED) $(MKOCTFILE) -o $@ gpr_predict.o $(OBJS_GPR_PRED) $(LIBS) + +pdist2_mw.oct: pdist2_mw.cc + $(MKOCTFILE) -o $@ $< + clean: rm -f *.o *.oct distclean: clean Modified: trunk/octave-forge/main/octgpr/src/configure.in =================================================================== --- trunk/octave-forge/main/octgpr/src/configure.in 2008-04-17 08:09:27 UTC (rev 4948) +++ trunk/octave-forge/main/octgpr/src/configure.in 2008-04-17 13:53:44 UTC (rev 4949) @@ -1,4 +1,4 @@ -s# Copyright (C) 2008 VZLU Prague, a.s., Czech Republic +# Copyright (C) 2008 VZLU Prague, a.s., Czech Republic # # Author: Jaroslav Hajek <hi...@gm...> # Added: trunk/octave-forge/main/octgpr/src/pdist2_mw.cc =================================================================== --- trunk/octave-forge/main/octgpr/src/pdist2_mw.cc (rev 0) +++ trunk/octave-forge/main/octgpr/src/pdist2_mw.cc 2008-04-17 13:53:44 UTC (rev 4949) @@ -0,0 +1,287 @@ +/* Copyright (C) 2008 VZLU Prague, a.s., Czech Republic + * + * Author: Jaroslav Hajek <hi...@gm...> + * + * This file is part of OctGPR. + * + * OctGPR is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this software; see the file COPYING. If not, see + * <http://www.gnu.org/licenses/>. */ + +#include <cmath> +#include <oct.h> +#include <oct-cmplx.h> + +inline double abs2(double x) +{ + return x*x; +} +inline double abs2(Complex x) +{ + return abs2(x. real()) + abs2(x.imag ()); +} + +/* distance functors */ +template<typename T> +struct distfun_eu +{ + double operator() (octave_idx_type dim, const T *x, const T* y) + { + double d = 1, scl = 0; + for (octave_idx_type i = 0; i < dim; i++) + { + double t = std::abs (x[i]-y[i]); + if (scl < t) + { + d *= abs2 (scl/t); + d += 1; + scl = t; + } + else if (t != 0) + d += abs2 (t/scl); + } + return scl * std::sqrt (d); + } +}; + +template<typename T> +struct distfun_sqeu +{ + double operator() (octave_idx_type dim, const T *x, const T* y) + { + double d = 0; + for (octave_idx_type i = 0; i < dim; i++) + d += abs2(x[i]-y[i]); + return d; + } +}; + +template<typename T> +struct distfun_l1 +{ + double operator() (octave_idx_type dim, const T *x, const T* y) + { + double d = 0; + for (octave_idx_type i = 0; i < dim; i++) + d += std::abs(x[i]-y[i]); + return d; + } +}; + +template<typename T> +struct distfun_max +{ + double operator() (octave_idx_type dim, const T *x, const T* y) + { + double d = 0; + for (octave_idx_type i = 0; i < dim; i++) + { + double t = std::abs(x[i]-y[i]); + if (t > d) d = t; + } + return d; + } +}; + +template<typename T> +struct distfun_mw +{ + double p; + distfun_mw (double _p) : p(_p) {} + double operator() (octave_idx_type dim, const T *x, const T* y) + { + double d = 1, scl = 0; + for (octave_idx_type i = 0; i < dim; i++) + { + double t = std::abs (x[i]-y[i]); + if (scl < t) + { + d *= std::pow (scl/t, p); + d += 1; + scl = t; + } + else if (t != 0) + d += std::pow (t/scl, p); + } + return scl * std::pow (d, 1/p); + } +}; + +template<typename T, class distfun> +void +fill_dist_matrix (octave_idx_type dim, octave_idx_type nx, octave_idx_type ny, + const T *X, const T* Y, double *D, distfun df) +{ + octave_idx_type i,j; + for (j = 0; j < ny; j++) + { + const T *PX = X; + for (i = 0; i < nx; i++) + { + *(D++) = df (dim, PX, Y); + PX += dim; + } + Y += dim; + } +} + +template<typename T, class distfun> +void +fill_dist_matrix (octave_idx_type dim, octave_idx_type nx, + const T *X, double *D, distfun df) +{ + octave_idx_type i,j; + for (j = 0; j < nx; j++) + { + D += j; + for (i = -j; i < 0; i++) + D[i] = *(D + i*nx); + const T *PX = X; + for (i = j; i < nx; i++) + { + *(D++) = df (dim, PX, X); + PX += dim; + } + X += dim; + } +} + +template<typename T> +Matrix +get_dist_matrix (const MArray2<T>& X, bool ssq, double p = 0) +{ + Matrix D(X.rows (), X.rows ()); + + MArray2<T> XT = X.transpose (); + + if (ssq) + fill_dist_matrix (XT.rows (), XT.cols (), + XT.data (), D.fortran_vec (), distfun_sqeu<T> ()); + else if (p == 2) + fill_dist_matrix (XT.rows (), XT.cols (), + XT.data (), D.fortran_vec (), distfun_eu<T> ()); + else if (p == 1) + fill_dist_matrix (XT.rows (), XT.cols (), + XT.data (), D.fortran_vec (), distfun_l1<T> ()); + else if (xisinf (p)) + fill_dist_matrix (XT.rows (), XT.cols (), + XT.data (), D.fortran_vec (), distfun_max<T> ()); + else + fill_dist_matrix (XT.rows (), XT.cols (), + XT.data (), D.fortran_vec (), distfun_mw<T> (p)); + return D; + +} + +template<typename T> +Matrix +get_dist_matrix (const MArray2<T>& X, const MArray2<T>& Y, bool ssq, double p = 0) +{ + Matrix D(X.rows (), Y.rows ()); + + MArray2<T> XT = X.transpose (), YT = Y.transpose (); + + if (ssq) + fill_dist_matrix (XT.rows (), XT.cols (), YT.cols (), + XT.data (), YT.data (), D.fortran_vec (), distfun_sqeu<T> ()); + else if (p == 2) + fill_dist_matrix (XT.rows (), XT.cols (), YT.cols (), + XT.data (), YT.data (), D.fortran_vec (), distfun_eu<T> ()); + else if (p == 1) + fill_dist_matrix (XT.rows (), XT.cols (), YT.cols (), + XT.data (), YT.data (), D.fortran_vec (), distfun_l1<T> ()); + else if (xisinf (p)) + fill_dist_matrix (XT.rows (), XT.cols (), YT.cols (), + XT.data (), YT.data (), D.fortran_vec (), distfun_max<T> ()); + else + fill_dist_matrix (XT.rows (), XT.cols (), YT.cols (), + XT.data (), YT.data (), D.fortran_vec (), distfun_mw<T> (p)); + return D; + +} + +DEFUN_DLD(pdist2_mw,args,, +"-*- texinfo -*-\n\ +@deftypefn {Loadable Function} @var{D} = pdist2_mw (@var{X}, @var{Y}, @var{p})\n\ +@deftypefnx {Loadable Function} @var{D} = pdist2_mw (@var{X}, @var{p})\n\ +Assembles a pairwise minkowski-distance matrix for two given sets of points.\n\ +@var{X} and @var{Y} should be real or complex matrices with a point per row,\n\ +so numbers of columns must match. The matrix contains the pairwise\n\ +distances @code{D(i,j) = norm(X(i,:)-Y(j,:),P)}.\n\ +@var{p} can also be the string \'ssq\' requesting squared euclidean distance.\n\ +(not a metric, but often useful and faster than @code{@var{p}=2})\n\ +If @var{Y} is not given, a symmetric distance matrix is calculated efficiently.\n\ +@seealso{norm}\n\ +@end deftypefn") +{ + int nargin = args.length(); + octave_value_list retval; + + if (nargin < 2 || nargin > 3) + { + print_usage (); + return retval; + } + + octave_value argx = args(0), argy, argp; + bool sym = false; + + if (nargin > 2) + argy = args(1); + else + { + argy = argx; + sym = true; + } + + if (nargin > 2) + argp = args(2); + else + argp = args(1); + + bool ssq = (argp.is_string () && argp.string_value () == "ssq"); + + if (argx.is_matrix_type () && argy.is_matrix_type () && (ssq || argp.is_real_scalar ())) + { + double p = ssq ? 0 : argp.scalar_value (); + + if (argx.columns () == argy.columns ()) + { + if (argx.is_real_matrix () && argy.is_real_matrix ()) + { + if (sym) + retval(0) = get_dist_matrix (argx.matrix_value (), + ssq, p); + else + retval(0) = get_dist_matrix (argx.matrix_value (), + argy.matrix_value (), + ssq, p); + } + else + { + if (sym) + retval(0) = get_dist_matrix (argx.complex_matrix_value (), + ssq, p); + else + retval(0) = get_dist_matrix (argx.complex_matrix_value (), + argy.complex_matrix_value (), + ssq, p); + } + } + else + error ("pmwdmat: dimension mismatch"); + } + else + error ("pmwdmat: X and Y should be matrices, p a real scalar"); + +} This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <hi...@us...> - 2008-04-18 12:56:00
|
Revision: 4957 http://octave.svn.sourceforge.net/octave/?rev=4957&view=rev Author: highegg Date: 2008-04-18 05:56:05 -0700 (Fri, 18 Apr 2008) Log Message: ----------- demo correction & cosmetic changes Modified Paths: -------------- trunk/octave-forge/main/octgpr/DESCRIPTION trunk/octave-forge/main/octgpr/TODO trunk/octave-forge/main/octgpr/inst/demo_octgpr.m Modified: trunk/octave-forge/main/octgpr/DESCRIPTION =================================================================== --- trunk/octave-forge/main/octgpr/DESCRIPTION 2008-04-18 09:11:18 UTC (rev 4956) +++ trunk/octave-forge/main/octgpr/DESCRIPTION 2008-04-18 12:56:05 UTC (rev 4957) @@ -6,8 +6,7 @@ Maintainer: Jaroslav Hajek (hi...@gm...) Description: The package allows interpolating and smoothing scattered multidimensional data using Gaussian Process Regression (also known - as Kriging). It features two functions: one for training (building - statistical model) and one for prediction (inference). + as Kriging). License: GPL v2 Depends: octave (>= 2.9.7) Modified: trunk/octave-forge/main/octgpr/TODO =================================================================== --- trunk/octave-forge/main/octgpr/TODO 2008-04-18 09:11:18 UTC (rev 4956) +++ trunk/octave-forge/main/octgpr/TODO 2008-04-18 12:56:05 UTC (rev 4957) @@ -2,5 +2,5 @@ o add documentation where missing o add more demos -o add RBF models with GCV training (in progress) +o add sparse gaussian process regression Modified: trunk/octave-forge/main/octgpr/inst/demo_octgpr.m =================================================================== --- trunk/octave-forge/main/octgpr/inst/demo_octgpr.m 2008-04-18 09:11:18 UTC (rev 4956) +++ trunk/octave-forge/main/octgpr/inst/demo_octgpr.m 2008-04-18 12:56:05 UTC (rev 4957) @@ -73,10 +73,6 @@ contourf (xi, yi, zi, 20); pause; - if (!exist ("nsamp", "var") || !isnumeric (nsamp)) - nsamp = 150; - endif - tit = sprintf ("sampled at %d random points", nsamp); disp (tit); % create random samples This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <hi...@us...> - 2008-05-16 04:54:33
|
Revision: 5054 http://octave.svn.sourceforge.net/octave/?rev=5054&view=rev Author: highegg Date: 2008-05-15 21:54:40 -0700 (Thu, 15 May 2008) Log Message: ----------- allow non-interactive demo Modified Paths: -------------- trunk/octave-forge/main/octgpr/ChangeLog trunk/octave-forge/main/octgpr/inst/demo_octgpr.m Modified: trunk/octave-forge/main/octgpr/ChangeLog =================================================================== --- trunk/octave-forge/main/octgpr/ChangeLog 2008-05-15 15:34:09 UTC (rev 5053) +++ trunk/octave-forge/main/octgpr/ChangeLog 2008-05-16 04:54:40 UTC (rev 5054) @@ -1,3 +1,7 @@ +2008-05-15 Jaroslav Hajek <hi...@gm...> + + * inst/demo_octgpr.m: add support for non-interactive running. + 2008-04-17 Jaroslav Hajek <hi...@gm...> * inst/demo_octgpr.m: refactor demo Modified: trunk/octave-forge/main/octgpr/inst/demo_octgpr.m =================================================================== --- trunk/octave-forge/main/octgpr/inst/demo_octgpr.m 2008-05-15 15:34:09 UTC (rev 5053) +++ trunk/octave-forge/main/octgpr/inst/demo_octgpr.m 2008-05-16 04:54:40 UTC (rev 5054) @@ -37,16 +37,43 @@ % @end itemize % @end deftypefn function demo_octgpr (number, varargin) - switch (number) - case 1 - demo_octgpr1 (varargin{:}) - case 2 - demo_octgpr2 (varargin{:}) - otherwise - error ("demo_octgpr: invalid demo number") - endswitch + + global prntfmt = ''; + + if (nargin < 1) + print_usage (); + elseif (ischar (number)) + prntfmt = number; + elseif (isscalar (number)) + figure (); + if (! isempty (prntfmt)) + figure (gcf, "visible", "off"); + endif + + switch (number) + case 1 + demo_octgpr1 (varargin{:}) + case 2 + demo_octgpr2 (varargin{:}) + otherwise + error ("demo_octgpr: invalid demo number") + endswitch + else + print_usage (); + endif + endfunction +function demo_octgpr_pause (idemo, iplot) + global prntfmt; + if (isempty (prntfmt)) + pause; + else + print (sprintf (prntfmt, idemo, iplot)); + fflush (stdout); + endif +endfunction + % define the test function (the well-known matlab "peaks" plus some sines) function z = testfun1 (x, y) z = 4 + 3 * (1-x).^2 .* exp(-(x.^2) - (y+1).^2) ... @@ -56,6 +83,7 @@ endfunction function demo_octgpr1 (nsamp = 150) + global prntfmt; tit = "a peaked surface"; disp (tit); @@ -71,7 +99,7 @@ title (tit); subplot (2, 2, 3); contourf (xi, yi, zi, 20); - pause; + demo_octgpr_pause (1, 1); tit = sprintf ("sampled at %d random points", nsamp); disp (tit); @@ -87,7 +115,7 @@ title (tit); subplot (2, 2, 4); plot (xs, ys, ".+"); - pause + demo_octgpr_pause (1, 2); tit = "GPR model with heuristic hypers"; disp (tit); @@ -104,7 +132,7 @@ contourf (xi, yi, zm, 20); plot (xs, ys, "+6"); hold off - pause + demo_octgpr_pause (1, 3); tit = "GPR model with MLE training"; disp (tit); @@ -121,20 +149,21 @@ contourf (xi, yi, zm, 20); plot (xs, ys, "+6"); hold off - pause + demo_octgpr_pause (1, 4); + close - close endfunction function demo_octgpr2 (ncnt = 50, npt = 500) + global prntfmt; npt = ncnt*ceil (npt/ncnt); U = rand (ncnt, 2); cs = min (pdist2_mw (U, 2) + diag (Inf (ncnt, 1))); X = repmat (U, npt/ncnt, 1) + repmat (cs', npt/ncnt, 2) .* randn (npt, 2); disp ("slightly clustered random points") plot (X(:,1), X(:,2), "+"); - pause + demo_octgpr_pause (2, 1); [U, ur] = rbf_centers(X, ncnt); @@ -147,7 +176,7 @@ line (xc, yc); endfor hold off - pause + demo_octgpr_pause (2, 2); close endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <hi...@us...> - 2008-05-30 06:07:44
|
Revision: 5085 http://octave.svn.sourceforge.net/octave/?rev=5085&view=rev Author: highegg Date: 2008-05-29 23:07:51 -0700 (Thu, 29 May 2008) Log Message: ----------- add normalizing constant to likelihoods Modified Paths: -------------- trunk/octave-forge/main/octgpr/ChangeLog trunk/octave-forge/main/octgpr/src/nl0gpr.f trunk/octave-forge/main/octgpr/src/nllgpr.f Modified: trunk/octave-forge/main/octgpr/ChangeLog =================================================================== --- trunk/octave-forge/main/octgpr/ChangeLog 2008-05-27 16:23:05 UTC (rev 5084) +++ trunk/octave-forge/main/octgpr/ChangeLog 2008-05-30 06:07:51 UTC (rev 5085) @@ -1,3 +1,8 @@ +2008-05-30 Jaroslav Hajek <hi...@gm...> + + * src/nllgpr.f, src/nl0gpr.f: add normalizing constant to obtain the + true negative log likelihood. + 2008-05-15 Jaroslav Hajek <hi...@gm...> * inst/demo_octgpr.m: add support for non-interactive running. Modified: trunk/octave-forge/main/octgpr/src/nl0gpr.f =================================================================== --- trunk/octave-forge/main/octgpr/src/nl0gpr.f 2008-05-27 16:23:05 UTC (rev 5084) +++ trunk/octave-forge/main/octgpr/src/nl0gpr.f 2008-05-30 06:07:51 UTC (rev 5085) @@ -32,6 +32,7 @@ integer nx double precision y(nx),nu,nll0,nllinf double precision mu,ssq + parameter (l2pi = 1.83787706640935d0) integer i mu = 0 @@ -47,7 +48,7 @@ end do c set values ssq = ssq / nx - nllinf = 0.5d0 * nx * log(ssq) + nllinf = 0.5d0 * nx * (log(ssq) + l2pi) nll0 = nllinf + 0.5d0 * log(1 + nx/nu**2) end subroutine Modified: trunk/octave-forge/main/octgpr/src/nllgpr.f =================================================================== --- trunk/octave-forge/main/octgpr/src/nllgpr.f 2008-05-27 16:23:05 UTC (rev 5084) +++ trunk/octave-forge/main/octgpr/src/nllgpr.f 2008-05-30 06:07:51 UTC (rev 5085) @@ -66,8 +66,10 @@ external dwdis2,dcopy,daxpy,dscal,dtrsv,dtrsm,xerbla, +dpotrf,dtrtrs,dgeqr2,dorm2r integer i,j,nl1 - real*8 sums,sum,dwdis2 + real*8 sums,sum,dwdis2,l2pi2 + parameter (l2pi = 1.83787706640935d0) + c argument checks info = 0 if (ndim < 0) then @@ -166,7 +168,7 @@ end do c final negative log likelihood - nll = sum + 0.5d0 * nx*log(var) + nll = sum + 0.5d0 * nx*(log(var) + l2pi) c normal return info = 0 return This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <hi...@us...> - 2008-06-03 06:51:29
|
Revision: 5086 http://octave.svn.sourceforge.net/octave/?rev=5086&view=rev Author: highegg Date: 2008-06-02 23:51:35 -0700 (Mon, 02 Jun 2008) Log Message: ----------- improve demo Modified Paths: -------------- trunk/octave-forge/main/octgpr/ChangeLog trunk/octave-forge/main/octgpr/inst/demo_octgpr.m Modified: trunk/octave-forge/main/octgpr/ChangeLog =================================================================== --- trunk/octave-forge/main/octgpr/ChangeLog 2008-05-30 06:07:51 UTC (rev 5085) +++ trunk/octave-forge/main/octgpr/ChangeLog 2008-06-03 06:51:35 UTC (rev 5086) @@ -1,3 +1,8 @@ +2008-06-03 Jaroslav Hajek <hi...@gm...> + + * inst/demo_octgpr.m: display points also on bottom left subplot, + cosmetic fixes. + 2008-05-30 Jaroslav Hajek <hi...@gm...> * src/nllgpr.f, src/nl0gpr.f: add normalizing constant to obtain the Modified: trunk/octave-forge/main/octgpr/inst/demo_octgpr.m =================================================================== --- trunk/octave-forge/main/octgpr/inst/demo_octgpr.m 2008-05-30 06:07:51 UTC (rev 5085) +++ trunk/octave-forge/main/octgpr/inst/demo_octgpr.m 2008-06-03 06:51:35 UTC (rev 5086) @@ -113,6 +113,10 @@ subplot (2, 2, 2); plot3 (xs, ys, zs, ".+"); title (tit); + subplot (2, 2, 3); + hold on + plot (xs, ys, "+6"); + hold off subplot (2, 2, 4); plot (xs, ys, ".+"); demo_octgpr_pause (1, 2); @@ -127,7 +131,7 @@ subplot (2, 2, 2); mesh (xi, yi, zm); title (tit); - subplot(2, 2, 4) + subplot(2, 2, 4); hold on contourf (xi, yi, zm, 20); plot (xs, ys, "+6"); @@ -144,7 +148,7 @@ subplot (2, 2, 2); mesh (xi, yi, zm); title (tit); - subplot(2, 2, 4) + subplot(2, 2, 4); hold on contourf (xi, yi, zm, 20); plot (xs, ys, "+6"); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <hi...@us...> - 2008-06-10 13:24:54
|
Revision: 5099 http://octave.svn.sourceforge.net/octave/?rev=5099&view=rev Author: highegg Date: 2008-06-10 06:24:22 -0700 (Tue, 10 Jun 2008) Log Message: ----------- fix invalid flag in configure.in Modified Paths: -------------- trunk/octave-forge/main/octgpr/ChangeLog trunk/octave-forge/main/octgpr/src/configure trunk/octave-forge/main/octgpr/src/configure.in Modified: trunk/octave-forge/main/octgpr/ChangeLog =================================================================== --- trunk/octave-forge/main/octgpr/ChangeLog 2008-06-10 07:16:42 UTC (rev 5098) +++ trunk/octave-forge/main/octgpr/ChangeLog 2008-06-10 13:24:22 UTC (rev 5099) @@ -1,3 +1,9 @@ +2008-06-10 Jaroslav Hajek <hi...@gm...> + + * src/configure.in: correct FC to F77 (inferred variable from + mkoctfile). + * src/configure: regenerate. + 2008-06-03 Jaroslav Hajek <hi...@gm...> * inst/demo_octgpr.m: display points also on bottom left subplot, Modified: trunk/octave-forge/main/octgpr/src/configure =================================================================== --- trunk/octave-forge/main/octgpr/src/configure 2008-06-10 07:16:42 UTC (rev 5098) +++ trunk/octave-forge/main/octgpr/src/configure 2008-06-10 13:24:22 UTC (rev 5099) @@ -1736,7 +1736,7 @@ fi # get preferred compilers -F77=`$MKOCTFILE -p FC` +F77=`$MKOCTFILE -p F77` CC=`$MKOCTFILE -p CC` # get preferred flags FFLAGS=`$MKOCTFILE -p FFLAGS` Modified: trunk/octave-forge/main/octgpr/src/configure.in =================================================================== --- trunk/octave-forge/main/octgpr/src/configure.in 2008-06-10 07:16:42 UTC (rev 5098) +++ trunk/octave-forge/main/octgpr/src/configure.in 2008-06-10 13:24:22 UTC (rev 5099) @@ -32,7 +32,7 @@ fi # get preferred compilers -F77=`$MKOCTFILE -p FC` +F77=`$MKOCTFILE -p F77` CC=`$MKOCTFILE -p CC` # get preferred flags FFLAGS=`$MKOCTFILE -p FFLAGS` This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <hi...@us...> - 2008-06-19 06:16:48
|
Revision: 5120 http://octave.svn.sourceforge.net/octave/?rev=5120&view=rev Author: highegg Date: 2008-06-18 23:16:48 -0700 (Wed, 18 Jun 2008) Log Message: ----------- add missing return Modified Paths: -------------- trunk/octave-forge/main/octgpr/ChangeLog trunk/octave-forge/main/octgpr/src/pdist2_mw.cc Modified: trunk/octave-forge/main/octgpr/ChangeLog =================================================================== --- trunk/octave-forge/main/octgpr/ChangeLog 2008-06-18 08:33:11 UTC (rev 5119) +++ trunk/octave-forge/main/octgpr/ChangeLog 2008-06-19 06:16:48 UTC (rev 5120) @@ -1,3 +1,7 @@ +2008-06-19 Jaroslav Hajek <hi...@gm...> + + * src/pdist2_mw.cc: add missing return statement. + 2008-06-10 Jaroslav Hajek <hi...@gm...> * src/configure.in: correct FC to F77 (inferred variable from Modified: trunk/octave-forge/main/octgpr/src/pdist2_mw.cc =================================================================== --- trunk/octave-forge/main/octgpr/src/pdist2_mw.cc 2008-06-18 08:33:11 UTC (rev 5119) +++ trunk/octave-forge/main/octgpr/src/pdist2_mw.cc 2008-06-19 06:16:48 UTC (rev 5120) @@ -284,4 +284,5 @@ else error ("pmwdmat: X and Y should be matrices, p a real scalar"); + return retval; } This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <hi...@us...> - 2008-07-03 11:43:49
|
Revision: 5156 http://octave.svn.sourceforge.net/octave/?rev=5156&view=rev Author: highegg Date: 2008-07-03 04:43:56 -0700 (Thu, 03 Jul 2008) Log Message: ----------- fix invalid parameter declaration (corrects the nll normalization) Modified Paths: -------------- trunk/octave-forge/main/octgpr/ChangeLog trunk/octave-forge/main/octgpr/DESCRIPTION trunk/octave-forge/main/octgpr/src/nllgpr.f Modified: trunk/octave-forge/main/octgpr/ChangeLog =================================================================== --- trunk/octave-forge/main/octgpr/ChangeLog 2008-07-02 19:02:04 UTC (rev 5155) +++ trunk/octave-forge/main/octgpr/ChangeLog 2008-07-03 11:43:56 UTC (rev 5156) @@ -1,3 +1,7 @@ +2008-07-03 Jaroslav Hajek <hi...@gm...> + + * src/nllgpr.f: fix invalid parameter declaration + 2008-06-19 Jaroslav Hajek <hi...@gm...> * src/pdist2_mw.cc: add missing return statement. Modified: trunk/octave-forge/main/octgpr/DESCRIPTION =================================================================== --- trunk/octave-forge/main/octgpr/DESCRIPTION 2008-07-02 19:02:04 UTC (rev 5155) +++ trunk/octave-forge/main/octgpr/DESCRIPTION 2008-07-03 11:43:56 UTC (rev 5156) @@ -1,5 +1,5 @@ Name: OctGPR -Version: 1.1.3 +Version: 1.1.4 Date: 2008-04-29 Author: Jaroslav Hajek (hi...@gm...) Title: Package for full dense Gaussian Process Regression Modified: trunk/octave-forge/main/octgpr/src/nllgpr.f =================================================================== --- trunk/octave-forge/main/octgpr/src/nllgpr.f 2008-07-02 19:02:04 UTC (rev 5155) +++ trunk/octave-forge/main/octgpr/src/nllgpr.f 2008-07-03 11:43:56 UTC (rev 5156) @@ -66,7 +66,7 @@ external dwdis2,dcopy,daxpy,dscal,dtrsv,dtrsm,xerbla, +dpotrf,dtrtrs,dgeqr2,dorm2r integer i,j,nl1 - real*8 sums,sum,dwdis2,l2pi2 + real*8 sums,sum,dwdis2,l2pi parameter (l2pi = 1.83787706640935d0) This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |