//  autodif.sci by Shin'ichi OISHI
//  automatic gradient calculator
//  1999 11.21
//  1999 12.20 modified by OISHI
//  v should be an [m,1] tlist mat type
//  tlist dif has the form
//     tlist(['dif','dim','val','grad'],m,val,grad)
//        val is a scaler and grad an [1,m]-mat-tlist.

function x=init_dif(a)
   n=a('dim');
   m=n(1);
   s='tlist(["'dif'","'dim'","'val'","'grad'"],m,0,0)';
   ss=evstr(s);
   x=mat(m,1);
   for i=1:m,
      x(i+2)=ss;
      x(i+2)('val')=a(i+2);
      x(i+2)('grad')=mat(1,m);
      ii=i+2;
      x(i+2)('grad')(ii)=1;
   end;
endfunction  


function x=%dif_a_dif(a,b)
   x=a;
   x('val')=a('val')+b('val');
   x('grad')=a('grad')+b('grad'); 
endfunction 


function x=%dif_s_dif(a,b)
   x=a;
   x('val')=a('val')-b('val');
   x('grad')=a('grad')-b('grad'); 
endfunction 


function x=%s_a_dif(a,b)
   x=b;
   x('val')=a+b('val');
endfunction 


function x=%dif_a_s(a,b)
   x=a;
   x('val')=a('val')+b;
endfunction 


function x=%s_s_dif(a,b)
   x=b;
   x('val')=a-b('val');
   x('grad')=(-1)*b('grad');
endfunction 


function x=%dif_s_s(a,b)
   x=a;
   x('val')=a('val')-b;
endfunction 


function x=%s_m_dif(a,b)
   x=b;
   x('val')=a*b('val');
   x('grad')=a*b('grad');
endfunction 


function x=%dif_m_s(a,b)
   x=a;
   x('val')=b*a('val');
   x('grad')=b*a('grad');
endfunction 


function x=%dif_m_dif(a,b)
   x=a;
   x('val')=a('val')*b('val');
   x('grad')=b('val')*a('grad')+a('val')*b('grad'); 
endfunction 


function x=%s_r_dif(a,b)
   x=b; 
   x('val')=a/b('val');
   d=-a*b('grad'); 
   x('grad')=d/(b('val')*b('val'));
endfunction 


function x=%dif_r_s(a,b)
   x=a;
   x('val')=a('val')/b;
   x('grad')=a('grad')/b; 
endfunction 


function x=%dif_r_dif(a,b)
   x=a;
   x('val')=a('val')/b('val');
   d=a('grad')*b('val')-a('val')*b('grad'); 
   x('grad')=d/(b('val')*b('val'));
endfunction 


function x=%dif_sin(a)
   x=a;
   x('val')=sin(a('val'));
   x('grad')=cos(a('val'))*a('grad');
endfunction


function x=%dif_cos(a)
   x=a;
   x('val') = cos(a('val'));
   x('grad') = (-1)*sin(a('val'))*a('grad');
endfunction
  

function x=%dif_exp(a)
   x=a;
   x('val')=exp(a('val'));
   x('grad')=exp(a('val'))*a('grad');
endfunction


function x=%dif_log(a)
   x=a;
   x('val')=log(a('val'));
   x('grad')=a('grad')/a('val');
endfunction


function x=jacobi(a)
   s=a('dim');
   ss=s(1);
   x=eye(ss,ss);
   for i=1:ss,
      for j=1:ss,
         x(i,j)=a(i+2)(4)(j+2);
      end;
   end;
endfunction


function [x,y]=int_jacobi(a)
   s=a('dim');
   ss=s(1);
   x=eye(ss,ss);
   y=eye(ss,ss);
   for i=1:ss,
      for j=1:ss,
         x(i,j)=a(i+2)(4)(j+2)(2);
         y(i,j)=a(i+2)(4)(j+2)(3);
      end;
   end;
endfunction


function x=val(a)
   s=a('dim');
   ss=s(1);
   x=eye(ss,1);
   for i=1:ss;
      x(i)=a(i+2)(3);
   end;
endfunction