(* $Id: SymbolicRegression.m, v 1.0 1995/4/11 RBN $ *)

(*:Version: Mathematica 2.2 *)

(*:Name: SymbolicRegression *)

(*:Title: SymbolicRegression *)

(*:Summary:
Functions for genetic programming of symbolic regression problems.
*)

(*:Keywords: *)

(*:Requirements: No special system requirements. *)

(*:Warning: None. *)

(*:Sources:
John R. Koza, "Genetic Programming: On the Programming of 
Computers by Means of Natural Selection" MIT Press, 1992.
*)

Print["SymbolicRegression.m, version 1.0."] ;

(* set up the package context, included any imports *)

BeginPackage["SymbolicRegression`"]


(* usage messages for the exported functions and the context itself *)

SymbolicRegression::usage = "SymbolicRegression.m is a package that provides
functions for symbolic regression under the genetic programming paradigm."

plus::usage = "plus[a,b,...] represents addition of a, b, ...."

times::usage = "times[a,b,...] represents multiplication of a, b, ...."

subtract::usage = "subtract[a,b] represents the subtraction of b from a."

divide::usage = "divide[a,b] represents the division of a by b."

sqrt::usage = "sqrt[a] represents the square root of a."

log::usage = "log[a] represents the logarithm of a."

exp::usage = "exp[a] represents the exponential of a."

sin::usage = "sin[a] represents the sine of a."

cos::usage = "cos[a] represents the cosine of a."

aif::usage = "aif[t,n,z,p] evaluates t and returns n if
t<0, z if t==0, or p if t>0."

ltz::usage = "ltz[t,lt,ge] evaluates t and returns lt if
t<0 or ge if t=>0."

eqz::usage = "eqz[t,eq,ne] evaluates t and returns eq if
t==0 or ne if t!=0."

SetBig::usage = "SetBig[n] sets the finite representation of infinity for symbolic regression."

PDivide::usage = "PDivide[n,d] is the protected division function which
returns BigNumber if d is zero, n/d otherwise."

PSqrt::usage = "PSqrt[x] is the protected square root function
which returns the square root of the absolute value of x."

PLog::usage = "PLog[x] is the protected natural logarithm function
which returns -BigNumber if x is zero, Log[Abs[x]] otherwise."

PExp::usage = "PExp[x] is the protected exponential function
which returns BigNumber if x is greater than Log[BigNumber],
0 if x is less than -BigNumber, and Exp[x] otherwise."

Eval::usage= "Eval[expr] evaluates expr and returns an expression that can be used
to produce a numerical value."

SumAbs::usage = "SumAbs[{x,y},vars,exprs] computes the standardized fitness for
the list of expressions exprs as a sum of absolute errors for the data in
{x,y}; vars is the list of symbols for the columns of x."

SumSqr::usage = "SumSq[{x,y},vars,exprs] computes the standardized fitness for
the list of expressions exprs as a sum of squared errors for the data in
{x,y}; vars is the list of symbols for the columns of x."

GetData::usage= "GetData[cases, vars] checks that the data in cases is of the
appropriate form and consistant with the variables in vars. It returns a data
object that is suitable for input to a fitness function after any necessary
preprocessing."

(* usage messages for input options *)

(* usage messages for output options *)


(* options for the exported objects *)


Begin["`Private`"]    (* begin the private context *)


(* unprotect any system functions for which rules will be defined *)


(* definition of auxiliary functions and local (static) variables *)

bigInteger = 2^63-1

bigReal = N[bigInteger]

maxExp=Log[bigReal] ;


(* error messages for the exported objects *)

GetData::fitc = "Number of coordinates (`1`) is not equal to
the number of paramters (`2`)."

GetData::fitd = "First argument `1` to GetData is not a list or
rectangular array of data."


(* definition of the exported functions *)

(* associative property *)
plus[a___, b_plus, c___] := 
	Flatten[Unevaluated[plus[a, b, c]], 1, plus]
plus[a_] := a

times[a___, b_times, c___] := 
	Flatten[Unevaluated[times[a, b, c]], 1, times]
times[a_] := a

(* identity property *)
plus[a___, 0, b___] := plus[a, b]
plus[a___, 0., b___] := plus[a, b]
subtract[a___, 0] := a
subtract[a___, 0.] := a

