/php/rk4.php

<?php
	function cmp_key_len($a,$b) {
		if(strlen($a) > strlen($b)) {
			return -1;
		} else if(strlen($a) < strlen($b)) {
			return 1;	
		} else {
			return 0;	
		}
	}
	
	class rk4 {
		
		var $h; //step size
		var $start_t; //start time
		var $end_t; //end time
		
		var $odes; //simultaneous odes
		var $params; //parameter values
		var $ics; //initial conditions
		var $outputs; //outputs
		
		var $inputs; //input array, heavily modified from js for processing purposes
		
		var $inputs_raw; //unmodified inputs, for possible future use, might be needed for optimization

		
		var $ode_system = array(); //odes with param values replaced and initial conditions set
		var $data;
		
		function rk4($step_size, $end_time, $start_time = 0) {
			$this->h = $step_size;
			$this->start_T = $start_time;
			$this->end_t = $end_time;
		}
		
		function set_system($odeArr) {
			$this->odes = $odeArr;
		}
		
		function set_params($paramArr){
			$this->params = $paramArr;
		}
		
		function set_ics($icArr) {
			$this->ics = $icArr;
		}
		
		function set_outputs($outArr) {
			$this->outputs = $outArr;	
		}
		
		/**
		*	Input array format:
		*	$inputs['q1'][0]['type'] 		//bolus or infusion, currently only bolus implemented
		*					['start'] 		//start time
		*					['end']			//end time (ignored for bolus)
		*					['repeat']		//# of times to repeat (0 for none). repeats will only be processed until sim end tiem
		*					['interval']	//interval between repeats.  same unit of time as simulation
		*					['magnitude']	//magnitude of event
		*
		*
		*	Output array format:
		*	$inputs[<event time step>][0]['magnitude']	//magnitude of event
		*								['sv']			//state var
		*	
		*	repeat events will have multiple entries at the closest time step point.  constant infusion will have entries in range
		*/
		
		function set_inputs($inputArr) {
			$this->inputs_raw = $inputArr;
			foreach($inputArr as $sv => $inputs) {
				foreach($inputs as $input) {
					switch($input['type']) {
						case "bolus":
							$start = floor(1.0 * $input['start'] / $this->h);
							
							
							if($input['repeat'] > 0) {
								$end_time = min(($input['repeat'] * $input['interval'] * 1.0),$this->end_t);
								$num_entries = floor(1.0 * $end_time / $input['interval']);
								$interval = floor(1.0 * $input['interval'] / $this->h);
								$end = ceil(1.0 * $end_time / $this->h);
								
								for($i = $start; $i < $end; $i += $interval) {
									$this->inputs[$i][$sv][] = $input['magnitude'];
								}
							} else {
								$this->inputs[$start][$sv][] = $input['magnitude'];
							}
						break;
						case "infusion":
						
						break;
					}
				}
			}
		}
		
		function init_system() {
			uksort($this->params,'cmp_key_len');
			//print_r($this->params);
			//echo 'after sort';
			//print_r($this->params);
			foreach($this->odes as $sv=>$ode) {
				foreach($this->params as $param => $val) {
//					echo $param . ' ' . $val . ' ' . $ode . ' ';
					$ode = str_replace($param,$val,$ode);	
					
				}
				
				$this->ode_system[$sv] = $ode;
				
			}
			foreach($this->ics as $sv => $val) {
				$this->data[0][$sv] = ($val * 1.0) + array_sum($this->inputs[0][$sv]);	
			}
			
			//print_r($this->ode_system);
			//$this->_a($this->ode_system['q1'],1);
		}
		
		function calculate_outputs() {
			
			foreach($this->data as $i => $y_vec) { //for each time point
				foreach($this->outputs as $out => $eqn) {  //for each output
					$this->data[$i]['y' . $out] = $this->f($eqn,$this->i2t($i),$this->data[$i]);
				}
			}
		}
		
		function format_data($format = "js") {

			if($format == "js" ) {
				$js = array();			
				foreach($this->data as $x => $y_vec) {
					foreach($y_vec as $sv => $y) {
						if($sv != "t") {
							$js[$sv][] = array($this->i2t($x),$y);
						}
					}
				}
				return $js;
			}else if($format == "matlab") {
				//matlab array output
				$matlab = "";
				return $matlab;
			}
			else if($format == "gnu") {
				$gnu = "";
//print_r($this->data);
				foreach($this->data as $x =>$y_vec) {
					$gnu .= $this->data[$x]['t'] . "\t";
					foreach($y_vec as $y => $val) {
						$gnu .= $this->data[$x][$y] . "\t";
					}
					$gnu .="\n";
					 
				}
				//print_r($y_vec);
				return $gnu;
			}

		}
		
		function simulate() {
			$this->init_system();
			if(count($this->ode_system) != count($this->ics) ) {
				return false;	
			}
			
			$test = array();
			
			$num_steps = ($this->end_t - $this->start_t) / $this->h;
			
			//data already contains our initial conditions, so we start at $i = 1
			for($i = 0; $i < $num_steps; $i ++) { //main step loop
				$this->rk_iter($i);
				$this->data[$i]['t'] = $this->i2t($i);
			}
			$this->calculate_outputs();
		//	print_r($this->data);
			
		}
		
		//convert from step to time
		function i2t($i) {
			return ($i * $this->h * 1.0);
		}
		
		//ode, t, x (x = data[sv][t])
		//t is previous step, current t = t + h
		function f($ode,$t,$x_vec) {
			foreach(array_keys($this->ics) as $sv) {
				$ode = str_replace($sv,$x_vec[$sv],$ode);	
				//$ode = str_replace("TIME",$t,$ode);
			}	
			//echo 'echo '. $t . " " .  $ode . ";\n\n\n";
  	 	 	eval("\$result = $ode;");
			//echo $t . " " . $ode . "=" . $result . "\n";
			return $result;
		}

		
		//i = previous iteration
		function _k1($i) {
			$t = $this->i2t($i);
			$k1 = array();
			foreach($this->ode_system as $sv=>$ode) {
				/*foreach($this->ics as $sv=>$ic) {
					$ode = str_replace($sv,$this->data[$sv][$t],$ode);
				}*/
				$k1[$sv] = $this->f($ode,$t,$this->data[$i]);
			}
			return $k1;
		}
		
		function _k2($k1,$i) {
			$t = $this->i2t($i);
			//print_r($k1);
			foreach($this->ode_system as $sv=>$ode) {
				
				//echo $i . " ".  $this->data[$i][$sv] . " ";
				$x_a[$sv] = $this->data[$i][$sv] + (($this->h/2.0) * $k1[$sv]);
				//echo "x_a[$sv] " . $x_a[$sv] . " "; 
			}
			//echo "\n";
			foreach($this->ode_system as $sv=>$ode) {
				$k2[$sv] = $this->f($ode,($t+($this->h/2.0)),$x_a);
			}
			return $k2;
		}
		
		function _k3($k2,$i) {
			$t = $this->i2t($i);
			//print_r($k2);
			foreach($this->ode_system as $sv=>$ode) {
				$x_b[$sv] = $this->data[$i][$sv]	+ (($this->h/2.0) * $k2[$sv]);
			}
			//print_r($x_b);
			foreach($this->ode_system as $sv=>$ode) {
				$k3[$sv] = $this->f($ode,($t+($this->h/2.0)),$x_b);
			}
			return $k3;
		}
		function _k4($k3,$i) {
			$t = $this->i2t($i);
			//print_r($k3);
			foreach($this->ode_system as $sv=>$ode) {
				$x_c[$sv] = $this->data[$i][$sv]	+ ($this->h * $k3[$sv]);
			}
			//print_r($x_c);
			foreach($this->ode_system as $sv=>$ode) {
				$k4[$sv] = $this->f($ode,($t+$this->h),$x_c);
			}
			return $k4;	
		}
		
		function rk_iter($i) {
			//print_r($k4);
			$t = $this->i2t($i);
			$k1 = $this->_k1($i);
			$k2 = $this->_k2($k1,$i);
			$k3 = $this->_k3($k2,$i);
			$k4 = $this->_k4($k3,$i);
			/*echo 'k1: ';
			print_r($k1);
			echo ' k2: ';
			print_r($k2);
			echo ' k3: ';
			print_r($k3);
			echo ' k4: ';
			print_r($k4);*/
			foreach($this->ode_system as $sv=>$ode) {
				$this->data[$i+1][$sv] = $this->data[$i][$sv] + ($this->h/6.0)*($k1[$sv] + (2.0 * $k2[$sv]) + (2.0 * $k3[$sv]) + $k4[$sv]);
				//foreach($this->inputs[$i+1][$sv] as $inp) {
				$this->data[$i+1][$sv] += array_sum($this->inputs[$i+1][$sv]);
				//}
			}
			//return $x1;
		}
	}

?>