# Copyright (C) 2001 by Peter J. Gawthrop

###############################################################
## Version control history
###############################################################
## $Header$
## $Log$
## Revision 1.309.2.1  2001/05/04 04:07:24  geraint
## Numerical solution of algebraic equations.
## sys_ae.cc written for unsolved inputs.
## Solution of equations using hybrd from MINPACK (as used by Octave fsolve).
##
## Revision 1.309  2001/04/28 03:15:03  geraint
## Fixed comment (interfered with "mtt help representations").
##
## Revision 1.308  2001/04/25 22:17:45  geraint
## Fixed icd.txt2 dependency.
##
## Revision 1.307  2001/04/15 21:15:41  geraint
## Added interface definition rep: _ICD.(txt|c|cc|m).
##
## Revision 1.305  2001/04/11 09:44:26  gawthrop

# Default not verbose
verbose=' -s'

# Default integration method
integration_method=implicit;

# Default algebraic equation solver
algebraic_solver=Reduce_Solver;

# Default no info
info_switch=''

# Default use m, not oct files
m='m';

# Default use ps files

			integration_method=rk4;
			mtt_switches="$mtt_switches rk4";
			;;
		    *)
			echo $1 is an unknown integration method - use euler, rk4 or implicit;
                        exit;;
		esac;;
	-ae )
                mtt_switches="$mtt_switches $1";
		case $2 in
		    fsolve | hybrd)
			mtt_switches="$mtt_switches $2";
			algebraic_solver=Hybrd_Solver;
			shift;
			;;
		    hooke)
			mtt_switches="$mtt_switches $2";
			algebraic_solver=HJ_Solver;
			shift;
			;;
		    reduce)
			mtt_switches="$mtt_switches $2";
			algebraic_solver=Reduce_Solver;			
			Solve='-A';
			shift;
			;;
		    *)
			echo $1 is an unknown solver - use hybrd, hooke or reduce;
			exit;;
	        esac;;
	-s )
                mtt_switches="$mtt_switches $1";
		sensitivity=sensitivity ;;
	-ss )
                mtt_switches="$mtt_switches $1";
		steadystate_computation=yes ;;
	-d )

    echo '         -I  prints more information'
    echo '         -abg start at abg.m representation'
    echo '         -c  c-code generation'
    echo '         -d  <dir>  use directory <dir>'
    echo '         -dc Maximise derivative (not integral) causality'
    echo '         -dc Maximise derivative (not integral) causality'
    echo '         -i <implicit|euler|rk4>   Use implicit, euler or rk4 integration'
    echo '         -ae <reduce|hybrd|hooke>   Solve algebraic equations with Reduce, hybrd (fsolve) or "Hooke and Jeeves"'
    echo '         -o ode is same as dae'
    echo '         -oct use oct files in place of m files where appropriate'
    echo '         -opt optimise code generation'
    echo '         -p  print environment variables'
    echo '         -partition partition hierachical system'
    echo '         -pdf generate pdf in place of ps'
    echo '         -r  reset time stamp on representation'

  rm -f *_csex.cc *_numpar.cc *_smxa.cc *_state.cc 
  rm -f *_ode.cc *_odeo.cc
  rm -f *_logic.m *_logic.cc *_logic.oct
  rm -f *_state.m *_state.cc *_state.oct
  rm -f *_ode2odes.* *.dat2 MTT.core
  rm -f *_modpar.txt *_modpar.r
  rm -f *_ICD.txt *_ICD.c *_ICD.cc *_ICD.m
  rm -f *_ae.r *_ae.m *_ae.cc *_ae.oct
  rm -fR *_rep MTT_work
  exit
fi

# Clean up named system
if [ "$2" = "Clean" ] && [ "$3" = "" ]; then
  echo 'Removing all generated files for system ' $1

  rm -f $1_csex.cc $1_numpar.cc $1_smxa.cc $1_state.cc 
  rm -f $1_ode.cc $1_odeo.cc
  rm -f $1_logic.m $1_logic.cc $1_logic.oct 
  rm -f $1_state.m $1_state.cc $1_state.oct
  rm -f $1_ode2odes.* $1.dat2
  rm -f $1_modpar.txt $1_modpar.r
  rm -f $1_ICD.txt $1_ICD.c $1_ICD.cc $1_ICD.m
  rm -f $1_ae.r $1_ae.m $1_ae.cc $1_ae.oct
  rm -fR $1_rep MTT_work
  exit
fi

if [ "$2" = "rep" ]; then
     documenttype=article
      # See if we are making a book -- ie representation rep on a directory

.PRECIOUS: %.cc # Don't let mtt delete them
$1_%.cc:  $1_%.m
	mtt_m2cc.sh  $1 \$* cat 

mtt_%.cc: ${MTT_LIB}/cc/mtt_%.cc
	cp ${MTT_LIB}/cc/mtt_\$*.cc mtt_\$*.cc

mtt_%.hh: ${MTT_LIB}/cc/mtt_%.hh
	cp ${MTT_LIB}/cc/mtt_\$*.hh mtt_\$*.hh

## .o files
.PRECIOUS: $1_%.o
$1_%.o: $1_%.cc $1_def.h $1_sympar.h $1_cr.h
	echo Compiling $1_\$*.cc
	${MTT_CXX} ${MTT_CXXFLAGS} ${MTT_CXXINCS} -c $< -DSTANDALONE

.PRECIOUS: mtt_%.o

else
	mtt_update $1 logic
endif

# Dummy target
FORCE:

$1_ode2odes_common.m    : $1_ae.m   $1_input.m   $1_logic.m   $1_numpar.m   $1_simpar.m   $1_state.m
	@echo > /dev/null
$1_ode2odes_common.cc   : $1_ae.cc  $1_input.cc  $1_logic.cc  $1_numpar.cc  $1_simpar.cc  $1_state.cc
	@echo > /dev/null
$1_ode2odes_common.o    : $1_ae.o   $1_input.o   $1_logic.o   $1_numpar.o   $1_simpar.o   $1_state.o
	@echo "Creating $1_ode2odes_common.o"
	ar -cr \$@ \$^
$1_ode2odes_common.oct  : $1_ae.oct $1_input.oct $1_logic.oct $1_numpar.oct $1_simpar.oct $1_state.oct
	@echo > /dev/null

$1_ode2odes_euler.m  $1_ode2odes_rk4.m    : $1_ode.m     $1_odeo.m
	@echo > /dev/null
$1_ode2odes_euler.cc $1_ode2odes_rk4.cc   : $1_ode.cc    $1_odeo.cc
	@echo > /dev/null
$1_ode2odes_euler.o  $1_ode2odes_rk4.o    : $1_ode.o     $1_odeo.o   mtt_euler.o

	@echo > /dev/null
$1_ode2odes_implicit.o  : $1_cseo.o    $1_csex.o   $1_smxa.o   $1_smxax.o   mtt_implicit.o
	@echo "Creating $1_ode2odes_implicit.o"
	ar -cr \$@ \$^
$1_ode2odes_implicit.oct: $1_cseo.oct  $1_csex.oct $1_smxa.oct $1_smxax.oct mtt_implicit.oct
	@echo > /dev/null

mtt_Solver.cc:	mtt_Solver.hh
$1_ode2odes_${algebraic_solver}.cc:	mtt_Solver.cc mtt_${algebraic_solver}.hh mtt_${algebraic_solver}.cc
$1_ode2odes_${algebraic_solver}.o:	mtt_Solver.o mtt_${algebraic_solver}.o
	@echo "Creating $1_ode2odes_${algebraic_solver}.o"
	ar -cr \$@ \$^

#SUMMARY numpar	numerical parameter declaration (m) 
$1_numpar.m:  $1_numpar.txt $1_sympars.txt
	mtt_txt2m $1 numpar

#SUMMARY numpar	numerical parameter declaration (c) 
#SUMMARY numpar	numerical parameter declaration (view) 

$1_dae.m: $1_def.r $1_dae.r  $1_sympars.txt
	dae_r2m $1; matlab_tidy $1_dae.m; matlab_tidy $1_daeo.m
$1_dae.c: $1_def.r $1_dae.r  $1_sympar.r
	dae_r2c $1; c_tidy $1_dae.c
$1_dae.tex: $1_dae.r $1_simp.r
	dae_r2tex $partition $1; latex_tidy $1_dae.tex

#SUMMARY ae	algebraic equations - unknown inputs (r)
#SUMMARY ae	algebraic equations - unknown inputs (m)
#SUMMARY ae	algebraic equations - unknown inputs (cc)

$1_ae.r: $1_cse.r
	touch $1_ae.r
$1_ae.m: $1_ae.r
	mtt_r2m $1 ae

#SUMMARY cse	constrained-state equations (r)
#SUMMARY cse*	constrained-state equations (m)
#SUMMARY cse*	constrained-state equations (oct)
#SUMMARY cse	constrained-state equations (tex)
#SUMMARY cse*	constrained-state equations (view)
#SUMMARY cse	constrained-state equations (ps)

                $1_smxa.m $1_smxax.m\
                $1_simpar.m $1_numpar.m $1_state.m $1_input.m \
                $1_csex.m $1_cseo.m  $1_logic.m
ifeq ($using_oct,yes)
	touch $1_ode2odes.m # Create a dummy which wont' be used
	mtt $mtt_switches -q -u $1 ode2odes oct
else
	make_ode2odes $1 m $integration_method $algebraic_solver
endif
endif
ifneq ($integration_method,implicit)
$1_ode2odes.m : $1_def.r $1_sympars.txt\
		$1_simpar.m $1_numpar.m $1_state.m $1_input.m \
		$1_ode.m $1_odeo.m  $1_logic.m
ifeq ($using_oct,yes)
	echo "*** Warning: Shouldn't be here! Creating dummy $1_ode2odes.m"
	touch $1_ode2odes.m # Create a dummy which wont' be used
	mtt $mtt_switches -q -u $1 ode2odes oct
else
	make_ode2odes $1 m $integration_method $algebraic_solver
endif

endif

#SUMMARY ode2odes Simulation function (m)
#SUMMARY ode2odes Simulation function (cc)
#SUMMARY ode2odes Simulation function (oct)
#SUMMARY ode2odes Simulation function (exe)
$1_ode2odes.exe: $1_def.h $1_sympar.h\
		 $1_ode2odes.o $1_ode2odes_common.o $1_ode2odes_${integration_method}.o $1_ode2odes_${algebraic_solver}.o
	echo Creating $1_ode2odes.exe
	${MTT_CXX} ${MTT_CXXFLAGS} -o $1_ode2odes.exe\
	$1_ode2odes.o $1_ode2odes_common.o $1_ode2odes_${integration_method}.o $1_ode2odes_${algebraic_solver}.o\
	${MTT_LDFLAGS} ${MTT_CXXLIBS}

$1_ode2odes.o: $1_ode2odes.cc $1_ode2odes_common.o $1_ode2odes_${integration_method}.o $1_ode2odes_${algebraic_solver}.o
	echo Creating $1_ode2odes.o
	${MTT_CXX} ${MTT_CXXFLAGS} ${MTT_CXXINCS} -c $1_ode2odes.cc -DSTANDALONE

$1_ode2odes.oct: $1_ode2odes.cc $1_ode2odes_common.oct $1_ode2odes_${integration_method}.oct
	touch $1_ode2odes.m
	echo Creating $1_ode2odes.oct
	$MKOCTFILE $1_ode2odes.cc

$1_ode2odes.cc: $1_def.r $1_sympars.txt\
		$1_ode2odes_common.m $1_ode2odes_common.cc\
		$1_ode2odes_${integration_method}.m $1_ode2odes_${integration_method}.cc mtt_Solver.cc mtt_${algebraic_solver}.cc mtt_${algebraic_solver}.hh
	touch $1_ode2odes.m
	make_ode2odes $1 cc $integration_method $algebraic_solver

#Conversion of m to p to c
#SUMMARY ode	ordinary differential equations (c)
#SUMMARY ode	ordinary differential equations (p)
#SUMMARY state	state declaration (c) 
#SUMMARY state	state declaration (p) 
$1_simpar.p : $1_def.r $1_simpar.m

# Copyright (c) P.J.Gawthrop, 1991, 1994.

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.1  2000/12/28 12:21:31  peterg
## Put under RCS
##
## Revision 1.2  1997/01/06 21:17:10  peterg
## Removed:  OFF Exp; OFF GCD;
##
## Revision 1.1  1996/08/25 10:05:45  peter
## Initial revision
##
###############################################################

in "$1_def.r";

%Read the substitution file
in "$1_subs.r";

%Read the constrained-state equations file
in "$1_cse.r";

%Read the algebraic equations file
in "$1_ae.r";

IF MTTNx>0 THEN
 IF MTTNz>0 THEN
   MTTdXX := MTTE^(-1)*MTTEdX
 ELSE
   MTTdXX := MTTEdX;

# Copyright (c) P.J.Gawthrop 1991, 1992, 1994.

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.15.2.2  2001/05/05 20:50:16  geraint
## Fixed errors when MTTNx=0.
##
## Revision 1.15.2.1  2001/05/04 04:07:24  geraint
## Numerical solution of algebraic equations.
## sys_ae.cc written for unsolved inputs.
## Solution of equations using hybrd from MINPACK (as used by Octave fsolve).
##
## Revision 1.15  2001/03/19 02:28:52  geraint
## Branch merge: merging-ode2odes-exe back to MAIN.
##
## Revision 1.14.2.1  2001/03/19 00:29:08  geraint
## Parse switches (-A) before calling def2write_r.
## Update $1_def.* instead of removing.
##
## Revision 1.14  2000/12/28 12:24:35  peterg
## *** empty log message ***
##

  BEGIN
    MTT_sol_i := first(MTT_sol); MTT_sol := rest(MTT_sol);
    set(lhs(MTT_sol_i),rhs(MTT_sol_i));
  END;

  % No algebraic variables left!
  MTTNYz := 0;
END; % IF MTTNyz>0 and $solve

OUT "$1_ae.r";
IF (MTTNyz>0) THEN % not $solve (or perhaps solution failed?)
BEGIN
    WRITE "MATRIX MTTyz(",MTTNyz,",1)";
    WRITE "%File: $1_ae.r";
    FOR i := 1:MTTNyz DO
	WRITE "MTTyz(",i,",1) := ",MTTyz(i,1);
END; % if MTTNyz>0 (and !$solve)
WRITE ";END;";
SHUT "$1_ae.r";

% Create the matrix declarations
OUT "$1_cse.r1";
write "%";
IF (MTTNx > 0) THEN
BEGIN
    write "MATRIX MTTEdx(", MTTNx, ",", 1, ")$"; 
    write "MATRIX MTTE(", MTTNx, ",", MTTNx, ")$"; 
END;
SHUT "$1_cse.r1";

OUT "$1_csex.r1";
write "%File:$1_csex.r1";
IF (MTTNx > 0) THEN
    write "MATRIX MTTEdx(", MTTNx, ",", 1, ")$"; 
SHUT "$1_csex.r1";

IF MTTNy>0 THEN
BEGIN
  OUT "$1_cseo.r1";
  write "MATRIX MTTY(", MTTNy, ",", 1, ")$"; 
  SHUT "$1_cseo.r1";
END;

%%Create the _cse.r file
OUT "$1_cse.r2";
write "%File: $1_cse.r";
in ("$1_cse_write.r");

# Copyright (C) 2000 by Peter J. Gawthrop

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.7  2001/04/11 09:44:26  gawthrop
## Fixed cc and c problems to do with pow(x,y) and integers
## mtt/lib/reduce/fix_c.r is included in rdae2dae and cse2smx_lang for
## -c, -cc and -oct options
##
## Revision 1.6  2000/11/29 20:48:53  peterg
## Zapped unnecessary Npar creation
##
## Revision 1.5  2000/11/09 10:12:24  peterg
## Removed debugging line
##
## Revision 1.4  2000/10/14 16:19:54  peterg

	;;
    ode)
        matrices='dX'
        ns="$Nx"
        ms="1"
	;;
    odeo)
        matrices='Y'
        ns="$Ny"
        ms="1 1"
	;;
    sm)
        matrices='A B C D'
        ns="$Nx $Nx $Ny $Ny"
        ms="$Nx $Nu $Nx $Nu"
	;;

% %%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Field separator
if nargin<3
  FS = ';';
end;

arg = '';
if strcmp(args, '')==0
  L = length(args);
  args_count = 0;
  for i=1:n
    arg_count = 0;
    arg = '';
    if args_count == L
      break;
    end;  
    while args_count < L
      args_count = args_count+1;
      arg_count = arg_count+1;
      ch = str2ch(args,args_count);
      if ch==FS
	break;
      end;
      arg = [arg ch];
    end;
  end;
end;

#! /bin/sh

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

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.57.2.1  2001/05/04 04:07:24  geraint
## Numerical solution of algebraic equations.
## sys_ae.cc written for unsolved inputs.
## Solution of equations using hybrd from MINPACK (as used by Octave fsolve).
##
## Revision 1.57  2001/04/01 03:38:54  geraint
## Reset row to zero after write to file, ready for subsequent runs.
## Eliminates SIGSEGV in Octave when _ode2odes called multiple times.
##
## Revision 1.56  2001/03/30 15:13:58  gawthrop
## Rationalised simulation modes to each return mtt_data
##
## Revision 1.55  2001/03/27 13:21:59  geraint
## Octave version compatibility for save_ascii_data(_for_plotting).
##
## Revision 1.54  2001/03/27 01:14:27  geraint

filename=${sys}_ode2odes.${lang}

if [ -n "$3" ]; then
  method=$3    
else
  method=implicit  
fi

if [ -n "$4" ]; then
    algebraic_solver=$4
else
  algebraic_solver="Reduce_Solver"
fi

echo Creating $filename with $method integration method

# Find system constants
Nx=`mtt_getsize $sys x` # States
Nu=`mtt_getsize $sys u` # Inputs 
Ny=`mtt_getsize $sys y` # Inputs  

	feval_header=parse.h
	save_ascii_data_function=save_ascii_data_for_plotting
	;;
esac

cat <<EOF  > $filename
#include <octave/oct.h>
#include <octave/load-save.h>
#include <octave/lo-mappers.h>
#include <octave/ov-struct.h>
#include <octave/variables.h>

#ifndef STANDALONE
#include <octave/${feval_header}>
#endif

#include "${sys}_def.h"
#include "${sys}_sympar.h"

#ifdef STANDALONE
#include <csignal>
#include <fstream>

#include "mtt_${algebraic_solver}.hh"

extern ColumnVector F${sys}_ae (
	ColumnVector &x,
	ColumnVector &u,
	const double &t,
	ColumnVector &par);

extern ColumnVector F${sys}_input (
	ColumnVector &x,
	ColumnVector &y,
	const double &t,
	ColumnVector &par);

	ColumnVector &u,
	const double &t,
	ColumnVector &par);
 
EOF
fi
cat <<EOF >> $filename

void set_signal_handlers (void);

#endif // STANDALONE

inline ColumnVector
mtt_input (ColumnVector &x,
	   ColumnVector &y,
	   const double &t,
	   ColumnVector &par)
{
#ifdef STANDALONE
  static ColumnVector u  (MTTNU);
  static ColumnVector ui (MTTNYZ);
  static ColumnVector U  (MTTNU+MTTNYZ);
 
  static ${algebraic_solver} ae(F${sys}_ae,MTTNPAR,MTTNU,MTTNX,MTTNY,MTTNYZ);

  u = F${sys}_input (x, y, t, par);
  if (MTTNYZ == 0)
    {
      return u;
    }
  else
    {
      return ae.solve(x,u,t,par);
    }
#else
  static octave_value_list args, f;
  args (0) = octave_value (x);
  args (1) = octave_value (y);
  args (2) = octave_value (t);
  args (3) = octave_value (par);
  f = feval ("${sys}_input", args, 1);

	   ostream &file = cout)
{
  static Matrix data;
  static int row;

  if (dump_data)
    {
       if (row > 0)
         {
           Matrix written_data = data.extract (0, 0, row-1, data.cols ()-1);
           $save_ascii_data_function (file, written_data, "MTT_data");
	 }
       return;
    }

  const int nx = x.length (), ny = y.length ();
  register int col = 0;

  if (0 == row)

#ifdef STANDALONE
void dump_data (ostream &file)
{
  ColumnVector null (0.0);
  mtt_write (0.0, null, null, 0, true, file);
}



void handle_signal (int signum)
{
  // handle some signals to ensure data is written.
  cerr << "# Writing data to MTT.core (signal " << signum << ")" << endl;
  ofstream corefile ("MTT.core");
  dump_data (corefile);
  switch (signum)
    {
      case SIGFPE:
        // Intel chips do not raise SIGFPE for DIVZERO :-(
//        raise (SIGABRT);
        break;
      case SIGINT:
	break;
      case SIGQUIT:
	signal (SIGQUIT, SIG_DFL);
	raise (SIGQUIT);
	break;

# Copyright (C) 2000 by Peter J. Gawthrop

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.31  2001/04/03 14:49:42  gawthrop
## Revised to incorporate new ssim (sensitivity simulation)
## representation (m only just now).
##
## Revision 1.30  2001/03/30 15:13:58  gawthrop
## Rationalised simulation modes to each return mtt_data
##
## Revision 1.29  2001/03/27 13:10:23  geraint
## Improved determination of Octave version.
##
## Revision 1.28  2001/02/17 03:48:17  geraint

# Representation-specific stuff

eqnargs='mttx,mttu,mttt,mttpar'
inputeqnargs='mttx,mtty,mttt,mttpar'
case $rep in
    ae)
	states=yes;
	inputs=yes;
	parameters=yes;
	output=mttyz;
	args=$eqnargs;
	;;
    cse)
	states=yes;
        inputs=yes;
	parameters=yes;
        output='mttdx,mtte'
        args=$eqnargs
	;;

	  states=no;
	  parameters=no;
          declareinputs=no;
        else
	  states=yes;
	  parameters=yes;
          declareinputs=yes
	  declarestates=yes
        fi
	;;
    logic)
	states=no;
	inputs=no;
	parameters=yes;
        output=mttopen

