###
### Re-implementation of simple shooting in Maple 6 (DBM, 8 Sep 2000)
### 
### Based on implementation originally prepared for Maple V, Release 3,
### and updated to Release 4. The reimplementation is necessitated by
### changes to the limitations on using the results from dsolve,numeric
### in fsolve.

shoot := proc( )
  local ode, ic, bc, fn, init_param, v, endpts, bc_terms, bc_term1, bc_term2,
        optI, optR, optG,
        method,
        fn_name, fn_aux, n_fn, n_param, _ode_aux,
        aux_fn, aux_ode, aux_ic,
        full_fn, full_ode, full_ic,
        convert_aux_fn, convert_aux_f,
        F_new, F_old,
        param, param_name, _rename, L, new_var,
        i, j, l, n, r, R, FN, ODE, IC, F, DF, H, DH;
#  option `Copyright 1995-2000 by Douglas B. Meade, University of South Carolina`;

### Process arguments

  ode, ic, bc, fn, init_param, v, endpts, bc_terms
   := `shoot/process_required_args`(args[1..5]);
  optG, optI, optR
   := `shoot/process_optional_args`(args[6..-1]);

#### Extract commonly used information

  fn_name := map2(op,0,fn);
  param_name := v[2..-1];
  convert_aux_fn := [ seq( op([
                       f(op(v))=f(v[1]),
                       seq( op([
                        Diff(f(op(v)),p)
                         =cat(f,`_`,p)(v[1]),
                        Diff(Diff(f(op(v)),v[1]),p)
                         =Diff(cat(f,`_`,p)(v[1]),v[1]) ]),
                        p=param_name )]),
                       f=fn_name ) ];
  convert_aux_f := [ seq(
                        seq( op([
                         f(op(subs(v[1]=e,v))) = f(e),
                         seq(
                          diff(f(op(subs(v[1]=e,v))),p) = cat(f,`_`,p)(e),
                          p=param_name )]),
                         f=fn_name ),
                        e=endpts ) ];
  aux_fn := [ seq(
               seq(
                cat(f,`_`,p)(v[1]),
                p=param_name ),
               f=fn_name ) ];
  aux_ode := map(convert,ode,Diff);
  aux_ode := [ seq(
                seq(
                 convert(map(diff,eqn,p),Diff),
                 p=param_name),
                eqn=subs(convert_var_list(fn,[v[1]],v), ode) ) ];
  aux_ode := subs( convert_aux_fn, aux_ode );
  aux_ode := subs( convert_aux_fn, convert(aux_ode,diff) );
  aux_ic := [ seq(
               seq(
                cat(op(0,lhs(eqn)),`_`,p)(endpts[1]) = diff(rhs(eqn),p),
                p=param_name),
               eqn=ic ) ];

  full_fn  := [ op(fn),  op(aux_fn)  ];
  full_ode := [ op(ode), op(aux_ode) ];
  full_ic  := [ op(ic), op(aux_ic) ];

  method := subs( optG, `iter_method` );
  bc_term1,bc_term2 := selectremove(t->evalb(op(t)=endpts[1]),bc_terms);
  new_var := [op(convert_var_list(bc_term1,[endpts[1]],subsop(1=endpts[1],v))),
              op(convert_var_list(bc_term2,[endpts[2]],subsop(1=endpts[2],v)))];

  F := bc;
  DF := linalg[jacobian] (subs(new_var,F), param_name );
  DF := subs( convert_aux_f, evalm(DF) );

  if method=`newton` then
    param := `shoot/newton`(F,DF,init_param,
                            full_ode,full_ic,full_fn,
                            endpts,optI,optG);
  elif method=`modified` then
    H := linalg[innerprod](F,F);
    DH := linalg[grad] (subs(new_var,H), param_name );
    DH := subs( convert_aux_f, evalm(DH) );
    param := `shoot/modified_newton`(F,DF,DH,init_param,
                                     full_ode,full_ic,full_fn,
                                     endpts,optI,optG);
  elif method=`broyden` then
    param := `shoot/broyden`(F,DF,init_param,
                             full_ode,full_ic,full_fn,
                             endpts,optI,optG);
  elif method=`BFGS` then
    error "BFGS method not yet implemented"
  fi;

  if param=FAIL then 
    error sprintf("no convergence detected in %d iterations", max_iter)
  end if;

  return dsolve( { op(ode), op(subs(param,ic)) }, fn, op(optR) );
