#! /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` ## Number of parameters
## Parameter information for reduce.
parameters=`echo $2 | sed 's/,/ /g' |\
gawk '{
for (i=1; i<=NF; i++) {
printf("mttpar(%i,1) := %s;\n", i, $i);
printf("mttcoef(%i,1) := %ss;\n", i, $i);
}
}'`
## Update sympar list.
echo Recreating $1_sympar.txt
# Zap any sensitivity coeficients
mv $1_sympar.txt mtt_junk
grep -v MTT_Sensitivity_Coefficients mtt_junk > $1_sympar.txt
# Create the new sens coeffs at end of list.
echo $2 | sed 's/,/ /g' |\
gawk '{
for (i=1; i<=NF; i++) {
printf("%ss\tMTT_Sensitivity_Coefficients\n",$i);
}
}' >> $1_sympar.txt
touch $1_sympar.txt
matrix="matrix mttpar("$n",1); matrix mttcoef("$n",1);"
# Number of states
Nx=`grep "MTTNx " <$system_def | gawk '{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";
% Constants
write "% New constants";
write "MTTNx := ", MTTNx2, ";";
write "MTTNy := ", MTTNy2, ";";
% E matrix
MTT_Matrix := MTTsE$
MTT_Matrix_name := "MTTE"$
MTT_Matrix_n := MTTNx2$
MTT_Matrix_m := MTTNx2$
Reduce_Matrix()$
% State derivative
MTT_Matrix := MTTsEdX$
MTT_Matrix_name := "MTTEdX"$
MTT_Matrix_n := MTTNx2$
MTT_Matrix_m := 1$
Reduce_Matrix()$
% Output
MTT_Matrix := MTTsY$
MTT_Matrix_name := "MTTY"$
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"
gawk '{
## 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