$Lc Inputs $Rc
EOF

N=`n2m 1 $Nu`
for i in $N; do
  echo $constant_declaration 'mttu'$i' = mttu('$i$minusone');'
done
cat <<EOF

$Lc Unknown Inputs $Rc
EOF

Ni=`n2m 1 $Nyz`
for i in $Ni; do
  echo $constant_declaration 'mttui'$i' = mttu('$Nu+$i$minusone');'
done
fi

if [ "$declareinputs" = "yes" ]; then
  declare_vars input
fi

if [ "$declaredummies" = "yes" ]; then

	    arg_type="Octave_map"
	    arg_size="";
	    ;;
	mttu)
	    arg_type="ColumnVector"
	    arg_size="(MTTNU)"
	    ;;
	mttyz)
	    arg_type="ColumnVector"
	    arg_size="(MTTNYZ)"
	    ;;
	mtty)
	    arg_type="ColumnVector"
	    arg_size="(MTTNY)"
	    ;;
	mttt)
	    arg_type="const double"
	    arg_size=""

# Copyright (C) 2000 by Peter J. Gawthrop

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.27  2001/02/08 00:39:56  geraint
## Removed clear from code generated by: mtt -c -i euler -o
##
## Revision 1.26  2000/11/27 11:52:10  peterg
## Changed some greps to be case insensitive
##
## Revision 1.25  2000/10/14 08:50:07  peterg
## Use new mtt_header
##
## Revision 1.24  2000/10/11 09:07:53  peterg

echo Creating $out

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

#Header
case $rep in
    ae|csex|cseo|ode|odeo)
    vectorise=yes
       ;;
    *)
esac

    mtt_header  $1 $rep $ext  > mtt_junk

#define pi 3.14159 // Predefine pi

#ifndef HAVE_USEFUL_FUNCTIONS_HH
#define HAVE_USEFUL_FUNCTIONS_HH


#ifndef __cplusplus
#define inline			/* strip */



#endif // ! __cplusplus


static inline double
max (const double x1, const double x2)
{
  return static_cast<double>((x1 >= x2) ? x1 : (x1 < x2) ? x2 : 0);
}

static inline double
min (const double x1, const double x2)
{
  return static_cast<double>((x1 <= x2) ? x1 : (x1 > x2) ? x2 : 0);
}

static inline double
nonsingular (const double x)
{
  return static_cast<double>((x == 0) ? 1.0e-30 : x);
}

static inline double
sign (const double x)
{
  return static_cast<double>((x > 0) ? +1 : (x < 0) ? -1 : 0);
}


// Octave functions
#ifdef __cplusplus

#include "mtt_HJ_Solver.hh"

// http://www.netlib.org/opt/hooke.c
// Hooke and Jeeves solution
  
HJ_Solver *HJ_Solver::static_ptr;

double
HJ_Solver::f (double tryUi[], int nyz)
{
  for (int i = 0; i < nyz; i++)
    HJ_Solver::static_ptr->_ui(i) = tryUi[i];

  HJ_Solver::static_ptr->_yz = HJ_Solver::static_ptr->eval(HJ_Solver::static_ptr->_ui);

  double ans = 0.0;
  for (int i = 0; i < nyz; i++)
    ans += HJ_Solver::static_ptr->_yz(i)*HJ_Solver::static_ptr->_yz(i);
  return ans;
}

void
HJ_Solver::Solve (void)
{
  double startpt[_nyz];
  double endpt[_nyz];
  for (int i = 0; i < _nyz; i++)
    startpt[i] = _ui(i);
  hooke(_nyz,startpt,endpt);
  for (int i = 0; i < _nyz; i++)
    _ui(i) = endpt[i];
}

  // adapted from http://www.netlib.org/opt/hooke.c:

/* Nonlinear Optimization using the algorithm of Hooke and Jeeves  */
/*      12 February 1994        author: Mark G. Johnson            */


/* Find a point X where the nonlinear function f(X) has a local    */
/* minimum.  X is an n-vector and f(X) is a scalar.  In mathe-     */
/* matical notation  f: R^n -> R^1.  The objective function f()    */
/* is not required to be continuous.  Nor does f() need to be      */
/* differentiable.  The program does not use or require            */
/* derivatives of f().                                             */

/* The software user supplies three things: a subroutine that      */
/* computes f(X), an initial "starting guess" of the minimum point */
/* X, and values for the algorithm convergence parameters.  Then   */
/* the program searches for a local minimum, beginning from the    */
/* starting guess, using the Direct Search algorithm of Hooke and  */
/* Jeeves.                                                         */

/* This C program is adapted from the Algol pseudocode found in    */
/* "Algorithm 178: Direct Search" by Arthur F. Kaupe Jr., Commun-  */
/* ications of the ACM, Vol 6. p.313 (June 1963).  It includes the */
/* improvements suggested by Bell and Pike (CACM v.9, p. 684, Sept */
/* 1966) and those of Tomlin and Smith, "Remark on Algorithm 178"  */
/* (CACM v.12).  The original paper, which I don't recommend as    */
/* highly as the one by A. Kaupe, is:  R. Hooke and T. A. Jeeves,  */
/* "Direct Search Solution of Numerical and Statistical Problems", */
/* Journal of the ACM, Vol. 8, April 1961, pp. 212-229.            */

/* Calling sequence:                                               */
/*  int hooke(nvars, startpt, endpt, rho, epsilon, itermax)        */
/*                                                                 */
/*     nvars       {an integer}  This is the number of dimensions  */
/*                 in the domain of f().  It is the number of      */
/*                 coordinates of the starting point (and the      */
/*                 minimum point.)                                 */
/*     startpt     {an array of doubles}  This is the user-        */
/*                 supplied guess at the minimum.                  */
/*     endpt       {an array of doubles}  This is the location of  */
/*                 the local minimum, calculated by the program    */
/*     rho         {a double}  This is a user-supplied convergence */
/*                 parameter (more detail below), which should be  */
/*                 set to a value between 0.0 and 1.0.  Larger     */
/*                 values of rho give greater probability of       */
/*                 convergence on highly nonlinear functions, at a */
/*                 cost of more function evaluations.  Smaller     */
/*                 values of rho reduces the number of evaluations */
/*                 (and the program running time), but increases   */
/*                 the risk of nonconvergence.  See below.         */
/*     epsilon     {a double}  This is the criterion for halting   */
/*                 the search for a minimum.  When the algorithm   */
/*                 begins to make less and less progress on each   */
/*                 iteration, it checks the halting criterion: if  */
/*                 the stepsize is below epsilon, terminate the    */
/*                 iteration and return the current best estimate  */
/*                 of the minimum.  Larger values of epsilon (such */
/*                 as 1.0e-4) give quicker running time, but a     */
/*                 less accurate estimate of the minimum.  Smaller */
/*                 values of epsilon (such as 1.0e-7) give longer  */
/*                 running time, but a more accurate estimate of   */
/*                 the minimum.                                    */
/*     itermax     {an integer}  A second, rarely used, halting    */
/*                 criterion.  If the algorithm uses >= itermax    */
/*                 iterations, halt.                               */


/* The user-supplied objective function f(x,n) should return a C   */
/* "double".  Its  arguments are  x -- an array of doubles, and    */
/* n -- an integer.  x is the point at which f(x) should be        */
/* evaluated, and n is the number of coordinates of x.  That is,   */
/* n is the number of coefficients being fitted.                   */

/* rho, the algorithm convergence control                          */
/*      The algorithm works by taking "steps" from one estimate of */
/*    a minimum, to another (hopefully better) estimate.  Taking   */
/*    big steps gets to the minimum more quickly, at the risk of   */
/*    "stepping right over" an excellent point.  The stepsize is   */
/*    controlled by a user supplied parameter called rho.  At each */
/*    iteration, the stepsize is multiplied by rho  (0 < rho < 1), */
/*    so the stepsize is successively reduced.                     */
/*      Small values of rho correspond to big stepsize changes,    */
/*    which make the algorithm run more quickly.  However, there   */
/*    is a chance (especially with highly nonlinear functions)     */
/*    that these big changes will accidentally overlook a          */
/*    promising search vector, leading to nonconvergence.          */
/*      Large values of rho correspond to small stepsize changes,  */
/*    which force the algorithm to carefully examine nearby points */
/*    instead of optimistically forging ahead.  This improves the  */
/*    probability of convergence.                                  */
/*      The stepsize is reduced until it is equal to (or smaller   */
/*    than) epsilon.  So the number of iterations performed by     */
/*    Hooke-Jeeves is determined by rho and epsilon:               */
/*          rho**(number_of_iterations) = epsilon                  */
/*      In general it is a good idea to set rho to an aggressively */
/*    small value like 0.5 (hoping for fast convergence).  Then,   */

  
/*    if the user suspects that the reported minimum is incorrect  */
/*    (or perhaps not accurate enough), the program can be run     */
/*    again with a larger value of rho such as 0.85, using the     */
/*    result of the first minimization as the starting guess to    */
/*    begin the second minimization.                               */

/* Normal use: (1) Code your function f() in the C language        */
/*             (2) Install your starting guess {or read it in}     */
/*             (3) Run the program                                 */
/*             (4) {for the skeptical}: Use the computed minimum   */
/*                    as the starting point for another run        */

/* Data Fitting:                                                   */
/*      Code your function f() to be the sum of the squares of the */
/*      errors (differences) between the computed values and the   */
/*      measured values.  Then minimize f() using Hooke-Jeeves.    */
/*      EXAMPLE: you have 20 datapoints (ti, yi) and you want to   */
/*      find A,B,C such that  (A*t*t) + (B*exp(t)) + (C*tan(t))    */
/*      fits the data as closely as possible.  Then f() is just    */
/*      f(x) = SUM (measured_y[i] - ((A*t[i]*t[i]) + (B*exp(t[i])) */
/*                                + (C*tan(t[i]))))^2              */
/*      where x[] is a 3-vector consisting of {A, B, C}.           */

