#COPYRIGHT NOTICE:  Copyright (c) 1995-98 by Sasha Cyganowski.  Permission is 
 #granted to anyone to use, modify, or redistribute this software freely, 
 #subject to the following restrictions:  1. The author accepts no
 #responsibility for any consequences of this software and makes no guarantee
 #that this software is free of defects.  2.  The origin of this software must 
 #not be misrepresented, either by explicit claim or by omission. 
 #3.  This notice and the copyright must be included with all copies or
 #modified versions of this software.  4.  This software may not be included or
 #redistributed as part of any package to be sold for profit unlesss the author
 #has given explicit written permission to do so.  
 # Sasha Cyganowski, School of Computing and Mathematics, Deakin University
 #
 # This file contains the function definitions for the package
 #"stochastic", developed by Sasha Cyganowski.  The package contains functions
 #which can be used to find explicit solutions of Stochastic Differential
 #Equations (SDEs), construct numerical schemes up to strong order 2 and weak
 #order 3, check for commutative noise of the first and second kind, and
 #convert SDEs into their coloured noise form.  Other useful routines are also
 #included.        
 #The general schemes, solutions, and conditions from which the code
 #was developed can be found in Kloeden, P.E, Platen, E.:  Numerical
 #Solution of Stochastic Differential Equations, (Springer-Verlag, 1992).
 #
 #
 #A manual for this package is available from 
 #
 #    http://www.math.uni-frankfurt.de/~numerik/maplestoch/
 #
 #This version of the software is known to work for Maple 9 and 9.5
 #
 #Additional Routines by Thomas Pohl, updates and bugfixes by Lars Gruene
 #and Stefanie Schoenbrodt
 #
 #For questions, comments etc. please contact Lars Gruene 
 #(http://www.uni-bayreuth.de/departments/math/~lgruene/)
 #
 #Last Change: June 21, 2005



stochastic[linearsde]:=proc(a::algebraic,b::algebraic)
local temp1,alpha,beta,gamma,delta,fundsoln,fundsoln2,soln,default1,default2,
default3;
    if diff(a,x,x) <> 0 or diff(b,x,x) <> 0 then 
    ERROR(`SDE not linearsde, try a reducible procedure`)
    else
        alpha := diff(a,x);
        alpha := subs(t = s,alpha);
        beta := diff(b,x);
        beta := subs(t = s,beta);
        if diff(beta,s) = 0 then temp1 := beta*W;
	else temp1:=Int(beta,W = 0 .. t);
	fi;
        gamma := coeff(a,x,0);
        gamma := subs(t = s,gamma);
        delta := coeff(b,x,0);
        delta := subs(t = s,delta);
        fundsoln := exp(int(alpha-1/2*beta^2,s = 0 .. t)+temp1);
        fundsoln2 := subs(t = s,fundsoln);
        if beta = 0 then
            soln := fundsoln*(X[0]+int(1/fundsoln2*(gamma-beta*delta),
	    s = 0 .. t)+Int(1/fundsoln2*delta,W = 0 .. t))
        else soln := fundsoln*(X[0]+Int(1/fundsoln2*(gamma-beta*delta),
	s = 0 .. t)+Int(1/fundsoln2*delta,W = 0 .. t))
        fi;
        default1 := Int(0,W = 0 .. t) = 0;
        default2 := Int(0,W = 0 .. s) = 0;
        default3 := Int(0,s = 0 .. t) = 0;
        soln := X[t] = subs(default1,default2,default3,soln);
    fi;
end:




stochastic[reducible]:=proc(a::algebraic,b::algebraic)
 local beta,temp1,h,temp3,alpha,soln,soln1,II;
    h := int(1/b,x);
    temp1 := alpha*b*h+1/2*b*simplify(diff(b,x));
    temp1 = a;
    alpha := simplify(solve(%,alpha));
    beta := alpha*h;
    if diff(alpha,x) = 0 then
	if alpha=0 then
	    soln:=h=subs(x=X[0],h)+W;
	    X[t]=simplify(solve(soln,x));
        else
        soln1 := h = exp(alpha*t)*subs(x = X[0],h)+exp(alpha*t)*II;
        X[t] = subs(II=Int(exp(-alpha*s),W=0..t),solve(soln1,x));
        fi
    elif diff(beta,x) = 0 then
    X[t]=simplify(solve(h = beta*t+W+subs(x = X[0],h),x));
    else ERROR(`non-linearsde SDE not reducible`)
    fi
end:


stochastic[explicit]:=
proc(a::algebraic,b::algebraic) if diff(a,x,x) = 0 and diff(b,x,x) = 0 then linearsde
(a,b) else reducible(a,b) fi end:


stochastic[LJ]:=
proc(X::algebraic,b::list(list(algebraic)),j::integer)  add(op(j,op(k,b))*
diff(X,x[k]),k = 1 .. nops(b)) end:


stochastic[L0]:=proc(X::algebraic,a::list(algebraic),b::list(list(algebraic)))
local part1,part2,part3;
    part1 := diff(X,t);
    part2 := add(a[k]*diff(X,x[k]),k = 1 .. nops(a));
    part3 := 1/2*add(
    add(add(op(j,op(k,b))*op(j,op(l,b))*diff(X,x[k],x[l]),
    j = 1 .. nops(op(1,b))),k = 1 .. nops(a)),l = 1 .. nops(a)
    );
    part1+part2+part3;
end:


stochastic[Euler]:=proc(a::list(algebraic),b::list(list(algebraic)))
local i,u,soln;
    for i to nops(a) do
        soln[i] := Y||i[n+1] = Y||i[n]+L0(x[i],a,b)*Delta[n]+add(LJ(x[i],b,j)*
	DeltaW||j[n],j = 1 .. nops(op(1,b)));
        for u to nops(a) do  soln[i] := subs(x[u] = Y||u[n],soln[i]) od
    od;
    RETURN(eval(soln))
end:


stochastic[Milstein]:=proc(a::list(algebraic),b::list(list(algebraic)))
local u,i,soln;
    for i to nops(a) do
        soln[i] := Y||i[n+1] = Y||i[n]+L0(x[i],a,b)*Delta[n]+add(LJ(x[i],b,j)*
	DeltaW||j[n],j = 1 .. nops(op(1,b)))+
            add(add(LJ(op(j2,op(i,b)),b,j1)*I[j1,j2],
	    j1 = 1 .. nops(op(1,b))),j2 = 1 .. nops(op(1,b)));
        for u to nops(a) do  soln[i] := subs(x[u] = Y||u[n],soln[i]) od
    od;
    RETURN(eval(soln))
end:

 
stochastic[Taylor1hlf]:=proc(a::list(algebraic),b::list(list(algebraic)))
local u,i,soln;
    for i to nops(a) do
        soln[i] := Y||i[n+1] = Y||i[n]+a[i]*Delta[n]+1/2*L0(a[i],a,b)*Delta[n]^2+
            add(op(j,op(i,b))*DeltaW||j[n]+L0(op(j,op(i,b)),a,b)*I[0,j]+
	    LJ(a[i],b,j)*I[j,0],j = 1 .. nops(op(1,b)))+
            add(add(LJ(op(j2,op(i,b)),b,j1)*I[j1,j2],
	    j1 = 1 .. nops(op(1,b))),j2 = 1 .. nops(op(1,b)))+add(
            add(add(LJ(LJ(op(p3,op(i,b)),b,p2),b,p1)*I[p1,p2,p3],
	    p1 = 1 .. nops(op(1,b))),p2 = 1 .. nops(op(1,b))),
            p3 = 1 .. nops(op(1,b)));
        for u to nops(a) do  soln[i] := subs(x[u] = Y||u[n],soln[i]) od
    od;
    RETURN(eval(soln))
