Commit 6c5510bb authored by Andrey Filippov's avatar Andrey Filippov

debugging 1d

parent 71f850b0
...@@ -13,7 +13,7 @@ class Matrix ...@@ -13,7 +13,7 @@ class Matrix
$this->M = $M; $this->M = $M;
if (!is_array($M[0])){ // single-column array if (!is_array($M[0])){ // single-column array
for ($i = 0; $i< sizeof($M); $i++){ for ($i = 0; $i< sizeof($M); $i++){
$this->M[i] = array($M[i]); $this->M[$i] = array($M[$i]);
} }
} }
} }
...@@ -24,9 +24,31 @@ class Matrix ...@@ -24,9 +24,31 @@ class Matrix
} }
public function getColumnPackedCopy(){ public function getColumnPackedCopy(){
$rows = sizeof($this->M);
$cols = sizeof($this->M[0]);
$R = array_fill(0, $rows * $cols, 0.0);
$indx = 0;
for ($i = 0; $i < $rows; $i++){
for ($j = 0; $j < $cols; $j++){
$R[$indx++] = $this->M[$i][$j];
}
}
return $R;
} }
public function getRowPackedCopy(){
$rows = sizeof($this->M);
$cols = sizeof($this->M[0]);
$R = array_fill(0, $rows * $cols, 0.0);
$indx = 0;
for ($i = 0; $i < $rows; $i++){
for ($j = 0; $j < $cols; $j++){
$R[$indx++] = $this->M[$j][$i];
}
}
return $R;
}
public function times($B){ public function times($B){
if ($B instanceof Matrix){ if ($B instanceof Matrix){
$B = $B->get(); $B = $B->get();
...@@ -160,20 +182,19 @@ class Matrix ...@@ -160,20 +182,19 @@ class Matrix
$this->LUPDecompose(); $this->LUPDecompose();
} }
if ($b instanceof Matrix){ if ($b instanceof Matrix){
$b = $b->get(); $b = $b->getColumnPackedCopy();
} }
$N = sizeof($this->M); $N = sizeof($this->M);
$x = array_fill(0,$N,0.0); $x = array_fill(0,$N,0.0);
for ($i = 0; $i < $N; $i++) { for ($i = 0; $i < $N; $i++) {
$x[$i] = $b[$this->P[$i]]; $x[$i] = $b[$this->P[$i]];
for ($k = 0; $k < $i; $k++) for ($k = 0; $k < $i; $k++)
$x[$i] -= $this->M[$i][$k] * $x[$k]; $x[$i] -= $this->A[$i][$k] * $x[$k];
} }
for ($i = $N - 1; $i >= 0; $i--) { for ($i = $N - 1; $i >= 0; $i--) {
for ($k = $i + 1; $k < $N; $k++) for ($k = $i + 1; $k < $N; $k++)
$x[$i] -= $this->M[$i][$k] * $x[$k]; $x[$i] -= $this->A[$i][$k] * $x[$k];
$x[$i] = $x[$i] / $this->M[$i][$i]; $x[$i] = $x[$i] / $this->A[$i][$i];
} }
return new Matrix ($x, $this->Tol); return new Matrix ($x, $this->Tol);
} }
...@@ -257,7 +278,7 @@ class Matrix ...@@ -257,7 +278,7 @@ class Matrix
} }
} }
/*
// test code here // test code here
$A = array( $A = array(
array( 10, -15, 30, 6, -8 ), array( 10, -15, 30, 6, -8 ),
...@@ -298,5 +319,5 @@ print("\nMS nonsingular=".$MS->isNonsingular()."\n"); ...@@ -298,5 +319,5 @@ print("\nMS nonsingular=".$MS->isNonsingular()."\n");
print("\nM1 * IM1=\n"); print("\nM1 * IM1=\n");
$M1->times($IM1)->print(); $M1->times($IM1)->print();
*/
?> ?>
...@@ -10,7 +10,7 @@ class PolynomialApproximation ...@@ -10,7 +10,7 @@ class PolynomialApproximation
public function polynomialApproximation1d($data, $N){ public function polynomialApproximation1d($data, $N){
//$my_array = array_fill(0, $size_of_the_array, $some_data); //$my_array = array_fill(0, $size_of_the_array, $some_data);
if ($this->debugFile === null){ if ($this->debugFile === null){
$this->debugLevel = 0; $this->debugLevel = 1;
} }
$S=array_fill(0,2*$N +1,0); // new double [2*N+1]; $S=array_fill(0,2*$N +1,0); // new double [2*N+1];
$SF=array_fill(0,$N+1,0); // new double [N+1]; $SF=array_fill(0,$N+1,0); // new double [N+1];
...@@ -29,13 +29,11 @@ class PolynomialApproximation ...@@ -29,13 +29,11 @@ class PolynomialApproximation
$wxn *= $x; $wxn *= $x;
} }
$S[2*$N] += $wxn; $S[2*$N] += $wxn;
if ($this->debugLevel > 1){ if ($this->debugLevel > 1) {
if (sizeof($data[$i]) > 2) if (sizeof($data[$i]) > 2)
fprintf($this->debugFile, fprintf($this->debugFile, "polynomialApproximation1d() |" . $i . "|: x=|" . $data[$i][0] . "| f(x)=|" . $data[$i][1] . "| (w=\t|" . $data[$i][2] . "|\n)");
"polynomialApproximation1d() |".$i."|: x=|".$data[$i][0]."| f(x)=|".$data[$i][1]."| (w=\t|".$data[$i][2]."|\t)"); else
else fprintf($this->debugFile, "polynomialApproximation1d() |" . $i, "|: x=|" . $data[$i][0] . "| f(x)=|" . $data[$i][1] . "|\n)");
fprintf($this->debugFile,
"polynomialApproximation1d() |".$i,"|: x=|".$data[$i][0]."| f(x)=|".$data[$i][1]."|\t)");
} }
} }
} }
...@@ -51,18 +49,16 @@ class PolynomialApproximation ...@@ -51,18 +49,16 @@ class PolynomialApproximation
if ($this->debugLevel > 1){ if ($this->debugLevel > 1){
fprintf($this->debugFile, fprintf($this->debugFile,
"polynomialApproximation1d(data,".$N.") M:\n"); "\npolynomialApproximation1d(data,".$N.") M:\n");
$M->print($this->debugFile); $M->print($this->debugFile);
fprintf($this->debugFile, fprintf($this->debugFile,
"polynomialApproximation1d() B:\n"); "\npolynomialApproximation1d() B:\n");
$B->print($this->debugFile); $B->print($this->debugFile);
} }
// while (!(new LUDecomposition(M)).isNonsingular() && (N1>0)){ // while (!(new LUDecomposition(M)).isNonsingular() && (N1>0)){
// $df = ($this->debugLevel > 2) ? $this->debugFile : null; // $df = ($this->debugLevel > 2) ? $this->debugFile : null;
while ($N1 >= 0) { // make N=0 legal ? while ($N1 >= 0) { // make N=0 legal ?
try { try {
// $M_inv = $M.invert();
// $R = $M_inv->times($B);
$R = $M->solve($B); $R = $M->solve($B);
if ($this->debugLevel > 1){ if ($this->debugLevel > 1){
...@@ -70,10 +66,7 @@ class PolynomialApproximation ...@@ -70,10 +66,7 @@ class PolynomialApproximation
"polynomialApproximation1d() solution=\n"); "polynomialApproximation1d() solution=\n");
$R->print($this->debugFile); $R->print($this->debugFile);
} }
$result=array_fill(0,$N+1,0); // new double [N+1]; $result = $R->getColumnPackedCopy();
for ($i = 0; $i <= $N1; $i++){
$result[$i] = $R[$i][0];
}
return $result; return $result;
} catch (Exception $e) { } catch (Exception $e) {
......
...@@ -246,7 +246,7 @@ function findMax1d($data, $fract, $pow, $frac_outside = 0.2){ ...@@ -246,7 +246,7 @@ function findMax1d($data, $fract, $pow, $frac_outside = 0.2){
} }
$pa = new PolynomialApproximation(); $pa = new PolynomialApproximation();
$pa->debugLevel = 3; $pa->debugLevel = 1;
$pa->debugFile = $dbg_file; $pa->debugFile = $dbg_file;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment