#!/usr/bin/perl
#				Problem Statement
# Analysis of an amplifier circuit involves taking many simple equations and
# solving them simultaneously.  For small signal amplitudes, the model can be
# linearized and solved using matrix algebra with matlab.  This generally
# involves taking N equations and forming a NxN matrix, inverting it, and
# acting on a column vector of input values and constants.  This matrix is
# generally very sparse, so writing it out in row/column form is a very poor
# way to represent the model.  A much more compact and readable expression
# of the model is the original linear equations in plain text.  These can
# then be automatically scanned and converted into matrix form for solving.
#
#				Synopsis
# mpp.pl - a primitive preprocessor for matlab .m files
#
#	Usage: mpp.pl [-d] {one or more mpp prototypes}
#
# The mpp preprocessor reads an m-file prototype like the example below,
# and outputs a valid m-function containing valid matlab statements that
# solve for the unknowns in terms of the inputs and known constants.  The
# option -d generates verbose debug output to standard output.
#
#    % example mpp prototype m-file
#    R2 = 50; % Ohms
#    Ie0 = 0.011; % A
#    Vtr = 0.025; % V
#    Q0 = Ie0/Vtr; % /Ohms
#    Vcc = 5 + ?; % V
#    Vb0 = 1.9 + ?; % V
#    0 = Vcc - Ie*R2 - Vb;
#    0 = Ie - Ie0 - Vb*Q0 + Vb0*Q0;
#
# The first 4 lines in the above program define the values of some constants
# in the model, and are valid matlab m-file statements.  However, the last
# four lines are not valid m-file statements, and must be rewritten by mpp.
# The two lines that assign Vcc and Vb0 contain the additive factor "?".
# This identifies these variables as inputs that drive the circuit.  The
# quiescent (dc) value of the inputs is obtained by replacing the "?" with 0.
# The last two lines express the constraints that the model must satisfy.
# The constraints must be written as shown, with 0 as the assignment target.
# The only operators allowed on the right-hand side are +, -, and *.  Any
# variables that are not assigned constant values or declared as inputs are
# treated as unknowns.
#
# author: richard.t.jones at uconn.edu
# version: July 26, 2011

use FileHandle;

our $inputs;
our %inputs;
our %knowns;
our $knowns;
our $unknowns;
our %unknowns;
our $constraints;
our @constraints;
our %keywords = { function => 0,
                  if => 0,
                  else => 0,
                  end => 0,
                  classdef => 0,
                  methods => 0
                };
our %matrix;
our $debug=0;

if (@ARGV < 1) {
   warn "Usage: mpp [-d] {yourfile.mpp}\n";
   exit 1;
}

while (@ARGV) {
   my $infile = shift @ARGV;
   if ($infile =~ /^-d/) {
      $debug = 1;
      next;
   }
   if (!-r $infile) {
      warn "Error: $infile not found or not readable\n";
      next;
   }
   my $outfile = $infile;
   $outfile =~ s/\.mpp$/\.m/ || next;
   my $ifd = FileHandle->new;
   my $ofd = FileHandle->new;
   my $afd = FileHandle->new;
   $ifd->open("<$infile");
   $ofd->open(">$outfile");
   $afd->open(">varnames.m");
   $inputs = 0;
   %inputs = ();
   %knowns = ();
   %matrix = ();
   $unknowns = 0;
   %unknowns = ();
   $constraints = 0;
   @constraints = ();
   $start_of_code = 1;
   while (my $line = <$ifd>) {
      my $statement = $line;
      chop $statement;
      if ($start_of_code && ($statement =~ /^\s*$/)) {
         $afd->print("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
                     "% This file was automatically generated by mpp.pl\n",
                     "% based on $infile so do not make changes here!\n",
                     "% Edit $infile to make changes, and regenerate\n",
                     "% this file using the mpp.pl tool.\n",
                     "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n");
         $ofd->print("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
                     "% This file was automatically generated by mpp.pl\n",
                     "% based on $infile so do not make changes here!\n",
                     "% Edit $infile to make changes, and regenerate\n",
                     "% this file using the mpp.pl tool.\n",
                     "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n");
         $start_of_code = 0;
      }
      while ($statement =~ s/[^%]*\.\.\..*//) {
         my $cline = <$ifd>;
         $line .= $cline;
         chop $cline;
         $statement .= $cline;
      }
      if ($statement !~ /\s*([^% =]+)\s*=\s*([^%=]+)(.*)/) {
         $ofd->print("$line");
         next;
      }
      my ($target,$value,$comment) = ($1,$2,$3);
      if (exists $keywords{$target}) {
         $ofd->print($line);
         next;
      }
      elsif ($target !~ /^0$/) {
         if ($value =~ s/\?/0/) {
            $inputs{$target} = $inputs++;
            $line =~ s/\?/0/;
         }
         $knowns{$target} = $value;
         $ofd->print($line);
         next;
      }
      else {
         if ($debug) {
            print "found constraint $target = $value\n";
         }
         parseformula($value,$constraints++);
         push @constraints, $line;
      }
   }
   $ifd->close;

   my $unknowns = scalar (keys %unknowns);

   if ($debug) {
      print "knowns are: ",join(", ", sort keys %knowns),"\n";
      print "inputs are: ",join(", ", sort keys %inputs),"\n";
      print "unknowns are: ",join(", ", sort keys %unknowns),"\n";
      print "There are $constraints constraints and $unknowns unknowns.\n";
   }

   if ($unknowns != $constraints) {
      warn "Error - ill-posed problem has $unknowns unknowns",
           " but $constraints constraints, cannot continue.\n";
      exit 6;
   }
   generate_solver($ofd,$afd);
   $afd->close;
   $ofd->close;
}

sub parseformula
{
   my $formula = shift;
   my $constraint = shift;
   $formula =~ s/([1-9][.0-9]*)[eEdD]\+([0-9]+)/$1%$2/g;
   $formula =~ s/([1-9][.0-9]*)[eEdD]-([0-9]+)/$1%%$2/g;
   $formula =~ s/\s//g;
   $formula =~ s/;$//g;
   foreach my $sumterm (split /\+/, $formula) {
      my @terms = split /-/, $sumterm;
      for (my $s=0; $s<@terms; ++$s) {
         my $term = $terms[$s];
         my $input = undef;
         my $unknown = undef;
         my $coefficient = undef;
         foreach my $factor (split /\*/, $term) {
            if ($factor =~ /^\s*([A-Za-z_][A-Za-z0-9_]*)\s*$/) {
               if (!exists $knowns{$factor}) {
                  if (defined $unknown || defined $input) {
                     warn "Error - more than one unknown factor found",
                          " in term \"$term\" of formula \"$formula\"\n";
                     exit 3;
                  }
                  else {
                     $unknowns{$factor} = ++$unknowns;
                     $unknown = $factor;
                     next;
                  }
               }
               elsif (exists $inputs{$factor}) {
                  if (defined $input || defined $unknown) {
                     warn "Error - more than one variable factor found",
                          " in term \"$term\" of formula \"$formula\"\n";
                     exit 4;
                  }
                  $input = $factor;
               }
            }
            if (defined $coefficient) {
               $coefficient .= "*$factor";
            }
            else {
               $coefficient = $factor;
            }
         }
         my $sign = ($s > 0)? "-" : "+";
         $unknown = (defined $unknown)? $unknown : "constant";
         $coefficient = ($coefficient)? $coefficient : "1";
         $coefficient =~ s/%%/e-/g;
         $coefficient =~ s/%/e/g;
         if (exists $matrix{$unknown}) {
            $matrix{$unknown}->[$constraint] .= $sign.$coefficient;
         }
         else {
            $matrix{$unknown} = [];
            $matrix{$unknown}->[$constraint] = $sign.$coefficient;
         }
         if (defined $input) {
            unless ($coefficient =~ s/\*$input// ||
                    $coefficient =~ s/$input\*// ||
                    $coefficient =~ s/^$input$/1/) {
               warn "Error - unable to factor symbol $input",
                    " out of expression $coefficient, unable to continue.\n";
               exit 5;
            }
            if (exists $matrix{$input}) {
               $matrix{$input}->[$constraint] .= $sign.$coefficient;
            }
            else {
               $matrix{$input} = [];
               $matrix{$input}->[$constraint] = $sign.$coefficient;
            }
         }
         if ($debug && $unknown ne "constant") {
            print "  coefficient of $unknown is ",
                  $matrix{$unknown}->[$constraint],"\n";
         }
         if ($debug && defined $input) {
            print "  coefficient of $input is ",
                  $matrix{$input}->[$constraint],"\n";
         }
      }
   }
   if ($debug && exists $matrix{"constant"}->[$constraint]) {
      print "  constant term is ", $matrix{"constant"}->[$constraint],"\n";
   }
   elsif ($debug) {
      print "  constant term is 0\n";
   }
}

sub generate_solver
{
   my $ofd = shift;
   my $afd = shift;
   my $columns = 0;
   my @columns = ();
   foreach my $column (sort {$unknowns{$a} <=> $unknowns{$b}} keys %unknowns) {
      $column[$columns++] = $column;
   }
   $column[$columns++] = "constant";
   foreach my $input (sort {$inputs{$a} <=> $inputs{$b}} keys %inputs) {
      $column[$columns++] = $input;
   }
   $afd->print("Nvariables = $columns\n");
   $afd->print("Nconstraints = $constraints\n");
   $ofd->print("  Nvariables = $columns;\n");
   for (my $var=0; $var < $columns; ++$var) {
      if ($var == 0) {
         $afd->print("% these are the unknowns\n");
      }
      elsif ($var == $constraints) {
         $afd->print("% these are the external inputs\n");
      }
      $afd->print("i$column[$var] = ",$var+1,";\n");
   }
   $ofd->print("% these are the constraints\n");
   $ofd->print("  Nconstraints = $constraints;\n");
   $ofd->print("  constraint = zeros(Nconstraints,Nvariables);\n");
   for (my $c=0; $c < $constraints; ++$c) {
      $ofd->print("% ",$constraints[$c]);
      for (my $u=0; $u < $columns; ++$u) {
         if (exists $matrix{$column[$u]} &&
             defined $matrix{$column[$u]}->[$c]) {
            $ofd->print("  constraint(",$c+1,",",$u+1,") = ",
                        $matrix{$column[$u]}->[$c],"; % $column[$u]\n");
         }
      }
   }
   $ofd->print("\n","% now solve for the desired variable\n",
                    "% in response to the designated input signal\n",
         "  ivar = index_lookup(var);\n",
         "  iinp = index_lookup(inp);\n",
         "  solution = constraint(:,1:Nconstraints)\\constraint(:,iinp);\n",
         "  if ivar == 0\n",
         "    response = -solution;\n",
         "  else\n",
         "    response = -solution(ivar);\n",
         "  end\n");
  $ofd->print("\n",
         "  function ivar = index_lookup(var)\n",
         "    if isnumeric(var)\n",
         "      ivar = var;\n");
  for ($col=0; $col<$columns; ++$col) {
    $ofd->print("    elseif strcmp(var,'$column[$col]')\n",
                "      ivar = ",$col+1,";\n");
  }
  $ofd->print("    end\n","  end\n","end\n");
}