times[a___, 1, b___] := times[a, b]
times[a___, 1., b___] := times[a, b]
divide[a___, 1] := a
divide[a___, 1.] := a

(* inverse property *)
subtract[a_, a_] := 0
divide[a_, a_] := 1

(* annihilation property *)
(*
times[a___, 0, b___] := 0
times[a___, 0., b___] := 0.
divide[0, _] := 0
divide[0., _] := 0.
*)

SetBig[big_Integer]:=
	(bigInteger=big ;
	bigReal=N[big] ;
	maxExp=Log[bigReal] ;
	)
SetBig[big_Real]:=
	(bigInteger=Round[big] ;
	bigReal=big ;
	maxExp=Log[bigReal] ;
	)

PDivide[n_Integer, 0] := bigInteger
PDivide[n_Integer, 0.] := bigInteger
PDivide[n_, 0] := bigInteger
PDivide[n_Real, 0] := bigReal
PDivide[n_Real, 0.] := bigReal
PDivide[n_, 0.] := bigReal
PDivide[n_, d_?NumberQ] := Divide[n, d]

PSqrt[x_Real]:=Sqrt[x]
PSqrt[x_Integer]:=Sqrt[x]
PSqrt[x_?(#<0&)]:=-Sqrt[-x]

PLog[b_,x_]:=PLog[x]/PLog[b]
PLog[x_?(#<0&)]:=Log[-x]
PLog[x_?NumberQ]:=Log[x]
PLog[0]:=-bigInteger
PLog[0.]:=-bigReal

PExp[x_?(#>maxExp&)]:=bigReal
PExp[x_?(#<-maxExp&)]:=0.
PExp[x_?NumberQ]:=Exp[x]

Attributes[aif]={HoldRest} ;

Attributes[ltz]={HoldRest} ;

Attributes[eqz]={HoldRest} ;

exprRules = {
	plus -> Plus, 
	times -> Times, 
	subtract -> Subtract, 
	divide -> PDivide, 
	sqrt -> PSqrt, 
	log -> PLog, 
	exp -> PExp, 
	sin -> Sin, 
	cos -> Cos, 
	aif[t_/;t<0, n_, z_, p_] :> n, 
	aif[0, n_, z_, p_] :> z, 
	aif[0., n_, z_, p_] :> z, 
	aif[t_/;t>0, n_, z_, p_] :> p, 
	ltz[t_/;t<0, lt_, ge_] :> lt, 
	ltz[t_/;t>=0, lt_, ge_] :> ge, 
	eqz[0, eq_, ne_] :> eq, 
	eqz[0., eq_, ne_] :> eq, 
	eqz[t_/;t!=0, eq_, ne_] :> ne
	}


Eval[expr_] := expr /. exprRules

GetData[data_?VectorQ, varlist_] := 
	GetData[Transpose[{Range[Length[data]], data}], varlist]

GetData[data_?MatrixQ, varlist_] := 
	Module[{nc=Length[First@data]-1, nvar=Length[varlist]}, 
		If[nc==nvar,
			{Drop[#,-1]& /@ data, Last /@ data}, 
			Message[GetData::fitc, nc, nvar]; $Failed
			]
		]

GetData[data_, varlist_] := (Message[GetData::fitd,data]; $Failed)

SumAbs[{x_, y_}, vars_, pop_, comp_:True] := 
	Module[{expr, f},
		Check[expr=Eval[#] ;
			f=If[comp, Compile[Evaluate[vars],Evaluate[expr]], 
				Function[Evaluate[vars],Evaluate[expr]]] ;
			Plus @@ Abs[y - Apply[f, x, {1}]] // N, 
			bigReal, General::ovfl]& /@ pop
		]		

SumSqr[{x_, y_}, vars_, pop_, comp_:True] := 
	Module[{f, e},
		Check[f=If[comp, Compile@@{vars,Eval[#]}, 
				Function@@{vars,Eval[#]}] ;
			e = (y - Apply[f, x, {1}]) // N ;
			Dot[e,e],
			bigReal, General::ovfl]& /@ pop
		]		


(* rules for system functions *)

Protect[  ]	(* restore protection of system symbols *)

End[]	(* end the private context *)

Protect[]	(* protect exported symbols *)

EndPackage[]	(* end the package context *)