SQLITE_NOTICE(283): recovered 5 frames from WAL file /data/mtt.fossil-wal

File mttroot/mtt/bin/trans/ode2obs_r artifact 117332d406 part of check-in 5d463df513

#! /bin/sh

     ###################################### 
     ##### Model Transformation Tools #####
     ######################################

# Bourne shell script: ode2obs_r
# Odrinary differential equations to observer function equations
# P.J.Gawthrop 14 June 1991, 8 Aug 1991, 2 April 1992, 14 April 1994, 28 Dec 94,
#               12th July 1995, April 1996
# Copyright (c) P.J.Gawthrop 1991, 1992, 1994, 1995, 1996.

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.5  1998/04/07 08:12:12  peterg
## Added affine form.
##
## Revision 1.4  1998/04/07 05:45:12  peterg
## Reverted to an older version due to untraceable bug in new version
##
## Revision 1.3  1996/08/25 10:07:05  peter
## Remove a du state ment causaing touble
##      - but needs more work.
##
## Revision 1.2  1996/08/25 08:38:14  peter
## Error handling added.
##
## Revision 1.1  1996/08/25 08:37:44  peter
## Initial revision
##
###############################################################

#Inform user
echo Creating $1_obs.r

# Remove the old log file
rm -f ode2obs_r.log

# Use reduce to accomplish the transformation
reduce >ode2obs_r.log << EOF

%Read the formatting function
in "$MTTPATH/trans/reduce_matrix.r";


OFF Echo;
OFF Nat;
ON NERO;

in "$1_def.r";
MTTdxs := MTTdX;  %Save the symbolic form of dX

%Set default values - reset by obspar file.
MTTGPCNy := 2;
MTTGPCNu := 0;

%%in "$1_sympar.r";
in "$1_ode.r";
%%in "$1_simp.r";

% Read the obs form parameters.
in "$1_obspar.r";

%Create the U vector of input derivatives.
MTTNuu := (MTTGPCNu+1)*MTTNu;
MTTNuu1 := MTTGPCNu*MTTNu;


%FOR i := 1:MTTNu DO
%  MTTUU(i,1) := MTTU(i,1);
%END;
%
%IF MTTGPCNu>0 THEN
%BEGIN
%  FOR i := 1:MTTNu DO
%    MTTUU(i+MTTNu,1) := MTTdU(i,1);
%  END;
%END;

MTTU := MTTU;
MTTdU := MTTdU; % ---- removed temporarily, needs def file change? 

%Create the Y vector of output derivatives.
MTTNyy := (MTTGPCNy+1)*MTTNy;
Matrix MTTYY(MTTNyy,1);

FOR i := 1:MTTNy DO
  MTTYY(i,1) := MTTY(i,1);
END;

l := MTTNy;
FOR i := 1:MTTGPCNy DO
    FOR j := 1:MTTNy DO
    BEGIN
      l := l+1;
      MTTYY(l,1) := 0;
      FOR k := 1:MTTNx DO %Derivatives wrt x
      BEGIN
        xk := MTTX(k,1);
        MTTYY(l,1) := MTTYY(l,1) + df(MTTYY(l-MTTNy,1), xk, 1)*MTTdX(k,1);
      END;
      IF MTTGPCNu>0 THEN
      FOR k := 1:MTTGPCNu DO %Non-zero derivatives of u
      BEGIN
        FOR kk := 1:MTTNu DO
        BEGIN
          uk := MTTUU(kk,k);
          MTTYY(l,1) := MTTYY(l,1) + df(MTTYY(l,1), uk, 1)*MTTUU(1,k+1);
        END;
      END;
    END;
END;

%%Create O_x - derivative of YY wrt x
%Matrix MTTO_x(MTTNyy,MTTNx);
%FOR j := 1:MTTNx DO
%  BEGIN
%  xj := MTTX(j,1);
%  FOR i := 1:MTTNyy DO
%    BEGIN
%    MTTO_x(i,j) := df(MTTYY(i,1), xj);
%    END;
%  END;