end:

 
stochastic[SL0]:=proc(X::algebraic,a::list(algebraic),b::list(list(algebraic)))
local part1,part2;
    part1 := diff(X,t); part2 := add(a[k]*diff(X,x[k]),k = 1 .. nops(a));
    part1+part2;
end:

 
stochastic[correct]:=proc(a::list(algebraic),b::list(list(algebraic)),i) a[i]-
 1/2*add(LJ(op(j,op(i,b)),b,j),j = 1 .. nops(op(1,b))); end:


stochastic[Taylor2]:=
proc(a::list(algebraic),b::list(list(algebraic)))
local u,i,soln;
    for i to nops(a) do
        soln[i] := Y||i[n+1] = Y||i[n]+correct(a,b,i)*Delta[n]+
	1/2*SL0(correct(a,b,i),a,b)*Delta[n]^2+
          add(op(j,op(i,b))*DeltaW||j[n]+SL0(op(j,op(i,b)),a,b)*J[0,j]+
	  LJ(correct(a,b,i),b,j)*J[j,0],j = 1 .. nops(op(1,b)))
          +add(add(LJ(op(j2,op(i,b)),b,j1)*J[j1,j2]+
	  SL0(LJ(op(j2,op(i,b)),b,j1),a,b)*J[0,j1,j2]+
          LJ(SL0(op(j2,op(i,b)),a,b),b,j1)*J[j1,0,j2]+
	  LJ(LJ(correct(a,b,i),b,j2),b,j1)*J[j1,j2,0],
	  j1 = 1 .. nops(op(1,b))),
          j2 = 1 .. nops(op(1,b)))+add(
          add(add(LJ(LJ(op(p3,op(i,b)),b,p2),b,p1)*J[p1,p2,p3],
	  p1 = 1 .. nops(op(1,b))),p2 = 1 .. nops(op(1,b))),
          p3 = 1 .. nops(op(1,b)))+add(add(add(
          add(LJ(LJ(LJ(op(m4,op(i,b)),b,m3),b,m2),b,m1)*J[m1,m2,m3,m4],
	  m1 = 1 .. nops(op(1,b))),m2 = 1 .. nops(op(1,b)))
          ,m3 = 1 .. nops(op(1,b))),m4 = 1 .. nops(op(1,b)));
        for u to nops(a) do  soln[i] := subs(x[u] = Y||u[n],soln[i]) od
    od;
    RETURN(eval(soln))
end:

 
stochastic[wkeuler]:=
proc(a::list(algebraic),b::list(list(algebraic)))
local u,i,soln;
    for i to nops(a) do
        soln[i] := Y||i[n+1] = Y||i[n]+L0(x[i],a,b)*Delta[n]+
	add(LJ(x[i],b,j)*DeltaWs||j[n],j = 1 .. nops(op(1,b)));
        for u to nops(a) do  soln[i] := subs(x[u] = Y||u[n],soln[i]) od
    od;
    RETURN(eval(soln))
end:


stochastic[wktay2]:=
proc(a::list(algebraic),b::list(list(algebraic)))
local u,i,soln;
    for i to nops(a) do
        soln[i] := Y||i[n+1] = Y||i[n]+a[i]*Delta[n]+1/2*L0(a[i],a,b)*Delta[n]^2+
            add((op(j,op(i,b))+1/2*Delta[n]*(L0(op(j,op(i,b)),a,b)+
	    LJ(a[i],b,j)))*DeltaWs||j[n],j = 1 .. nops(op(1,b)))+1/2*
            add(add(LJ(op(j2,op(i,b)),b,j1)*(DeltaWs||j1[n]*DeltaWs||j2[n]+
	    V[j1,j2]),j1 = 1 .. nops(op(1,b))),
            j2 = 1 .. nops(op(1,b)));
        for u to nops(a) do  soln[i] := subs(x[u] = Y||u[n],soln[i]) od
    od;
    RETURN(eval(soln))
end:

 
stochastic[MLJ]:=
proc(X::algebraic,a::list(algebraic),b::list(list(algebraic)),j::integer)
local flag;
    flag := 0;
    if j = 0 then flag := L0(X,a,b) fi;
    if flag = 0 then flag := add(op(j,op(k,b))*diff(X,x[k]),
    k = 1 .. nops(b)) fi;
    RETURN(flag)
end:

 
stochastic[wktay3]:=proc(a::list(algebraic),b::list(list(algebraic)))
local u,i,soln;
    for i to nops(a) do
        soln[i] := Y||i[n+1] = Y||i[n]+a[i]*Delta[n]+
	add(op(j,op(i,b))*DeltaW||j[n],j = 1 .. nops(op(1,b)))+
         add(MLJ(a[i],a,b,j0)*I[j0,0],j0 = 0 .. nops(op(1,b)))+
         add(add(MLJ(op(j2,op(i,b)),a,b,j1)*I[j1,j2],
	 j2 = 1 .. nops(op(1,b))),j1 = 0 .. nops(op(1,b)))+
         add(add(MLJ(MLJ(a[i],a,b,k2),a,b,k1)*I[k1,k2,0],
	 k1 = 0 .. nops(op(1,b))),k2 = 0 .. nops(op(1,b)))+add(
         add(add(MLJ(MLJ(op(m3,op(i,b)),a,b,m2),a,b,m1)*I[m1,m2,m3],
	 m3 = 1 .. nops(op(1,b))),m2 = 0 .. nops(op(1,b))),
         m1 = 0 .. nops(op(1,b)));
        for u to nops(a) do  soln[i] := subs(x[u] = Y||u[n],soln[i]) od
    od;
    RETURN(eval(soln))
end:

 
stochastic[colour]:=
proc(a::list(algebraic),b::list(algebraic))
local temp1,i;
    for i to nops(a) do  temp1[i] := dx[i][t] = a[i]*dt+b[i]*z[t]*dt od;
    temp1[i] := dz[t] = -gamma*z[t]*dt+beta*dW[t];
    RETURN(eval(temp1))
end:


stochastic[comm1]:=proc()
local LJ1,LJ2,k,j1,j2,flag,p;
    for p to nargs do
        if type(args[p],list) <> true then 
	ERROR(`Expecting input to be an expression sequence of lists`) fi
    od;
    for k to nargs do
        for j1 to nops(args[1]) do
            for j2 to nops(args[1]) do
                LJ1 := add(op(j1,args[l])*diff(op(j2,args[k]),x[l]),
		l = 1 .. nargs);
                LJ2 := add(op(j2,args[l])*diff(op(j1,args[k]),x[l]),
		l = 1 .. nargs);
                if LJ1 <> LJ2 then flag := 1 fi
            od
        od
    od;
    if flag = 1 then 
    RETURN("Commutative noise of the first kind doesn't exist for this system")
    else RETURN("This system exhibits commutative noise of the first kind")
    fi;
end:


stochastic[comm2]:=
proc()
local LJ1LJ2,LJ2LJ1,k,p,j1,j2,j3,flag;
    for p to nargs do
        if type(args[p],list) <> true then
	ERROR(`Expecting input to be an expression sequence of lists`) fi;
    od;
    for k to nargs do
        for j1 to nops(args[1]) do
            for j2 to nops(args[1]) do
                for j3 to nops(args[1]) do
                    LJ1LJ2 := add(op(j1,args[m])*diff(add(op(j2,args[l])*
		    diff(op(j3,args[k]),x[l]),l = 1 .. nargs),x[m]),
                        m = 1 .. nargs);
                    LJ2LJ1 := add(op(j2,args[m])*diff(add(op(j1,args[l])*
		    diff(op(j3,args[k]),x[l]),l = 1 .. nargs),x[m]),
                        m = 1 .. nargs);
                    if LJ1LJ2 <> LJ2LJ1 then flag := 1 fi;
                od;
            od;
        od;
    od;
    if flag = 1 then 
    RETURN("Commutative noise of the second kind doesn't exist for this system")
    else RETURN("This system exhibits commutative noise of the second kind")
    fi;
