#! /bin/sh

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

# Bourne shell script: cse2scse_r
# Reduce constrained-state equations to sensitivity version
# P.J.Gawthrop 10th May 199, 8th August 1991, April 1994, Dec 1994
# Copyright (c) P.J.Gawthrop, 1999

# 	$Id$	

# Arguments
system=$1; 
system_def=$1_def.r
system_cse=$1_cse.r
system_scse=$1_scse.r

# Parameters
n=`echo $2 | sed 's/,/ /g' |wc -w`
echo $n_parameters

parameters=`echo $2 | sed 's/,/ /g' |\
 awk '{
   for (i=1; i<=NF; i++) {
      printf("mttpar(%i,1) := %s;\n", i, $i);
      printf("mttcoef(%i,1) := %ss;\n", i, $i);
   }
  }'`

matrix="matrix mttpar("$n",1); matrix mttcoef("$n",1);"

echo $parameters
echo $matrix

# Number of states
Nx=`grep "MTTNx " <$system_def | awk '{print $3}' | sed 's/;//'`

#Inform user
echo Creating $system_scse "(for parameters $2, $Nx states)"

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

# Use reduce to accomplish the transformation
$SYMBOLIC  << EOF >cse2scse_r.log

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

%Read the definitions file
in "$system_def";

%Read the constrained-state equations file
in "$system_cse";

% Declare the parameter matrix and fill it
$matrix
$parameters
mtt_n_par := $n;

MTTNx2 := 2*MTTNx;
MTTNy2 := 2*MTTNy;

% Compute the sensitivity E matrix
matrix MTTssE(MTTNx,MTTNx);
clear MTTx; % Dont use - mkid is better

FOR ii := 1:MTTNx DO
  FOR jj := 1:MTTNx DO
  IF MTTE(ii,jj) NEQ 0 THEN IF MTTE(ii,jj) NEQ 1 THEN
  BEGIN
    % First with respect to theta...
    FOR kk := 1:mtt_n_par DO
    BEGIN
      mttpar_k := mttpar(kk,1);
      mttcoef_k := mttcoef(kk,1);
      MTTssE(ii,jj) := MTTssE(ii,jj) + df(MTTE(ii,jj), mttpar_k)*mttcoef_k;
    END;
    % then with respect to x
    FOR i := 1:MTTNx DO
    BEGIN
      xi := mkid(MTTx,i);
      sxi := mkid(MTTsx,i);
      MTTssE(ii,jj) := MTTssE(ii,jj) + df(MTTE(ii,jj), xi)*sxi;
    END;
  END;

% Compute the sensitivity of the RHS of the cse
matrix MTTssEdx(MTTNx,1);


FOR ii := 1:MTTNx DO
BEGIN 
  % First with respect to theta
    FOR kk := 1:mtt_n_par DO
    BEGIN
      mttpar_k := mttpar(kk,1);
      mttcoef_k := mttcoef(kk,1);
      MTTssEdx(ii,1) := MTTssEdx(ii,1) + df(MTTEdx(ii,1), mttpar_k)*mttcoef_k;
    END;
% Then with respect to x
  FOR i := 1:MTTNx DO
  BEGIN
    xi := mkid(MTTx,i);
    sxi := mkid(MTTsx,i);
    MTTssEdx(ii,1) := MTTssEdx(ii,1) + df(MTTEdx(ii,1), xi)*sxi;
  END; 
END; 

% Sensitivity output function
matrix MTTssY(MTTNy,1);
FOR ii := 1:MTTNy DO
BEGIN 
  % First with respect to theta
    FOR kk := 1:mtt_n_par DO
    BEGIN
      mttpar_k := mttpar(kk,1);
      mttcoef_k := mttcoef(kk,1);
      MTTssY(ii,1) := MTTssY(ii,1) + df(MTTY(ii,1), mttpar_k)*mttcoef_k;
    END;
% Then with respect to x
  FOR i := 1:MTTNx DO
  BEGIN
    xi := mkid(MTTx,i);
    sxi := mkid(MTTsx,i);
    MTTssY(ii,1) := MTTssY(ii,1) + df(MTTY(ii,1), xi)*sxi;
  END; 
END; 


% Now reorganise everything into composite system
%  - odd rows are the system 
%  - even rows are the sensitivity system
% NB at this stage, the states are numbered incorrectly - sorted out below.

%E matrix 

matrix MTTsE(MTTNx2,MTTNx2);
FOR i := 1:MTTNx DO
  FOR j := 1:MTTNx DO
  BEGIN
    MTTsE(2*i-1,2*j-1) := MTTE(i,j);  % System
    MTTsE(2*i,2*j)     := MTTE(i,j);  % Sensitivity system
    MTTsE(2*i,2*j-1)   := MTTssE(i,j);  % Sensitivity system
  END; 

%dX matrix
matrix MTTsEdX(MTTNx2,1);

FOR i := 1:MTTNx DO
BEGIN
  MTTsEdX(2*i-1,1) := MTTEdx(i,1);
  MTTsEdX(2*i,1)   := MTTssEdx(i,1);
END; 

%Y matrix
matrix MTTsY(MTTNy2,1);%dX matrix
matrix MTTsEdX(MTTNx2,1);

FOR i := 1:MTTNx DO
BEGIN
  MTTsEdX(2*i-1,1) := MTTEdx(i,1);
  MTTsEdX(2*i,1)   := MTTssEdx(i,1);
END; 


FOR i := 1:MTTNy DO
BEGIN
  MTTsY(2*i-1,1) := MTTY(i,1);
  MTTsY(2*i,1)   := MTTssY(i,1);
END; 




OFF nat;

OUT "$system_scse";

write "%File: $system_scse";


% E matrix
MTT_Matrix := MTTsE$ 
MTT_Matrix_name := "MTTsE"$
MTT_Matrix_n := MTTNx2$
MTT_Matrix_m := MTTNx2$
Reduce_Matrix()$

% State derivative
MTT_Matrix := MTTsEdX$ 
MTT_Matrix_name := "MTTsEdX"$
MTT_Matrix_n := MTTNx2$
MTT_Matrix_m := 1$
Reduce_Matrix()$

% Output
MTT_Matrix := MTTsY$ 
MTT_Matrix_name := "MTTsY"$
MTT_Matrix_n := MTTNy2$
MTT_Matrix_m := 1$
Reduce_Matrix()$

write "END;";

SHUT "$system_scse";
quit;
EOF

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


## Now reorganise the  states
mv -f $system_scse mtt_junk


##echo "Nx = $Nx"
awk '{
  ## Make sure all MTTn variables are followed by a space
  gsub(/mttx[0-9]*/, "& ");

  for (i=Nx;i>0;i--) {

    state = sprintf("mttx%i ",i);
    newstate = sprintf("mttx%i ",2*i-1);
    sstate = sprintf("mttsx%i",i);
    newsstate = sprintf("mttx%i",2*i);

    gsub(state,newstate);
    gsub(sstate,newsstate);

  }
  print $0
}' Nx=$Nx <mtt_junk > $system_scse