%%Create O_u - derivative of YY wrt u (Assumes GPC Nu = 0)
%MTTNNu := (MTTGPCNu+1)*MTTNu;
%Matrix MTTO_u(MTTNyy,MTTNNu);
%FOR j := 1:MTTNNu DO
%  BEGIN
%  uj := MTTu(j,1);     
%  FOR i := 1:MTTNyy DO
%    BEGIN
%    MTTO_u(i,j) := df(MTTYY(i,1), uj);
%    END;
%  END;


%%Create O_uu - derivative of O_u wrt u (Assumes GPC Nu = 0)
%%This is a multi-dimensional matrix kth elements stacked sideways.
%Matrix MTTO_uu(MTTNyy,MTTNNu*MTTNNu);
%FOR k := 1:MTTNNu DO
%  BEGIN
%  uk := MTTu(k,1);     
%  FOR j := 1:MTTNNu DO
%    BEGIN
%    FOR i := 1:MTTNyy DO
%      BEGIN
%      jk := j+(k-1)*MTTNu;
%      MTTO_uu(i,jk) := df(MTTO_u(i,j), uk);
%      END;
%    END;
%  END;

%%Create O_ux - derivative of O_u wrt x 
%%This is a multi-dimensional matrix kth elements stacked sideways.
%Matrix MTTO_ux(MTTNyy,MTTNu*MTTNx);
%FOR k := 1:MTTNx DO
%  BEGIN
%  xk := MTTx(k,1);     
%  FOR j := 1:MTTNu DO
%    BEGIN
%    FOR i := 1:MTTNyy DO
%      BEGIN
%      jk := j+(k-1)*MTTNu;
%      MTTO_ux(i,jk) := df(MTTO_u(i,j), xk);
%      END;
%    END;
%  END;


% Try and split the obs function into affine form (O(x,u) = o(x) + h(x)u)
affine := 1;
MTTNuu := (MTTGPCNu+1)*MTTNu;
write "1 affine := ", affine, ";";

%IF (affine=1) THEN
%BEGIN
  MATRIX MTTObs_o(MTTNyy,1);
  MATRIX MTTObs_h(MTTNyy,MTTNuu);
write "2 affine := ", affine, ";";
  FOR i := 1:MTTNyy DO
  BEGIN
write "3 affine := ", affine, ";";
    MTTObs_o(i,1) := MTTYY(i,1);
    k := 0;
    FOR j := 1:MTTNu DO
    FOR jj := 0:MTTGPCNu DO
    BEGIN
      k := k+1;
      %Expand as polynomial in u_j^[jj]
      coeffs := coeff(MTTObs_o(i,1), MTTUU(j,jj+1));
      MTTObs_o(i,1) := first(coeffs);
      IF length(coeffs)>1 THEN
        MTTObs_h(i,k) := second(coeffs);
      IF length(coeffs)>2 THEN 
        affine := -1;
    END;
  END;
%END;


%%Create the _obs.r file
OUT "$1_obs.r";

IF affine=1 THEN
  write "affine := 1; % The O function is affine in u"
ELSE
  write "affine := 0; % The O function is not affine in u";


write "% The matrix sizes";
write "MTTNyy := ", MTTNyy, ";";
write "MTTNuu := ", MTTNuu, ";";

mtt_matrix := MTTYY$
mtt_matrix_n := MTTNyy$
mtt_matrix_m := 1$
mtt_matrix_name := "MTTYY"$
reduce_matrix();

IF affine=1 THEN
BEGIN
  mtt_matrix := MTTObs_o$
  mtt_matrix_n := MTTNyy$
  mtt_matrix_m := 1$
  mtt_matrix_name := "MTTObs_o"$
  reduce_matrix();

  mtt_matrix := MTTObs_h$
  mtt_matrix_n := MTTNyy$
  mtt_matrix_m := MTTNuu$
  mtt_matrix_name := "MTTObs_h"$
  reduce_matrix();
END;

%MTTO_x := MTTO_x;
%MTTO_u := MTTO_u;
%MTTO_uu := MTTO_uu;
%MTTO_ux := MTTO_ux;

write ";END;";

SHUT "$1_obs.r";
quit;

EOF

# Now invoke the standard error handling.
mtt_error_r ode2obs_r.log