end proc:

convert_var_list := proc( L::list, curr_var::list, new_var::list )
  zip( (a,b)->a=b, L, map2(subs,op(curr_var)=op(new_var), L) );
end proc:

reverse := (e::equation) -> rhs(e)=lhs(e):
move_to_lhs := (e::equation) -> lhs(e)-rhs(e):

`shoot/G` := proc( F, ivp, fn )
  local i, SOL, SOLb;

  SOL := dsolve( ivp, fn, args[4..-1] );
  SOLb := [ op(subs( SOL[1,1][1]=SOL[2,1][1,1],
                     [seq( SOL[1,1][i]=SOL[2,1][1,i],
                           i=2..linalg[vectdim](SOL[1,1]) )]
                   )),
            op(subs( SOL[1,1][1]=SOL[2,1][2,1],
                     [seq( SOL[1,1][i]=SOL[2,1][2,i],
                           i=2..linalg[vectdim](SOL[1,1]) )]
                   ))];
  userinfo(3,`shoot`,`solution to IVP:`,map(evalm,SOLb));
  return subs( SOLb, map(evalm,F) )
end proc:

`shoot/process_required_args` := proc(ODE::{list,set,equation,algebraic},
                                      IC::{list,set,equation},
                                      BC::{list,set,equation},
                                      FN::{list,set,equation},
                                      INIT_PARAM::{list,set,equation})
  local ode, ic, bc, fn, param, _v, _z, z0, z1, bc_terms, endpts;

  if   type(ODE,list)      then ode:=args[1]
  elif type(ODE,set)       then ode:=convert(args[1],list)
  elif type(ODE,equation)  then ode:=[args[1]]
  elif type(ODE,algebraic) then ode:=[args[1]=0]
  else
    error "argument 1 must be the ODE(s) as a set, list, single equation, or single algebraic expression"
  end if;
  ode := map( isolate, ode, diff );
  userinfo(5,`shoot`,`system of ODEs          : `,ode);

  if   type(IC,list)       then ic:=args[2]
  elif type(IC,set)        then ic:=convert(args[2],list)
  elif type(IC,equation)   then ic:=[args[2]]
  else
    error "argument 2 must be the ICs as a set, list, or single equation"
  end if;
  userinfo(5,`shoot`,`initial conditions      : `,ic);

  if   type(BC,list)       then bc:=args[3]
  elif type(BC,set)        then bc:=convert(args[3],list)
  elif type(BC,equation)   then bc:=[args[3]]
  else
    error "argument 3 must be the BCs as a set, list, or single equation"
  end if;
  bc := map(move_to_lhs,bc);
  userinfo(5,`shoot`,`boundary conditions     : `,bc);

  if   type(FN,list)       then fn:=args[4]
  elif type(FN,set)        then fn:=convert(args[4],list)
  elif type(FN,function)   then fn:=[args[4]]
  else
    error "argument 4 must be the functions as a set, list, or single function"
  end if;
  userinfo(5,`shoot`,`function names          : `,fn);

  if   type(INIT_PARAM,list)     then param:=args[5]
  elif type(INIT_PARAM,set)      then param:=convert(args[5],list)
  elif type(INIT_PARAM,equation) then param:=[args[5]]
  else
    error "argument 5 must be the initial shooting parameters as a set, list, or single equation"
  end if;
  userinfo(5,`shoot`,`initial parameter values: `,param);

### Check for same number of functions, ODEs, and ICs

  if nops(map(nops,{fn,ode,ic}))<>1
    then error "invalid arguments, different number of ODEs, ICs, and functions"
  end if;

### Identify names and initial values of independent variable and parameter(s)

  _z := op(convert(map(op,fn),set)):
  if nops({_z})<>1 then
    error "invalid function list, expecting only one independent variable"
  end if;
  _v := [ _z, op(map(lhs,param)) ];
  userinfo(5,`shoot`,`auxiliary variables     : `, _v);

### Get first boundary point (for initial conditions)

  z0 := op(convert( map(op@lhs,ic), set ));
  if nops({z0}) <> 1
    then error "invalid initial conditions, must be given at one point"
  end if;

### Extract constraints (i.e., boundary conditions)

  bc_terms := convert( indets(bc,function(realcons)), list );
  z1 := op(map2(op,1,bc_terms));
  if nops({z1} minus {z0}) = 1 then
    z1 := op( {z1} minus {z0} );
  else
    error "invalid boundary conditions, must be given at one point"
  end if;
  endpts := [z0,z1];
  userinfo(5,`shoot`,`endpoints               : `, endpts);
  userinfo(5,`shoot`,`terms in BCs            : `, bc_terms);

  return ode, ic, bc, fn, param, _v, endpts, bc_terms;

end proc:

`shoot/process_optional_args` := proc()
  local i, l, opt, r, optG, optI, optR;

### - optG contains general parameters for the iterative solver
### - optI contains optional dsolve arguments for iterations
### - optR contains optional dsolve arguments for returned value

  optG := [
           abserr = Float(1,2-Digits),
           relerr = Float(1,2-Digits),
           max_iter = 10,
           iter_method = 'newton',        # modified, Broyden, BFGS
           broyden_init = 'identity'      # grad
          ];

  opt:=NULL;
  for i in args do
    if not(type(i,`=`)) then
      WARNING("optional argument, %1, ignored, must be an equation",i);
      next
    end if;
    l := lhs(i);
    r := rhs(i);
    if l='abserr' then
      optG := subsop(1=i,optG);
      next
    end if;
    if l='ic_err_tol' then
      optG := subsop(2=i,optG);
      next
    end if;
    if l='max_iter' then
      optG := subsop(3=i,optG);
      next
    end if;
    if l='iter_method' then
      if member(r,{'newton','modified','broyden','BFGS'}) then
        optG := subsop(4=i,optG);
      next
      else
        WARNING("unknown iterative solver, %1, using default (%2)",
                r, rhs(optG[4]));
      end if;
    end if;
    if l='broyden_init' then
      if member(r,{'identity','gradient'}) then
        optG := subsop(5=i,optG);
      next
      else
        WARNING("unknown iterative solver, %1, using default (%2)",
                r, rhs(optG[5]));
      end if;
    end if;
    if l='type' and r<>'numeric' then
      WARNING("optional argument, type=... ignored; type=numeric is required");
      next
    end if;
    opt := opt, i;
  end do;

  optI:=remove(has,[opt],{'type','value','output'});
  if member('output=procedurelist',optI) then
    WARNING("optional argument, output=procedurelist, ignored; array of values is returned");
    optI := remove(has,optI,'output');
  end if;
  optI := [
           type=numeric,
#           output=listprocedure,
           op(optI)
          ];

  optR:=remove(has,[opt],{'type'});
  optR := [
           type=numeric,
           op(optR)
          ];

  userinfo(5,`shoot`,`optI  : `,optI);
  userinfo(5,`shoot`,`optR  : `,optR);
  userinfo(5,`shoot`,`optG  : `,optG);

  return optG, optI, optR;

end proc:

`shoot/newton` := proc(F,DF,init_pt,ode,ic,fn,endpts,optI,optG)
  local s, ds, Fs, normFs, norms, normds, DFs,
        ivp, param, n, N, opt, abs_tol, rel_tol;

  s := linalg[vector](map(rhs,[op(init_pt)]));
  ds := linalg[vector](nops(init_pt),0);

### Iterate until sufficiently close, or too many iterations

  N := subs( optG, `max_iter` );
  abs_tol := subs( optG, `abserr` );
  rel_tol := subs( optG, `relerr` );
  for n from 1 to N do
    param := zip( (l,r)->l=r, map(lhs,init_pt), convert(s, list) );
    userinfo(1,`shoot`,`Step # `,n);
    userinfo(1,`shoot`,`Parameter values : `,op(param));

    ## numerical solution of IVP, evaluated at boundary point

    ivp := { op(ode), op(subs(param,ic)) };
    opt := value=convert(endpts,vector), op(optI);
    Fs,DFs := op(`shoot/G`( map(evalm,[F,DF]), ivp, fn, opt ));
    userinfo(3,`shoot`,`value of Fs :`,map(evalm,Fs));
    userinfo(3,`shoot`,`value of DFs:`,map(evalm,DFs));

    ## check BC for convergence

    normFs := linalg[norm](Fs,2);
    norms := linalg[norm](s,2);
    normds := linalg[norm](ds,2);
    userinfo(2,`shoot`,`Absolute BC Error: `,normFs);

    if normds<>0 then
      userinfo(2,`shoot`,`Relative BC Error: `,norms/normds)
    end if;
    if normFs<abs_tol or (normds<>0 and norms/normds<rel_tol) then
      return param;    ## return numerical solution
    end if;

### one more step of the Newton-Raphson method

    ds  := evalm( -linalg[inverse](DFs) &* Fs );
    s   := evalm( s + ds );
  end do;

  return FAIL;

end proc:

`shoot/modified_newton` := proc(F,DF,DH,init_pt,ode,ic,fn,endpts,optI,optG)
  local s, ds, Fs, normFs, norms, normds, DFs, DHs,
        ivp, param, n, N, opt, abs_tol, rel_tol,
        C, gamma, j, lambda, hk_min, h_k, hk_zero;

  s := linalg[vector](map(rhs,[op(init_pt)]));
  ds := linalg[vector](nops(init_pt),0);

### Iterate until sufficiently close, or too many iterations

  N := subs( optG, `max_iter` );
  abs_tol := subs( optG, `abserr` );
  rel_tol := subs( optG, `relerr` );
  for n from 1 to N do
    param := zip( (l,r)->l=r, map(lhs,init_pt), convert(s, list) );
    userinfo(1,`shoot`,`Step # `,n);
    userinfo(1,`shoot`,`Parameter values : `,op(param));

    ## numerical solution of IVP, evaluated at boundary point

    ivp := { op(ode), op(subs(param,ic)) };
    opt := value=convert(endpts,vector), op(optI);
    Fs,DFs,DHs := op(`shoot/G`( map(evalm,[F,DF,DH]), ivp, fn, opt ));
    userinfo(3,`shoot`,`value of Fs :`,map(evalm,Fs));
    userinfo(3,`shoot`,`value of DFs:`,map(evalm,DFs));
    userinfo(3,`shoot`,`value of DHs:`,map(evalm,DHs));

    ## check BC for convergence

    normFs := linalg[norm](Fs,2);
    norms := linalg[norm](s,2);
    normds := linalg[norm](ds,2);
    userinfo(2,`shoot`,`Absolute BC Error: `,normFs);

    if normds<>0 then
      userinfo(2,`shoot`,`Relative BC Error: `,norms/normds)
    end if;
    if normFs<abs_tol or (normds<>0 and norms/normds<rel_tol) then
      return param;    ## return numerical solution
    end if;

### auxiliary numerical function: h = innerprod( f, f )

    h_k := proc(tau::numeric,x::vector,d::vector)
      local f, S, X;

      userinfo(3,`h_k`,`calling arguments:`);
      userinfo(3,`h_k`,` tau =`,tau);
      userinfo(3,`h_k`,` x   =`,evalm(x));
      userinfo(3,`h_k`,` d   =`,evalm(d));
      X := zip( (l,r)->l=r, map(lhs,init_pt), convert(evalm(x-tau*d),list) );
      userinfo(3,`h_k`,` x-tau d =`,evalm(X));
      ivp := { op(ode), op(subs(X,ic)) };
      f := `shoot/G`( F, ivp, fn, opt );
      userinfo(3,`h_k`,`value of f:     `,map(evalm,f));

      return linalg[innerprod](f,f)

    end proc:
#    debug( h_k ):

### one more step of the modified Newton method

    ds  := evalm( linalg[inverse](DFs) &* Fs );
    userinfo(3,`shoot`,`d_k      = `,evalm(ds));
    gamma := 1/linalg[cond](DFs,2);
    userinfo(3,`shoot`,`gamma_k  = `,gamma);
    C := gamma*linalg[norm](ds,2)*linalg[norm](DHs,2);
    userinfo(3,`shoot`,`constant = `,C);
    hk_zero := h_k(0,s,ds);
    hk_min := hk_zero;
    lambda := 1;
    for j from 0 by 1 do
      hk||j := h_k(2^(-j),s,ds);
      userinfo(3,`shoot`,` j=`,j,`: is `, hk||j,`<=`,hk_zero-C*2^(-j),`?`);
      if hk||j <= hk_zero-C*2^(-j) then break end if;
      if hk||j < hk_min then
        hk_min := hk||j;
        lambda := 2^(-j)
      end if;
    end do;
    userinfo(3,`shoot`,` lambda_k = `,lambda);
    userinfo(3,`shoot`,` hk_min =   `,hk_min);
    s   := evalm( s - lambda * ds );
    userinfo(3,`shoot`,`x_(k+1) = `,evalm(s));
  end do;

  return FAIL;

end proc:

`shoot/broyden` := proc(F,DF,init_pt,ode,ic,fn,endpts,optI,optG)
  local s, d, normFs, norms, normds, DFs, DHs,
        ivp, param, n, N, opt, abs_tol, rel_tol,
        DD, DF0, F_curr, F_prev, y, init_D;

  s := linalg[vector](map(rhs,[op(init_pt)]));
  d := linalg[vector](nops(init_pt),0);

  param := zip( (l,r)->l=r, map(lhs,init_pt), convert(s, list) );
  userinfo(1,`shoot`,`Step # `,0);
  userinfo(1,`shoot`,`Parameter values : `,op(param));

  ## numerical solution of IVP, evaluated at boundary point

  ivp := { op(ode), op(subs(param,ic)) };
  opt := value=convert(endpts,vector), op(optI);
  F_curr,DF0 := op(`shoot/G`( map(evalm,[F,DF]), ivp, fn, opt ));
  userinfo(3,`shoot`,`value of F_curr:`,map(evalm,F_curr));

  init_D := subs( optG, `broyden_init` );
  if init_D=`identity` then
    DD := linalg[matrix](nops(init_pt),nops(init_pt),
                        (i,j)->if i=j then 1 else 0 end if);
  elif init_D=`gradient` then
    DD := evalm(DF0);
  end if;

  userinfo(3,`shoot`,`initial matrix, D0 = `,evalm(DD));

### Iterate until sufficiently close, or too many iterations

  N := subs( optG, `max_iter` );
  abs_tol := subs( optG, `abserr` );
  rel_tol := subs( optG, `relerr` );
  for n from 1 to N do
    param := zip( (l,r)->l=r, map(lhs,init_pt), convert(s, list) );
    userinfo(1,`shoot`,`Step # `,n);
    userinfo(1,`shoot`,`Parameter values : `,op(param));

    ## check BC for convergence

    normFs := linalg[norm](F_curr,2);
    norms := linalg[norm](s,2);
    normds := linalg[norm](d,2);
    userinfo(2,`shoot`,`Absolute BC Error: `,normFs);

    if normds<>0 then
      userinfo(2,`shoot`,`Relative BC Error: `,norms/normds)
    end if;
    if normFs<abs_tol or (normds<>0 and norms/normds<rel_tol) then
      return param;    ## return numerical solution
    end if;

### one more step of Broyden's method

    d := evalm( -linalg[inverse](DD) &* F_curr );
    s := evalm( s + d );

    F_prev := evalm( F_curr );
    param := zip( (l,r)->l=r, map(lhs,init_pt), convert(s, list) );
    ivp := { op(ode), op(subs(param,ic)) };
    opt := value=convert(endpts,vector), op(optI);
    F_curr := `shoot/G`( evalm(F), ivp, fn, opt );
    userinfo(3,`shoot`,`value of F_curr:`,map(evalm,F_curr));

    y := evalm( F_curr - F_prev );
    DD := evalm( DD + (F_curr &* linalg[transpose](d))/linalg[innerprod](d,d) );

  end do;

  return FAIL;

end proc: