diff --git a/.gitattributes b/.gitattributes
index 69e400e494cf50de55f77b250d4bd1d7b6accbf6..1c65e282edbd35be2ee64a2a94967c420a0df675 100644
--- a/.gitattributes
+++ b/.gitattributes
@@ -186,10 +186,11 @@ CEP/Calibration/ExpIon/CMakeLists.txt -text
CEP/Calibration/ExpIon/src/CMakeLists.txt -text
CEP/Calibration/ExpIon/src/__init__.py -text
CEP/Calibration/ExpIon/src/acalc.py -text
+CEP/Calibration/ExpIon/src/baselinefitting.cc -text
+CEP/Calibration/ExpIon/src/baselinefitting.py -text
CEP/Calibration/ExpIon/src/calibrate-ion -text
CEP/Calibration/ExpIon/src/client.py -text
CEP/Calibration/ExpIon/src/error.py -text
-CEP/Calibration/ExpIon/src/fitting.cc -text
CEP/Calibration/ExpIon/src/format.py -text
CEP/Calibration/ExpIon/src/io.py -text
CEP/Calibration/ExpIon/src/ionosphere.py -text
diff --git a/CEP/Calibration/ExpIon/src/CMakeLists.txt b/CEP/Calibration/ExpIon/src/CMakeLists.txt
index 8f8090bc197c5ec79be66f4ffa606551ef800019..109a238bcc7224acb4e3620424011d8ac88aac9c 100644
--- a/CEP/Calibration/ExpIon/src/CMakeLists.txt
+++ b/CEP/Calibration/ExpIon/src/CMakeLists.txt
@@ -2,6 +2,20 @@
include(PythonInstall)
+install(PROGRAMS calibrate-ion readms-part.sh readms-part.py parmdbwriter.py DESTINATION bin)
+
+lofar_add_library(_baselinefitting MODULE baselinefitting.cc)
+
+set_target_properties(_baselinefitting PROPERTIES
+ PREFIX ""
+ LIBRARY_OUTPUT_DIRECTORY ${PYTHON_BUILD_DIR}/lofar/expion)
+
+# This is a quick-and-dirty fix to install the Python binding module in the
+# right place. It will now be installed twice, because lofar_add_library()
+# will install it in $prefix/$libdir
+install(TARGETS _baselinefitting
+ DESTINATION ${PYTHON_INSTALL_DIR}/lofar/expion)
+
# Python modules.
python_install(
__init__.py
@@ -15,19 +29,7 @@ python_install(
mpfit.py
readms.py
parmdbmain.py
+ baselinefitting.py
repairGlobaldb.py
DESTINATION lofar/expion)
-install(PROGRAMS calibrate-ion readms-part.sh readms-part.py parmdbwriter.py DESTINATION bin)
-
-lofar_add_library(fitting MODULE fitting.cc)
-
-set_target_properties(fitting PROPERTIES
- PREFIX "")
-# LIBRARY_OUTPUT_DIRECTORY ${PYTHON_BUILD_DIR}/lofar/expion)
-
-# This is a quick-and-dirty fix to install the Python binding module in the
-# right place. It will now be installed twice, because lofar_add_library()
-# will install it in $prefix/$libdir
-install(TARGETS fitting
- DESTINATION ${PYTHON_INSTALL_DIR}/lofar/expion)
diff --git a/CEP/Calibration/ExpIon/src/fitting.cc b/CEP/Calibration/ExpIon/src/baselinefitting.cc
similarity index 72%
rename from CEP/Calibration/ExpIon/src/fitting.cc
rename to CEP/Calibration/ExpIon/src/baselinefitting.cc
index 679d8800b4ecef2d3fbcd4218e9d8ae9488baf40..fc6dede8b5c41e7dcd8f2b761acdf454b5f4761a 100644
--- a/CEP/Calibration/ExpIon/src/fitting.cc
+++ b/CEP/Calibration/ExpIon/src/baselinefitting.cc
@@ -23,7 +23,7 @@
#include <lofar_config.h>
#include <Common/OpenMP.h>
-// #include <omp.h>
+
#include <casa/Containers/ValueHolder.h>
#include <casa/Containers/Record.h>
#include <casa/Utilities/CountedPtr.h>
@@ -46,47 +46,71 @@ namespace LOFAR
namespace ExpIon
{
-ValueHolder fit(const ValueHolder &phases_vh, const ValueHolder &A_vh, const ValueHolder &init_vh = ValueHolder(), const ValueHolder &flags_vh = ValueHolder())
+ValueHolder fit(const ValueHolder &phases_vh, const ValueHolder &A_vh, const ValueHolder &init_vh,
+ const ValueHolder &flags_vh, const ValueHolder &constant_parm_vh )
{
+ // Arrays are passed as ValueHolder
+ // They are Array is extracted by the asArray<Type> methods
+ // They pointer to the actual data is obtained by the Array.data() method.
+
+ // Get the phases
Array<Float> phases = phases_vh.asArrayFloat();
Float *phases_data = phases.data();
- Array<Bool> flags;
- Bool noflags = flags_vh.isNull();
+
+ // Get the flags
+ Array<Bool> flags(flags_vh.asArrayBool());
+ Bool noflags = (flags.ndim() == 0);
+ Bool *flags_data;
if (!noflags)
{
- flags = flags_vh.asArrayBool();
+ flags_data = flags.data();
}
- Bool *flags_data = flags.data();
+
+ IPosition s = phases.shape();
+ int N_station = s[0];
+ int N_freq = s[1];
+
+ // Get the matrix with basis functions
Array<Float> A = A_vh.asArrayFloat();
Float *A_data = A.data();
- IPosition s = phases.shape();
+
IPosition s1 = A.shape();
- int N_station = s[0];
- int N_freq = s[1];
int N_coeff = s1[0];
int N_parm = N_station * N_coeff;
Float sol[N_parm];
- int N_thread = OpenMP::maxThreads();
- LSQFit lnl[N_thread];
- lnl[0] = LSQFit(N_parm);
- lnl[0].setMaxIter(10000);
- uInt nr = 0;
-
- if (!init_vh.isNull())
+ Array<Float> init(init_vh.asArrayFloat());
+ if (init.ndim() == 0)
{
- Array<Float> init = init_vh.asArrayFloat();
- for(int i = 0; i < N_parm; i++) sol[i] = init.data()[i];
+ for(int i = 0; i < N_parm; i++) sol[i] = 0.0;
}
else
{
- for(int i = 0; i < N_parm; i++) sol[i] = 0.0;
+ for(int i = 0; i < N_parm; i++) sol[i] = init.data()[i];
}
- int iter = 0;
- while (!lnl[0].isReady())
+
+ // Get the flags indicating which parameters are constant_parm
+ // i.e. are not a free parameter in the minimization problem
+ Array<Bool> constant_parm(constant_parm_vh.asArrayBool());
+ Bool no_constant_parm = (constant_parm.ndim() == 0);
+ Bool *constant_parm_data;
+ if (!no_constant_parm)
{
- for(int i = 1; i<N_thread; i++) {
+ constant_parm_data = constant_parm.data();
+ }
+
+ Float cEq[N_parm];
+ for(int i=0; i<N_parm; ++i) cEq[i] = 0.0;
+
+ int N_thread = OpenMP::maxThreads();
+ LSQFit lnl[N_thread];
+
+ uInt nr = 0;
+
+ for (int iter = 0; iter<1000; iter++)
+ {
+ for(int i = 0; i<N_thread; i++) {
lnl[i] = LSQFit(N_parm);
}
#pragma omp parallel
@@ -112,13 +136,10 @@ ValueHolder fit(const ValueHolder &phases_vh, const ValueHolder &A_vh, const Val
Float phase_ij_obs = phases_data_k_i - phases_data_k[j];
Float phase_ij_model = 0.0;
- int idx_i = i*N_coeff;
- int idx_j = j*N_coeff;
-
for (int l = 0; l<N_coeff; l++)
{
Float coeff = A_data_k[l];
- phase_ij_model += (sol[idx_i + l] - sol[idx_j + l]) * coeff ;
+ phase_ij_model += (sol[i + l*N_station] - sol[j + l*N_station]) * coeff ;
}
Float sin_dphase, cos_dphase;
sincosf(phase_ij_obs - phase_ij_model, &sin_dphase, &cos_dphase);
@@ -137,8 +158,8 @@ ValueHolder fit(const ValueHolder &phases_vh, const ValueHolder &A_vh, const Val
derivatives_re[N_coeff + l] = -a;
derivatives_im[l] = b;
derivatives_im[N_coeff + l] = -b;
- idx[l] = idx_i+l;
- idx[N_coeff + l] = idx_j+l;
+ idx[l] = i+l*N_station;
+ idx[N_coeff + l] = j+l*N_station;
}
lnl[threadNum].makeNorm(uInt(2*N_coeff), (uInt*) idx, (Float*) derivatives_re, Float(1.0), residual_re);
lnl[threadNum].makeNorm(uInt(2*N_coeff), (uInt*) idx, (Float*) derivatives_im, Float(1.0), residual_im);
@@ -146,44 +167,53 @@ ValueHolder fit(const ValueHolder &phases_vh, const ValueHolder &A_vh, const Val
}
}
}
-// #pragma omp sections
-// {
-// { lnl[0].merge(lnl[1]); }
-// #pragma omp section
-// { lnl[2].merge(lnl[3]); }
-// }
}
-// lnl[0].merge(lnl[2]);
for(int i = 1; i<N_thread; i++)
{
lnl[0].merge(lnl[i]);
}
+
+ Float eq[2];
+ eq[0] = 1.0;
+ eq[1] = -1.0;
+ uInt idx[2];
+
+ if ((!no_constant_parm) )
+ {
+ for (int i = 0; i<N_parm; i++)
+ {
+ if (constant_parm_data[i])
+ {
+ cEq[i] = 1.0;
+ lnl[0].addConstraint( (Float*) cEq, 0.0);
+ cEq[i] = 0.0;
+ }
+ }
+ }
+
if (!lnl[0].solveLoop(nr, sol, True))
{
cout << "Error in loop: " << nr << endl;
break;
}
- iter++;
+ if (lnl[0].isReady())
+ {
+ break;
+ }
}
- cout << iter << endl;
- Array<Float> solutions(IPosition(2, N_coeff, N_station ), sol);
+ Array<Float> solutions(IPosition(2, N_station, N_coeff), sol);
ValueHolder result(solutions);
return result;
};
-ValueHolder fit2(const ValueHolder &phases_vh, const ValueHolder &A_vh)
-{
- return fit(phases_vh, A_vh);
-}
-
} // namespace ExpIon
} // namespace LOFAR
-BOOST_PYTHON_MODULE(fitting)
+BOOST_PYTHON_MODULE(_baselinefitting)
{
casa::pyrap::register_convert_excp();
casa::pyrap::register_convert_basicdata();
@@ -191,5 +221,4 @@ BOOST_PYTHON_MODULE(fitting)
casa::pyrap::register_convert_casa_record();
def("fit", LOFAR::ExpIon::fit);
- def("fit", LOFAR::ExpIon::fit2);
}
diff --git a/CEP/Calibration/ExpIon/src/baselinefitting.py b/CEP/Calibration/ExpIon/src/baselinefitting.py
new file mode 100644
index 0000000000000000000000000000000000000000..4708eea126511836feeeec40920344c173cebd05
--- /dev/null
+++ b/CEP/Calibration/ExpIon/src/baselinefitting.py
@@ -0,0 +1,35 @@
+"""Fit basefunctions to phases using all baselines
+
+Application: Clock-TEC separation
+
+
+The phase can be described by a linear model
+
+.. math::
+ phases = A p
+
+where the columns of matrix A contain the basefunctions and p are the parameters
+
+The same equation holds for the phases of multiple stations
+The dimension of phases is then N_freq x N_station and
+the dimension of p is N_freq x N_param
+
+The cost function that will be minimized is
+
+.. math::
+ \\sum_{i=1}^{N_{station}-1}\\sum_{j=0}^{i-1} \\sum_{k=0}^{N_{freq}} \\| \\exp(\imath(\\Delta phase_{ijk} - \\Delta modelphase_{ijk})) \\|^2
+where
+.. math::
+ \\Delta phase_{ijk} =
+ \\Delta modelphase_{ijk} =
+
+
+
+"""
+
+
+import _baselinefitting
+
+def fit(phases, A, p_0 = None, flags = None, constant_parms = None):
+ """see module description for detailed info"""
+ return _baselinefitting.fit(phases, A, p_0, flags, constant_parms)
\ No newline at end of file