Re: nonconformant arguments in leasqr.m

[Olaf Till]

On Wed, May 04, 2011 at 03:29:53PM +0200, Olga Kruglova wrote:
> On Wed, May 4, 2011 at 3:04 PM, Olaf Till <address@hidden> wrote:
> 
>     On Wed, May 04, 2011 at 01:57:59PM +0200, Olga Kruglova wrote:
>     > On Wed, May 4, 2011 at 1:55 PM, Andy Buckle <address@hidden> wrote:
>     >
>     > > On Wed, May 4, 2011 at 12:50 PM, olga_kruglova <address@hidden>
>     wrote:
>     > > > Hello everybody,
>     > > >
>     > > > I try to use leasqr.m to fit data with function convoluted with 
> known
>     > > > resolution function but it gives me an error:
>     > > >>> Convolution
>     > > > Variables in the current scope:
>     > > >
>     > > >  Attr Name        Size                     Bytes  Class
>     > > >  ==== ====        ====                     =====  =====
>     > > >   f   y        1024x1                       8192  double
>     > > >
>     > > > Total is 1024 elements using 8192 bytes
>     > > >
>     > > > Variables in the current scope:
>     > > >
>     > > >  Attr Name        Size                     Bytes  Class
>     > > >  ==== ====        ====                     =====  =====
>     > > >       data     1024x1                       8192  double
>     > > >
>     > > > Total is 1024 elements using 8192 bytes
>     > > >
>     > > > Variables in the current scope:
>     > > >
>     > > >  Attr Name        Size                     Bytes  Class
>     > > >  ==== ====        ====                     =====  =====
>     > > >   f   wt          4x1                         32  double
>     > > >
>     > > > Total is 4 elements using 32 bytes
>     > > >
>     > > > Variables in the current scope:
>     > > >
>     > > >  Attr Name        Size                     Bytes  Class
>     > > >  ==== ====        ====                     =====  =====
>     > > >   f   y        1024x1                       8192  double
>     > > >
>     > > > Total is 1024 elements using 8192 bytes
>     > > >
>     > > > error: product: nonconformant arguments (op1 is 4x1, op2 is 1024x1)
>     > > > error: called from:
>     > > > error:   C:\OCTAVE~2\leasqr.m at line 184, column 3
>     > > > error:   C:\OCTAVE~2\Convolution.m at line 25, column 67
>     > > >  I can not understand why. Could you please help me out with it?
>     Thank
>     > > you
>     > > > very much in advance!
>     > > >
>     > > > Kind regards, Olga
>     > >
>     > > Is Convolution.m your script?
>     > >
>     > > I think that in order to help, we have to see line 25 of 
> Convolution.m.
>     > >
>     > > --
>     > > /* andy buckle */
>     > >
>     >
>     > Yes, it is
>     >                  global verbose;
>     >                   verbose=1;
>     >  line 25:           [f1, p1, kvg1, iter1, corp1, covp1, covr1, 
> stdresid1,
>     > Z1, r21] = ...
>     >                         leasqr (t, data, pin, F, stol, niter, options);
> 
>     In order to reproduce the failure, we have to know the values of all
>     these variables: (t, data, pin, F, stol, niter, options), and probably
>     have to know the code of your function F.
> 
>     But first: It seems you are using Matlab, and as far as I know the
>     current version of our 'leasqr' does not run in Matlab at all. Where
>     do you have your 'leasqr' from? If from Octave-Forge, from which
>     version of the optim package?
> 
>     Olaf
> 
> 
> I use octave 3.2.4 i686-pc-mingw32 gcc-4.4.0 on WindowsXP SP3, not Matlab. 
> Here
> is Convolution.m:
> function Convolution
> s = load ('vesbufIV.txt');
> t = s(:,1); t = 0.55e-9.*t;
> z = s(:,3);
> data = leasqrfunc(t, z);
> whos data;
> q = s (:,2);
> 
> 
>   F = @leasqrfunc;
>   %dFdp = @leasqrdfdp; % exact derivative
>   %dFdp = @dfdp;     % estimated derivative
>   %dp = [0.001; 0.001];
>   pin = [0.5; 1.0e9];
>   stol=0.001; niter=50;
>   minstep = [0.01; 0.01e9];
>   maxstep = [0.8; 0.8e9];
>   options = [minstep, maxstep];
> 
> 
>   global verbose;
>   verbose=1;
>   [f1, p1, kvg1, iter1, corp1, covp1, covr1, stdresid1, Z1, r21] = ...
>     leasqr (t, data, pin, F, stol, niter, options);
>    
> function y = leasqrfunc(omega,p);
> %Period
>    s = load ('vesbufIV.txt');
>    t = s(:,1); t = 0.55e-9.*t;
>    omega=zeros(1024,1);
>    omega = 2*pi./ t;
> 
> %Convolution
> 
>  
> %nh = length(h);
> 
>   ML = zeros(1024,1);
>    for i=1:1024, ML = (2*sin(pi*p(1)/2)./(omega.*p(2)).^p(1)+2*cos(pi*p(1)/2)+
> (omega.*p(2)).^(-p(1))); end;
>    
> 
> %nYf = length(ML);
> %nfft = 2.^nextpow2(nYf+nh-1);
> 
> %Yfzp = [ML, zeros(nfft-nYf,1)];
> %hzp = [h, zeros(nfft-nh,1)];
>   s = load ('vesbufIV.txt');
>   q = s (:,2);
>   Hf = zeros(1024,1);
>   Hf = fft(q);
> 
> % Restored signal
>   eXf = Hf.*ML; 
> 
>   format bank;
> 
>   Xt =real( ifft(eXf));
> 
>   y =  abs(Xt(1:1024,1));
>   whos y;
>   save("-ascii","test.txt","y");
> 
> endfunction
> 
> Thank you very much for help!
> 
> --
> Kind regards, Olga

This would be not enough to reproduce the failure since we do not know
what is in 'vesbufIV.txt'.

But I see only now --- sorry --- that the call to leasqr is
wrong. Leave off the argument 'options' (do you really already know
that you need to specify a sufficient precision and a maximum
fractional step?), or if you provide it, it must be the 10th argument
(not the 7th), so try:

[f1, p1, kvg1, iter1, corp1, covp1, covr1, stdresid1, Z1, r21] = \
leasqr (t, data, pin, F, stol, niter, ones (2, 1), .001 * ones (2, 1), \
"dfdp", options);

I guess this was the reason for the failure.

By the way, according to your code, you should note a few things:

It's no use first assigning zeros(..) and then a different value to a
variable. First assigning zeros(..) to a variable is only more
efficient if then its single elements are assigned to in a loop (which
in many cases can be avoided).

In 'leasqrfunc' you have a 'for loop' which does the same thing 1024
times; the index 'i' is never used in it.

Using 'load' in 'leasqrfunc' is inefficient and can be avoided, but
let's not make things too complicated for now. But at least you should
not load the same file twice in 'leasqrfunc'.

Olaf