/*                                                                 */
/*  The author of this software is M.G. Johnson.                   */
/*  Permission to use, copy, modify, and distribute this software  */
/*  for any purpose without fee is hereby granted, provided that   */
/*  this entire notice is included in all copies of any software   */
/*  which is or includes a copy or modification of this software   */
/*  and in all copies of the supporting documentation for such     */
/*  software.  THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT    */
/*  ANY EXPRESS OR IMPLIED WARRANTY.  IN PARTICULAR, NEITHER THE   */
/*  AUTHOR NOR AT&T MAKE ANY REPRESENTATION OR WARRANTY OF ANY     */
/*  KIND CONCERNING THE MERCHANTABILITY OF THIS SOFTWARE OR ITS    */
/*  FITNESS FOR ANY PARTICULAR PURPOSE.                            */
/*                                                                 */

double
HJ_Solver::best_nearby(double *delta,
		       double *point,
		       double prevbest,
		       int    nvars)
{
  double          z[VARS];
  double          minf, ftmp;
  int             i;
  minf = prevbest;
  for (i = 0; i < nvars; i++)
    z[i] = point[i];
  for (i = 0; i < nvars; i++) {
    z[i] = point[i] + delta[i];
    ftmp = f(z, nvars);
    if (ftmp < minf)
      minf = ftmp;
    else {
      delta[i] = 0.0 - delta[i];
      z[i] = point[i] + delta[i];
      ftmp = f(z, nvars);
      if (ftmp < minf)
	minf = ftmp;
      else
	z[i] = point[i];
    }
  }
  for (i = 0; i < nvars; i++)
    point[i] = z[i];
  return (minf);
}


int
HJ_Solver::hooke (int	nvars,			// MTTNYZ
		  double	startpt[],	// user's initial guess
		  double	endpt[],	// result
		  double	rho,		// geometric shrink factor
		  double	epsilon,	// end value stepsize
		  int	itermax)		// max # iterations 
{
  const int VARS = _nyz;
  double          delta[VARS];
  double          newf, fbefore, steplength, tmp;
  double          xbefore[VARS], newx[VARS];
  int             i, j, keep;
  int             iters, iadj;
  for (i = 0; i < nvars; i++) {
    newx[i] = xbefore[i] = startpt[i];
    delta[i] = fabs(startpt[i] * rho);
    if (delta[i] == 0.0)
      delta[i] = rho;
  }
  iadj = 0;
  steplength = rho;
  iters = 0;
  fbefore = f(newx, nvars);
  newf = fbefore;
  while ((iters < itermax) && (steplength > epsilon)) {
    iters++;
    iadj++;
  /* find best new point, one coord at a time */
  for (i = 0; i < nvars; i++) {
    newx[i] = xbefore[i];
  }
  newf = best_nearby(delta, newx, fbefore, nvars);
  /* if we made some improvements, pursue that direction */
  keep = 1;
  while ((newf < fbefore) && (keep == 1)) {
    iadj = 0;
    for (i = 0; i < nvars; i++) {
      /* firstly, arrange the sign of delta[] */
      if (newx[i] <= xbefore[i])
	delta[i] = 0.0 - fabs(delta[i]);
      else
	delta[i] = fabs(delta[i]);                                   /* now, move further in this direction */
      tmp = xbefore[i];
      xbefore[i] = newx[i];
      newx[i] = newx[i] + newx[i] - tmp;
    }
    fbefore = newf;
    newf = best_nearby(delta, newx, fbefore, nvars);
    /* if the further (optimistic) move was bad.... */
    if (newf >= fbefore)
      break;
    /* make sure that the differences between the new */
    /* and the old points are due to actual */
    /* displacements; beware of roundoff errors that */
    /* might cause newf < fbefore */
    keep = 0;
    for (i = 0; i < nvars; i++) {
      keep = 1;
      if (fabs(newx[i] - xbefore[i]) >
	  (0.5 * fabs(delta[i])))
	break;
      else
	keep = 0;
    }
  }
  if ((steplength >= epsilon) && (newf >= fbefore)) {
    steplength = steplength * rho;
    for (i = 0; i < nvars; i++) {
      delta[i] *= rho;
    }
  }
}
for (i = 0; i < nvars; i++)
  endpt[i] = xbefore[i];
return (iters);
}

#include "mtt_Solver.hh"

class HJ_Solver : public Solver {

  // http://www.netlib.org/opt/hooke.c
  // Hooke and Jeeves solution
  
public:
  
  HJ_Solver (sys_ae ae,
	     const int npar,
	     const int nu,
	     const int nx,
	     const int ny,
	     const int nyz)
    : Solver (ae,npar,nu,nx,ny,nyz)
  { static_ptr = this; VARS = nyz; };
  
  static double
  f (double tryUi[], int nyz);
  
protected:

  void
  Solve (void);

  double
  best_nearby (double    	delta[],
	       double   	point[],
	       double   	prevbest,
	       int      	nvars);
  
  int
  hooke (int	nvars,			  // MTTNYZ
	 double	startpt[],		  // user's initial guess
	 double	endpt[],		  // result
	 double	rho		= 0.05,	  // geometric shrink factor
	 double	epsilon		= 1e-3,	  // end value stepsize
	 int	itermax		= 5000);  // max # iterations 

private:

  int VARS;

public:

  static HJ_Solver *static_ptr;

};

#include "mtt_Hybrd_Solver.hh"

// http://www.netlib.org/minpack/hybrd.f
// used by Octave's fsolve
  
Hybrd_Solver *Hybrd_Solver::static_ptr;

ColumnVector
Hybrd_Solver::f_hybrd (const ColumnVector &tryUi)
{
  Hybrd_Solver::static_ptr->_yz = Hybrd_Solver::static_ptr->eval(tryUi);
  return Hybrd_Solver::static_ptr->_yz;
}

void
Hybrd_Solver::Solve (void)
{    
  int info;
  static int input_errors;
  static int user_errors;
  static int convergences;
  static int progress_errors;
  static int limit_errors;
  static int unknown_errors;
  
  NLFunc fcn(&Hybrd_Solver::f_hybrd);
  NLEqn	 eqn(Solver::_ui,fcn);
  //  eqn.set_tolerance(0.01);
  Solver::_ui = eqn.solve(info);

  switch (info)
    {
    case 1:
      convergences++;
      break;
    case -2:
      input_errors++;
      break;
    case -1:
      user_errors++;
      break;
    case 3:
      progress_errors++;
      break;
    case 4:
      //      if (abs(eval(_ui).max()) > 1.0e-6)
      limit_errors++;
      //      else
      //	convergences++;
      break;
    default:
      unknown_errors++;
      break;
    }
  if (1 != info)
    {
      cerr << "input (" << input_errors << ") "
	   << "  user (" << user_errors << ") "
	   << "  converge (" << convergences << ") "
	   << "  progress (" << progress_errors << ") "
	   << "  limit (" << limit_errors << ")"
	   << "  unknown (" << unknown_errors << ")"
	   << "  (max error = " << abs(eval(_ui).max()) << ")" << endl;
    }
}

#include "mtt_Solver.hh"
#include <octave/NLEqn.h>

class Hybrd_Solver : public Solver {

  // http://www.netlib.org/minpack/hybrd.f
  // used by Octave's fsolve
  
public:

  Hybrd_Solver (sys_ae ae,
		const int npar,
		const int nu,
		const int nx,
		const int ny,
		const int nyz)
  : Solver (ae,npar,nu,nx,ny,nyz)
  {
    static_ptr = this;
  }

  static ColumnVector
  f_hybrd (const ColumnVector &tryUi);

protected:

  void
  Solve (void);

public:

  static Hybrd_Solver *static_ptr;

};

#include "mtt_Reduce_Solver.hh"


void
Reduce_Solver::Solve (void)
{
  cerr << "Error:"
       << " Symbolic solution of equations failed during model build" << endl
       << "       Try using one of the other algebraic solution methods" << endl;
}

#include "mtt_Solver.hh"

class Reduce_Solver : public Solver {

  // Dummy class
  // This will not be used unless the Reduce solver has failed earlier
  // in the model build process

public:

  Reduce_Solver (sys_ae ae,
		 const int npar,
		 const int nu,
		 const int nx,
		 const int ny,
		 const int nyz)
    : Solver (ae,npar,nu,nx,ny,nyz)
  { ; };
	
  void
  Solve (void);

};
   

#include "mtt_Solver.hh"

Solver::Solver (sys_ae ae,
		const int npar,
		const int nu,
		const int nx,
		const int ny,
		const int nyz)
{
  _ae  = ae;
  _np  = npar; 
  _nu  = nu;
  _nx  = nx;
  _ny  = ny;
  _nyz = nyz;

  _uui = ColumnVector (_nu+_nyz);
};

ColumnVector
Solver::solve (const ColumnVector	&x,
	       const ColumnVector	&u,
	       const double		&t,
	       const ColumnVector	&par)
{
  _x = x;
  _uui.insert(u,0);
  _t = t;
  _par = par;
  _ui = ColumnVector(_nyz,0.2);
  Solve ();
  _uui.insert(_ui,_nu);
  return _uui;
}    

ColumnVector
Solver::eval (const ColumnVector &ui)
{
  _uui.insert(ui,_nu);
  return _ae (_x, _uui, _t, _par);
}

#ifndef MTT_SOLVER
#define MTT_SOLVER

#include <cmath>
#include <cstdlib>
#include <iostream>

#include <octave/oct.h>

class Solver {

  typedef ColumnVector (*sys_ae) // pointer to F${sys}_ae function
    (ColumnVector &,ColumnVector &,const double &t,ColumnVector &);

public:

  Solver (sys_ae ae,
	  const int npar,
	  const int nu,
	  const int nx,
	  const int ny,
	  const int nyz);

  ColumnVector
  solve (const ColumnVector	&x,
	 const ColumnVector	&u,
	 const double		&t,
	 const ColumnVector	&par);
  
  ColumnVector
  eval (const ColumnVector	&ui);
  
protected:
  
  virtual void
  Solve (void) = 0;

protected:

  ColumnVector			_x;
  ColumnVector			_uui;
  double			_t;
  ColumnVector			_par;

  ColumnVector  		_ui;
  ColumnVector          	_yz;

  int   			_nu;
  int				_np;
  int				_nx;
  int				_ny;
  int				_nyz;
  
  sys_ae			_ae;

};

#endif // MTT_SOLVER