User:Karthik0112358/Coefficient Calculator

/*	ASCEND modelling environment
	Copyright (C) 1997, 2006 Carnegie Mellon University
	Copyright (C) 1993, 1994 Joseph James Zaher, Benjamin Andrew Allan
	Copyright (C) 1993 Joseph James Zaher
	Copyright (C) 1990 Thomas Guthrie Epperly, Karl Michael Westerberg

	This program is free software; you can redistribute it and/or modify
	it under the terms of the GNU General Public License as published by
	the Free Software Foundation; either version 2, or (at your option)
	any later version.

	This program is distributed in the hope that it will be useful,
	but WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
	GNU General Public License for more details.

	You should have received a copy of the GNU General Public License
	along with this program; if not, write to the Free Software
	Foundation, Inc., 59 Temple Place - Suite 330,
	Boston, MA 02111-1307, USA.
*/
/*
KARTHIK C S, June 2011
*/
#include <math.h>
#include "relation_util.h"
#include <ascend/general/ascMalloc.h>
#include <ascend/general/panic.h>
#include "relation_type.h"
#include "expr_types.h"
#include "instance_enum.h"
#include "functype.h"
#include "mathinst.c"
/*------------------------------------------------------------------------------
  FINDING COEFFICIENTS OF THE DERIVATIVE TERMS IN A RELATION
*/

/**
	Compute the coefficients of the derivative terms in the relation (compiler-side routine)

	@param i Instance pointer through which the coefficients of the derivative terms of the relation are found
	@param coefficient pointer to a double in which the coefficients shall be stored (must have already been allocated)
	@param Info of structure Info_abt_deriv which contains the real value of all the derivative variables

 	@return 0 on success, and otherwise on failure

	Here we compute the coefficients of the derivative terms by first extracting each term. 
	We then check if its a derivative variable and if so we compute its coefficient by considering the terms previous to it and after it. 
	In this process we also check if the equation is linear wrt its derivatives. Suitable error messages are displayed if its not.
*/
/* Started reorganizing the entire Inspect_Equation(). Providing modularity for efficiency and flexibility.*/
/*
static int Inspect_Equation(struct Info_abt_deriv Info,int *order,double *var,double *func,int *ope,unsigned long length_rel){
	int i,j,k;
	int varc=-1,opec=-1,funcc=-1,numc=-1;
	int flag,count,c,denom_flag,point;
	unsigned long bin_num,pass,new_bin_num;
	for (i=0;i<Info.no_of_deriv;i++){
		flag=0;
		for(j=0;j<length_rel;j++){
			switch(order[j]){
				case 1:
					++numc;
					break;
				case 2:
					++varc;
					//Write Error catching analysis code here
					if(Info.deriv_values[i]==var[varc]) {
						if(order[j-1]==2) {
							count=1;
							c=1;
							k=0;
							while(count!=0&&(j+c)<=length_rel) {
								count=(order[c+j]==3)?--count:++count;
								k=(order[c+j]==3)?++k:k;
								++c;
							}
							--c;
							if(ope[opec+k]!=1||ope[opec+k]!=2){
								CONSOLE_DEBUG("The coefficient of the derivative term contains a non constant term");
							}
						}else if(order[j+1]==2) {
							count=1;
							c=2;
							k=0;
							while(count!=0&&(j+c)<=length_rel){	
								count=(order[c+j]==3)?--count:++count;
								k=(order[c+j]==3)?++k:k;
								++c;
							}
							--c;
							if(ope[opec+k]!=1||ope[opec+k]!=2){
								CONSOLE_DEBUG("The coefficient of the derivative term contains a non constant term");
							}	
						}
						c=1;
						count=0;
						denom_flag=0;
						while(order[j+c]!=4&&(j+c)<=length_rel) {
							if(order[c+j]==3) {
								++count;
								if(ope[opec+count]==4) {
									if(count>c-count) {
										denom_flag=(denom_flag==0)?1:0;
									}
								}					
							}							
							++c;						
						}
						if(denom_flag==1){
							CONSOLE_DEBUG("The derivative is in denominator");
						}
						c=1;
						count=0;
						while(order[j+c]!=4&&(j+c)<=length_rel) {
							if(order[c+j]==3) {
								++count;
								if(ope[opec+count]==6) {
									if(count>c-count) {
										CONSOLE_DEBUG("The derivative is in index");
									}
								}					
							}							
							++c;						
						}
						/* We implement here an algorithm which finds the span of an operator. The span of an operator 
						is defined as the number of terms the operator links in an expression. We see if the derivative 
						is in the upper span (the set of terms to the left of an operator in infix form) of the operator, 
						by using a variable to point at the location of the derivative in a binary representation of the 
						relation where 0 means a number or variable and 1 means an operator. Then we display an error is 
						the derivative is raised to a higher/lower power.*/
						c=0;
						count=-1;
						bin_num=0;
						pass=0;
						/* Need to write code here to check if the derivative is being raised to power 1. 
						The idea is to locate every ^ operator and see if the item immediately following 
						it in the infix notation is 1 which does not have any prior connection with any other
						operator. */
						/* As a second thought, I think I will segregate each Error case as a function, considering the recursive 							algorithms used in some of the cases. */
						if(flag==0){						
							while(order[c]!=4){
								++c;		
							}
							for(k=c-1;k>=0;k--) {
								if(order[k]==3) {
									bin_num+=pow(2,c-k-1);
								}
							}
							point=j;
							for(k=c-1;k>=0;k--) {
								if(bin_num-pass>0) {	
									++count;									
									pass+=pow(2,k);							
									if(point+k==c-3) {
										if(ope[count]==6) {
											CONSOLE_DEBUG("The derivative is raised to a power other than one");
										} else {
											
										}
									} else {
	
									}
								}
							}	
						}else {
							while(order[c]!=4){
								++c;
							}
							++c;
							while(c<=length_rel) {
								++c;						
							}
						}
					}
					break;
				case 3:
					++opec;
					break;
				case 4:
					flag=1;
					break;
				default:
					/* Give Error Message */	
					break;		
			}
		}		
	}
	return 1;
}
*/

/* So here it begins */

static int Inspect_Equation(struct Info_abt_deriv Info,int *order,double *var,double *func,int *ope,unsigned long length_rel){
	double var_index;
	char *left_bin_num,*right_bin_num;
	left_bin_num=find_left_bin_number(order,length_rel);
	right_bin_num=find_right_bin_number(order,length_rel);
/* We write here an inspection for special cases, like : 3=4+Q(x,y)/x, where is Q(x,y) is a linear function in x */

/* End of special cases */

/* We start inspecting general cases */

	flag=0;
	for(j=0;j<length_rel;j++){
		switch(order[j]){
			case 1:
				++numc;
				break;
			case 2:
				++varc;
				/* Write Error catching analysis code here */
				var_index=find_variable(Info,var[varc]);
				if(var_index!=0){
					
				}
				break;
			case 3:
				++opec;
				break;
			case 4:
				flag=1;
				break;
			default:
				/* Give Error Message */	
				break;		
		}
	}
	return 1;
}

void coefficient_match(double *coefficient,struct Info_abt_deriv Info,int *order,double *num,double *var,double *func,int *ope, unsigned long length_rel){
	int i,j,k;
	double value;
	int flag;
	int varc=-1,opec=-1,funcc=-1,numc=-1;
	int c;
	if(Inspect_Equation(Info,order,var,func,ope,length_rel)==1){
		for (i=0;i<Info.no_of_deriv;i++){
			value=1.0;
			flag=0;
			for(j=0;j<length_rel;j++){
				switch(order[j]){
					case 1:
						++numc;
						break;
					case 2:
						++varc;
						if(Info.deriv_values[i]==var[varc]) {
							if(order[j-1]==1&&order[j+1]==3) {
								if(ope[opec+1]==3) {
									value=num[numc];
								}
							}else if(order[j+1]==1&&order[j+2]==3) {
								if(ope[opec+1]==3) {
									value=num[numc+1];
								}
							}else {							
								while(order[c+j]!=4){
									c=1;
									while(order[c+j]==2||order[c+j]==1) {
										c++;
									}		
								}								
							}
						}
						break;
					case 3:
						++opec;
						break;
					case 4:
						flag=1;
						break;
					default:
					/* Give Error Message */	
						break;		
				}
			}	
		}
	}else {
		CONSOLE_DEBUG("Equation cannot be accepted by solver");
	}	
}

static int coefficient_calculator(struct Instance *i,double *coefficient,struct Info_abt_deriv Info){
	struct relation *r;
	enum Expr_enum reltype;
	unsigned long t;       			/* the current term in the relation r */
	unsigned long num_var; 			/* the number of variables in the relation r */
	unsigned long v;      			/* loop variable */
	int lhs;               			/* looking at left(=1) or right(=0) hand side of r */
	unsigned long length_lhs, length_rhs;
	double *num;				/* pointer for storing all constants in relation */
	double *var;				/* pointer for storing all real value of variables in relation */
	int *ope;				/* pointer for storing all operators in relation */
	double *func;				/* pointer for storing all functions in relation */
	int *order;				/* pointer for storing order of all constants, variables and operators in postfix form in relation */
	int varc=0,numc=0,orderc=0,opec=0,funcc=0;	/* counters for various arrays */
	CONST struct relation_term *term;
	CONST struct Func *fxnptr;
	r = (struct relation *)GetInstanceRelation(i, &reltype);
	if( r == NULL ) {
		CONSOLE_DEBUG("null relation\n");
		return -1;
	}
	num_var = NumberVariables(r);
	length_lhs = RelationLength(r, 1);
	length_rhs = RelationLength(r, 0);
	num=ASC_NEW_ARRAY(double,length_lhs+length_rhs);
	var=ASC_NEW_ARRAY(double,length_lhs+length_rhs);
	ope=ASC_NEW_ARRAY(int,length_lhs+length_rhs);
	func=ASC_NEW_ARRAY(double,length_lhs+length_rhs);
	order = ASC_NEW_ARRAY(int,(length_lhs+length_rhs+1));
	lhs = 1;
	t = 0;

	for (v=0; v<length_lhs+length_rhs; v++) {
/* Initializing all entries to 0 */
		order[v]=0;
		func[v]=0.0;
		ope[v]=0;
		num[v]=0.0;
		var[v]=0.0;
	}
	order[length_lhs+length_rhs]=0;
	while(1) {
		if( lhs && (t >= length_lhs) ) {
/* need to switch to the right hand side--if it exists */
			if( length_rhs ) {
				lhs = t = 0;
				order[orderc++]=4;
			} else {
/* 
	We know that (length_lhs+length_rhs>0) and that (length_rhs==0),
	the length_lhs must be > 0
*/
				coefficient_match(coefficient,Info,order,num,var,func,ope,length_lhs+length_rhs+1);
			        return 0;
      			}
    		} else if( (!lhs) && (t >= length_rhs) ) {
/* we have processed both sides, quit */
				coefficient_match(coefficient,Info,order,num,var,func,ope,length_lhs+length_rhs+1);
			        return 0;
		}
		term = NewRelationTerm(r, t++, lhs);
		switch( RelationTermType(term) ) {
			case e_zero:
				num[numc++]=0.0;
				order[orderc++]=1;
				break;
			case e_real:
				num[numc++]=(double)TermReal(term);
				order[orderc++]=1;
				break;
			case e_int:
				num[numc++]=(double)TermInteger(term);
				order[orderc++]=1;
				break;
			case e_var:
				var[varc++]= TermVariable(r, term);
				order[orderc++]=2;
				break;
			case e_plus:
				ope[opec++]=1;
				order[orderc++]=3;
				break;
			case e_minus:
				ope[opec++]=2;
				order[orderc++]=3;
				break;
			case e_times:
				ope[opec++]=3;
				order[orderc++]=3;
				break;
			case e_divide:
				ope[opec++]=4;
				order[orderc++]=3;
				break;
			case e_uminus:
				ope[opec++]=5;
				order[orderc++]=3;
				break;
			case e_power:
				ope[opec++]=6;
				order[orderc++]=3;
				break;
			case e_ipower:
				ope[opec++]=7;
				order[orderc++]=3;
				break;
			case e_func:
				ope[opec++]=8;
				order[orderc++]=3;
				fxnptr = TermFunc(term);
				func[funcc] = FuncEval( fxnptr, func[funcc] );
				++funcc;
				break;
			default:
				ASC_PANIC("Unknown relation term type");
				break;
		}

	}
}