end:

 
stochastic[itoformula]:=proc(U::list(algebraic),a::list(algebraic),
                 b::list(list(algebraic)))
  local i,k,l0,lj,soln;
  for i from 1 to nops(U) do
    l0:=L0(U[i],a,b)*dt;
    lj:=0;
    for k from 1 to nops(op(1,b)) do
      lj:=lj+LJ(U[i],b,k)*dW||k;
    od;
    soln[i]:=dX||i=l0 +lj;
  od;
  RETURN(eval(soln));
end:


stochastic[LFP]:=proc(P::algebraic,a::list(algebraic),b::list(list(algebraic)))
local part1,part2,part3;
  part1 := diff(P,t);
  part2 := add(diff(a[k]*P,y[k]),k = 1 .. nops(a));
  part3 := 1/2*add(
  add(add(diff(op(j,op(k,b))*op(j,op(l,b))*P,y[k],y[l]),
  j = 1 .. nops(op(1,b))),k = 1 .. nops(a)),l = 1 .. nops(a)
  );
  part1+part2-part3;
end:


stochastic[conv]:=proc(a::list(algebraic),b::list(list(algebraic)),c::algebraic)
local temp,i;
if c=ito then
      for i from 1 to nops(a) do temp[i]:=op(i,a)-1/2*add(add(op(k,op(j,b))
          *diff(op(k,op(i,b)),x[j]),k=1..nops(op(1,b))),j=1..nops(a));
     od;
elif c=strat then
       for i from 1 to nops(a) do
temp[i]:=op(i,a)+1/2*add(add(op(k,op(j,b))
           *diff(op(k,op(i,b)),x[j]),k=1..nops(op(1,b))),j=1..nops(a));
     od;
else
     ERROR(`Must enter either ito or strat for the 3rd argument`)
fi;
RETURN(map(simplify,eval(temp)))
end:


stochastic[chainrule]:=proc(U::list(algebraic),a::list(algebraic),
                            b::list(list(algebraic)))
  local i,k,l0,lj,soln;
  for i from 1 to nops(U) do
    l0:=SL0(U[i],a,b)*dt;
    lj:=0;
    for k from 1 to nops(b) do
      lj:=lj+LJ(U[i],b,k)*odW||i;
    od;
    soln[i]:=dX||i=l0 +lj;
  od;
  RETURN(eval(soln));
end:


stochastic[linearize]:=proc(a::list(algebraic),b::list(list(algebraic)),c::list(algebraic))
local i,tempA,tempB,j,k,l;
tempA:=array(1..nops(a),1..nops(a));
for i from 1 to nops(a) do for j from 1 to nops(a)
                        do tempA[i,j]:=diff(op(i,a),x[j]);
  od; od;
for i from 1 to nops(a) do for j from 1 to nops(a) do for l from 1 to nops(c)
                        do tempA[i,j]:=subs(x[l]=op(l,c),tempA[i,j]);
 od; od; od;
for k from 1 to nops(op(1,b)) do tempB[k]:=array(1..nops(a),1..nops(a));
for i from 1 to nops(a) do for j from 1 to nops(a)
                        do tempB[k][i,j]:=diff(op(k,op(i,b)),x[j]);
 od; od;
for i from 1 to nops(a) do for j from 1 to nops(a) do for l from 1 to nops(c)
                        do tempB[k][i,j]:=subs(x[l]=op(l,c),tempB[k][i,j]);
 od; od; od;

od;
RETURN(A=map(simplify,convert(eval(tempA),matrix)),B=eval(tempB))
end:


stochastic[sphere]:=proc(a,b)
  global q,q0,qk,h,hk;
  local i,j,k,tempa,tempb,stempbs,N,tmp;

  if type(a,array) then tmp:=convert(a,listlist);else tmp:=a; fi;
  N:=nops(tmp);
  hk:=evaln(hk); h:=evaln(h); qk:=evaln(qk);
  q:=evaln(q); q0:=evaln(q0);tempa:=evaln(tempa);tempb:=evaln(tempb);

  q0:=add(add(s[i]*a[i,j],i=1..N)*s[j],j=1..N);
  for k from 1 to nops(b) do
    qk[k]:=add(add(s[i]*b[k][i,j],i=1..N)*s[j],j=1..N);
  od;

  for k from 1 to nops(b) do

    stempbs[k]:=add(add(s[i]*(b[k][i,j]+b[k][j,i]),i=1..N)*s[j],j=1..N);
  od;

  q:=q0+ add(0.5*stempbs[k]-qk[k]^2,k=1..nops(b));

  for i from 1 to N do
    for j from 1 to N do
      if (i=j) then tempa[i,i]:=a[i,i]-q0;
               else tempa[i,j]:=a[i,j];
      fi;
    od;
  od;

  for i from 1 to N do
    h[i]:=add(tempa[i,j],j=1..N);
  od;

  for k from 1 to nops(b) do
    for i from 1 to N do
      for j from 1 to N do
        if (i=j) then tempb[k][i,i]:=b[k][i,i]-qk[k];
                 else tempb[k][i,j]:=b[k][i,j];
        fi;
      od;
    od;
  od;

  for k from 1 to nops(b) do
    for i from 1 to N do
      hk[k][i]:=add(tempb[k][i,j],j=1..N);
    od;
  od;

end:


stochastic[position]:=proc(N,i,j)
  global stelle;
  stelle:=add(N-k+1,k=1..i-1)+j-i+1;
end:


stochastic[ap]:=proc(A)
  local i,j,k,Atmp,N,counter;
  global B1;
  if type(A,array) then Atmp:=convert(A,listlist);else Atmp:=A; fi;
  N:=nops(Atmp);
  B1:=array(1..N*(N+1)/2,1..N*(N+1)/2);
  for i from 1 to N*(N+1)/2 do
    for j from 1 to N*(N+1)/2 do
      B1[i,j]:=0;
    od;
  od;
  counter:=0:
  for i from 1 to N do
    for j from i to N do
      counter:=counter+1;
      for k from 1 to N do
        if (j<=k)  then B1[counter,position(N,j,k)]:=A[i,k];
                   else B1[counter,position(N,k,j)]:=A[i,k];
        fi;
      od;
    od;
  od;
  RETURN(evalm(B1));
end:


stochastic[pa]:=proc(A)
  local i,j,k,Atmp,N,counter;
  global B2;
  if type(A,array) then Atmp:=convert(A,listlist);else Atmp:=A; fi;
  N:=nops(Atmp);
  B2:=array(1..N*(N+1)/2,1..N*(N+1)/2);
  for i from 1 to N*(N+1)/2 do
    for j from 1 to N*(N+1)/2 do
      B2[i,j]:=0;
    od;
  od;
  counter:=0:
  for i from 1 to N do
    for j from i to N do
      counter:=counter+1;
      for k from 1 to N do
        if (i<=k)  then B2[counter,position(N,i,k)]:=A[j,k];
                   else B2[counter,position(N,k,i)]:=A[j,k];
        fi;
      od;
    od;
  od;
  RETURN(evalm(B2));
end:


stochastic[bpb]:=proc(B)
  local i,j,k,l,Btmp,N,counter;
  global B3;
  if type(B,array) then Btmp:=convert(B,listlist);else Btmp:=B; fi;
  N:=nops(Btmp);
  B3:=array(1..N*(N+1)/2,1..N*(N+1)/2);
  for i from 1 to N*(N+1)/2 do
    for j from 1 to N*(N+1)/2 do
      B3[i,j]:=0;
    od;
  od;
  counter:=0:
  for i from 1 to N do
    for j from i to N do
      counter:=counter+1;
      for l from 1 to N do
        for k from 1 to N do
          if (k<=l)  then
            B3[counter,position(N,k,l)]:=B3[counter,position(N,k,l)]
                          +B[i,k]*B[j,l];
                     else
            B3[counter,position(N,l,k)]:=B3[counter,position(N,l,k)]
                          +B[i,k]*B[j,l];
          fi;
        od;
      od;
    od;
  od;
  RETURN(evalm(B3));
end:


stochastic[momenteqn]:=proc(A,B)
  local i,j,k,N,Btmp,Ctmp;
  global New_A;
  if type(A,array) then Btmp:=convert(A,listlist);else Btmp:=A; fi;
  N:=nops(Btmp);

  New_A:=array(1..N*(N+1)/2,1..N*(N+1)/2);
  Ctmp:=array(1..N*(N+1)/2,1..N*(N+1)/2);
  stochastic[ap](A);
  stochastic[pa](A);
  for i from 1 to N*(N+1)/2 do
    for j from 1 to N*(N+1)/2 do
      Ctmp[i,j]:=0;
    od;
  od;
  for k from 1 to nops(B) do
    stochastic[bpb](B[k]);
    for i from 1 to N*(N+1)/2 do
      for j from 1 to N*(N+1)/2 do
        Ctmp[i,j]:=Ctmp[i,j]+B3[i,j];
      od;
    od;
  od;
  for i from 1 to N*(N+1)/2 do
    for j from 1 to N*(N+1)/2 do
      New_A[i,j]:=B1[i,j]+B2[i,j]+Ctmp[i,j];
    od;
  od;
  RETURN(evalm(New_A));
end:


stochastic[pmatrix2pvector]:=proc(p)
  local i,j,k,ptmp;
  global pvector;
  if type(p,array) then ptmp:=convert(p,listlist);else ptmp:=p; fi;
  pvector:=array(1..nops(ptmp)*(nops(ptmp)+1)/2);
  k:=0;
  for i from 1 to nops(ptmp) do
    if (i>1) then k:=k+(nops(ptmp)-i+2); fi;
    for j from i to nops(ptmp) do
       pvector[k+j-i+1]:=ptmp[i,j];
    od;
  od;
  RETURN(eval(pvector));
end:


stochastic[pvector2pmatrix]:=proc(pvector)
  local i,j,k,ptmp,N;
  global p;
  if type(pvector,array) then ptmp:=convert(pvector,list);else
  ptmp:=pvector; fi;
  N:=-1/2+sqrt(1/4+2*nops(ptmp));
  p:=array(1..N,1..N);
  k:=0;
  for i from 1 to N do
    if (i>1) then k:=k+(N-i+2); fi;
    for j from i to N do
      p[i,j]:=ptmp[k+j-i+1];
      if (i<>j) then p[j,i]:=p[i,j]; fi;
    od;
  od;
  RETURN(eval(p));
end:


stochastic[milcomm]:=proc(a::list(algebraic),b::list(list(algebraic)))
  local u,i,l,soln;
  for i to nops(a) do 
    soln[i]:=Y||i[n+1]=Y||i[n]+L0(x[i],a,b)*Delta[n]
           +add(LJ(x[i],b,j)*DeltaW||j[n],j=1..nops(op(1,b)))
           +1/2*add(add(LJ(op(j2,op(i,b)),b,j1)*
                 (DeltaW||j1[n])*(DeltaW||j2[n]),
                 j1=1..nops(op(1,b))),j2=1..nops(op(1,b)));
    for l to nops(op(1,b)) do 	
      soln[i]:=subs((DeltaW||l[n])^2=((DeltaW||l[n])^2-Delta[n]),
                                       soln[i]) od;
    for u to nops(a) do
      soln[i]:=subs(x[u]=Y||u[n],soln[i]);
  od; od;
  RETURN(eval(soln));
end: