Skip to content

M
motosat
Commit 3ed6ebc5 authored by Andrey Filippov's avatar Andrey Filippov
Browse files
Options

debugging 2D peaking

parent 56ca0c83
Hide whitespace changes
Inline Side-by-side
Showing with 236 additions and 46 deletions
Matrix.php
View file @ 3ed6ebc5
......@@ -24,7 +24,7 @@ class Matrix
}
public function set($i,$j,$v){
$M[$i][$j] = $v;
$this->M[$i][$j] = $v;
}
......@@ -46,9 +46,10 @@ class Matrix
$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];
// var_dump($this->M);
for ($j = 0; $j < $cols; $j++){
for ($i = 0; $i < $rows; $i++){
$R[$indx++] = $this->M[$i][$j];
}
}
return $R;
......
PolynomialApproximation.php
View file @ 3ed6ebc5
......@@ -89,6 +89,42 @@ class PolynomialApproximation
}
return null;
}
public function quadraticMax2d($data, $thresholdQuad = 1.0E-15, $debugLevel = 1)
{
$coeff = $this->quadraticApproximation($data, false, null, 1.0E-20, $thresholdQuad, $debugLevel);
if ($coeff == null)
return null;
if (sizeof($coeff[0]) < 6)
return null;
// double [][] aM={
// {2*coeff[0][0], coeff[0][2]}, // | 2A, C |
// { coeff[0][2],2*coeff[0][1]}}; // | C, 2B |
$aM = array(
array(2 * $coeff[0][0], $coeff[0][2]), // | 2A, C |)
array($coeff[0][2], 2 * $coeff[0][1]) // | C, 2B |
);
$M = new Matrix($aM);
$nmQ = PolynomialApproximation::normMatix($aM);
// if (debugLevel>3) System.out.println("M.det()="+M.det()+" PolynomialApproximation::normMatix(aM)="+nmQ+" data.length="+data.length);
if (($nmQ == 0.0) || (abs($M -> det()) / $nmQ < $thresholdQuad)) {
// if (debugLevel>3) System.out.println("quadraticMax2d() failed: M.det()="+M.det()+" PolynomialApproximation::normMatix(aM)="+PolynomialApproximation::normMatix(aM));
return null;
}
// double [][] aB={
// {-coeff[0][3]}, // | - D |
// {-coeff[0][4]}}; // | - E |
// $xy= $M.solve(new Matrix(aB)).getColumnPackedCopy();
$aB = array(
- $coeff[0][3], // | - D |
- $coeff[0][4] // | - E |
);
$mxy = $M -> solve($aB);
$xy = $mxy->getColumnPackedCopy();
return $xy;
}
public function quadraticApproximation(
$data,
......@@ -102,7 +138,7 @@ class PolynomialApproximation
$this->debugLevel = 0;
}
if ((data == null) || (data.length == 0)) {
if (($data == null) || (sizeof($data) == 0)) {
return null;
}
/* ix, iy - the location of the point with maximal value. We'll approximate the vicinity of that maximum using a
......@@ -257,16 +293,16 @@ class PolynomialApproximation
$mLin=new Matrix ($mAarrayL);
if ($mDampingLin !== null){
if (isset($mDampingLin)){
$mLin->plusEquals($mDampingLin);
}
// TODO Maybe bypass determinant checks for damped ?
// if (debugLevel>3) System.out.println(">>> n="+n+" det_lin="+mLin.det()+" norm_lin="+normMatix(mAarrayL));
$nmL = normMatix($mAarrayL);
if (($nmL == 0.0) || (abs($mLin.det()) / $nmL < $thresholdLin)){
// if (debugLevel>3) System.out.println(">>> n="+n+" det_lin="+mLin.det()+" norm_lin="+PolynomialApproximation::normMatix(mAarrayL));
$nmL = PolynomialApproximation::normMatix($mAarrayL);
if (($nmL == 0.0) || (abs($mLin->det()) / $nmL < $thresholdLin)){
// return average value for each channel
if ($S00 == 0.0) return null; // not even average
$ABCDEF = Matrix::ZeroMatrix(zDim, 3);
$ABCDEF = Matrix::ZeroMatrix($zDim, 3);
for ($i = 0; $i < $zDim; $i++) {
$ABCDEF[$i][0] = 0.0;
$ABCDEF[$i][1] = 0.0;
......@@ -281,9 +317,9 @@ class PolynomialApproximation
$zAarrayL[1]=$SZ01[$i];
$zAarrayL[2]=$SZ00[$i];
$Z = new Matrix ($zAarrayL); // ,3);
$ABCDEF[$i]= $mLin.solve($Z).getRowPackedCopy();
$ABCDEF[$i]= $mLin->solve($Z)->getRowPackedCopy();
}
if (forceLinear) return ABCDEF;
if ($forceLinear) return $ABCDEF;
// quote try quadratic approximation
$mAarrayQ = array(
array($S40, $S22, $S31, $S30, $S21, $S20),
......@@ -298,16 +334,16 @@ class PolynomialApproximation
}
// if (debugLevel>3) {
// System.out.println(" n="+n+" det_quad="+mQuad.det()+" norm_quad="+normMatix(mAarrayQ)+" data.length="+data.length);
// System.out.println(" n="+n+" det_quad="+mQuad.det()+" norm_quad="+PolynomialApproximation::normMatix(mAarrayQ)+" data.length="+data.length);
// mQuad.print(10,5);
// }
$nmQ = normMatix($mAarrayQ);
if (($nmQ == 0.0) || (abs($mQuad.det())/normMatix($mAarrayQ) < $thresholdQuad)) {
$nmQ = PolynomialApproximation::normMatix($mAarrayQ);
if (($nmQ == 0.0) || (abs($mQuad->det())/$nmQ < $thresholdQuad)) {
// if (debugLevel>0) System.out.println("Using linear approximation, M.det()="+mQuad.det()+
// " normMatix(mAarrayQ)="+normMatix(mAarrayQ)+
// " PolynomialApproximation::normMatix(mAarrayQ)="+PolynomialApproximation::normMatix(mAarrayQ)+
// ", thresholdQuad="+thresholdQuad+
// ", nmQ="+nmQ+
// ", Math.abs(M.det())/normMatix(mAarrayQ)="+(Math.abs(mQuad.det())/normMatix(mAarrayQ))); //did not happen
// ", Math.abs(M.det())/PolynomialApproximation::normMatix(mAarrayQ)="+(Math.abs(mQuad.det())/PolynomialApproximation::normMatix(mAarrayQ))); //did not happen
return $ABCDEF; // not enough data for the quadratic approximation, return linear
}
// double [] zAarrayQ={SZ20,SZ02,SZ11,SZ10,SZ01,SZ00};
......@@ -320,7 +356,7 @@ class PolynomialApproximation
$zAarrayQ[4] = $SZ01[$i];
$zAarrayQ[5] = $SZ00[$i];
$Z = new Matrix ($zAarrayQ); // ,6);
$ABCDEF[i]= $mQuad.solve($Z).getRowPackedCopy();
$ABCDEF[$i]= $mQuad->solve($Z)->getRowPackedCopy();
}
return $ABCDEF;
}
......@@ -330,7 +366,7 @@ class PolynomialApproximation
// calcualte "volume" made of the matrix row-vectors, placed orthogonally
// to be compared to determinant
public function normMatix($a) {
public static function normMatix($a) {
$norm=1.0;
for ($i=0; $i<sizeof($a); $i++) {
$d=0;
......
test_01.php
View file @ 3ed6ebc5
......@@ -55,10 +55,29 @@ $dbg_file = fopen("/home/eyesis/git/motosat/attic/logs/test04.log","w");
fprintf($dbg_file,"test log\n");
if (true) {
$data = Array(5,6,7,9,9,10,11,11,12,12,11,11,10,10,10,36,40,45,47,51,51,51,49,48,43,38,30,30,30,10,10,10,10);
$rslt = findMax1d($data, $USABLE_FRACT, $USABLE_POW, $ARGMAX_OURSIDE);
$arr = Array(5,6,7,9,9,10,11,11,12,12,11,11,10,10,10,36,40,45,47,51,51,51,49,48,43,38,30,30,30,10,10,10,10);
$rslt = findMax1d($arr, $USABLE_FRACT, $USABLE_POW, $ARGMAX_OURSIDE);
print_r($rslt);
print("\n");
// Try 2D:
$arr2d = array();
$rows = sizeof($arr);
$cols = $rows;
for ($i= 0; $i < $rows; $i++){
$line = $arr;
for ($j= 0; $j < $cols; $j++){
$line[$j] *= $arr[$i]/51;
$line[$j] += rand(-10,10);
if ($line[$j] <1){
$line[$j] = 1;
}
}
$arr2d [] = $line;
}
// $rslt = findMax2d($arr2d, $USABLE_FRACT, $USABLE_POW, $ARGMAX_OURSIDE);
$rslt = findMax2d($arr2d, 0.55, $USABLE_POW, $ARGMAX_OURSIDE);
var_dump($rslt);
exit (0);
}
......@@ -212,6 +231,165 @@ function peakAzimuthOrElevation($xml_state, $el_not_az, $azel_play, $azel_range,
}
}
function findMax2d($data, $fract, $pow, $frac_outside = 0.2){ // $pow not yet used
global $dbg_file;
$min= $data[0][0];
$max = $data[0][0];
$rows = sizeof($data);
$cols = sizeof($data[0]); // should be rectangular
$argmax = array(0.0, 0.0);
foreach ($data as $row => $line) {
foreach ($line as $col => $v) {
if ($v > $max) {
$max = $v;
$argmax = array($col, $row,);
} else if ($v < $min)
$min = $v;
}
}
$threshold = $min + ($max-$min)*$fract;
// Find a 2D cluster around $argmax exceeding $threshold
$cluster = $data;
for ($i = 0; $i < $rows; $i++){
for ($j = 0; $j < $cols; $j++){
$cluster[$i][$j] = false;
}
}
// Wve method (1 directions)
$poly_data = array();
$queue = array();
$xy = $argmax; // x,y pair
// put a cell into a queue, in $poly_data and mark a cell as used
$queue[] = $xy;
$poly_data[] = array(
array(0.0,0.0), // $xy,
array(
$data[$xy[1]][$xy[0]]
)
);
$cluster[$xy[1]][$xy[0]] = true;
while (sizeof($queue) > 0) {
$xy0 = $queue[0];
$queue = array_slice($queue, 1);
// print("<=" . print_r($xy0, 1));
for ($dir = 0; $dir < 4; $dir ++) {
$xy = $xy0;
switch ($dir) {
case 0:
$xy[0] ++;
if ($xy[0] >= $cols) {
continue;
}
// print($dir . ":" . print_r($xy, 1));
break;
case 1:
$xy[1] ++;
if ($xy[1] >= $rows) {
continue;
}
// print($dir . ":" . print_r($xy, 1));
break;
case 2:
$xy[0] --;
if ($xy[0] < 0) {
continue;
}
// print($dir . ":" . print_r($xy, 1));
break;
case 3:
$xy[1] --;
if ($xy[1] < 0) {
continue;
}
// print($dir . ":" . print_r($xy, 1));
break;
}
if ($cluster[$xy[1]][$xy[0]])
continue; // already used
if ($data[$xy[1]][$xy[0]] < $threshold) // signal too weak
continue; // already used
$queue[] = $xy;
$poly_data[] = array(
array($xy[0]-$argmax[0], $xy[1]-$argmax[1]),
array(
$data[$xy[1]][$xy[0]]
)
);
$cluster[$xy[1]][$xy[0]] = true;
}
}
print ("*** sizeof(poly_data)=".sizeof($poly_data)."\n");
if (true) {
printf("argmax= (x=%d ,y=%d), max = %f, threshold = %f\n",
$argmax[0], $argmax[1], $max, $threshold);
for ($i = 0; $i < $rows; $i ++) {
for ($j = 0; $j < $cols; $j ++) {
if (($i==$argmax[1]) && ($j==$argmax[0])) {
printf("[%02d]", $data[$i][$j]);
} else if ($cluster[$i][$j] > 0) {
printf("<%02d>", $data[$i][$j]);
} else {
printf(" %02d ", $data[$i][$j]);
}
}
printf("\n");
}
}
$pa = new PolynomialApproximation();
$pa->debugLevel = 1;
$pa->debugFile = $dbg_file;
$rslt = $pa->quadraticMax2d($poly_data);
if ($rslt == null){
return null;
}
$rslt[0] += $argmax[0];
$rslt[1] += $argmax[1];
if ($dbg_file) {
fprintf($dbg_file, "argmax= (x=%d ,y=%d), max = %f, threshold = %f, rslt:\n",
$argmax[0], $argmax[1], $max, $threshold);
fprintf($dbg_file, print_r($rslt,1));
}
if (is_nan($rslt[0]) || is_nan($rslt[1])){
if ($dbg_file) {
fprintf($dbg_file, "Could not find maximum (NaN): coeff = \n".print_r($rslt,1)."\n");
}
return null; // no maximum
}
if ($rslt[0] < -($frac_outside * $cols)){
if ($dbg_file) {
fprintf($dbg_file, "argmax[0] is too low = %f, range=%d\n",$rslt[0],$cols );
}
return null; // az too low
}
if ($rslt[0] > ((1.0 + $frac_outside) * $cols)){
if ($dbg_file) {
fprintf($dbg_file, "argmax[0] is too high = %f, range=%d\n",$rslt[0],$cols );
}
return null; // too low
}
if ($rslt[1] < -($frac_outside * $rows)){
if ($dbg_file) {
fprintf($dbg_file, "argmax[1] is too low = %f, range=%d\n",$rslt[1],$rows );
}
return null; // az too low
}
if ($rslt[1] > ((1.0 + $frac_outside) * $rows)){
if ($dbg_file) {
fprintf($dbg_file, "argmax[1] is too high = %f, range=%d\n",$rslt[1],$cols );
}
return null; // too low
}
return $rslt; // $rslt;
}
function findMax1d($data, $fract, $pow, $frac_outside = 0.2){
......@@ -287,31 +465,6 @@ function findMax1d($data, $fract, $pow, $frac_outside = 0.2){
}
return null; // no maximum
}
/*
// debugging
$rslt[] = $argmax;
$rslt[] = $argmax - $rslt[1] / (2 * $rslt[2]);
$pre = array_fill(0, sizeof($data),0);
$diff= array_fill(0, sizeof($data),0);
for ($k = $k_min; $k <= $k_max; $k++){
$x = $k-$argmax;
$pre[$k]=$rslt[0]+$rslt[1]*$x+$rslt[2]*$x*$x;
$diff[$k]=$pre[$k]-$datap[$k];
}
for ($k = $k_min; $k <= $k_max; $k++){
printf("[%02d]: %8.3f %8.3f %8.3f %8.3f\n", $k, $data[$k], $datap[$k], $pre[$k], $diff[$k]);
}
print("data\n");
print_r($data);
print("datap\n");
print_r($datap);
print("---datap\n");
// print_r($pre);
// print_r($diff);
*/
// $frac_outside = 0.2;
$rel_argmax = $argmax - $rslt[1] / (2 * $rslt[2]);
if ($rel_argmax < -($frac_outside * sizeof($data))){
if ($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