Open Build Service

Changes of Revision 9

atmos.spec Changed

atmos-2.96.9.tar.gz/configure -> atmos-2.97.1.tar.gz/configure Changed

@@ -1,6 +1,6 @@
 #! /bin/sh
 # Guess values for system-dependent variables and create Makefiles.
-# Generated by GNU Autoconf 2.63 for atmos 2.96.9.
+# Generated by GNU Autoconf 2.63 for atmos 2.97.1.
 #
 # Report bugs to <mass-general@curl.sai.msu.ru>.
 #
@@ -596,8 +596,8 @@
 # Identity of this package.
 PACKAGE_NAME='atmos'
 PACKAGE_TARNAME='atmos'
-PACKAGE_VERSION='2.96.9'
-PACKAGE_STRING='atmos 2.96.9'
+PACKAGE_VERSION='2.97.1'
+PACKAGE_STRING='atmos 2.97.1'
 PACKAGE_BUGREPORT='mass-general@curl.sai.msu.ru'
 
 ac_unique_file="src/atmos.h"
@@ -1293,7 +1293,7 @@
   # Omit some internal or obsolete options to make the list less imposing.
   # This message is too long to be a string in the A/UX 3.1 sh.
   cat <<_ACEOF
-\`configure' configures atmos 2.96.9 to adapt to many kinds of systems.
+\`configure' configures atmos 2.97.1 to adapt to many kinds of systems.
 
 Usage: $0 OPTION... VAR=VALUE...
 
@@ -1359,7 +1359,7 @@
 
 if test -n "$ac_init_help"; then
   case $ac_init_help in
-     short | recursive ) echo "Configuration of atmos 2.96.9:";;
+     short | recursive ) echo "Configuration of atmos 2.97.1:";;
    esac
   cat <<\_ACEOF
 
@@ -1449,7 +1449,7 @@
 test -n "$ac_init_help" && exit $ac_status
 if $ac_init_version; then
   cat <<\_ACEOF
-atmos configure 2.96.9
+atmos configure 2.97.1
 generated by GNU Autoconf 2.63
 
 Copyright (C) 1992, 1993, 1994, 1995, 1996, 1998, 1999, 2000, 2001,
@@ -1463,7 +1463,7 @@
 This file contains any messages produced by compilers while
 running configure, to aid debugging if configure makes a mistake.
 
-It was created by atmos $as_me 2.96.9, which was
+It was created by atmos $as_me 2.97.1, which was
 generated by GNU Autoconf 2.63.  Invocation command line was
 
   $ $0 $@
@@ -2312,7 +2312,7 @@
 
 # Define the identity of the package.
  PACKAGE='atmos'
- VERSION='2.96.9'
+ VERSION='2.97.1'
 
 
 cat >>confdefs.h <<_ACEOF
@@ -3978,11 +3978,15 @@
 
 
 
-for ac_header in boost/numeric/ublas/io.hpp \
+
+
+for ac_header in boost/numeric/ublas/banded.hpp \
+	boost/numeric/ublas/matrix_expression.hpp \
 	boost/numeric/ublas/matrix.hpp \
 	boost/numeric/ublas/matrix_proxy.hpp \
 	boost/numeric/ublas/symmetric.hpp \
 	boost/numeric/ublas/triangular.hpp \
+	boost/numeric/ublas/vector_expression.hpp \
 	boost/numeric/ublas/vector.hpp \
 	boost/numeric/ublas/vector_proxy.hpp
 do
@@ -4152,6 +4156,10 @@
 
 
 
+
+
+
+
 for ac_header in boost/algorithm/string.hpp \
 	boost/bind.hpp \
 	boost/cstdint.hpp \
@@ -4161,10 +4169,14 @@
 	boost/format.hpp \
 	boost/function.hpp \
 	boost/io/ios_state.hpp \
-	boost/iterator.hpp \
+	boost/iterator/transform_iterator.hpp \
 	boost/lexical_cast.hpp \
+	boost/math/tools/roots.hpp \
+	boost/mem_fn.hpp \
+	boost/mpl/if.hpp \
 	boost/program_options.hpp \
-	boost/tuple/tuple.hpp
+	boost/tuple/tuple.hpp \
+	boost/type_traits/is_complex.hpp
 do
 as_ac_Header=`$as_echo "ac_cv_header_$ac_header" | $as_tr_sh`
 if { as_var=$as_ac_Header; eval "test \"\${$as_var+set}\" = set"; }; then
@@ -4321,160 +4333,6 @@
 
 
 
-for ac_header in boost/math/tools/roots.hpp \
-	tr1/tuple
-do
-as_ac_Header=`$as_echo "ac_cv_header_$ac_header" | $as_tr_sh`
-if { as_var=$as_ac_Header; eval "test \"\${$as_var+set}\" = set"; }; then
-  { $as_echo "$as_me:$LINENO: checking for $ac_header" >&5
-$as_echo_n "checking for $ac_header... " >&6; }
-if { as_var=$as_ac_Header; eval "test \"\${$as_var+set}\" = set"; }; then
-  $as_echo_n "(cached) " >&6
-fi
-ac_res=`eval 'as_val=${'$as_ac_Header'}
-		 $as_echo "$as_val"'`
-	       { $as_echo "$as_me:$LINENO: result: $ac_res" >&5
-$as_echo "$ac_res" >&6; }
-else
-  # Is the header compilable?
-{ $as_echo "$as_me:$LINENO: checking $ac_header usability" >&5
-$as_echo_n "checking $ac_header usability... " >&6; }
-cat >conftest.$ac_ext <<_ACEOF
-/* confdefs.h.  */
-_ACEOF
-cat confdefs.h >>conftest.$ac_ext
-cat >>conftest.$ac_ext <<_ACEOF
-/* end confdefs.h.  */
-$ac_includes_default
-#include <$ac_header>
-_ACEOF
-rm -f conftest.$ac_objext
-if { (ac_try="$ac_compile"
-case "(($ac_try" in
-  *\"* | *\`* | *\\*) ac_try_echo=\$ac_try;;
-  *) ac_try_echo=$ac_try;;
-esac
-eval ac_try_echo="\"\$as_me:$LINENO: $ac_try_echo\""
-$as_echo "$ac_try_echo") >&5
-  (eval "$ac_compile") 2>conftest.er1
-  ac_status=$?
-  grep -v '^ *+' conftest.er1 >conftest.err
-  rm -f conftest.er1
-  cat conftest.err >&5
-  $as_echo "$as_me:$LINENO: \$? = $ac_status" >&5
-  (exit $ac_status); } && {
-	 test -z "$ac_cxx_werror_flag" ||
-	 test ! -s conftest.err
-       } && test -s conftest.$ac_objext; then
-  ac_header_compiler=yes
-else
-  $as_echo "$as_me: failed program was:" >&5
-sed 's/^/| /' conftest.$ac_ext >&5
-
-	ac_header_compiler=no
-fi
-
-rm -f core conftest.err conftest.$ac_objext conftest.$ac_ext
-{ $as_echo "$as_me:$LINENO: result: $ac_header_compiler" >&5
-$as_echo "$ac_header_compiler" >&6; }
-
-# Is the header present?
-{ $as_echo "$as_me:$LINENO: checking $ac_header presence" >&5
-$as_echo_n "checking $ac_header presence... " >&6; }
-cat >conftest.$ac_ext <<_ACEOF
-/* confdefs.h.  */
-_ACEOF
-cat confdefs.h >>conftest.$ac_ext
-cat >>conftest.$ac_ext <<_ACEOF
-/* end confdefs.h.  */
-#include <$ac_header>
-_ACEOF
-if { (ac_try="$ac_cpp conftest.$ac_ext"
-case "(($ac_try" in
-  *\"* | *\`* | *\\*) ac_try_echo=\$ac_try;;
-  *) ac_try_echo=$ac_try;;
-esac
-eval ac_try_echo="\"\$as_me:$LINENO: $ac_try_echo\""
-$as_echo "$ac_try_echo") >&5
-  (eval "$ac_cpp conftest.$ac_ext") 2>conftest.er1
-  ac_status=$?
-  grep -v '^ *+' conftest.er1 >conftest.err
-  rm -f conftest.er1
-  cat conftest.err >&5
-  $as_echo "$as_me:$LINENO: \$? = $ac_status" >&5
-  (exit $ac_status); } >/dev/null && {
-	 test -z "$ac_cxx_preproc_warn_flag$ac_cxx_werror_flag" ||
-	 test ! -s conftest.err
-       }; then
-  ac_header_preproc=yes
-else
-  $as_echo "$as_me: failed program was:" >&5
-sed 's/^/| /' conftest.$ac_ext >&5
-
-  ac_header_preproc=no
-fi
-
-rm -f conftest.err conftest.$ac_ext
-{ $as_echo "$as_me:$LINENO: result: $ac_header_preproc" >&5
-$as_echo "$ac_header_preproc" >&6; }
-
-# So?  What about this header?
-case $ac_header_compiler:$ac_header_preproc:$ac_cxx_preproc_warn_flag in
-  yes:no: )
-    { $as_echo "$as_me:$LINENO: WARNING: $ac_header: accepted by the compiler, rejected by the preprocessor!" >&5
-$as_echo "$as_me: WARNING: $ac_header: accepted by the compiler, rejected by the preprocessor!" >&2;}
-    { $as_echo "$as_me:$LINENO: WARNING: $ac_header: proceeding with the compiler's result" >&5
-$as_echo "$as_me: WARNING: $ac_header: proceeding with the compiler's result" >&2;}
-    ac_header_preproc=yes
-    ;;
-  no:yes:* )
-    { $as_echo "$as_me:$LINENO: WARNING: $ac_header: present but cannot be compiled" >&5
-$as_echo "$as_me: WARNING: $ac_header: present but cannot be compiled" >&2;}
-    { $as_echo "$as_me:$LINENO: WARNING: $ac_header:     check for missing prerequisite headers?" >&5
-$as_echo "$as_me: WARNING: $ac_header:     check for missing prerequisite headers?" >&2;}
-    { $as_echo "$as_me:$LINENO: WARNING: $ac_header: see the Autoconf documentation" >&5
-$as_echo "$as_me: WARNING: $ac_header: see the Autoconf documentation" >&2;}
-    { $as_echo "$as_me:$LINENO: WARNING: $ac_header:     section \"Present But Cannot Be Compiled\"" >&5
-$as_echo "$as_me: WARNING: $ac_header:     section \"Present But Cannot Be Compiled\"" >&2;}
-    { $as_echo "$as_me:$LINENO: WARNING: $ac_header: proceeding with the preprocessor's result" >&5
-$as_echo "$as_me: WARNING: $ac_header: proceeding with the preprocessor's result" >&2;}
-    { $as_echo "$as_me:$LINENO: WARNING: $ac_header: in the future, the compiler will take precedence" >&5
-$as_echo "$as_me: WARNING: $ac_header: in the future, the compiler will take precedence" >&2;}
-    ( cat <<\_ASBOX
-## ------------------------------------------- ##
-## Report this to mass-general@curl.sai.msu.ru ##
-## ------------------------------------------- ##
-_ASBOX
-     ) | sed "s/^/$as_me: WARNING:     /" >&2
-    ;;
-esac
-{ $as_echo "$as_me:$LINENO: checking for $ac_header" >&5
-$as_echo_n "checking for $ac_header... " >&6; }
-if { as_var=$as_ac_Header; eval "test \"\${$as_var+set}\" = set"; }; then
-  $as_echo_n "(cached) " >&6
-else
-  eval "$as_ac_Header=\$ac_header_preproc"
-fi
-ac_res=`eval 'as_val=${'$as_ac_Header'}
-		 $as_echo "$as_val"'`
-	       { $as_echo "$as_me:$LINENO: result: $ac_res" >&5
-$as_echo "$ac_res" >&6; }
-
-fi
-as_val=`eval 'as_val=${'$as_ac_Header'}
-		 $as_echo "$as_val"'`
-   if test "x$as_val" = x""yes; then
-  cat >>confdefs.h <<_ACEOF
-#define `$as_echo "HAVE_$ac_header" | $as_tr_cpp` 1
-_ACEOF
-
-fi
-
-done
-
-
-
-
 for ac_header in gsl/gsl_errno.h \
 	gsl/gsl_integration.h
 do
@@ -5445,7 +5303,7 @@
 # report actual input values of CONFIG_FILES etc. instead of their
 # values after options handling.
 ac_log="
-This file was extended by atmos $as_me 2.96.9, which was
+This file was extended by atmos $as_me 2.97.1, which was
 generated by GNU Autoconf 2.63.  Invocation command line was
 
   CONFIG_FILES    = $CONFIG_FILES
@@ -5508,7 +5366,7 @@
 _ACEOF
 cat >>$CONFIG_STATUS <<_ACEOF || ac_write_fail=1
 ac_cs_version="\\
-atmos config.status 2.96.9
+atmos config.status 2.97.1
 configured by $0, generated by GNU Autoconf 2.63,
   with options \\"`$as_echo "$ac_configure_args" | sed 's/^ //; s/\\""\`\$/\\\\&/g'`\\"

atmos-2.96.9.tar.gz/configure.in -> atmos-2.97.1.tar.gz/configure.in Changed

@@ -2,7 +2,7 @@
 # Process this file with autoconf to produce a configure script.
 
 AC_PREREQ(2.59)
-AC_INIT(atmos, 2.96.9, mass-general@curl.sai.msu.ru)
+AC_INIT(atmos, 2.97.1, mass-general@curl.sai.msu.ru)
 AM_INIT_AUTOMAKE()
 AC_CONFIG_HEADERS(src/config.h)
 AC_CONFIG_SRCDIR(src/atmos.h)
@@ -18,11 +18,13 @@
 
 # Checks for header files.
 AC_CHECK_HEADERS(stdlib.h string.h unistd.h)
-AC_CHECK_HEADERS(boost/numeric/ublas/io.hpp \
+AC_CHECK_HEADERS(boost/numeric/ublas/banded.hpp \
+	boost/numeric/ublas/matrix_expression.hpp \
 	boost/numeric/ublas/matrix.hpp \
 	boost/numeric/ublas/matrix_proxy.hpp \
 	boost/numeric/ublas/symmetric.hpp \
 	boost/numeric/ublas/triangular.hpp \
+	boost/numeric/ublas/vector_expression.hpp \
 	boost/numeric/ublas/vector.hpp \
 	boost/numeric/ublas/vector_proxy.hpp,,AC_MSG_ERROR(can't find boost ublas headers))
 
@@ -35,13 +37,14 @@
 	boost/format.hpp \
 	boost/function.hpp \
 	boost/io/ios_state.hpp \
-	boost/iterator.hpp \
+	boost/iterator/transform_iterator.hpp \
 	boost/lexical_cast.hpp \
+	boost/math/tools/roots.hpp \
+	boost/mem_fn.hpp \
+	boost/mpl/if.hpp \
 	boost/program_options.hpp \
-	boost/tuple/tuple.hpp,,AC_MSG_ERROR(can't find boost headers))
-
-AC_CHECK_HEADERS(boost/math/tools/roots.hpp \
-	tr1/tuple)
+	boost/tuple/tuple.hpp \
+	boost/type_traits/is_complex.hpp,,AC_MSG_ERROR(can't find boost headers))
 
 AC_CHECK_HEADERS(gsl/gsl_errno.h \
 	gsl/gsl_integration.h ,,AC_MSG_ERROR(can't find gsl headers))

atmos-2.96.9.tar.gz/src/Makefile.am -> atmos-2.97.1.tar.gz/src/Makefile.am Changed

@@ -2,15 +2,48 @@
 
 atmos_CPPFLAGS = -I @abs_top_srcdir@/third-party/lsp/include -I @abs_top_srcdir@/third-party/yalpa/include
 
-atmos_SOURCES  = apply_to_all.h atmos.h average.h background_factory.h background_factory.cpp \
-		dimm.h dimm.cpp \
-		file_input.cpp file_input.h file_output.cpp file_output.h fft.h \
-		indices.cpp indices.h integrals.cpp integrals.h \
-		io.cpp io.h gsl_wrap.h lambert.h main.cpp nocase_less.h params.h \
-		photometry.h profile.cpp profile.h spectral_factory.cpp \
-		spectral_factory.h spline.h table_storage.h typesdef.h unstable_remove_if.h utils.h \
-		vector_generator.h worksheet.h worksheet.cpp \
-		weif.cpp weif.h
+atmos_SOURCES  = apply_to_all.h \
+	atmos.h \
+	average.h \
+	background_factory.h \
+	covariance_decomposition.h \
+	dimm.h \
+	fft.h \
+	file_input.h \
+	file_output.h \
+	gsl_wrap.h \
+	indices.h \
+	integrals.h \
+	io.h \
+	lambert.h \
+	local_context.h \
+	nocase_less.h \
+	params.h \
+	photometry.h \
+	profile.h \
+	spectral_factory.h \
+	spline.h \
+	table_storage.h \
+	typesdef.h \
+	unstable_remove_if.h \
+	utils.h \
+	vector_generator.h \
+	weif.h \
+	worksheet.h \
+	background_factory.cpp \
+	dimm.cpp \
+	file_input.cpp \
+	file_output.cpp \
+	indices.cpp \
+	integrals.cpp \
+	io.cpp \
+	main.cpp \
+	params.cpp \
+	profile.cpp \
+	spectral_factory.cpp \
+	utils.cpp \
+	weif.cpp \
+	worksheet.cpp
 
 atmos_LDADD    = @abs_top_builddir@/third-party/yalpa/src/.libs/libyalpa.a
 atmos_DEPENDENCIES = @abs_top_builddir@/third-party/yalpa/src/.libs/libyalpa.a

atmos-2.96.9.tar.gz/src/Makefile.in -> atmos-2.97.1.tar.gz/src/Makefile.in Changed

@@ -50,9 +50,10 @@
 	atmos-dimm.$(OBJEXT) atmos-file_input.$(OBJEXT) \
 	atmos-file_output.$(OBJEXT) atmos-indices.$(OBJEXT) \
 	atmos-integrals.$(OBJEXT) atmos-io.$(OBJEXT) \
-	atmos-main.$(OBJEXT) atmos-profile.$(OBJEXT) \
-	atmos-spectral_factory.$(OBJEXT) atmos-worksheet.$(OBJEXT) \
-	atmos-weif.$(OBJEXT)
+	atmos-main.$(OBJEXT) atmos-params.$(OBJEXT) \
+	atmos-profile.$(OBJEXT) atmos-spectral_factory.$(OBJEXT) \
+	atmos-utils.$(OBJEXT) atmos-weif.$(OBJEXT) \
+	atmos-worksheet.$(OBJEXT)
 atmos_OBJECTS = $(am_atmos_OBJECTS)
 DEFAULT_INCLUDES = -I.@am__isrc@
 depcomp = $(SHELL) $(top_srcdir)/depcomp
@@ -159,15 +160,48 @@
 top_builddir = @top_builddir@
 top_srcdir = @top_srcdir@
 atmos_CPPFLAGS = -I @abs_top_srcdir@/third-party/lsp/include -I @abs_top_srcdir@/third-party/yalpa/include
-atmos_SOURCES = apply_to_all.h atmos.h average.h background_factory.h background_factory.cpp \
-		dimm.h dimm.cpp \
-		file_input.cpp file_input.h file_output.cpp file_output.h fft.h \
-		indices.cpp indices.h integrals.cpp integrals.h \
-		io.cpp io.h gsl_wrap.h lambert.h main.cpp nocase_less.h params.h \
-		photometry.h profile.cpp profile.h spectral_factory.cpp \
-		spectral_factory.h spline.h table_storage.h typesdef.h unstable_remove_if.h utils.h \
-		vector_generator.h worksheet.h worksheet.cpp \
-		weif.cpp weif.h
+atmos_SOURCES = apply_to_all.h \
+	atmos.h \
+	average.h \
+	background_factory.h \
+	covariance_decomposition.h \
+	dimm.h \
+	fft.h \
+	file_input.h \
+	file_output.h \
+	gsl_wrap.h \
+	indices.h \
+	integrals.h \
+	io.h \
+	lambert.h \
+	local_context.h \
+	nocase_less.h \
+	params.h \
+	photometry.h \
+	profile.h \
+	spectral_factory.h \
+	spline.h \
+	table_storage.h \
+	typesdef.h \
+	unstable_remove_if.h \
+	utils.h \
+	vector_generator.h \
+	weif.h \
+	worksheet.h \
+	background_factory.cpp \
+	dimm.cpp \
+	file_input.cpp \
+	file_output.cpp \
+	indices.cpp \
+	integrals.cpp \
+	io.cpp \
+	main.cpp \
+	params.cpp \
+	profile.cpp \
+	spectral_factory.cpp \
+	utils.cpp \
+	weif.cpp \
+	worksheet.cpp
 
 atmos_LDADD = @abs_top_builddir@/third-party/yalpa/src/.libs/libyalpa.a
 atmos_DEPENDENCIES = @abs_top_builddir@/third-party/yalpa/src/.libs/libyalpa.a
@@ -278,8 +312,10 @@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/atmos-integrals.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/atmos-io.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/atmos-main.Po@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/atmos-params.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/atmos-profile.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/atmos-spectral_factory.Po@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/atmos-utils.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/atmos-weif.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/atmos-worksheet.Po@am__quote@
 
@@ -409,6 +445,20 @@
 @AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCXX_FALSE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o atmos-main.obj `if test -f 'main.cpp'; then $(CYGPATH_W) 'main.cpp'; else $(CYGPATH_W) '$(srcdir)/main.cpp'; fi`
 
+atmos-params.o: params.cpp
+@am__fastdepCXX_TRUE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT atmos-params.o -MD -MP -MF $(DEPDIR)/atmos-params.Tpo -c -o atmos-params.o `test -f 'params.cpp' || echo '$(srcdir)/'`params.cpp
+@am__fastdepCXX_TRUE@	$(am__mv) $(DEPDIR)/atmos-params.Tpo $(DEPDIR)/atmos-params.Po
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	source='params.cpp' object='atmos-params.o' libtool=no @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o atmos-params.o `test -f 'params.cpp' || echo '$(srcdir)/'`params.cpp
+
+atmos-params.obj: params.cpp
+@am__fastdepCXX_TRUE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT atmos-params.obj -MD -MP -MF $(DEPDIR)/atmos-params.Tpo -c -o atmos-params.obj `if test -f 'params.cpp'; then $(CYGPATH_W) 'params.cpp'; else $(CYGPATH_W) '$(srcdir)/params.cpp'; fi`
+@am__fastdepCXX_TRUE@	$(am__mv) $(DEPDIR)/atmos-params.Tpo $(DEPDIR)/atmos-params.Po
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	source='params.cpp' object='atmos-params.obj' libtool=no @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o atmos-params.obj `if test -f 'params.cpp'; then $(CYGPATH_W) 'params.cpp'; else $(CYGPATH_W) '$(srcdir)/params.cpp'; fi`
+
 atmos-profile.o: profile.cpp
 @am__fastdepCXX_TRUE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT atmos-profile.o -MD -MP -MF $(DEPDIR)/atmos-profile.Tpo -c -o atmos-profile.o `test -f 'profile.cpp' || echo '$(srcdir)/'`profile.cpp
 @am__fastdepCXX_TRUE@	$(am__mv) $(DEPDIR)/atmos-profile.Tpo $(DEPDIR)/atmos-profile.Po
@@ -437,19 +487,19 @@
 @AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCXX_FALSE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o atmos-spectral_factory.obj `if test -f 'spectral_factory.cpp'; then $(CYGPATH_W) 'spectral_factory.cpp'; else $(CYGPATH_W) '$(srcdir)/spectral_factory.cpp'; fi`
 
-atmos-worksheet.o: worksheet.cpp
-@am__fastdepCXX_TRUE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT atmos-worksheet.o -MD -MP -MF $(DEPDIR)/atmos-worksheet.Tpo -c -o atmos-worksheet.o `test -f 'worksheet.cpp' || echo '$(srcdir)/'`worksheet.cpp
-@am__fastdepCXX_TRUE@	$(am__mv) $(DEPDIR)/atmos-worksheet.Tpo $(DEPDIR)/atmos-worksheet.Po
-@AMDEP_TRUE@@am__fastdepCXX_FALSE@	source='worksheet.cpp' object='atmos-worksheet.o' libtool=no @AMDEPBACKSLASH@
-@AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
-@am__fastdepCXX_FALSE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o atmos-worksheet.o `test -f 'worksheet.cpp' || echo '$(srcdir)/'`worksheet.cpp
-
-atmos-worksheet.obj: worksheet.cpp
-@am__fastdepCXX_TRUE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT atmos-worksheet.obj -MD -MP -MF $(DEPDIR)/atmos-worksheet.Tpo -c -o atmos-worksheet.obj `if test -f 'worksheet.cpp'; then $(CYGPATH_W) 'worksheet.cpp'; else $(CYGPATH_W) '$(srcdir)/worksheet.cpp'; fi`
-@am__fastdepCXX_TRUE@	$(am__mv) $(DEPDIR)/atmos-worksheet.Tpo $(DEPDIR)/atmos-worksheet.Po
-@AMDEP_TRUE@@am__fastdepCXX_FALSE@	source='worksheet.cpp' object='atmos-worksheet.obj' libtool=no @AMDEPBACKSLASH@
+atmos-utils.o: utils.cpp
+@am__fastdepCXX_TRUE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT atmos-utils.o -MD -MP -MF $(DEPDIR)/atmos-utils.Tpo -c -o atmos-utils.o `test -f 'utils.cpp' || echo '$(srcdir)/'`utils.cpp
+@am__fastdepCXX_TRUE@	$(am__mv) $(DEPDIR)/atmos-utils.Tpo $(DEPDIR)/atmos-utils.Po
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	source='utils.cpp' object='atmos-utils.o' libtool=no @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o atmos-utils.o `test -f 'utils.cpp' || echo '$(srcdir)/'`utils.cpp
+
+atmos-utils.obj: utils.cpp
+@am__fastdepCXX_TRUE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT atmos-utils.obj -MD -MP -MF $(DEPDIR)/atmos-utils.Tpo -c -o atmos-utils.obj `if test -f 'utils.cpp'; then $(CYGPATH_W) 'utils.cpp'; else $(CYGPATH_W) '$(srcdir)/utils.cpp'; fi`
+@am__fastdepCXX_TRUE@	$(am__mv) $(DEPDIR)/atmos-utils.Tpo $(DEPDIR)/atmos-utils.Po
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	source='utils.cpp' object='atmos-utils.obj' libtool=no @AMDEPBACKSLASH@
 @AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
-@am__fastdepCXX_FALSE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o atmos-worksheet.obj `if test -f 'worksheet.cpp'; then $(CYGPATH_W) 'worksheet.cpp'; else $(CYGPATH_W) '$(srcdir)/worksheet.cpp'; fi`
+@am__fastdepCXX_FALSE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o atmos-utils.obj `if test -f 'utils.cpp'; then $(CYGPATH_W) 'utils.cpp'; else $(CYGPATH_W) '$(srcdir)/utils.cpp'; fi`
 
 atmos-weif.o: weif.cpp
 @am__fastdepCXX_TRUE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT atmos-weif.o -MD -MP -MF $(DEPDIR)/atmos-weif.Tpo -c -o atmos-weif.o `test -f 'weif.cpp' || echo '$(srcdir)/'`weif.cpp
@@ -465,6 +515,20 @@
 @AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCXX_FALSE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o atmos-weif.obj `if test -f 'weif.cpp'; then $(CYGPATH_W) 'weif.cpp'; else $(CYGPATH_W) '$(srcdir)/weif.cpp'; fi`
 
+atmos-worksheet.o: worksheet.cpp
+@am__fastdepCXX_TRUE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT atmos-worksheet.o -MD -MP -MF $(DEPDIR)/atmos-worksheet.Tpo -c -o atmos-worksheet.o `test -f 'worksheet.cpp' || echo '$(srcdir)/'`worksheet.cpp
+@am__fastdepCXX_TRUE@	$(am__mv) $(DEPDIR)/atmos-worksheet.Tpo $(DEPDIR)/atmos-worksheet.Po
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	source='worksheet.cpp' object='atmos-worksheet.o' libtool=no @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o atmos-worksheet.o `test -f 'worksheet.cpp' || echo '$(srcdir)/'`worksheet.cpp
+
+atmos-worksheet.obj: worksheet.cpp
+@am__fastdepCXX_TRUE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT atmos-worksheet.obj -MD -MP -MF $(DEPDIR)/atmos-worksheet.Tpo -c -o atmos-worksheet.obj `if test -f 'worksheet.cpp'; then $(CYGPATH_W) 'worksheet.cpp'; else $(CYGPATH_W) '$(srcdir)/worksheet.cpp'; fi`
+@am__fastdepCXX_TRUE@	$(am__mv) $(DEPDIR)/atmos-worksheet.Tpo $(DEPDIR)/atmos-worksheet.Po
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	source='worksheet.cpp' object='atmos-worksheet.obj' libtool=no @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o atmos-worksheet.obj `if test -f 'worksheet.cpp'; then $(CYGPATH_W) 'worksheet.cpp'; else $(CYGPATH_W) '$(srcdir)/worksheet.cpp'; fi`
+
 ID: $(HEADERS) $(SOURCES) $(LISP) $(TAGS_FILES)
 	list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \
 	unique=`for i in $$list; do \

atmos-2.96.9.tar.gz/src/apply_to_all.h -> atmos-2.97.1.tar.gz/src/apply_to_all.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: apply_to_all.h 155 2010-09-19 08:30:20Z matwey $
+    $Id: apply_to_all.h 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -21,10 +21,10 @@
 
 #include <cmath>
 
-#include "lambert.h"
+#include <boost/numeric/ublas/matrix_expression.hpp>
+#include <boost/numeric/ublas/vector_expression.hpp>
 
-#include <boost/numeric/ublas/matrix.hpp>
-#include <boost/numeric/ublas/vector.hpp>
+#include "lambert.h"
 
 namespace utils {
 namespace ublas {
@@ -78,6 +78,14 @@
 			return x + value_type(1);
 		}
 	};
+	template <class T> class abs {
+	public:
+		typedef T value_type;
+		typedef T result_type;
+		static result_type apply( const value_type& x) {
+			return std::abs(x);
+		}
+	};
 }
 
 using namespace boost::numeric::ublas;
@@ -146,6 +154,15 @@
 	return apply_to_all( e, detail::lambert_W0< value_type >() );
 }
 
+template<class E> BOOST_UBLAS_INLINE typename vector_unary_traits<E, detail::abs< typename E::value_type > >::result_type abs(const vector_expression<E> &e ) {
+	typedef typename E::value_type value_type;
+	return apply_to_all( e, detail::abs< value_type >() );
+}
+template<class E> BOOST_UBLAS_INLINE typename matrix_unary1_traits<E, detail::abs< typename E::value_type > >::result_type abs(const matrix_expression<E> &e ) {
+	typedef typename E::value_type value_type;
+	return apply_to_all( e, detail::abs< value_type >() );
+}
+
 } //ublas
 } //utils

atmos-2.96.9.tar.gz/src/atmos.h -> atmos-2.97.1.tar.gz/src/atmos.h Changed

atmos-2.96.9.tar.gz/src/average.h -> atmos-2.97.1.tar.gz/src/average.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: average.h 143 2010-04-04 08:42:34Z matwey $
+    $Id: average.h 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2010  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -20,153 +20,156 @@
 #define _AVERAGE_H
 
 #include <cmath>
-#include <iterator>
+#include <iterator> // std::iterator_traits
 
-#include <boost/numeric/ublas/matrix.hpp>
-#include <boost/numeric/ublas/vector.hpp>
-#include <boost/numeric/ublas/io.hpp>
+#include <boost/numeric/ublas/matrix_expression.hpp>
+#include <boost/numeric/ublas/symmetric.hpp>
+#include <boost/numeric/ublas/vector_expression.hpp>
 
 #include "apply_to_all.h"
 
 namespace utils {
 namespace algo {
 
-	template<class T> struct mean;
-	template<class T> struct variance;
+namespace detail {
 
-	template<class T> struct mean {
-	public:
-		typedef T value_type;
-		typedef T result_type;
-
-		void operator() ( result_type& result, const value_type& x ) const {
-			result += x;
-		}
-		result_type initial( const value_type& x ) const {
-			return result_type( x );
-		}
-		template<class S> void result( result_type& result, S size ) const {
-			result /= size;
-		}
-	private:
-	};
-
-	struct variance_traits {
-		static double prod( double x, double y ){
-			return x*y;
-		}
-		static float prod( float x, float y ){
-			return x*y;
-		}
-		template< class V1, class V2 > static
-		typename boost::numeric::ublas::vector_binary_traits< V1, V2, boost::numeric::ublas::scalar_multiplies< typename V1::value_type, typename V2::value_type > >::result_type
-		prod( const boost::numeric::ublas::vector_expression< V1 >& e1, const boost::numeric::ublas::vector_expression< V2 >& e2 ) {
-			return element_prod( e1, e2 );
-		}
-		template< class M1, class M2 > static
-		typename boost::numeric::ublas::matrix_binary_traits< M1, M2, boost::numeric::ublas::scalar_multiplies< typename M1::value_type, typename M2::value_type > >::result_type
-		prod( const boost::numeric::ublas::matrix_expression< M1 >& e1, const boost::numeric::ublas::matrix_expression< M2 >& e2 ) {
-			return element_prod( e1, e2 );
-		}
-
-		static double sqrt( double x ){
-			return std::sqrt(x);
-		}
-		static float sqrt( float x ){
-			return std::sqrt(x);
-		}
-		template< class V1 > static
-		typename boost::numeric::ublas::vector_unary_traits< V1, utils::ublas::detail::sqrt< typename V1::value_type > >::result_type
-		sqrt( const boost::numeric::ublas::vector_expression< V1 >& e1) {
-			return utils::ublas::sqrt( e1 );
-		}
-		template< class M1 > static
-		typename boost::numeric::ublas::matrix_unary1_traits< M1, utils::ublas::detail::sqrt< typename M1::value_type > >::result_type
-		sqrt( const boost::numeric::ublas::matrix_expression< M1 >& e1) {
-			return utils::ublas::sqrt( e1 );
-		}
-	};
-
-	struct covariance_traits {
-		template< class V1, class V2 > static
-		boost::numeric::ublas::symmetric_adaptor< const typename boost::numeric::ublas::vector_matrix_binary_traits<V1, V2, boost::numeric::ublas::scalar_multiplies<typename V1::value_type, typename V2::value_type> >::result_type, boost::numeric::ublas::upper >
-		prod( const boost::numeric::ublas::vector_expression< V1 >& e1, const boost::numeric::ublas::vector_expression< V2 >& e2 ) {
-			typedef typename boost::numeric::ublas::vector_matrix_binary_traits<V1, V2, boost::numeric::ublas::scalar_multiplies<typename V1::value_type, typename V2::value_type> >::result_type outer_prod_type;
-			return boost::numeric::ublas::symmetric_adaptor< const outer_prod_type, boost::numeric::ublas::upper >( outer_prod( e1, e2 ) );
-		}
-	};
-	
-	template<class T> struct variance {
-	public:
-		typedef T value_type;
-		typedef T result_type;
-		
-		variance( const value_type& mean ): mean_(mean) {}
-		
-		void operator() ( result_type& result, const value_type& x ) const {
-			result += variance_traits::prod( x - mean_, x - mean_ );
-		}
-		result_type initial( const value_type& x ) const {
-			return result_type( variance_traits::prod( x - mean_, x - mean_ ) );
-		}
-		template<class S> void result( result_type& result, S size ) const {
-			result /= size-1;
-		}
-	private:
-		value_type mean_;
-	};
-
-	template<class T> struct covariance {
-	public:
-		typedef T value_type;
-		typedef boost::numeric::ublas::symmetric_matrix< typename value_type::value_type > result_type;
-
-		covariance( const value_type& mean ): mean_(mean) {}
-
-		void operator() ( result_type& result, const value_type& x ) const {
-			result += covariance_traits::prod( x - mean_, x - mean_ );
-		}
-		result_type initial( const value_type& x ) const {
-			return result_type( covariance_traits::prod( x - mean_, x - mean_ ) );
-		}
-		template<class S> void result( result_type& result, S size ) const {
-			result /= size-1;
-		}
-	private:
-		value_type mean_;
-	};
-
-	template<class T> struct error_of_mean: variance<T> {
-		error_of_mean( const typename variance<T>::value_type& mean ): variance<T>(mean) {}
-		template<class S> void result( typename variance<T>::result_type& result, S size ) const {
-			variance<T>::result( result, size );
-			result = variance_traits::sqrt( result / size );
-		}
-	};
+struct variance_traits {
+	static double prod( double x, double y ){
+		return x*y;
+	}
+	static float prod( float x, float y ){
+		return x*y;
+	}
+	template< class V1, class V2 > static
+	typename boost::numeric::ublas::vector_binary_traits< V1, V2, boost::numeric::ublas::scalar_multiplies< typename V1::value_type, typename V2::value_type > >::result_type
+	prod( const boost::numeric::ublas::vector_expression< V1 >& e1, const boost::numeric::ublas::vector_expression< V2 >& e2 ) {
+		return element_prod( e1, e2 );
+	}
+	template< class M1, class M2 > static
+	typename boost::numeric::ublas::matrix_binary_traits< M1, M2, boost::numeric::ublas::scalar_multiplies< typename M1::value_type, typename M2::value_type > >::result_type
+	prod( const boost::numeric::ublas::matrix_expression< M1 >& e1, const boost::numeric::ublas::matrix_expression< M2 >& e2 ) {
+		return element_prod( e1, e2 );
+	}
+};
+
+struct error_of_mean_traits {
+	static double sqrt( double x ){
+		return std::sqrt(x);
+	}
+	static float sqrt( float x ){
+		return std::sqrt(x);
+	}
+	template< class V1 > static
+	typename boost::numeric::ublas::vector_unary_traits< V1, utils::ublas::detail::sqrt< typename V1::value_type > >::result_type
+	sqrt( const boost::numeric::ublas::vector_expression< V1 >& e1) {
+		return utils::ublas::sqrt( e1 );
+	}
+	template< class M1 > static
+	typename boost::numeric::ublas::matrix_unary1_traits< M1, utils::ublas::detail::sqrt< typename M1::value_type > >::result_type
+	sqrt( const boost::numeric::ublas::matrix_expression< M1 >& e1) {
+		return utils::ublas::sqrt( e1 );
+	}
+};
+
+struct covariance_traits {
+	template< class V1, class V2 > static
+	boost::numeric::ublas::symmetric_adaptor< const typename boost::numeric::ublas::vector_matrix_binary_traits<V1, V2, boost::numeric::ublas::scalar_multiplies<typename V1::value_type, typename V2::value_type> >::result_type, boost::numeric::ublas::upper >
+	prod( const boost::numeric::ublas::vector_expression< V1 >& e1, const boost::numeric::ublas::vector_expression< V2 >& e2 ) {
+		typedef typename boost::numeric::ublas::vector_matrix_binary_traits<V1, V2, boost::numeric::ublas::scalar_multiplies<typename V1::value_type, typename V2::value_type> >::result_type outer_prod_type;
+		return boost::numeric::ublas::symmetric_adaptor< const outer_prod_type, boost::numeric::ublas::upper >( outer_prod( e1, e2 ) );
+	}
+};
 	
-	template< class Iterator, template<class T> class F >
-	typename F< typename std::iterator_traits<Iterator>::value_type >::result_type
-	average( Iterator begin, Iterator end, F< typename std::iterator_traits<Iterator>::value_type > op ) {
-		typedef Iterator iterator;
-		typedef typename std::iterator_traits< iterator >::value_type value_type;
-		typedef typename std::iterator_traits< iterator >::difference_type difference_type;
-		typedef difference_type size_type;
-		typedef F< value_type > operator_type;
-		typedef typename operator_type::result_type result_type;
+}
+
+template<class T> struct mean {
+public:
+	typedef T value_type;
+	typedef T result_type;
 
-		if( begin == end )
-			return result_type();
+	void operator() ( result_type& result, const value_type& x ) const {
+		result += x;
+	}
+	result_type initial( const value_type& x ) const {
+		return result_type( x );
+	}
+	template<class S> void result( result_type& result, S size ) const {
+		result /= size;
+	}
+private:
+};
+
+template<class T> struct variance {
+public:
+	typedef T value_type;
+	typedef T result_type;
 		
-		result_type ret( op.initial( *begin ) );
-		size_type n = 1;
+	variance( const value_type& mean ): mean_(mean) {}
 		
-		for( iterator it = ++begin; it != end; ++it,++n ) {
-			op( ret, *it );
-		}
+	void operator() ( result_type& result, const value_type& x ) const {
+		result += detail::variance_traits::prod( x - mean_, x - mean_ );
+	}
+	result_type initial( const value_type& x ) const {
+		return result_type( detail::variance_traits::prod( x - mean_, x - mean_ ) );
+	}
+	template<class S> void result( result_type& result, S size ) const {
+		result /= size-1;
+	}
+private:
+	value_type mean_;
+};
+
+template<class T> struct covariance {
+public:
+	typedef T value_type;
+	typedef boost::numeric::ublas::symmetric_matrix< typename value_type::value_type > result_type;
 
-		op.result( ret, n );
-		return ret;
+	covariance( const value_type& mean ): mean_(mean) {}
+
+	void operator() ( result_type& result, const value_type& x ) const {
+		result += detail::covariance_traits::prod( x - mean_, x - mean_ );
+	}
+	result_type initial( const value_type& x ) const {
+		return result_type( detail::covariance_traits::prod( x - mean_, x - mean_ ) );
 	}
+	template<class S> void result( result_type& result, S size ) const {
+		result /= size-1;
+	}
+private:
+	value_type mean_;
+};
+
+template<class T> struct error_of_mean: variance<T> {
+	error_of_mean( const typename variance<T>::value_type& mean ): variance<T>(mean) {}
+	
+	template<class S> void result( typename variance<T>::result_type& result, S size ) const {
+		variance<T>::result( result, size );
+		result = detail::error_of_mean_traits::sqrt( result / size );
+	}
+};
+	
+template< class Iterator, template<class T> class F >
+typename F< typename std::iterator_traits<Iterator>::value_type >::result_type
+average( Iterator begin, Iterator end, F< typename std::iterator_traits<Iterator>::value_type > op ) {
+	typedef Iterator iterator;
+	typedef typename std::iterator_traits< iterator >::value_type value_type;
+	typedef typename std::iterator_traits< iterator >::difference_type difference_type;
+	typedef difference_type size_type;
+	typedef F< value_type > operator_type;
+	typedef typename operator_type::result_type result_type;
+
+	if( begin == end )
+		return result_type();
+
+	result_type ret( op.initial( *begin ) );
+	size_type n = 1;
+
+	for( iterator it = ++begin; it != end; ++it,++n )
+		op( ret, *it );
+
+	op.result( ret, n );
+	return ret;
+}
 
 } // algo
 } // utils

atmos-2.96.9.tar.gz/src/background_factory.cpp -> atmos-2.97.1.tar.gz/src/background_factory.cpp Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: background_factory.cpp 147 2010-04-17 09:25:25Z matwey $
+    $Id: background_factory.cpp 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -16,73 +16,79 @@
     along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
 
+#include <stdexcept>
+#include <boost/format.hpp>
+
 #include "background_factory.h"
 
-#include <boost/format.hpp>
+namespace atmos {
+namespace background {
 
-background_factory::background_factory(const select_model_function_type& select_model):
+factory::factory(const select_model_function_type& select_model):
 	m_curmodel(0),
 	m_select_model(select_model) {
 }
-background_factory::background_factory(const background_factory& factory) {
+factory::factory(const factory& factory) {
 	m_select_model = factory.m_select_model;
 	m_curmodel = std::auto_ptr<model_type>(factory.m_curmodel->clone());
 }
-void background_factory::swap(background_factory& factory) throw() {
-	factory.m_select_model.swap( m_select_model );
-	std::swap( factory.m_curmodel, m_curmodel );
+void factory::swap(factory& f) throw() {
+	f.m_select_model.swap( m_select_model );
+	std::swap( f.m_curmodel, m_curmodel );
 }
-background_factory& background_factory::operator=(const background_factory& factory){
-	background_factory tmp( factory );
+factory& factory::operator=(const factory& f){
+	factory tmp( f );
 	tmp.swap(*this );
 	return *this;
 }
-background_factory::~background_factory() {
+factory::~factory() {
 }
-background_factory::vector_type background_factory::background( const boost::posix_time::ptime& time ) const {
+factory::flux_type factory::background( const time_type& time ) const {
 	if( ! m_curmodel.get() )
 		throw std::logic_error("There is not any background model");
 	return m_curmodel->background( time );
 }
-background_factory::model_type* background_factory::create_new_model( const boost::posix_time::ptime& time ) const {
+factory::model_type* factory::create_new_model( const time_type& time ) const {
 	return m_select_model( time );
 }
-background_factory::model_type* background_factory::release_current_model() {
+factory::model_type* factory::release_current_model() {
 	return m_curmodel.release();
 }
-void background_factory::handle_new_model( model_type* model ) {
+void factory::handle_new_model( model_type* model ) {
 	m_curmodel = std::auto_ptr<model_type>( model );
 }
 
-background_model::background_model( const boost::posix_time::ptime& time ): m_inittime(time) {
+namespace model {
+
+abstract::abstract( const time_type& time ): m_inittime(time) {
 }
-background_model::~background_model() {
+abstract::~abstract() {
 }
 
-background_daytime_model::background_daytime_model( const boost::posix_time::ptime& time ):
-	background_model(time),
+daytime::daytime( const time_type& time ):
+	abstract(time),
 	m_size(0) {
 }
-background_daytime_model::size_type background_daytime_model::size() const {
+daytime::size_type daytime::size() const {
 	return m_size;
 }
-background_daytime_model::vector_type background_daytime_model::background( const boost::posix_time::ptime& time ) const {
+daytime::flux_type daytime::background( const time_type& time ) const {
 	throw std::logic_error("Daytime background observation");
 }
-void background_daytime_model::accumulate( const boost::posix_time::ptime& time, const vector_type& mean ) {
+void daytime::accumulate( const time_type& time, const flux_type& flux ) {
 	m_size++;
 }
-background_model* background_daytime_model::clone() const {
-	return new background_daytime_model(*this);
+abstract* daytime::clone() const {
+	return new daytime(*this);
 }
 
-background_twilight_model::background_twilight_model( const boost::posix_time::ptime& time, const sun_altitude_function_type& sun_altitude ):
-	background_model(time),
+twilight::twilight( const time_type& time, const sun_altitude_function_type& sun_altitude ):
+	abstract(time),
 	m_size(0),
 	m_background(0),
 	m_sun_altitude(sun_altitude) {
 }
-background_twilight_model::value_type background_twilight_model::extrapolate( const boost::posix_time::ptime& timef ) const {
+value_type twilight::extrapolate( const time_type& timef ) const {
 	value_type hi = m_sun_altitude( inittime() ) * 180.0 / M_PI; // degrees
 	value_type hf = m_sun_altitude( timef ) * 180.0 / M_PI; // degrees
 	if( hi < -12.0 ) hi = -12.0;
@@ -90,62 +96,66 @@
 	value_type scale = 0.52189 - 0.26242*hi - 0.0235464*hi*hi; // model
 	return std::exp( scale * ( hf - hi ) );
 }
-background_twilight_model::size_type background_twilight_model::size() const {
+twilight::size_type twilight::size() const {
 	return m_size;
 }
-background_twilight_model::vector_type background_twilight_model::background( const boost::posix_time::ptime& time ) const {
+twilight::flux_type twilight::background( const time_type& time ) const {
 	return m_background * extrapolate( time );
 }
-void background_twilight_model::accumulate( const boost::posix_time::ptime& time, const vector_type& mean ) {
+void twilight::accumulate( const time_type& time, const flux_type& flux ) {
 
 	if( m_size == 0 ) {
-		m_background = mean / extrapolate( time );
+		m_background = flux / extrapolate( time );
 		m_size++;
 		return;
 	}
 
-	m_background = m_background * ( value_type( m_size ) / value_type( m_size + 1 ) ) + mean / extrapolate( time ) * ( value_type(1) / value_type( m_size + 1 ) );
+	m_background = m_background * ( value_type( m_size ) / value_type( m_size + 1 ) ) + flux / extrapolate( time ) * ( value_type(1) / value_type( m_size + 1 ) );
 	m_size++;
 }
-background_model* background_twilight_model::clone() const {
-	return new background_twilight_model(*this);
+abstract* twilight::clone() const {
+	return new twilight(*this);
 }
 
-background_night_model::background_night_model( const boost::posix_time::ptime& time ):
-	background_model(time),
+night::night( const time_type& time ):
+	abstract(time),
 	m_size(0),
 	m_background(0) {
 }
-background_night_model::size_type background_night_model::size() const {
+night::size_type night::size() const {
 	return m_size;
 }
-background_night_model::vector_type background_night_model::background( const boost::posix_time::ptime& time ) const {
+night::flux_type night::background( const time_type& time ) const {
 	return m_background;
 }
-void background_night_model::accumulate( const boost::posix_time::ptime& time, const vector_type& mean ) {
+void night::accumulate( const time_type& time, const flux_type& flux ) {
 	if( m_size == 0 ) {
-		m_background = mean;
+		m_background = flux;
 		m_size++;
 		return;
 	}
 
-	assert( mean.size() == m_background.size() );
+	assert( flux.size() == m_background.size() );
 
-	m_background = m_background * ( value_type( m_size ) / value_type( m_size + 1 ) ) + mean * ( value_type(1) / value_type( m_size + 1 ) );
+	m_background = m_background * ( value_type( m_size ) / value_type( m_size + 1 ) ) + flux * ( value_type(1) / value_type( m_size + 1 ) );
 	m_size++;
 }
-background_model* background_night_model::clone() const {
-	return new background_night_model(*this);
+abstract* night::clone() const {
+	return new night(*this);
 }
 
-background_nonpoisson_filter::background_nonpoisson_filter( const value_type& threshold, const value_type& nonpoissonD, const value_type& numexp ):
+}
+
+namespace filter {
+
+nonpoisson::nonpoisson( const value_type& threshold, const value_type& nonpoissonD, const value_type& numexp ):
 	m_threshold( threshold ),
 	m_nonpoissonD( nonpoissonD ),
 	m_numexp( numexp ) {
 }
-background_nonpoisson_filter::~background_nonpoisson_filter() {
+nonpoisson::~nonpoisson() {
 }
-bool background_nonpoisson_filter::operator() ( const boost::posix_time::ptime& time, const flux_type& flux, const vars_type& vars ) const {
+bool nonpoisson::operator() ( const time_type& time, const flux_type& flux, const vars_type& vars ) const {
 	double poisson = ( vars(3,3)*flux(3) );
 	double relepsp = 2.0 * ( 1.0 + 1.0/flux(3) ) / ( m_numexp );
 	double pmax = m_nonpoissonD * ( 1.0 + m_threshold * std::sqrt(relepsp) );
@@ -155,3 +165,6 @@
 	}
 	return false;
 }
+
+}
+}}

atmos-2.96.9.tar.gz/src/background_factory.h -> atmos-2.97.1.tar.gz/src/background_factory.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: background_factory.h 147 2010-04-17 09:25:25Z matwey $
+    $Id: background_factory.h 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -19,8 +19,7 @@
 #ifndef _BACKGROUND_FACTORY_H
 #define _BACKGROUND_FACTORY_H
 
-#include <list>
-#include <memory>
+#include <memory> // std::auto_ptr
 
 #include <boost/function.hpp>
 #include <boost/date_time/posix_time/posix_time.hpp>
@@ -28,133 +27,137 @@
 #include <boost/numeric/ublas/vector.hpp>
 #include <boost/numeric/ublas/symmetric.hpp>
 
-using namespace boost::numeric;
+namespace atmos {
+namespace background {
 
-/* Forward declaration */
-class background_factory;
+using namespace boost::numeric::ublas;
 
-class background_model;
-class background_daytime_model;
-class background_twilight_model;
-class background_night_model;
-
-class background_abstract_filter;
-class background_nonpoisson_filter;
-
-/* Declaration */
-class background_factory {
-	public:
-		typedef background_abstract_filter filter_type;
-		typedef background_model  model_type;
-
-		typedef boost::function1<model_type*, const boost::posix_time::ptime&> select_model_function_type;
-
-		typedef double                     value_type;
-		typedef ublas::vector<value_type>  vector_type;
-		typedef vector_type::size_type     size_type;
-
-		typedef vector_type flux_type;
-		typedef ublas::symmetric_matrix<value_type> vars_type;
-	private:
-		std::auto_ptr< model_type > m_curmodel;
-		select_model_function_type  m_select_model;
-	public:
-		background_factory(const select_model_function_type& select_model);
-		background_factory(const background_factory& factory);
-		void swap(background_factory& factory) throw();
-		background_factory& operator=(const background_factory& factory);
-		~background_factory();
-
-		vector_type background( const boost::posix_time::ptime& time ) const;
-
-		model_type* create_new_model( const boost::posix_time::ptime& time ) const;
-		model_type* release_current_model();
-		void handle_new_model( model_type* model );
-};
-
-class background_model {
-	public:
-		typedef background_factory::value_type  value_type;
-		typedef background_factory::vector_type vector_type;
-		typedef background_factory::size_type   size_type;
-	private:
-		boost::posix_time::ptime m_inittime;
-	public:
-		background_model( const boost::posix_time::ptime& time );
-		virtual ~background_model() = 0;
-
-		const boost::posix_time::ptime& inittime() const { return m_inittime; }
-
-		virtual size_type size() const = 0;
-		virtual vector_type background( const boost::posix_time::ptime& time ) const = 0;
-		virtual void accumulate( const boost::posix_time::ptime& time, const vector_type& mean ) = 0;
-
-		virtual background_model* clone() const = 0;
-};
-
-class background_daytime_model: public background_model {
-	private:
-		size_type   m_size;
-	public:
-		background_daytime_model( const boost::posix_time::ptime& time );
-
-		size_type size() const;
-		vector_type background( const boost::posix_time::ptime& time ) const;
-		void accumulate( const boost::posix_time::ptime& time, const vector_type& mean );
-
-		background_model* clone() const;
-};
-
-class background_twilight_model: public background_model {
-	public:
-		typedef boost::function1<value_type, const boost::posix_time::ptime&> sun_altitude_function_type;
-	private:
-		size_type   m_size;
-		vector_type m_background;
-		sun_altitude_function_type m_sun_altitude;
-	private:
-		value_type extrapolate( const boost::posix_time::ptime& timef ) const;
-	public:
-		background_twilight_model( const boost::posix_time::ptime& time, const sun_altitude_function_type& sun_altitude );
-
-		size_type size() const;
-		vector_type background( const boost::posix_time::ptime& time ) const;
-		void accumulate( const boost::posix_time::ptime& time, const vector_type& mean );
-
-		background_model* clone() const;
-};
-
-class background_night_model: public background_model {
-	private:
-		size_type   m_size;
-		vector_type m_background;
-	public:
-		background_night_model( const boost::posix_time::ptime& time );
-
-		size_type size() const;
-		vector_type background( const boost::posix_time::ptime& time ) const;
-		void accumulate( const boost::posix_time::ptime& time, const vector_type& mean );
-
-		background_model* clone() const;
-};
-
-class background_nonpoisson_filter {
-	public:
-		typedef background_factory::value_type  value_type;
-		typedef background_factory::vector_type vector_type;
-		typedef background_factory::size_type   size_type;
-
-		typedef background_factory::flux_type   flux_type;
-		typedef background_factory::vars_type   vars_type;
-	private:
-		value_type m_threshold;
-		value_type m_nonpoissonD;
-		value_type m_numexp;
-	public:
-		background_nonpoisson_filter( const value_type& threshold, const value_type& nonpoissonD, const value_type& numexp );
-		~background_nonpoisson_filter();
+typedef boost::posix_time::ptime   time_type;
+typedef double value_type;
+typedef vector<value_type>         vector_type;
+typedef vector_type::size_type     size_type;
+typedef symmetric_matrix< value_type > symmetric_matrix_type;
+
+typedef vector_type           flux_type;
+typedef symmetric_matrix_type vars_type;
+
+namespace model {
+	class abstract;
+}
+
+class factory {
+public:
+	typedef model::abstract              model_type;
+	typedef atmos::background::time_type time_type;
+	typedef atmos::background::flux_type flux_type;
+	
+	typedef boost::function1<model_type*, const time_type&> select_model_function_type;
+private:
+	std::auto_ptr< model_type > m_curmodel;
+	select_model_function_type  m_select_model;
+public:
+	factory(const select_model_function_type& select_model);
+	factory(const factory& factory);
+	void swap(factory& factory) throw();
+	factory& operator=(const factory& factory);
+	~factory();
+
+	flux_type background( const time_type& time ) const;
+
+	model_type* create_new_model( const time_type& time ) const;
+	model_type* release_current_model();
+	void handle_new_model( model_type* model );
+};
+
+namespace model {
+	
+class abstract {
+public:
+	typedef atmos::background::size_type size_type;
+	typedef atmos::background::time_type time_type;
+	typedef atmos::background::flux_type flux_type;
+private:
+	time_type m_inittime;
+public:
+	abstract( const time_type& time );
+	virtual ~abstract() = 0;
+
+	const time_type& inittime() const { return m_inittime; }
+
+	virtual size_type size() const = 0;
+	virtual flux_type background( const time_type& time ) const = 0;
+	virtual void accumulate( const time_type& time, const flux_type& flux ) = 0;
+
+	virtual abstract* clone() const = 0;
+};
+
+class daytime: public abstract {
+private:
+	size_type   m_size;
+public:
+	daytime( const time_type& time );
+
+	size_type size() const;
+	flux_type background( const time_type& time ) const;
+	void accumulate( const time_type& time, const flux_type& flux );
+
+	abstract* clone() const;
+};
+
+class twilight: public abstract {
+public:
+	typedef boost::function1<value_type, const time_type&> sun_altitude_function_type;
+private:
+	size_type   m_size;
+	flux_type   m_background;
+	sun_altitude_function_type m_sun_altitude;
+private:
+	value_type extrapolate( const time_type& timef ) const;
+public:
+	twilight( const time_type& time, const sun_altitude_function_type& sun_altitude );
+
+	size_type size() const;
+	flux_type background( const time_type& time ) const;
+	void accumulate( const time_type& time, const flux_type& flux );
+
+	abstract* clone() const;
+};
+
+class night: public abstract {
+private:
+	size_type   m_size;
+	flux_type m_background;
+public:
+	night( const time_type& time );
+
+	size_type size() const;
+	flux_type background( const time_type& time ) const;
+	void accumulate( const time_type& time, const flux_type& mean );
+
+	abstract* clone() const;
+};
+}
+
+namespace filter {
+
+class nonpoisson {
+public:
+	typedef atmos::background::time_type  time_type;
+	typedef atmos::background::flux_type  flux_type;
+	typedef atmos::background::vars_type  vars_type;
+	typedef atmos::background::value_type value_type;
+private:
+	value_type m_threshold;
+	value_type m_nonpoissonD;
+	value_type m_numexp;
+public:
+	nonpoisson( const value_type& threshold, const value_type& nonpoissonD, const value_type& numexp );
+	~nonpoisson();
 
-		bool operator() ( const boost::posix_time::ptime& time, const flux_type& flux, const vars_type& vars ) const;
+	bool operator() ( const time_type& time, const flux_type& flux, const vars_type& vars ) const;
 };
 
+}
+}}
+
 #endif // _BACKGROUND_FACTORY_H

atmos-2.96.9.tar.gz/src/config.h.in -> atmos-2.97.1.tar.gz/src/config.h.in Changed

@@ -29,8 +29,9 @@
 /* Define to 1 if you have the <boost/io/ios_state.hpp> header file. */
 #undef HAVE_BOOST_IO_IOS_STATE_HPP
 
-/* Define to 1 if you have the <boost/iterator.hpp> header file. */
-#undef HAVE_BOOST_ITERATOR_HPP
+/* Define to 1 if you have the <boost/iterator/transform_iterator.hpp> header
+   file. */
+#undef HAVE_BOOST_ITERATOR_TRANSFORM_ITERATOR_HPP
 
 /* Define to 1 if you have the <boost/lexical_cast.hpp> header file. */
 #undef HAVE_BOOST_LEXICAL_CAST_HPP
@@ -38,8 +39,19 @@
 /* Define to 1 if you have the <boost/math/tools/roots.hpp> header file. */
 #undef HAVE_BOOST_MATH_TOOLS_ROOTS_HPP
 
-/* Define to 1 if you have the <boost/numeric/ublas/io.hpp> header file. */
-#undef HAVE_BOOST_NUMERIC_UBLAS_IO_HPP
+/* Define to 1 if you have the <boost/mem_fn.hpp> header file. */
+#undef HAVE_BOOST_MEM_FN_HPP
+
+/* Define to 1 if you have the <boost/mpl/if.hpp> header file. */
+#undef HAVE_BOOST_MPL_IF_HPP
+
+/* Define to 1 if you have the <boost/numeric/ublas/banded.hpp> header file.
+   */
+#undef HAVE_BOOST_NUMERIC_UBLAS_BANDED_HPP
+
+/* Define to 1 if you have the <boost/numeric/ublas/matrix_expression.hpp>
+   header file. */
+#undef HAVE_BOOST_NUMERIC_UBLAS_MATRIX_EXPRESSION_HPP
 
 /* Define to 1 if you have the <boost/numeric/ublas/matrix.hpp> header file.
    */
@@ -57,6 +69,10 @@
    file. */
 #undef HAVE_BOOST_NUMERIC_UBLAS_TRIANGULAR_HPP
 
+/* Define to 1 if you have the <boost/numeric/ublas/vector_expression.hpp>
+   header file. */
+#undef HAVE_BOOST_NUMERIC_UBLAS_VECTOR_EXPRESSION_HPP
+
 /* Define to 1 if you have the <boost/numeric/ublas/vector.hpp> header file.
    */
 #undef HAVE_BOOST_NUMERIC_UBLAS_VECTOR_HPP
@@ -71,6 +87,10 @@
 /* Define to 1 if you have the <boost/tuple/tuple.hpp> header file. */
 #undef HAVE_BOOST_TUPLE_TUPLE_HPP
 
+/* Define to 1 if you have the <boost/type_traits/is_complex.hpp> header file.
+   */
+#undef HAVE_BOOST_TYPE_TRAITS_IS_COMPLEX_HPP
+
 /* Define to 1 if you have the <gsl/gsl_errno.h> header file. */
 #undef HAVE_GSL_GSL_ERRNO_H
 
@@ -118,9 +138,6 @@
 /* Define to 1 if you have the <sys/types.h> header file. */
 #undef HAVE_SYS_TYPES_H
 
-/* Define to 1 if you have the <tr1/tuple> header file. */
-#undef HAVE_TR1_TUPLE
-
 /* Define to 1 if you have the <unistd.h> header file. */
 #undef HAVE_UNISTD_H

atmos-2.97.1.tar.gz/src/covariance_decomposition.h Added

@@ -0,0 +1,137 @@
+/*
+    $Id: covariance_decomposition.h 182 2011-01-31 11:26:03Z matwey $
+    Copyright (C) 2011  Sternberg Astronomical Institute, MSU
+
+    This program is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef _COVARIANCE_DECOMPOSITION_H
+#define _COVARIANCE_DECOMPOSITION_H
+
+#include <cmath>
+
+#include <boost/numeric/ublas/banded.hpp>
+#include <boost/numeric/ublas/matrix.hpp>
+#include <boost/numeric/ublas/symmetric.hpp>
+#include <boost/numeric/ublas/triangular.hpp>
+#include <boost/tuple/tuple.hpp>
+
+#include <lsp/cholesky_decomposition.h>
+
+#include "apply_to_all.h"
+
+namespace utils {
+namespace ublas {
+
+using namespace boost::numeric::ublas;
+
+namespace detail {
+	template<class M> boost::tuple< banded_matrix< typename M::value_type >, symmetric_matrix< typename M::value_type > > correlation_matrix( const M& cov ) {
+		typedef M matrix_type;
+		typedef typename M::value_type value_type;
+		typedef typename M::size_type  size_type;
+
+		symmetric_matrix< value_type > cor(cov);
+		banded_matrix< value_type > sigma = utils::ublas::sqrt( banded_adaptor< symmetric_matrix< value_type > >(cor) );
+		for( typename banded_matrix< value_type >::size_type i = 0; i < sigma.size1(); ++i ) {
+			row(cor,i)    /= sigma(i,i);
+			column(cor,i) /= sigma(i,i);
+		}
+
+		return boost::make_tuple( sigma, cor );
+	}
+	template<class M> triangular_matrix< typename M::value_type > sqrt_matrix_decomposition( const M& cov ) {
+		typedef M matrix_type;
+		typedef typename M::value_type value_type;
+		typedef typename M::size_type  size_type;
+
+		symmetric_matrix< value_type > decomp(cov);
+		lsp::cholesky_decomposition( decomp );
+		triangular_adaptor< symmetric_matrix< value_type > > result(decomp);
+		
+		return result;
+	}
+	template<class M> triangular_matrix< typename M::value_type > inv_sqrt_matrix_decomposition( const M& cov ) {
+		typedef M matrix_type;
+		typedef typename M::value_type value_type;
+		typedef typename M::size_type  size_type;
+
+		triangular_matrix< value_type > triang( sqrt_matrix_decomposition(cov) );
+		matrix< value_type > inv = solve( triang, identity_matrix< value_type >( triang.size1() ), lower_tag() );
+
+		return triangular_adaptor< matrix< value_type > >(inv);
+	}
+	template<class M> triangular_matrix< typename M::value_type > aetkin_weight_matrix( const M& cov, const typename M::value_type& cond_threshold = 10000 ) {
+		typedef M matrix_type;
+		typedef typename M::value_type value_type;
+		typedef typename M::size_type  size_type;
+
+		if( BOOST_UBLAS_SAME( cov.size1(), cov.size2() ) < 1 )
+			return triangular_matrix< value_type >(0,0);
+
+		for( size_type i = 0; i < BOOST_UBLAS_SAME( cov.size1(), cov.size2() ); ++i ) {
+			if( cov(i,i) <= value_type(0) )
+				throw std::logic_error("Non-positive diagonal value");
+		}
+
+		boost::tuple< banded_matrix< value_type >, symmetric_matrix< value_type > > sigma_cor = correlation_matrix(cov);
+		triangular_matrix< value_type > weight = sqrt_matrix_decomposition( sigma_cor.template get<1>() );
+
+		value_type min = std::abs( weight(0,0) );
+		value_type max = std::abs( weight(0,0) );
+		for( typename triangular_matrix< value_type >::size_type i = 1; i < weight.size1(); ++i ) {
+			if( std::abs( weight(i,i) ) > max ) max = std::abs( weight(i,i) );
+			if( std::abs( weight(i,i) ) < min ) min = std::abs( weight(i,i) );
+		}
+
+		value_type cond = max / min;
+		if( !std::isfinite( cond ) || cond > cond_threshold ) {
+			throw std::logic_error( (boost::format("The condition number of covariance matrix is %1%") % cond).str() );
+		}
+
+		matrix< value_type > inv = solve( weight, identity_matrix< value_type >( weight.size1() ), lower_tag() );
+		inv = prod( inv, utils::ublas::inv( sigma_cor.template get<0>() ) );
+		return triangular_adaptor< matrix< value_type > >(inv);
+	}
+}
+
+template<class T> boost::tuple< banded_matrix<T>, symmetric_matrix<T> > correlation_matrix( const symmetric_matrix<T>& cov ) {
+	return detail::correlation_matrix( cov );
+}
+template<class T> triangular_matrix<T> sqrt_matrix_decomposition( const symmetric_matrix<T>& cov ) {
+	return detail::sqrt_matrix_decomposition( cov );
+}
+template<class T> triangular_matrix<T> inv_sqrt_matrix_decomposition( const symmetric_matrix<T>& cov ) {
+	return detail::inv_sqrt_matrix_decomposition( cov );
+}
+template<class T> triangular_matrix<T> aetkin_weight_matrix( const symmetric_matrix<T>& cov ) {
+	return detail::aetkin_weight_matrix( cov );
+}
+
+template<class M> boost::tuple< banded_matrix< typename M::value_type >, symmetric_matrix< typename M::value_type > > correlation_matrix( const symmetric_adaptor<M>& cov ) {
+	return detail::correlation_matrix( cov );
+}
+template<class M> triangular_matrix< typename M::value_type > sqrt_matrix_decomposition( const symmetric_adaptor<M>& cov ) {
+	return detail::sqrt_matrix_decomposition( cov );
+}
+template<class M> triangular_matrix< typename M::value_type > inv_sqrt_matrix_decomposition( const symmetric_adaptor<M>& cov ) {
+	return detail::inv_sqrt_matrix_decomposition( cov );
+}
+template<class M> triangular_matrix< typename M::value_type >  aetkin_weight_matrix( const symmetric_adaptor<M>& cov ) {
+	return detail::aetkin_weight_matrix( cov );
+}
+
+}}
+
+#endif // _COVARIANCE_DECOMPOSITION_H

atmos-2.96.9.tar.gz/src/dimm.cpp -> atmos-2.97.1.tar.gz/src/dimm.cpp Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: dimm.cpp 115 2010-03-14 08:48:08Z matwey $
+    $Id: dimm.cpp 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2010  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -17,13 +17,15 @@
 */
 
 #include <cmath>
-
-#include <boost/bind.hpp>
+#include <boost/mem_fn.hpp>
 
 #include "dimm.h"
 
-dimm::dimm_vars_type dimm::correction( const dimm_vars_type& vars ) {
-	dimm_vars_type ret( vars );
+namespace atmos {
+namespace dimm {
+
+vars_type correction( const vars_type& vars ) {
+	vars_type ret( vars );
 
 	ret(1) = ret(1)*ret(1) - ret(3)*ret(3); // CCD noise removal
 	value_type rho = ret(2)/ret(1);
@@ -33,48 +35,24 @@
 	return ret;
 }
 
-dimm::value_type dimm::turbulence( value_type scale, value_type Knumber, value_type variance ) {
+value_type turbulence( value_type scale, value_type Knumber, value_type variance ) {
 	value_type ret = variance;
 	ret *= scale * scale / 16.7 / Knumber;
 	return ret;
 }
 
-dimm::value_type dimm::seeing( value_type turbulence_power ) {
+value_type seeing( value_type turbulence_power ) {
 	value_type ret = 2.0e7 * std::pow( turbulence_power, 0.6 );
 	return ret;
 }
 
-const dimm::dimm_flux_type& dimm::left_flux( const dimm_sample& x ) {
-	return x.left_flux;
-}
-const dimm::dimm_flux_type& dimm::right_flux( const dimm_sample& x ) {
-	return x.right_flux;
-}
-const dimm::dimm_vars_type& dimm::xvars( const dimm_sample& x ) {
-	return x.xvars;
-}
-const dimm::dimm_vars_type& dimm::yvars( const dimm_sample& x ) {
-	return x.yvars;
-}
-
-dimm_storage::dimm_storage( size_type reserve ):
-	table_storage< dimm::dimm_sample >( reserve ),
-	m_left_flux_column( *this, boost::bind(&dimm::left_flux,_1) ),
-	m_right_flux_column( *this, boost::bind(&dimm::right_flux,_1) ),
-	m_xvars_column( *this, boost::bind(&dimm::xvars,_1) ),
-	m_yvars_column( *this, boost::bind(&dimm::yvars,_1) )
+storage::storage( size_type reserve ):
+	table_storage< detail::sample >( reserve ),
+	m_left_flux_column( *this, boost::mem_fn(&detail::sample::left_flux) ),
+	m_right_flux_column( *this, boost::mem_fn(&detail::sample::right_flux) ),
+	m_xvars_column( *this, boost::mem_fn(&detail::sample::xvars) ),
+	m_yvars_column( *this, boost::mem_fn(&detail::sample::yvars) )
 	{
 }
 
-void dimm_storage::correct_low_vars( ) {
-	double meanx = 0;
-	double meany = 0;
-	for( const_iterator it = begin(); it != end(); ++it ) {
-		meanx += it->xvars(0) / size();
-		meany += it->yvars(0) / size();
-	}
-	for( iterator it = begin(); it != end(); ++it ) {
-		it->xvars(1) += ( it->xvars(0) - meanx ) * ( it->xvars(0) - meanx );
-		it->yvars(1) += ( it->yvars(0) - meany ) * ( it->yvars(0) - meany );
-	}
-}
+}}

atmos-2.96.9.tar.gz/src/dimm.h -> atmos-2.97.1.tar.gz/src/dimm.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: dimm.h 115 2010-03-14 08:48:08Z matwey $
+    $Id: dimm.h 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2010  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -19,71 +19,99 @@
 #ifndef _DIMM_H
 #define _DIMM_H
 
-#include "table_storage.h"
-
 #include <boost/date_time/posix_time/posix_time.hpp>
 #include <boost/numeric/ublas/vector.hpp>
 
+#include "table_storage.h"
+
+namespace atmos {
 namespace dimm {
-	using namespace boost::numeric::ublas;
 
-	typedef boost::posix_time::ptime time_type;
-	
-	typedef double                 value_type;
-	typedef vector<value_type>     vector_type;
-	typedef vector_type::size_type size_type;
-
-	typedef vector_type dimm_flux_type;
-	typedef vector_type dimm_vars_type;
-
-	dimm_vars_type correction( const dimm_vars_type& vars );
-	value_type turbulence( value_type scale, value_type Knumber, value_type variance );
-	value_type seeing( value_type turbulence_power );
-
-	struct dimm_sample {
-		time_type      timestamp;
-		unsigned int   framenum;
-		dimm_flux_type left_flux;
-		dimm_flux_type right_flux;
-		dimm_vars_type xvars;
-		dimm_vars_type yvars;
-
-		dimm_sample( const time_type& time, unsigned int framenum, const dimm_flux_type& left, const dimm_flux_type& right, const dimm_vars_type& xvars, const dimm_flux_type& yvars ):
-			timestamp(time), framenum(framenum), left_flux(left), right_flux(right), xvars(xvars), yvars(yvars) {}
-	};
-
-	const dimm_flux_type& left_flux( const dimm_sample& x );
-	const dimm_flux_type& right_flux( const dimm_sample& x );
-	const dimm_vars_type& xvars( const dimm_sample& x );
-	const dimm_vars_type& yvars( const dimm_sample& x );
+using namespace boost::numeric::ublas;
+
+typedef boost::posix_time::ptime time_type;
+typedef double                   value_type;
+typedef vector<value_type>       vector_type;
+typedef vector_type::size_type   size_type;
+
+typedef vector_type  flux_type;
+typedef vector_type  vars_type;
+typedef unsigned int franenum_type;
+
+vars_type  correction( const vars_type& vars );
+value_type turbulence( value_type scale, value_type Knumber, value_type variance );
+value_type seeing( value_type turbulence_power );
+
+template<class Iterator, class T> void correct_low_vars( Iterator begin, Iterator end, const T& mean ) {
+	typedef Iterator iterator_type;
+	typedef T value_type;
+
+	for( iterator_type it = begin; it != end; ++it ) {
+		(*it)(1) += ((*it)(0)-mean)*((*it)(0)-mean);
+	}
+}
+
+namespace detail {
+
+struct sample {
+	typedef atmos::dimm::time_type     timestamp_type;
+	typedef atmos::dimm::franenum_type framenum_type;
+	typedef atmos::dimm::flux_type     flux_type;
+	typedef atmos::dimm::vars_type     vars_type;
+		
+	timestamp_type timestamp;
+	framenum_type  framenum;
+	flux_type      left_flux;
+	flux_type      right_flux;
+	vars_type      xvars;
+	vars_type      yvars;
+
+	sample( const timestamp_type& time,
+		     const framenum_type& framenum,
+		     const flux_type& left,
+		     const flux_type& right,
+		     const vars_type& xvars,
+		     const vars_type& yvars ):
+		timestamp(time), framenum(framenum), left_flux(left), right_flux(right), xvars(xvars), yvars(yvars) {}
 };
+}
 
-class dimm_storage:
-	public table_storage< dimm::dimm_sample > {
+class storage:
+	public table_storage< detail::sample > {
 public:
-	typedef table_storage_column< dimm::dimm_sample, boost::function1< dimm::dimm_flux_type, dimm::dimm_sample > > left_flux_column_type;
-	typedef table_storage_column< dimm::dimm_sample, boost::function1< dimm::dimm_flux_type, dimm::dimm_sample > > right_flux_column_type;
-	typedef table_storage_column< dimm::dimm_sample, boost::function1< dimm::dimm_vars_type, dimm::dimm_sample > > xvars_column_type;
-	typedef table_storage_column< dimm::dimm_sample, boost::function1< dimm::dimm_vars_type, dimm::dimm_sample > > yvars_column_type;
+	typedef storage self_type;
+	typedef table_storage_column_proxy< self_type, detail::sample::flux_type > left_flux_column_type;
+	typedef table_storage_column_proxy< self_type, detail::sample::flux_type > right_flux_column_type;
+	typedef table_storage_column_proxy< self_type, detail::sample::vars_type > xvars_column_type;
+	typedef table_storage_column_proxy< self_type, detail::sample::vars_type > yvars_column_type;
 private:
 	left_flux_column_type  m_left_flux_column;
 	right_flux_column_type m_right_flux_column;
 	xvars_column_type      m_xvars_column;
 	yvars_column_type      m_yvars_column;
 public:
-	dimm_storage( size_type reserve = 60 );
+	storage( size_type reserve = 60 );
 
 	const left_flux_column_type&  left_flux()  const { return m_left_flux_column; }
 	const right_flux_column_type& right_flux() const { return m_right_flux_column; }
 	const xvars_column_type&      xvars()      const { return m_xvars_column; }
 	const yvars_column_type&      yvars()      const { return m_yvars_column; }
+	left_flux_column_type&  left_flux()  { return m_left_flux_column; }
+	right_flux_column_type& right_flux() { return m_right_flux_column; }
+	xvars_column_type&      xvars()      { return m_xvars_column; }
+	yvars_column_type&      yvars()      { return m_yvars_column; }
+
+	void push_back( const time_type& time,
+			const franenum_type& framenum,
+			const flux_type& left,
+			const flux_type& right,
+			const vars_type& xvars,
+			const vars_type& yvars ) {
 
-	void push_back( const dimm::time_type& time, unsigned int framenum, const dimm::dimm_flux_type& left, const dimm::dimm_flux_type& right, const dimm::dimm_vars_type& xvars, const dimm::dimm_flux_type& yvars ) {
-		// FIXME: dimm_sample_ref
-		table_storage< dimm::dimm_sample >::push_back( dimm::dimm_sample( time, framenum, left, right, xvars, yvars ) );
+		table_storage< detail::sample >::push_back( detail::sample( time, framenum, left, right, xvars, yvars ) );
 	}
-
-	void correct_low_vars();
 };
 
+}}
+
 #endif // _DIMM_H

atmos-2.96.9.tar.gz/src/fft.h -> atmos-2.97.1.tar.gz/src/fft.h Changed

atmos-2.96.9.tar.gz/src/file_input.cpp -> atmos-2.97.1.tar.gz/src/file_input.cpp Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: file_input.cpp 142 2010-04-03 13:50:39Z matwey $
+    $Id: file_input.cpp 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -16,34 +16,34 @@
     along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-#include <cassert>
-#include <stdexcept>
-#include <sstream>
 #include <iostream>
+#include <sstream>
+#include <stdexcept>
 
-#include <ctime>
-
-#include <boost/io/ios_state.hpp>
-#include <boost/numeric/ublas/vector.hpp>
-
+#include <boost/algorithm/string.hpp>
 #include <boost/date_time/posix_time/posix_time.hpp>
 
 #include "atmos.h"
 #include "file_input.h"
+#include "io.h"
+
+namespace atmos {
+namespace input {
 
-file_input::file_input( std::istream& stm ):
+file::file( std::istream& stm ):
 	m_linenum(0),
 	m_stm(stm),
 	m_cached( false ),
 	m_locale( m_stm.getloc(), new boost::posix_time::time_input_facet("%Y-%m-%d %H:%M:%S") ) {
 }
-file_input::~file_input() {
+file::~file() {
 }
-bool file_input::unshift() throw() {
+
+bool file::unshift() throw() {
 	m_cached = true;
 	return true;
 }
-bool file_input::shift() {
+bool file::shift() {
 	struct linenum_guard {
 			linenum_guard( size_type& linenum ): m_linenum(linenum) {}
 			~linenum_guard() { m_linenum++; }
@@ -65,16 +65,15 @@
 	
 	return false;
 }
-file_input::size_type file_input::line() const throw() {
+size_type file::line() const throw() {
 	return m_linenum;
 }
-//FIXME: use boost string_algo
-std::string file_input::strip_spaces( const std::string& str ) {
+std::string file::strip_spaces( const std::string& str ) {
 	std::string ret(str);
 	boost::trim_if( ret, boost::is_space() );
 	return ret;
 }
-void file_input::parse_line( size_type linenum, const std::string& line ) {
+void file::parse_line( size_type linenum, const std::string& line ) {
 	if( line.size() == 0 )
 		return; // skip empty line
 
@@ -96,10 +95,10 @@
 		throw failure("Bad line format");
 	}
 }
-void file_input::parse_comment_line( std::istringstream& stm ) {
+void file::parse_comment_line( std::istringstream& stm ) {
 	on_comment( stm.str() );
 }
-void file_input::parse_param_line( std::istringstream& stm ) {
+void file::parse_param_line( std::istringstream& stm ) {
 	static time_type timestamp; // NOTE: thread unsafe
 	std::string key,value;
 	
@@ -109,7 +108,7 @@
 
 	on_parameter(timestamp, strip_spaces(key), strip_spaces(value));
 }
-void file_input::parse_object_line( std::istringstream& stm ) {
+void file::parse_object_line( std::istringstream& stm ) {
 	static time_type   timestamp; // NOTE: thread unsafe
 	static object_type obj;       // NOTE: thread unsafe
 
@@ -117,7 +116,7 @@
 
 	on_object( timestamp, obj );
 }
-void file_input::parse_mode_line( std::istringstream& stm ) {
+void file::parse_mode_line( std::istringstream& stm ) {
 	static time_type timestamp; // NOTE: thread unsafe
 	
 	stm >> timestamp;
@@ -132,7 +131,7 @@
 		throw failure("Unknown mode");
 	}
 }
-void file_input::parse_measurement_line( std::istringstream& stm ) {
+void file::parse_measurement_line( std::istringstream& stm ) {
 	static time_type timestamp;               // NOTE: thread unsafe
 	static flux_type flux( atmos::channels ); // NOTE: thread unsafe
 	static vars_type vars( atmos::channels, atmos::channels ); // NOTE: thread unsafe
@@ -147,7 +146,7 @@
 
 	on_measurements( timestamp, flux, vars, rho1, rho2 );
 }
-void file_input::parse_dimm_line( std::istringstream& stm ) {
+void file::parse_dimm_line( std::istringstream& stm ) {
 	static time_type timestamp;  // NOTE: thread unsafe
 	static unsigned int framenum;
 	static dimm_flux_type left_flux( 3 );  // NOTE: thread unsafe
@@ -169,7 +168,10 @@
 	on_dimm( timestamp, framenum, left_flux, right_flux, xvars, yvars );
 }
 
-file_input::failure::failure( const std::string& what):
+file::failure::failure( const std::string& what):
 	std::runtime_error( what ) {
 }
 
+}}
+
+

atmos-2.96.9.tar.gz/src/file_input.h -> atmos-2.97.1.tar.gz/src/file_input.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: file_input.h 142 2010-04-03 13:50:39Z matwey $
+    $Id: file_input.h 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -23,37 +23,49 @@
 #include <stdexcept>
 #include <string>
 
-#include <boost/algorithm/string.hpp>
 #include <boost/date_time/posix_time/posix_time.hpp>
-
-#include <boost/numeric/ublas/vector.hpp>
 #include <boost/numeric/ublas/symmetric.hpp>
+#include <boost/numeric/ublas/vector.hpp>
 
-#include "io.h"
+#include "atmos.h" // size_type
+#include "typesdef.h" // object
 
-using namespace boost::numeric;
+namespace atmos {
+namespace input {
 
-class file_input {
-public:
-	typedef size_t size_type;
-	
-	typedef boost::posix_time::ptime time_type;
-	typedef ublas::vector<double>    flux_type;
-	typedef ublas::symmetric_matrix<double> vars_type;
-	typedef ublas::symmetric_matrix<double> rho1_type;
-	typedef ublas::vector<double>    rho2_type;
-	typedef object                   object_type;
+using namespace boost::numeric::ublas;
+
+typedef boost::posix_time::ptime time_type;
+typedef double               value_type;
+typedef vector< value_type > vector_type;
+typedef symmetric_matrix< value_type > symmetric_matrix_type;
+
+typedef vector_type flux_type;
+typedef symmetric_matrix_type vars_type;
+typedef symmetric_matrix_type rho1_type;
+typedef vector_type rho2_type;
+typedef vector_type dimm_flux_type;
+typedef vector_type dimm_vars_type;
+typedef object object_type;
 
-	typedef ublas::vector<double>    dimm_flux_type;
-	typedef ublas::vector<double>    dimm_vars_type;
+class file {
+public:
+	typedef atmos::size_type        size_type;
+	typedef atmos::input::time_type time_type;
+	typedef atmos::input::flux_type flux_type;
+	typedef atmos::input::vars_type vars_type;
+	typedef atmos::input::rho1_type rho1_type;
+	typedef atmos::input::rho2_type rho2_type;
+	typedef atmos::input::dimm_flux_type dimm_flux_type;
+	typedef atmos::input::dimm_vars_type dimm_vars_type;
+	typedef atmos::input::object_type object_type;
 private:
 	size_type m_linenum;
 	std::istream& m_stm;
 	std::string m_line;
 	bool m_cached;
 	std::locale m_locale;
-	private:
-	
+private:
 	static std::string strip_spaces( const std::string& str );
 
 	void parse_line( size_type linenum, const std::string& line );
@@ -63,7 +75,6 @@
 	void parse_mode_line( std::istringstream& stm );
 	void parse_measurement_line( std::istringstream& stm );
 	void parse_dimm_line( std::istringstream& stm );
-
 protected:
 	bool unshift() throw();
 public:
@@ -73,11 +84,10 @@
 		failure( const std::string& what );
 	};
 public:
-	file_input( std::istream& stm );
-	virtual ~file_input();
+	file( std::istream& stm );
+	virtual ~file();
 
 	bool shift();
-
 	size_type line() const throw() ;
 public:
 	virtual void on_comment( const std::string& line ) = 0;
@@ -90,4 +100,6 @@
 	virtual void on_dimm( const time_type& time, unsigned int framenum, const dimm_flux_type& left, const dimm_flux_type& right, const dimm_vars_type& xvars, const dimm_flux_type& yvars ) = 0;
 };
 
+}}
+
 #endif // _FILE_INPUT_H

atmos-2.96.9.tar.gz/src/file_output.cpp -> atmos-2.97.1.tar.gz/src/file_output.cpp Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: file_output.cpp 108 2010-03-07 11:40:13Z matwey $
+    $Id: file_output.cpp 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -16,30 +16,36 @@
     along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
 
+#include <boost/date_time/posix_time/posix_time.hpp>
+#include <boost/format.hpp>
 #include <boost/io/ios_state.hpp>
 
+#include "atmos.h"
 #include "file_output.h"
-#include "integrals.h"
+#include "io.h"
 
-file_output::file_output( std::ostream& stm ): m_stm(stm) {
+namespace atmos {
+namespace output {
+
+file::file( std::ostream& stm ): m_stm(stm) {
 	boost::posix_time::time_facet* facet(new boost::posix_time::time_facet("%Y-%m-%d %H:%M:%S"));
 	stm.imbue(std::locale(stm.getloc(), facet));
 }
-file_output::~file_output() {
+file::~file() {
 
 }
-void file_output::comment( const std::string& comment ) {
+void file::comment( const std::string& comment ) {
 	m_stm << "# " << comment << std::endl;
 }
-void file_output::normal( const time_type& time, double mz ) {
+void file::normal( const time_type& time, double mz ) {
 	boost::io::ios_all_saver asaver( m_stm );
 	m_stm << "M " << time << " Normal " << std::fixed << std::setprecision(3) << mz << std::endl;
 }
-void file_output::background( const time_type& time ) {
+void file::background( const time_type& time ) {
 	boost::io::ios_all_saver asaver( m_stm );
 	m_stm << "M " << time << " Background measurements" << std::endl;
 }
-void file_output::atmosphere( const time_type& time, double free_seeing, double free_seeing_err, double dimm_seeing, double dimm_seeing_err, double dimm_turbulence, double dimm_turbulence_err, double isoplanatic_angle, double isoplanatic_angle_err, double altitude_eff, double altitude_eff_err, double time_constant, double time_constant_err, double moment_0, double moment_0_err, double moment_2, double moment_2_err ) {
+void file::atmosphere( const time_type& time, double free_seeing, double free_seeing_err, double dimm_seeing, double dimm_seeing_err, double dimm_turbulence, double dimm_turbulence_err, double isoplanatic_angle, double isoplanatic_angle_err, double altitude_eff, double altitude_eff_err, double time_constant, double time_constant_err, double moment_0, double moment_0_err, double moment_2, double moment_2_err ) {
 	static datum skipped;
 	boost::io::ios_all_saver asaver( m_stm );
 	m_stm << "A " << time << " "
@@ -50,7 +56,7 @@
 		<< std::scientific << datum( moment_2, moment_2_err/moment_2 ) << " "
 		<< std::fixed << datum( time_constant, time_constant_err/time_constant ) << std::endl;
 }
-void file_output::profile( const time_type& time, const std::string& method, double cn2seeing, double chi2, const ublas::vector<double>& cn2, const ublas::vector<double>& ecn2, const ublas::vector<double>& grid ){
+void file::profile( const time_type& time, const std::string& method, double cn2seeing, double chi2, const vector_type& cn2, const vector_type& ecn2, const vector_type& grid ){
 	assert( cn2.size() == ecn2.size() );
 	assert( cn2.size() == grid.size() );
 
@@ -58,27 +64,27 @@
 	m_stm << "T " << time << " " << method << " "
 		<< std::setw(2) << cn2.size() << " " << std::setprecision(2)
 		<< std::fixed << chi2 << " " << cn2seeing;
-	for( ublas::vector<double>::size_type i = 0; i < cn2.size(); i++)
+	for( vector_type::size_type i = 0; i < cn2.size(); i++)
 		m_stm << " " << std::fixed << grid(i)
 			<< std::scientific  << " " << cn2(i) << " " << ecn2(i);
 	m_stm << std::endl;
 }
-void file_output::residual( const time_type& time, const std::string& method, const ublas::vector<double>& o_csi ) {
+void file::residual( const time_type& time, const std::string& method, const vector_type& o_csi ) {
 	boost::io::ios_all_saver asaver( m_stm );
 	m_stm << "R " << time << " " << method << " "
 		<< std::setprecision(3) << std::fixed;
-	for( ublas::vector<double>::size_type i = 0; i < o_csi.size(); i++ )
+	for( vector_type::size_type i = 0; i < o_csi.size(); i++ )
 		m_stm << "  " << o_csi(i);
 	m_stm << std::endl;
 }
-void file_output::fluxes( const time_type& time, const ublas::vector<double>& result_flux ) {
+void file::fluxes( const time_type& time, const vector_type& result_flux ) {
 	boost::io::ios_all_saver asaver( m_stm );
 	assert( result_flux.size() >= 4 );
 	m_stm << "F " << time << " " << std::fixed << datum( result_flux(0) ) << ' ' << datum( result_flux(1) ) << ' ' << datum( result_flux(2) ) << ' ' << datum( result_flux(3) ) << std::endl;
 }
-void file_output::fluxes( const time_type& time, const ublas::vector<double>& result_flux, const ublas::vector<double>& result_flux_err ) {
+void file::fluxes( const time_type& time, const vector_type& result_flux, const vector_type& result_flux_err ) {
 	boost::io::ios_all_saver asaver( m_stm );
-	assert( result_flux.size() >= 4 );
+	assert( result_flux.size() == atmos::channels );
 	assert( result_flux.size() == result_flux_err.size() );
 	m_stm << "F " << time << " " << std::fixed
 	<< datum( result_flux(0), result_flux_err(0)/result_flux(0) ) << ' '
@@ -86,12 +92,12 @@
 	<< datum( result_flux(2), result_flux_err(2)/result_flux(2) ) << ' '
 	<< datum( result_flux(3), result_flux_err(3)/result_flux(3) ) << std::endl;
 }
-void file_output::indices( const time_type& time, const ublas::symmetric_matrix<double>& result_indx, const ublas::symmetric_matrix<double>& result_indx_err, const ublas::vector<double>& result_desi, const ublas::vector<double>& result_desi_err ) {
+void file::indices( const time_type& time, const symmetric_matrix_type& result_indx, const symmetric_matrix_type& result_indx_err, const vector_type& result_desi, const vector_type& result_desi_err ) {
 	boost::io::ios_all_saver asaver( m_stm );
 	ios_indic_saver isaver( m_stm );
-	assert( result_desi.size() >= 4 );
-	assert( result_indx.size1() >= 4 );
-	assert( result_indx.size2() >= 4 );
+	assert( result_desi.size() == atmos::channels );
+	assert( result_indx.size1() == atmos::channels );
+	assert( result_indx.size2() == atmos::channels );
 	assert( result_desi.size() == result_desi_err.size() );
 	assert( result_indx.size1() == result_indx_err.size1() );
 	assert( result_indx.size1() == result_indx_err.size2() );
@@ -102,11 +108,13 @@
 		datum( result_indx(2,3), result_indx_err(2,3)/result_indx(2,3) ) <<
 		datum( result_desi(0), result_desi_err(0)/result_desi(0) ) << datum( result_desi(1), result_desi_err(1)/result_desi(1) ) << datum( result_desi(2), result_desi_err(2)/result_desi(2) ) << datum( result_desi(3), result_desi_err(3)/result_desi(3) ) << std::endl;
 }
-void file_output::parameter( const time_type& time, const std::string& key, const std::string& value ) {
+void file::parameter( const time_type& time, const std::string& key, const std::string& value ) {
 	boost::io::ios_all_saver asaver( m_stm );
 	m_stm << "P " << time << " " << key << " = " << value << std::endl;
 }
-void file_output::objectinfo( const time_type& time, const object_type& obj ){
+void file::objectinfo( const time_type& time, const object_type& obj ){
 	boost::io::ios_all_saver asaver( m_stm );
 	m_stm << "O " << time << " " << obj << std::endl;
 };
+
+}}

atmos-2.96.9.tar.gz/src/file_output.h -> atmos-2.97.1.tar.gz/src/file_output.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: file_output.h 140 2010-04-02 17:13:12Z matwey $
+    $Id: file_output.h 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -22,42 +22,46 @@
 #include <iosfwd>
 #include <string>
 
-#include <boost/format.hpp>
-#include <boost/io/ios_state.hpp>
+#include <boost/date_time/posix_time/posix_time.hpp>
 #include <boost/numeric/ublas/symmetric.hpp>
 #include <boost/numeric/ublas/vector.hpp>
 
-#include <boost/date_time/posix_time/posix_time.hpp>
+#include "typesdef.h" // object
 
-#include "io.h"
+namespace atmos {
+namespace output {
 
-using namespace boost::numeric;
+using namespace boost::numeric::ublas;
 
-class file_output {
-public:
-	typedef size_t size_type;
+typedef boost::posix_time::ptime time_type;
+typedef double               value_type;
+typedef vector< value_type > vector_type;
+typedef symmetric_matrix< value_type > symmetric_matrix_type;
 
-	typedef boost::posix_time::ptime time_type;
-	typedef object                   object_type;
+typedef object object_type;
+
+class file {
 private:
 	std::ostream& m_stm;
 public:
-	file_output( std::ostream& stm );
-	~file_output();
+	file( std::ostream& stm );
+	~file();
 
 	void comment( const std::string& comment );
 	void normal( const time_type& time, double mz );
 	void background( const time_type& time );
 	void atmosphere( const time_type& time, double free_seeing, double free_seeing_err, double dimm_seeing, double dimm_seeing_err, double dimm_turbulence, double dimm_turbulence_err, double isoplanatic_angle, double isoplanatic_angle_err, double altitude_eff, double altitude_eff_err, double time_constant, double time_constant_err, double moment_0, double moment_0_err, double moment_2, double moment_2_err );
-	void profile( const time_type& time, const std::string& method, double cn2seeing, double chi2, const ublas::vector<double>& cn2, const ublas::vector<double>& ecn2, const ublas::vector<double>& grid );
-	void residual( const time_type& time, const std::string& method, const ublas::vector<double>& o_csi );
-	void fluxes( const time_type& time, const ublas::vector<double>& result_flux );
-	void fluxes( const time_type& time, const ublas::vector<double>& result_flux, const ublas::vector<double>& result_flux_err );
-	void indices( const time_type& time, const ublas::symmetric_matrix<double>& result_indx, const ublas::symmetric_matrix<double>& result_indx_err, const ublas::vector<double>& result_desi, const ublas::vector<double>& result_desi_err );
+	void profile( const time_type& time, const std::string& method, double cn2seeing, double chi2, const vector_type& cn2, const vector_type& ecn2, const vector_type& grid );
+	void residual( const time_type& time, const std::string& method, const vector_type& o_csi );
+	void fluxes( const time_type& time, const vector_type& result_flux );
+	void fluxes( const time_type& time, const vector_type& result_flux, const vector_type& result_flux_err );
+	void indices( const time_type& time, const symmetric_matrix_type& result_indx, const symmetric_matrix_type& result_indx_err, const vector_type& result_desi, const vector_type& result_desi_err );
 	void parameter( const time_type& time, const std::string& key, const std::string& value );
 	void objectinfo( const time_type& time, const object_type& obj );
 
 	std::ostream& operator()() { return m_stm; }
 };
 
+}}
+
 #endif // _FILE_OUTPUT_H

atmos-2.96.9.tar.gz/src/gsl_wrap.h -> atmos-2.97.1.tar.gz/src/gsl_wrap.h Changed

atmos-2.96.9.tar.gz/src/indices.cpp -> atmos-2.97.1.tar.gz/src/indices.cpp Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: indices.cpp 147 2010-04-17 09:25:25Z matwey $
+    $Id: indices.cpp 181 2011-01-31 10:07:58Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -16,28 +16,27 @@
     along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-#include <algorithm>
-#include <numeric>
+#include <stdexcept>
 
 #include <boost/format.hpp>
 #include <boost/numeric/ublas/matrix_proxy.hpp>
 #include <boost/numeric/ublas/vector_proxy.hpp>
 
-#include "atmos.h"
 #include "indices.h"
 #include "io.h"
+#include "photometry.h"
 #include "utils.h"
 
-#include "photometry.h"
+namespace atmos {
+namespace indices {
 
-indices::flux_type indices::flux_nonlinearity_correction ( const flux_type& flux, const vector_type& deadtime ) {
+flux_type flux_nonlinearity_correction ( const flux_type& flux, const vector_type& deadtime ) {
 	return photometry::flux_nonlinearity_correction( flux, deadtime );
 }
-indices::flux_type indices::flux_background_correction( const flux_type& flux, const vector_type& scatter, const vector_type& background ) {
+flux_type flux_background_correction( const flux_type& flux, const vector_type& scatter, const vector_type& background ) {
 	return photometry::flux_scatter_correction( flux, scatter ) - background;
 }
-boost::tuple< indices::vars_type, indices::rho1_type, indices::rho2_type >
-	indices::vars_correction( const flux_type& flux, const vars_type& vars, const rho1_type& rho1, const rho2_type& rho2, const vector_type& deadtime, const vector_type& npoisson ) {
+boost::tuple< vars_type, rho1_type, rho2_type > vars_correction( const flux_type& flux, const vars_type& vars, const rho1_type& rho1, const rho2_type& rho2, const vector_type& deadtime, const vector_type& npoisson ) {
 	using namespace utils::ublas;
 
 	typedef matrix_vector_slice< vars_type > diagonal_type;
@@ -65,11 +64,10 @@
 	return ret;
 }
 
-boost::tuple< indices::indices_type, indices::indesis_type >
-	indices::calculate( const flux_type& flux, const vars_type& vars, const rho1_type& rho1, const rho2_type& rho2 ) {
+boost::tuple< indices_type, indesis_type > calculate( const flux_type& flux, const vars_type& vars, const rho1_type& rho1, const rho2_type& rho2 ) {
 	using namespace utils::ublas;
 
-	typedef matrix_vector_slice< indices_type >    indices_diagonal_type;
+	typedef matrix_vector_slice< const indices_type > indices_diagonal_type;
 	typedef matrix_vector_slice< const rho1_type > rho1_diagonal_type;
 	typedef boost::tuple< indices_type, indesis_type> result_type;
 
@@ -82,29 +80,27 @@
 	
 	rho1_diagonal_type rho1_diag = diag( rho1 );
 
-	ret.get<0>() = element_div( vars*1.25 - rho1*0.25, outer_prod( flux, flux ) );
-
-	indices_diagonal_type variances = diag( ret.get<0>() );
-
-	ret.get<1>() = ( variances*3.0 + element_div( rho2 - rho1_diag * 4.0, element_prod(flux,flux) ) ) / 4.5;
-
-	ret.get<0>() = utils::ublas::log( xp1( ret.get<0>() ) );
-
-	triangle( ret.get<0>(), unit_lower() ) = triangle(
-			outer_prod( variances, scalar_vector< value_type >( variances.size(), 1.0 ) )
-			+ outer_prod( scalar_vector< value_type >( variances.size(), 1.0 ), variances )
-			- triangle( ret.get<0>(), unit_lower() ) * 2.0,
-		unit_lower() );
+	ret.get<0>() = utils::ublas::log( xp1( element_div( vars*1.17 - rho1*0.17, outer_prod( flux, flux ) ) ) );
+	ret.get<1>() = ( diag( element_div( vars, outer_prod( flux, flux ) ) )*3.0 + element_div( rho2 - rho1_diag * 4.0, element_prod(flux,flux) ) ) / 4.5;
 
 	return ret;
 }
 
-indices_flux_threshold_filter::indices_flux_threshold_filter( value_type threshold ):
+storage::storage( size_type reserve ):
+	table_storage< detail::sample >( reserve ),
+	m_flux_column( *this, boost::mem_fn(&detail::sample::flux) ),
+	m_indices_column( *this, boost::mem_fn(&detail::sample::indices) ),
+	m_indesis_column( *this, boost::mem_fn(&detail::sample::indesis) ) {
+}
+
+namespace filter {
+	
+flux_threshold::flux_threshold( value_type threshold ):
 	m_threshold( threshold ) {
 }
-indices_flux_threshold_filter::~indices_flux_threshold_filter() {
+flux_threshold::~flux_threshold() {
 }
-bool indices_flux_threshold_filter::operator() ( const flux_type& flux ) const {
+bool flux_threshold::operator() ( const flux_type& flux ) const {
 	if( *std::min_element( flux.begin(), flux.end() ) < m_threshold ) {
 		throw std::logic_error( (boost::format("Object flux %.2f is under the threshold %.2f") % flux % m_threshold).str() );
 		return true;
@@ -112,19 +108,5 @@
 	return false;
 }
 
-const indices::flux_type& indices::flux( const indices::indices_sample& x ) {
-	return x.flux;
-}
-const indices::indices_type& indices::indexes( const indices::indices_sample& x ) {
-	return x.indices;
-}
-const indices::indesis_type& indices::indesis( const indices::indices_sample& x ) {
-	return x.indesis;
-}
-
-indices_storage::indices_storage( size_type reserve ):
-	table_storage< indices::indices_sample >( reserve ),
-	m_flux_column( *this, boost::bind(&indices::flux,_1) ),
-	m_indices_column( *this, boost::bind(&indices::indexes,_1) ),
-	m_indesis_column( *this, boost::bind(&indices::indesis,_1) ) {
 }
+}}

atmos-2.96.9.tar.gz/src/indices.h -> atmos-2.97.1.tar.gz/src/indices.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: indices.h 147 2010-04-17 09:25:25Z matwey $
+    $Id: indices.h 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -19,29 +19,23 @@
 #ifndef _INDICES_H
 #define _INDICES_H
 
-#include <stdexcept>
-#include <list>
-#include <vector>
-#include <iterator>
-
-#include <boost/bind.hpp>
 #include <boost/date_time/posix_time/posix_time.hpp>
-#include <boost/numeric/ublas/matrix.hpp>
 #include <boost/numeric/ublas/symmetric.hpp>
 #include <boost/numeric/ublas/vector.hpp>
 #include <boost/tuple/tuple.hpp>
 
 #include "table_storage.h"
 
+namespace atmos {
 namespace indices {
 	using namespace boost::numeric::ublas;
 
 	typedef boost::posix_time::ptime time_type;
-	
-	typedef double                value_type;
-	typedef vector< value_type >  vector_type;
-	typedef vector_type::size_type size_type;
+	typedef double                   value_type;
+	typedef vector< value_type >     vector_type;
+	typedef vector_type::size_type   size_type;
 	typedef symmetric_matrix< value_type > symmetric_matrix_type;
+	
 	typedef vector_type           flux_type;
 	typedef symmetric_matrix_type vars_type;
 	typedef symmetric_matrix_type rho1_type;
@@ -55,27 +49,33 @@
 		vars_correction( const flux_type& flux, const vars_type& vars, const rho1_type& rho1, const rho2_type& rho2, const vector_type& deadtime, const vector_type& npoisson );
 	boost::tuple< indices_type, indesis_type > calculate( const flux_type& flux, const vars_type& vars, const rho1_type& rho1, const rho2_type& rho2 );
 
-	struct indices_sample {
-		time_type    timestamp;
-		flux_type    flux;
-		indices_type indices;
-		indesis_type indesis;
-
-		indices_sample( const time_type& time, const flux_type& flux, const indices_type& indices, const indesis_type& indesis ):
+namespace detail {
+	struct sample {
+		typedef atmos::indices::time_type    timestamp_type;
+		typedef atmos::indices::flux_type    flux_type;
+		typedef atmos::indices::indices_type indices_type;
+		typedef atmos::indices::indesis_type indesis_type;
+
+		timestamp_type timestamp;
+		flux_type      flux;
+		indices_type   indices;
+		indesis_type   indesis;
+
+		sample( const timestamp_type& time,
+			const flux_type& flux,
+			const indices_type& indices,
+			const indesis_type& indesis ):
 			timestamp(time), flux(flux), indices(indices), indesis(indesis) {}
 	};
+}
 
-	const flux_type& flux( const indices_sample& x );
-	const indices_type& indexes( const indices_sample& x );
-	const indesis_type& indesis( const indices_sample& x );
-};
-
-class indices_storage:
-	public table_storage< indices::indices_sample > {
+class storage:
+	public table_storage< detail::sample > {
 public:
-	typedef table_storage_column< indices::indices_sample, boost::function1< indices::flux_type, indices::indices_sample > >    flux_column_type;
-	typedef table_storage_column< indices::indices_sample, boost::function1< indices::indices_type, indices::indices_sample > > indices_column_type;
-	typedef table_storage_column< indices::indices_sample, boost::function1< indices::indesis_type, indices::indices_sample > > indesis_column_type;
+	typedef storage self_type;
+	typedef table_storage_column_getter< self_type, detail::sample::flux_type >    flux_column_type;
+	typedef table_storage_column_getter< self_type, detail::sample::indices_type > indices_column_type;
+	typedef table_storage_column_getter< self_type, detail::sample::indesis_type > indesis_column_type;
 
 private:
 	flux_column_type    m_flux_column;
@@ -83,28 +83,35 @@
 	indesis_column_type m_indesis_column;
 
 public:
-	indices_storage( size_type reserve = 60 );
+	storage( size_type reserve = 60 );
 
 	const flux_column_type&    flux()    const { return m_flux_column; }
 	const indices_column_type& indices() const { return m_indices_column; }
 	const indesis_column_type& indesis() const { return m_indesis_column; }
 
-	void push_back( const indices::time_type& time, const indices::flux_type& flux, const indices::indices_type& indices, const indices::indesis_type& indesis ) {
-		table_storage< indices::indices_sample >::push_back( indices::indices_sample( time, flux, indices, indesis ) );
+	void push_back( const time_type& time,
+			const flux_type& flux,
+			const indices_type& indices,
+			const indesis_type& indesis ) {
+		
+		table_storage< detail::sample >::push_back( detail::sample( time, flux, indices, indesis ) );
 	}
 };
 
-class indices_flux_threshold_filter {
+namespace filter {
+	class flux_threshold {
 	public:
-		typedef indices::value_type value_type;
-		typedef indices::flux_type  flux_type;
+		typedef atmos::indices::flux_type  flux_type;
+		typedef atmos::indices::value_type value_type;
 	private:
 		value_type m_threshold;
 	public:
-		indices_flux_threshold_filter( value_type threshold );
-		~indices_flux_threshold_filter();
+		flux_threshold( value_type threshold );
+		~flux_threshold();
 
 		bool operator() ( const flux_type& flux ) const;
-};
+	};
+}
+}}
 
 #endif // _INDICES_H

atmos-2.96.9.tar.gz/src/integrals.cpp -> atmos-2.97.1.tar.gz/src/integrals.cpp Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: integrals.cpp 158 2010-09-19 08:38:25Z matwey $
+    $Id: integrals.cpp 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -17,18 +17,24 @@
 */
 
 #include <cmath>
-#include <algorithm>
-#include <numeric>
-#include <boost/bind.hpp>
+
 #include <boost/format.hpp>
+#include <boost/mem_fn.hpp>
 #include <boost/numeric/ublas/matrix_proxy.hpp>
 #include <boost/numeric/ublas/vector_proxy.hpp>
 #include <boost/numeric/ublas/symmetric.hpp>
+
 #include <lsp/singular_decomposition.h>
 
+#include "atmos.h"
 #include "apply_to_all.h"
-#include "indices.h"
 #include "integrals.h"
+#include "vector_generator.h"
+#include "utils.h"
+#include "weif.h"
+
+namespace atmos {
+namespace integrals {
 
 const double KM2M = 1000;
 const double POWSEE = 0 ;       // Power of Cn2 moment for seeing
@@ -64,11 +70,12 @@
 
 const double SPOW_TC  = -0.6;   // Power of DESI for atm. time constant
 
-integrals_factory::dcf_matrix_type integrals_factory::calculate( const std::string& spectral_type, const weight_type& awf ) {
+functor_factory::dcf_matrix_type functor_factory::calculate( const std::string& spectral_type, const weight_type& awf ) {
 	using namespace boost::numeric::ublas;
+	using namespace utils::algo;
 	using namespace utils::ublas;
 
-	const vector<double> grid( utils::ublas::vector_generator< lg_generator > ( lg_generator(0.0, 0.04, 0.25), 50 ) );
+	const vector<double> grid( vector_generator< lg_generator > ( lg_generator(0.0, 0.04, 0.25), 50 ) );
 
 	typedef dcf_matrix_type            result_type;
 	typedef dcf_matrix_type::size_type size_type;
@@ -81,7 +88,7 @@
 	assert( lower != grid.end() );
 	assert( higher != grid.rend() );
 
-	result_type res( npower, 12-2/*dimm*/ );
+	result_type res( npower, atmos::indices_size );
 
 	vector_range< const vector<double> > h( grid, range( lower.index(), higher.index()+1 ) );
 	matrix< double > wfv = vector_of_vector_to_matrix( project( grid, awf) );
@@ -94,7 +101,7 @@
 		vector<double> alt    = utils::ublas::exp( log( h )*moment_power ) * scale;
 		vector<double> weight = utils::ublas::exp( log( h )*weight_power );
 
-		matrix<double> A( trans( subrange( wfv, 0, 10, lower.index(), higher.index()+1 ) ) );
+		matrix<double> A( trans( subrange( wfv, 0, atmos::indices_size, lower.index(), higher.index()+1 ) ) );
 		matrix<double> left( identity_matrix<double>( A.size1() ) );
 		matrix<double> right( identity_matrix<double>( A.size2() ) );
 
@@ -121,22 +128,22 @@
 
 	return res;
 }
-integrals_factory::integrals_factory() {
+functor_factory::functor_factory() {
 }
-integrals_factory::integrals_factory( const integrals_factory& factory ) {
+functor_factory::functor_factory( const functor_factory& factory ) {
 	m_storage = factory.m_storage;
 }
-void integrals_factory::swap( integrals_factory& factory ) throw() {
+void functor_factory::swap( functor_factory& factory ) throw() {
 	m_storage.swap( factory.m_storage );
 }
-integrals_factory& integrals_factory::operator=( const integrals_factory& factory ) {
-	integrals_factory tmp( factory );
+functor_factory& functor_factory::operator=( const functor_factory& factory ) {
+	functor_factory tmp( factory );
 	tmp.swap( *this );
 	return *this;
 }
-integrals_factory::~integrals_factory() {
+functor_factory::~functor_factory() {
 }
-integrals_factory::functor_type integrals_factory::get_functor( const std::string& spectral_type, const weight_type& awf ) {
+functor_factory::functor_type functor_factory::operator() ( const std::string& spectral_type, const weight_type& awf ) {
 	storage_type::const_iterator it = m_storage.find( spectral_type );
 	if( it != m_storage.end() )
 		return functor_type( it->second );
@@ -146,12 +153,15 @@
 	return functor_type( ret.first->second );
 }
 
-integrals_functor::result_type integrals_functor::operator() ( const indices_type& indices, const indesis_type& indesis, double mz ) const {
-	assert( indices.size1() == 4 );
-	assert( indices.size2() == 4 );
-	assert( indesis.size() == 4 );
+functor::result_type functor::operator() ( const indices_type& indices, const indesis_type& indesis, double mz ) const {
+	typedef vector< double >       vector_type;
+	typedef vector_type::size_type size_type;
+	
+	assert( indices.size1() == atmos::channels );
+	assert( indices.size2() == atmos::channels );
+	assert( indesis.size()  == atmos::channels );
 
-	vector< double > aints(10);
+	vector_type aints( atmos::indices_size );
 	aints(0) = indices(0,0); aints(1) = indices(1,1); aints(2) = indices(2,2); aints(3) = indices(3,3);
 	aints(4) = indices(0,1); aints(5) = indices(0,2); aints(6) = indices(0,3);
 	aints(7) = indices(1,2); aints(8) = indices(1,3); aints(9) = indices(2,3);
@@ -159,52 +169,31 @@
 	aints = prod( m_dcf, aints );
 
 	assert( npower == aints.size() );
-	for( vector< double >::size_type i = 0; i < npower; ++i ){
+	for( size_type i = 0; i < npower; ++i ){
 		aints(i) = aints(i) / std::pow( mz, 1.0 + powersi );
 	}
 
-	result_type ret;
-
-	ret.free_seeing       = KS * LEFF_SEE * std::pow( aints(0), JPOW_SEE );
-	ret.isoplanatic_angle = KP * LEFF_ISP * std::pow( aints(1), JPOW_ISP );
-	ret.altitude_eff      = 0.001 * std::pow( aints(1) / aints(2), 1.0/( POWISP-POWEFF ) );
-	ret.time_constant     = KT * std::pow( indesis(0), SPOW_TC );
-	ret.moment_0          = aints(0);
-	ret.moment_2          = aints(3);
-
-	return ret;
-}
-
-double integrals::free_seeing( const integrals_sample& x ) {
-	return x.free_seeing;
-}
-double integrals::isoplanatic_angle( const integrals_sample& x ) {
-	return x.isoplanatic_angle;
-}
-double integrals::altitude_eff( const integrals_sample& x ) {
-	return x.altitude_eff;
-}
-double integrals::time_constant( const integrals_sample& x ) {
-	return x.time_constant;
-}
-double integrals::moment_0( const integrals_sample& x ) {
-	return x.moment_0;
-}
-double integrals::moment_2( const integrals_sample& x ) {
-	return x.moment_2;
-}
-
-integrals_storage::integrals_storage( size_type reserve ):
-	table_storage< integrals::integrals_sample >( reserve ),
-	m_free_seeing_column( *this, boost::bind(&integrals::free_seeing,_1) ),
-	m_isoplanatic_angle_column( *this, boost::bind(&integrals::isoplanatic_angle,_1) ),
-	m_altitude_eff_column( *this, boost::bind(&integrals::altitude_eff,_1) ),
-	m_time_constant_column( *this, boost::bind(&integrals::time_constant,_1) ),
-	m_moment_0_column( *this, boost::bind(&integrals::moment_0,_1) ),
-	m_moment_2_column( *this, boost::bind(&integrals::moment_2,_1) ) {
-}
-
-bool integrals_regular_filter::operator() ( const integrals::integrals_sample& x ) const {
+	return result_type( KS * LEFF_SEE * std::pow( aints(0), JPOW_SEE ),
+			    KP * LEFF_ISP * std::pow( aints(1), JPOW_ISP ),
+			    0.001 * std::pow( aints(1) / aints(2), 1.0/( POWISP-POWEFF ) ),
+			    KT * std::pow( indesis(0), SPOW_TC ),
+			    aints(0),
+			    aints(3) );
+}
+
+storage::storage( size_type reserve ):
+	table_storage< detail::sample >( reserve ),
+	m_free_seeing_column( *this, boost::mem_fn(&detail::sample::free_seeing) ),
+	m_isoplanatic_angle_column( *this, boost::mem_fn(&detail::sample::isoplanatic_angle) ),
+	m_altitude_eff_column( *this, boost::mem_fn(&detail::sample::altitude_eff) ),
+	m_time_constant_column( *this, boost::mem_fn(&detail::sample::time_constant) ),
+	m_moment_0_column( *this, boost::mem_fn(&detail::sample::moment_0) ),
+	m_moment_2_column( *this, boost::mem_fn(&detail::sample::moment_2) ) {
+}
+
+namespace filter {
+	
+bool regular::operator() ( const sample_type& x ) const {
 	return ! ( is_regular( x.free_seeing ) &&
 		is_regular( x.isoplanatic_angle ) &&
 		is_regular( x.altitude_eff ) &&
@@ -212,6 +201,10 @@
 		is_regular( x.moment_0 ) &&
 		is_regular( x.moment_2 ) );
 }
-bool integrals_regular_filter::is_regular( double x ) const {
+bool regular::is_regular( double x ) const {
 	return ! ( std::isnan( x ) || std::isinf( x ) || ( x < 0 ) );
 }
+
+}
+
+}}

atmos-2.96.9.tar.gz/src/integrals.h -> atmos-2.97.1.tar.gz/src/integrals.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: integrals.h 158 2010-09-19 08:38:25Z matwey $
+    $Id: integrals.h 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -19,86 +19,51 @@
 #ifndef _INTEGRALS_H
 #define _INTEGRALS_H
 
-#include <string>
-#include <stdexcept>
-
-#include <boost/tuple/tuple.hpp>
+#include <boost/function.hpp>
 #include <boost/numeric/ublas/matrix.hpp>
 #include <boost/numeric/ublas/vector.hpp>
+#include <boost/tuple/tuple.hpp>
 
-#include "table_storage.h"
+#include "indices.h"
 #include "nocase_less.h"
+#include "table_storage.h"
 
-#include "utils.h"
-#include "vector_generator.h"
-#include "weif.h"
-
+namespace atmos {
 namespace integrals {
-
-	struct integrals_sample {
-		double free_seeing;
-		double isoplanatic_angle;
-		double altitude_eff;
-		double time_constant;
-		double moment_0;
-		double moment_2;
+	using namespace boost::numeric::ublas;
+	
+	typedef double value_type;
+	typedef matrix< value_type > matrix_type;
+
+namespace detail {
+	struct sample {
+		typedef atmos::integrals::value_type value_type;
+
+		value_type free_seeing;
+		value_type isoplanatic_angle;
+		value_type altitude_eff;
+		value_type time_constant;
+		value_type moment_0;
+		value_type moment_2;
+
+		sample( const value_type& free_seeing,
+			const value_type& isoplanatic_angle,
+			const value_type& altitude_eff,
+			const value_type& time_constant,
+			const value_type& moment_0,
+			const value_type& moment_2 ):
+			free_seeing(free_seeing), isoplanatic_angle(isoplanatic_angle), altitude_eff(altitude_eff), time_constant(time_constant), moment_0(moment_0), moment_2(moment_2) {
+		}
 	};
+}
 
-	double free_seeing( const integrals_sample& x );
-	double isoplanatic_angle( const integrals_sample& x );
-	double altitude_eff( const integrals_sample& x );
-	double time_constant( const integrals_sample& x );
-	double moment_0( const integrals_sample& x );
-	double moment_2( const integrals_sample& x );
-};
-using namespace boost::numeric::ublas;
-
-class integrals_factory;
-class integrals_functor;
-
-class integrals_factory {
-	public:
-		typedef matrix< double > dcf_matrix_type;
-		typedef boost::function1< vector<double>, double > weight_type;
-		typedef std::map< std::string, dcf_matrix_type, utils::algo::nocase_less< std::string > > storage_type;
-		typedef integrals_functor functor_type;
-		typedef storage_type::const_iterator const_iterator;
-	private:
-		storage_type m_storage;
-
-		dcf_matrix_type calculate( const std::string& spectral_type, const weight_type& awf );
-	public:
-		integrals_factory();
-		integrals_factory( const integrals_factory& factory );
-		void swap( integrals_factory& factory ) throw();
-		integrals_factory& operator=( const integrals_factory& factory );
-		~integrals_factory();
-
-		functor_type get_functor( const std::string& spectral_type, const weight_type& awf );
-};
-
-class integrals_functor {
-	public:
-		typedef integrals_factory             factory_type;
-		typedef factory_type::dcf_matrix_type dcf_matrix_type;
-
-		typedef integrals::integrals_sample result_type;
-
-		typedef symmetric_matrix<double> indices_type;
-		typedef vector<double> indesis_type;
-	private:
-		const dcf_matrix_type& m_dcf;
-	public:
-		integrals_functor( const dcf_matrix_type& dcf ):
-			m_dcf( dcf ) {}
-
-		result_type operator() ( const indices_type& indices, const indesis_type& indesis, double mz ) const;
-};
-
-class integrals_storage:
-	public table_storage< integrals::integrals_sample > {
+class storage:
+	public table_storage< detail::sample > {
+public:
+	typedef storage self_type;
+private:
+	typedef table_storage_column_getter< self_type, detail::sample::value_type > common_column_type;
 public:
-	typedef table_storage_column< integrals::integrals_sample, boost::function1< double, integrals::integrals_sample > > common_column_type;
 	typedef common_column_type free_seeing_column_type;
 	typedef common_column_type isoplanatic_angle_column_type;
 	typedef common_column_type altitude_eff_column_type;
@@ -113,7 +78,7 @@
 	moment_0_column_type          m_moment_0_column;
 	moment_2_column_type          m_moment_2_column;
 public:
-	integrals_storage( size_type reserve = 60 );
+	storage( size_type reserve = 60 );
 
 	const free_seeing_column_type&       free_seeing()       const { return m_free_seeing_column; }
 	const isoplanatic_angle_column_type& isoplanatic_angle() const { return m_isoplanatic_angle_column; }
@@ -122,16 +87,61 @@
 	const moment_0_column_type&          moment_0()          const { return m_moment_0_column; }
 	const moment_2_column_type&          moment_2()          const { return m_moment_2_column; }
 
-	void push_back( const integrals::integrals_sample& integrals ) {
-		table_storage< integrals::integrals_sample >::push_back( integrals );
+	void push_back( const detail::sample& integrals ) {
+		table_storage< detail::sample >::push_back( integrals );
 	}
 };
 
-class integrals_regular_filter {
-public:
-	bool operator() ( const integrals::integrals_sample& x ) const;
-private:
-	bool is_regular( double x ) const;
+namespace filter {
+	class regular {
+	public:
+		typedef atmos::integrals::detail::sample sample_type;
+	public:
+		bool operator() ( const sample_type& x ) const;
+	private:
+		bool is_regular( double x ) const;
+	};
+}
+
+class functor_factory;
+class functor;
+
+class functor_factory {
+	public:
+		typedef matrix_type dcf_matrix_type;
+		typedef boost::function1< vector<double>, double > weight_type;
+		typedef std::map< std::string, dcf_matrix_type, utils::algo::nocase_less< std::string > > storage_type;
+		typedef functor functor_type;
+	private:
+		storage_type m_storage;
+
+		dcf_matrix_type calculate( const std::string& spectral_type, const weight_type& awf );
+	public:
+		functor_factory();
+		functor_factory( const functor_factory& factory );
+		void swap( functor_factory& factory ) throw();
+		functor_factory& operator=( const functor_factory& factory );
+		~functor_factory();
+
+		functor_type operator() ( const std::string& spectral_type, const weight_type& awf );
 };
 
+class functor {
+	public:
+		typedef matrix_type dcf_matrix_type;
+
+		typedef atmos::indices::indices_type indices_type;
+		typedef atmos::indices::indesis_type indesis_type;
+		typedef atmos::integrals::detail::sample result_type;
+	private:
+		const dcf_matrix_type& m_dcf;
+	public:
+		functor( const dcf_matrix_type& dcf ):
+			m_dcf( dcf ) {}
+
+		result_type operator() ( const indices_type& indices, const indesis_type& indesis, double mz ) const;
+};
+
+}}
+
 #endif // _INTEGRALS_H

atmos-2.96.9.tar.gz/src/io.cpp -> atmos-2.97.1.tar.gz/src/io.cpp Changed

@@ -32,7 +32,7 @@
 		stm.iword(iword_index) = true;
 		return stm;
 	}
-	std::ostream& operator<<( std::ostream& stm, const datum& m ){
+	std::ostream& operator<<( std::ostream& stm, const atmos::datum& m ){
 		boost::io::ios_flags_saver fsaver( stm );
 		if( std::ostream::sentry( stm ) ){
 			int w = ( stm.flags() & std::ios::scientific )? 8 : 5;
@@ -43,13 +43,13 @@
 		}
 		return stm;
 	}
-	std::ostream& operator<<( std::ostream& stm, const object& obj ){
+	std::ostream& operator<<( std::ostream& stm, const atmos::object& obj ){
 		if( std::ostream::sentry( stm ) ){
 			stm << obj.hd << stm.fill() << obj.name << stm.fill() << std::hours << std::sexagesimal << yalpa::make_angle( obj.coordinates.ra() ) << std::degrees << stm.fill() << yalpa::make_angle( obj.coordinates.dec() ) << stm.fill() << obj.vmag << stm.fill() << obj.bvind << stm.fill() << obj.sclass;
 		}
 		return stm;
 	}
-	std::istream& operator>>( std::istream& stm, object& obj ){
+	std::istream& operator>>( std::istream& stm, atmos::object& obj ){
 		if( std::istream::sentry( stm ) ){
 			//O 2009-07-19 18:04:29 5191 Eta_UMa 13 47 32 +49 18 48  1.86 -0.19 B55
 			double ra = 0,dec = 0;

atmos-2.96.9.tar.gz/src/io.h -> atmos-2.97.1.tar.gz/src/io.h Changed

@@ -50,23 +50,29 @@
 	template<class T> std::ostream& operator<<( std::ostream& stm, const std::vector<T>& x ) {
 		typedef typename std::vector<T>::const_iterator const_iterator;
 
-		for( const_iterator it = x.begin(); it != x.end(); ++it ){
-			stm << *it << stm.fill();
+		const_iterator it = x.begin();
+		if( it != x.end() )
+			stm << *(it++);
+		for( ; it != x.end(); ++it ){
+			stm << stm.fill() << *it;
 		}
 		return stm;
 	}
 	template<class T, class VA> std::ostream& operator<<( std::ostream& stm, const boost::numeric::ublas::vector<T,VA>& x ) {
 		typedef typename boost::numeric::ublas::vector<T,VA>::const_iterator const_iterator;
 
-		for( const_iterator it = x.begin(); it != x.end(); ++it ){
-			stm << *it << stm.fill();
+		const_iterator it = x.begin();
+		if( it != x.end() )
+			stm << *(it++);
+		for( ; it != x.end(); ++it ){
+			stm << stm.fill() << *it;
 		}
 		return stm;
 	}
 	std::ostream& indic( std::ostream& stm );
-	std::ostream& operator<<( std::ostream& stm, const datum& m );
-	std::ostream& operator<<( std::ostream& stm, const object& obj );
-	std::istream& operator>>( std::istream& stm, object& obj );
+	std::ostream& operator<<( std::ostream& stm, const atmos::datum& m );
+	std::ostream& operator<<( std::ostream& stm, const atmos::object& obj );
+	std::istream& operator>>( std::istream& stm, atmos::object& obj );
 }
 
 #endif // _IO_H

atmos-2.97.1.tar.gz/src/local_context.h Added

@@ -0,0 +1,41 @@
+/*
+    $Id: local_context.h 180 2011-01-30 15:56:33Z matwey $
+    Copyright (C) 2011  Sternberg Astronomical Institute, MSU
+
+    This program is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef _LOCAL_CONTEXT_H
+#define _LOCAL_CONTEXT_H
+
+namespace utils {
+
+template<class Context> class local_context {
+	public:
+		typedef Context                context_type;
+		typedef context_type&          context_reference_type;
+		typedef const context_type&    context_const_reference_type;
+		typedef local_context<Context> local_context_type;
+	private:
+		context_reference_type m_ctx;
+	protected:
+		inline context_reference_type       context()       throw() { return m_ctx; }
+		inline context_const_reference_type context() const throw() { return m_ctx; }
+	public:
+		local_context( context_reference_type ctx ) throw() : m_ctx(ctx) {}
+};
+
+}
+
+#endif // _LOCAL_CONTEXT_H

atmos-2.96.9.tar.gz/src/main.cpp -> atmos-2.97.1.tar.gz/src/main.cpp Changed

atmos-2.97.1.tar.gz/src/params.cpp Added

@@ -0,0 +1,132 @@
+/*
+    $Id: params.cpp 185 2011-01-31 12:01:43Z matwey $
+    Copyright (C) 2009  Sternberg Astronomical Institute, MSU
+
+    This program is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <stdexcept>
+
+#include <sstream>
+
+#include <yalpa/spheric.h>
+#include <yalpa/io.h>
+
+#include "atmos.h"
+#include "params.h"
+
+namespace atmos {
+namespace parameters {
+
+container::container() {}
+container::~container() {}
+bool container::key_equal( const key_type& k1, const key_type& k2 ) {
+	return !m_storage.key_comp()( k1, k2 ) && !m_storage.key_comp()( k2, k1 );
+}
+bool container::value_equal( const value_type& v1, const value_type& v2 ){
+	return !m_storage.value_comp()( v1, v2 ) && !m_storage.value_comp()( v2, v1 );
+}
+bool container::has(const storage_type::key_type& key) const {
+	storage_type::const_iterator it = m_storage.find( key );
+	return ( it != m_storage.end() );
+}
+const container::storage_type::mapped_type& container::at(const storage_type::key_type& key) const {
+	storage_type::const_iterator it = m_storage.find( key );
+	if( it != m_storage.end() )
+		return it->second;
+	throw std::range_error( std::string("Not found value for key \"") + key + std::string("\""));
+}
+
+proxy::proxy( const container& container ): m_container( container ) {}
+proxy::~proxy() {}
+
+value_type proxy::longitude() const {
+	value_type longitude = 0.0;
+	std::istringstream iss( m_container.at("general/site/longitude") );
+	iss.exceptions(std::ios_base::badbit | std::ios_base::failbit);
+	iss >> std::sexagesimal >> std::hours >> yalpa::make_angle( longitude );
+	return longitude;
+}
+value_type proxy::latitude() const {
+	value_type latitude = 0.0;
+	std::istringstream iss( m_container.at("general/site/latitude") );
+	iss.exceptions(std::ios_base::badbit | std::ios_base::failbit);
+	iss >> std::sexagesimal >> std::degrees >> yalpa::make_angle( latitude );
+	return latitude;
+}
+yalpa::geographic< value_type > proxy::site() const {
+	return yalpa::geographic< value_type >( longitude(), latitude() );
+}
+
+vector_type proxy::inner_diam()    const {
+	typedef vector_type::value_type value_type;
+	vector_type ret(atmos::channels);
+	ret(0) = m_container.as<value_type>("optics/segment a/inner");
+	ret(1) = m_container.as<value_type>("optics/segment b/inner");
+	ret(2) = m_container.as<value_type>("optics/segment c/inner");
+	ret(3) = m_container.as<value_type>("optics/segment d/inner");
+	return ret * ( magnification()*0.1 );
+}
+vector_type proxy::outer_diam()    const {
+	typedef vector_type::value_type value_type;
+	vector_type ret(atmos::channels);
+	ret(0) = m_container.as<value_type>("optics/segment a/outer");
+	ret(1) = m_container.as<value_type>("optics/segment b/outer");
+	ret(2) = m_container.as<value_type>("optics/segment c/outer");
+	ret(3) = m_container.as<value_type>("optics/segment d/outer");
+	return ret * ( magnification()*0.1 );
+}
+vector_type proxy::nonlinearity() const {
+	typedef vector_type::value_type value_type;
+	std::vector< std::string > ctr( counters() );
+	vector_type ret( ctr.size() );
+	for( std::vector< std::string >::size_type i = 0; i < ctr.size(); ++i ) {
+		ret(i) = m_container.as< value_type >( (boost::format("Devices/Bicounter %1%/NonLinearity%2%") % ctri0 % ctri1).str() );
+	}
+	return ret * 1e-6 / exposition();
+}
+vector_type proxy::nonpoisson() const {
+	typedef vector_type::value_type value_type;
+	std::vector< std::string > ctr( counters() );
+	vector_type ret( ctr.size() );
+	for( std::vector< std::string >::size_type i = 0; i < ctr.size(); ++i ) {
+		ret(i) = m_container.as< value_type >( (boost::format("Devices/Bicounter %1%/NonPoisson%2%") % ctri0 % ctri1).str() );
+	}
+	return ret;
+}
+vector_type proxy::scatter() const {
+	typedef vector_type::value_type value_type;
+	vector_type ret(atmos::channels);
+	ret(0) = m_container.as< value_type >("channels/scatter/channela",0.0);
+	ret(1) = m_container.as< value_type >("channels/scatter/channelb",0.0);
+	ret(2) = m_container.as< value_type >("channels/scatter/channelc",0.0);
+	ret(3) = m_container.as< value_type >("channels/scatter/channeld",0.0);
+	return ret;
+}
+std::vector< std::string > proxy::counters() const {
+	std::vector< std::string > ret;
+	ret.resize(atmos::channels);
+	ret0 = m_container.at("channels/configuration/channela");
+	ret1 = m_container.at("channels/configuration/channelb");
+	ret2 = m_container.at("channels/configuration/channelc");
+	ret3 = m_container.at("channels/configuration/channeld");
+	for( std::vector< std::string >::iterator it = ret.begin(); it != ret.end(); ++it ){
+		//FIXME: use regex
+		if( it->size() < 10 || it->substr(0,7) != "Counter" || ( (*it)8 != '1' && (*it)8 != '2' ) || ( (*it)9 != 'A' && (*it)9 != 'B' ) ) throw std::logic_error( std::string("Bad value in channel configuration: ") + *it);
+		*it = it->substr(8,2);
+	}
+	return ret;
+}
+	
+}}

atmos-2.96.9.tar.gz/src/params.h -> atmos-2.97.1.tar.gz/src/params.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: params.h 98 2010-03-01 12:55:46Z matwey $
+    $Id: params.h 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -19,61 +19,44 @@
 #ifndef _PARAMS_H
 #define _PARAMS_H
 
-#include <algorithm>
-#include <cmath>
+#include <map>
 #include <string>
-#include <exception>
-#include <stdexcept>
 
-#include <boost/numeric/ublas/vector.hpp>
-#include <boost/lexical_cast.hpp>
 #include <boost/format.hpp>
+#include <boost/lexical_cast.hpp>
+#include <boost/numeric/ublas/vector.hpp>
 
-#include <yalpa/spheric.h>
-#include <yalpa/io.h>
+#include "nocase_less.h"
 
-using namespace boost::numeric;
-// use ublas::vector for ublas vector and std::vector for std vector
+namespace atmos {
+namespace parameters {
 
-/**
-  @brief Container for handling parameters extracted from MASS file.
-*/
-template<class ContainerType> class parameters_container {
+using namespace boost::numeric::ublas;
+
+class container {
 public:
-	typedef ContainerType storage_type;
-	typedef typename storage_type::key_type   key_type;
-	typedef typename storage_type::value_type value_type;
-	typedef typename storage_type::const_iterator const_iterator;
+	typedef std::map< std::string, std::string, utils::algo::nocase_less< std::string > > storage_type;
+	typedef storage_type::key_type   key_type;
+	typedef storage_type::value_type value_type;
+	typedef storage_type::const_iterator const_iterator;
 private:
 	storage_type m_storage;
 public:
-	parameters_container() {}
-	~parameters_container() {}
+	container();
+	~container();
 
-	bool key_equal( const key_type& k1, const key_type& k2 ) {
-		return !m_storage.key_comp()( k1, k2 ) && !m_storage.key_comp()( k2, k1 );
-	}
-	bool value_equal( const value_type& v1, const value_type& v2 ){
-		return !m_storage.value_comp()( v1, v2 ) && !m_storage.value_comp()( v2, v1 ) ;
-	}
-
-	bool has(const typename storage_type::key_type& key) const {
-		typename storage_type::const_iterator it = m_storage.find( key );
-		return ( it != m_storage.end() );
-	}
-	const typename storage_type::mapped_type& at(const typename storage_type::key_type& key) const {
-		typename storage_type::const_iterator it = m_storage.find( key );
-		if( it != m_storage.end() )
-			return it->second;
-		throw std::range_error( std::string("Not found value for key \"") + key + std::string("\""));
-	}
-	typename storage_type::mapped_type& operator (const typename storage_type::key_type& key) {
+	bool key_equal( const key_type& k1, const key_type& k2 );
+	bool value_equal( const value_type& v1, const value_type& v2 );
+	bool has(const storage_type::key_type& key) const;
+	const storage_type::mapped_type& at(const storage_type::key_type& key) const;
+	
+	storage_type::mapped_type& operator (const storage_type::key_type& key) {
 		return m_storagekey;
 	}
-	template<class Y> Y as(const typename storage_type::key_type& key) const {
+	template<class Y> Y as(const storage_type::key_type& key) const {
 		return boost::lexical_cast<Y>( at(key) );
 	}
-	template<class Y> Y as(const typename storage_type::key_type& key, const Y& def ) const {
+	template<class Y> Y as(const storage_type::key_type& key, const Y& def ) const {
 		if( has( key ) )
 			return boost::lexical_cast<Y>( at(key) );
 		return def;
@@ -81,110 +64,43 @@
 
 	const_iterator begin() const { return m_storage.begin(); }
 	const_iterator end() const { return m_storage.end(); }
-
 };
 
+typedef double               value_type;
+typedef vector< value_type > vector_type;
 
-/**
-  @brief Proxy for access to parameters extracted from MASS file and stored in container.
-*/
-template<class ParamType> class parameters_proxy {
-public:
-	typedef ParamType params_type;
+class proxy {
 private:
-	const params_type& m_container;
+	const container& m_container;
 public:
-	parameters_proxy( const params_type& container ): m_container( container ) {}
-	~parameters_proxy() {}
+	proxy( const container& container );
+	~proxy();
 
 	std::string sitename() const { return m_container.at("general/site/sitename"); }
-	double longitude() const {
-		double longitude = 0.0;
-		std::istringstream iss( m_container.at("general/site/longitude") );
-		iss.exceptions(std::ios_base::badbit | std::ios_base::failbit);
-		iss >> std::sexagesimal >> std::hours >> yalpa::make_angle( longitude );
-		return longitude;
-	}
-	double latitude() const {
-		double latitude = 0.0;
-		std::istringstream iss( m_container.at("general/site/latitude") );
-		iss.exceptions(std::ios_base::badbit | std::ios_base::failbit);
-		iss >> std::sexagesimal >> std::degrees >> yalpa::make_angle( latitude );
-		return latitude;
-	}
-	yalpa::geographic<double> site() const {
-		return yalpa::geographic<double>( longitude(), latitude() );
-	}
-
-	double accumtime()     const { return m_container.template as<double>("operations/normal mode/accumtime"); }
-	double basetime()      const { return m_container.template as<double>("operations/normal mode/basetime"); }
-	double exposition()    const { return m_container.template as<double>("operations/normal mode/exposition"); }
-
-	double magnification() const { return m_container.template as<double>("optics/common/magnification"); }
-
-	//FIXME: make function with wildcards
-	ublas::vector<double> inner_diam()    const {
-		ublas::vector<double> ret(4);
-		ret(0) = m_container.template as<double>("optics/segment a/inner");
-		ret(1) = m_container.template as<double>("optics/segment b/inner");
-		ret(2) = m_container.template as<double>("optics/segment c/inner");
-		ret(3) = m_container.template as<double>("optics/segment d/inner");
-		return ret * ( magnification()*0.1 );
-	}
-	ublas::vector<double> outer_diam()    const {
-		ublas::vector<double> ret(4);
-		ret(0) = m_container.template as<double>("optics/segment a/outer");
-		ret(1) = m_container.template as<double>("optics/segment b/outer");
-		ret(2) = m_container.template as<double>("optics/segment c/outer");
-		ret(3) = m_container.template as<double>("optics/segment d/outer");
-		return ret * ( magnification()*0.1 );
-	}
-	ublas::vector<double> nonlinearity() const {
-		std::vector< std::string > ctr( counters() );
-		ublas::vector<double> ret( ctr.size() );
-		for( typename std::vector< std::string >::const_iterator it = ctr.begin(); it != ctr.end(); ++it ) {
-			ret( it - ctr.begin() ) = m_container.template as<double>( (boost::format("Devices/Bicounter %1%/NonLinearity%2%") % (*it)0 % (*it)1).str() );
-		}
-		return ret * 1e-6 / exposition();
-	}
-	ublas::vector<double> nonpoisson() const {
-		std::vector< std::string > ctr( counters() );
-		ublas::vector<double> ret( ctr.size() );
-		for( typename std::vector< std::string >::const_iterator it = ctr.begin(); it != ctr.end(); ++it ) {
-			ret( it - ctr.begin() ) = m_container.template as<double>( (boost::format("Devices/Bicounter %1%/NonPoisson%2%") % (*it)0 % (*it)1).str() );
-		}
-		return ret;
-	}
-	ublas::vector<double> scatter() const {
-		ublas::vector<double> ret( 4 );
-		ret(0) = m_container.template as<double>("channels/scatter/channela",0.0);
-		ret(1) = m_container.template as<double>("channels/scatter/channelb",0.0);
-		ret(2) = m_container.template as<double>("channels/scatter/channelc",0.0);
-		ret(3) = m_container.template as<double>("channels/scatter/channeld",0.0);
-		return ret;
-	}
-	//FIXME: use tuple
-	std::vector< std::string > counters() const {
-		std::vector< std::string > ret;
-		ret.resize(4);
-		ret0 = m_container.at("channels/configuration/channela");
-		ret1 = m_container.at("channels/configuration/channelb");
-		ret2 = m_container.at("channels/configuration/channelc");
-		ret3 = m_container.at("channels/configuration/channeld");
-		for( typename std::vector< std::string >::iterator it = ret.begin(); it != ret.end(); ++it ){
-			//FIXME: use regex
-			if( it->size() < 10 || it->substr(0,7) != "Counter" || ( (*it)8 != '1' && (*it)8 != '2' ) || ( (*it)9 != 'A' && (*it)9 != 'B' ) ) throw std::logic_error( std::string("Bad value in channel configuration: ") + *it);
-			*it = it->substr(8,2);
-		}
-		return ret;
-	}
-
-	double dimm_aperture_base() const { return 0.01 * m_container.template as<double>("general/dimm/aperturebase",20.0); }
-	double dimm_aperture_size() const { return 0.01 * m_container.template as<double>("general/dimm/aperturesize",10.0); }
-	double dimm_scale() const { return 4.8481368e-06 * m_container.template as<double>("camera/geometry/scale"); }
-	double dimm_Ktransversal() const { return 0.358*(1.0*std::pow( dimm_aperture_size(), -1.0/3.0 ) - 0.811*std::pow( dimm_aperture_base(), -1.0/3.0 )); }
-	double dimm_Klongitudinal() const { return 0.358*(1.0*std::pow( dimm_aperture_size(), -1.0/3.0 ) - 0.541*std::pow( dimm_aperture_base(), -1.0/3.0 )); }
+	value_type longitude() const;
+	value_type latitude() const;
+	yalpa::geographic< value_type > site() const;
+
+	value_type accumtime()     const { return m_container.as< value_type >("operations/normal mode/accumtime"); }
+	value_type basetime()      const { return m_container.as< value_type >("operations/normal mode/basetime"); }
+	value_type exposition()    const { return m_container.as< value_type >("operations/normal mode/exposition"); }
+	value_type magnification() const { return m_container.as< value_type >("optics/common/magnification"); }
+
+	vector_type inner_diam()   const;
+	vector_type outer_diam()   const;
+	vector_type nonlinearity() const;
+	vector_type nonpoisson() const;
+	vector_type scatter() const;
+	std::vector< std::string > counters() const;
+
+	value_type dimm_aperture_base() const { return 0.01 * m_container.as< value_type >("general/dimm/aperturebase",20.0); }
+	value_type dimm_aperture_size() const { return 0.01 * m_container.as< value_type >("general/dimm/aperturesize",10.0); }
+	value_type dimm_scale() const { return 4.8481368e-06 * m_container.as< value_type >("camera/geometry/scale"); }
+	value_type dimm_Ktransversal() const { return 0.358*(1.0*std::pow( dimm_aperture_size(), -1.0/3.0 ) - 0.811*std::pow( dimm_aperture_base(), -1.0/3.0 )); }
+	value_type dimm_Klongitudinal() const { return 0.358*(1.0*std::pow( dimm_aperture_size(), -1.0/3.0 ) - 0.541*std::pow( dimm_aperture_base(), -1.0/3.0 )); }
 
 };
 
+}}
+
 #endif // _PARAMS_H

atmos-2.96.9.tar.gz/src/photometry.h -> atmos-2.97.1.tar.gz/src/photometry.h Changed

atmos-2.96.9.tar.gz/src/profile.cpp -> atmos-2.97.1.tar.gz/src/profile.cpp Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: profile.cpp 158 2010-09-19 08:38:25Z matwey $
+    $Id: profile.cpp 182 2011-01-31 11:26:03Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -16,101 +16,137 @@
     along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-#include <lsp/nnls.h>
-#include <boost/bind.hpp>
+#include <algorithm>
+
 #include <boost/format.hpp>
+#include <boost/mem_fn.hpp>
+#include <boost/numeric/ublas/banded.hpp>
 #include <boost/numeric/ublas/matrix_proxy.hpp>
 #include <boost/numeric/ublas/vector_proxy.hpp>
 #include <boost/numeric/ublas/triangular.hpp>
-#include <boost/numeric/ublas/io.hpp>
-#include <iostream>
-#include <iomanip>
-#include <algorithm>
-#include <numeric>
 
-#include "io.h"
+#include <lsp/cholesky_decomposition.h>
+#include <lsp/nnls.h>
+
+#include "atmos.h"
+#include "apply_to_all.h"
+#include "covariance_decomposition.h"
 #include "profile.h"
 #include "utils.h"
 
-const double MAXCHI2 = 999.0;
-const double JPOW_SEE = 0.6 ;
-const double LPOW_SEE = -0.2;
-const double KS = 1.73588e+07 ;
-
-using namespace boost::numeric::ublas;
+namespace atmos {
+namespace mixture {
 
-restoration_functor::restoration_functor( const vector<double>& grid, const symmetric_matrix< double >& corr_matrix, double secz, const weight_type& basewf ):
-	m_grid( grid ), m_cmatrix(corr_matrix), m_secz(secz), m_lambda_eff( 0.5 ) {
-	using namespace utils::ublas;
+namespace {
+	
+vector_type reorder_flat( const atmos::indices::indices_type& indices ) {
+	assert( indices.size1() == atmos::channels );
+	assert( indices.size2() == atmos::channels );
+	
+	vector_type inds( atmos::indices_size );
+	inds(0) = indices(0,0); inds(1) = indices(1,1); inds(2) = indices(2,2); inds(3) = indices(3,3);
+	inds(4) = indices(0,1); inds(5) = indices(0,2); inds(6) = indices(0,3);
+	inds(7) = indices(1,2); inds(8) = indices(1,3); inds(9) = indices(2,3);
+	return inds;
+}
 
-	if( zero_elements( diag( m_cmatrix ) ) > 0 ) {
-		throw std::logic_error( (boost::format("Zero variances: %1% x %2% %3%") % m_cmatrix.size1() % m_cmatrix.size2()  % vector<double>( diag( m_cmatrix )) ).str().c_str() );
+}
+	
+storage::storage( size_type reserve ):
+	table_storage< detail::sample >( reserve ),
+	m_mixture_column( *this, boost::mem_fn(&detail::sample::mixture) ) {
+}
+storage::storage( const atmos::indices::storage& indices, const atmos::dimm::storage& dimm, double dimmscale ):
+	table_storage< detail::sample >( indices.size() ),
+	m_mixture_column( *this, boost::mem_fn(&detail::sample::mixture) ) {
+	using namespace boost::numeric::ublas;
+	
+	if( dimm.size() > 1 ) {
+		double dimmscale2 = dimmscale*dimmscale;
+		dimm::storage::const_iterator dimm_current = dimm.begin();
+		dimm::storage::const_iterator dimm_next = dimm_current + 1;
+		indices::storage::const_iterator it = indices.begin();
+
+		while( dimm_current != dimm.end() ) {
+			while( it != indices.end() ) {
+				if( dimm_next != dimm.end() && dimm_next->timestamp <= it->timestamp )
+					break;
+
+				mixture_type mixture( atmos::indices_size + atmos::dimm_indices_size );
+				project( mixture, range(0,atmos::indices_size) ) = reorder_flat(it->indices);
+				mixture(10) = dimmscale2 * dimm_current->xvars(1) * 1e+11;
+				mixture(11) = dimmscale2 * dimm_current->yvars(1) * 1e+11;
+				++it;
+				push_back( mixture );
+			}
+			dimm_current = dimm_next;
+			++dimm_next;
+		}
+	} else {
+		for( indices::storage::const_iterator it = indices.begin(); it != indices.end(); ++it ) {
+			push_back( reorder_flat(it->indices) );
+		}
 	}
+}
 
-	lsp::cholesky_decomposition( m_cmatrix, upper_tag() );
-	triangular_adaptor< symmetric_matrix< double >, lower > R( m_cmatrix );
-	m_wmatrix = identity_matrix< double >( corr_matrix.size1() );
-	inplace_solve( R, m_wmatrix, lower_tag() );
+}
+namespace profile {
 
-	if( zero_elements( diag( m_wmatrix ) ) > 0 ) {
-		throw std::logic_error( (boost::format("Zero weights: %1%") % vector<double>( diag( m_wmatrix )) ).str().c_str() );
-	}
+const double MAXCHI2 = 999.0;
+const double JPOW_SEE = 0.6 ;
+const double LPOW_SEE = -0.2;
+const double KL = 1.1486983549970351; // 0.5 ** LPOW_SEE
+const double KS = 1.73588e+07;
 
-	m_pmatrix = prod( m_wmatrix, subrange( vector_of_vector_to_matrix( project(grid*secz, basewf) ) ,0,m_wmatrix.size2(),0,grid.size()) );
-}
-vector< double > restoration_functor::flat_reorder( const matrix<double>& indices ) {
-	vector< double > inds(10);
-	inds(0) = indices(0,0); inds(1) = indices(1,1); inds(2) = indices(2,2); inds(3) = indices(3,3);
-	inds(4) = indices(0,1); inds(5) = indices(0,2); inds(6) = indices(0,3);
-	inds(7) = indices(1,2); inds(8) = indices(1,3); inds(9) = indices(2,3);
-	return inds;
+functor::functor( const grid_type& grid, const mixture_covariance_type& cov_matrix, double secz, const weight_type& basewf ):
+	m_grid( grid ),
+	m_weight_matrix( utils::ublas::aetkin_weight_matrix(cov_matrix) ),
+	m_invweight_matrix( solve( m_weight_matrix, identity_matrix< value_type >( m_weight_matrix.size1() ), lower_tag() ) ),
+	m_secz( secz ) {
+	using namespace utils::ublas;
+	using namespace std;
+
+	m_problem_matrix = prod( m_weight_matrix, subrange( vector_of_vector_to_matrix( project(grid*secz, basewf) ), 0, m_weight_matrix.size2(),0,grid.size()) );
 }
-restoration_functor::result_type restoration_functor::operator() ( const vector<double>& indices ) const {
-	vector< double > inds = prod( m_wmatrix, indices );
+
+functor::result_type functor::operator() ( const mixture_type& indices ) const {
+	vector_type weighted_indices = prod( m_weight_matrix, indices );
 
 	result_type ret;
 
-	matrix< double > cov = zero_matrix<double>( m_pmatrix.size2(), m_pmatrix.size2() );
+	matrix_type cov = zero_matrix<double>( m_problem_matrix.size2(), m_problem_matrix.size2() );
 
-	lsp::nnls< matrix< double >, vector< double > > nnls( m_pmatrix, inds );
+	lsp::nnls< matrix_type, vector_type > nnls( m_problem_matrix, weighted_indices );
 	nnls.solve( ret.cn2, cov );
 
-	ret.residual = prod( m_pmatrix, ret.cn2 ) - inds;
+	ret.residual = prod( m_problem_matrix, ret.cn2 ) - weighted_indices;
 
 	ret.chi2 = std::min( inner_prod( ret.residual, ret.residual ), MAXCHI2 );
 
 	ret.cn2 /= m_secz;
 	return ret;
 }
-double restoration_functor::seeing( const vector<double>& cn2 ) const {
-	double cn2_sum = sum( cn2 );
-	double secz = 1.0; // already corrected cn2 at input
-	return KS*pow( cn2_sum/secz, JPOW_SEE )*pow( m_lambda_eff, LPOW_SEE );
+double functor::seeing( const cn2_type& cn2 ) const {
+	return KS*KL*std::pow( sum( cn2 ), JPOW_SEE );
 }
-vector<double> restoration_functor::indices_residual( const vector<double>& residuals ) const {
-	return prod( triangular_adaptor< const matrix< double >, lower >(m_cmatrix), residuals );
+vector_type functor::indices_residual( const residual_type& residual ) const {
+	return prod( m_invweight_matrix, residual );
 }
 
-profile::chi2_type profile::chi2( const profile_sample& x ) {
-	return x.chi2;
-}
-const profile::cn2_type& profile::cn2( const profile_sample& x ) {
-	return x.cn2;
-}
-const profile::residual_type& profile::residual( const profile_sample& x ) {
-	return x.residual;
+storage::storage( size_type reserve ):
+	table_storage< detail::sample >( reserve ),
+	m_cn2_column( *this, boost::mem_fn(&detail::sample::cn2) ),
+	m_residual_column( *this, boost::mem_fn(&detail::sample::residual)  ),
+	m_chi2_column( *this, boost::mem_fn(&detail::sample::chi2) ) {
 }
 
-profile_storage::profile_storage( size_type reserve ):
-	table_storage< profile::profile_sample >( reserve ),
-	m_cn2_column( *this, boost::bind(&profile::cn2, _1) ),
-	m_residual_column( *this, boost::bind(&profile::residual, _1)  ),
-	m_chi2_column( *this, boost::bind(&profile::chi2, _1) )
-{
-}
+namespace filter {
 
-chi2_gt_than_filter::chi2_gt_than_filter( double chi2 ): chi2_(chi2) {
+chi2_gt_than::chi2_gt_than( const value_type& chi2 ): chi2_(chi2) {
 }
-bool chi2_gt_than_filter::operator()( const profile::profile_sample& x ) const {
+bool chi2_gt_than::operator()( const sample_type& x ) const {
 	return x.chi2 > chi2_;
 }
+
+}
+}}
\ No newline at end of file

atmos-2.96.9.tar.gz/src/profile.h -> atmos-2.97.1.tar.gz/src/profile.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: profile.h 158 2010-09-19 08:38:25Z matwey $
+    $Id: profile.h 182 2011-01-31 11:26:03Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -19,101 +19,146 @@
 #ifndef _PROFILE_H
 #define _PROFILE_H
 
-#include <string>
-#include <stdexcept>
-
-#include <lsp/cholesky_decomposition.h>
-
+#include <boost/function.hpp>
 #include <boost/numeric/ublas/matrix.hpp>
 #include <boost/numeric/ublas/vector.hpp>
 #include <boost/numeric/ublas/symmetric.hpp>
 
+#include "indices.h"
+#include "dimm.h"
 #include "table_storage.h"
-#include "weif.h"
 
-namespace profile {
+namespace atmos {
+namespace mixture {
 	using namespace boost::numeric::ublas;
 
 	typedef double                 value_type;
 	typedef vector<value_type>     vector_type;
 	typedef vector_type::size_type size_type;
+	typedef symmetric_matrix<value_type> symmetric_matrix_type;
 
-	typedef vector_type cn2_type;
-	typedef vector_type residual_type;
-	typedef value_type  chi2_type;
+	typedef vector_type           mixture_type;
+	typedef symmetric_matrix_type mixture_covariance_type;
 
-	struct profile_sample {
-		cn2_type       cn2;
-		residual_type  residual;
-		chi2_type      chi2;
-	};
+namespace detail {
+	struct sample {
+		typedef atmos::mixture::mixture_type mixture_type;
 
-	chi2_type chi2( const profile_sample& x );
-	const cn2_type& cn2( const profile_sample& x );
-	const residual_type& residual( const profile_sample& x );
-};
+		mixture_type mixture;
 
-template<class VE> typename VE::size_type zero_elements( const vector_expression<VE>&x ) {
-    typename VE::size_type n = 0;
-    for( typename VE::size_type j = 0; j < x().size(); j++ )
-	    n += ( x()(j) == 0 || std::isinf(x()(j)) || std::isnan(x()(j)) ) ? 1 : 0;
-    return n;
+		sample( const mixture_type& mixture ): mixture(mixture) {
+		}
+	};
 }
 
-class restoration_functor {
-	public:
-		typedef profile::profile_sample result_type;
+class storage:
+	public table_storage< detail::sample > {
+public:
+	typedef storage self_type;
+	typedef table_storage_column_getter< self_type, detail::sample::mixture_type > mixture_column_type;
+private:
+	mixture_column_type m_mixture_column;
+public:
+	storage( size_type reserve = 60 );
+	storage( const atmos::indices::storage& indices, const atmos::dimm::storage& dimm, double dimmscale );
 
-		typedef vector<double> grid_type;
-		typedef matrix<double> weight_matrix;
-		typedef boost::function1< vector<double>, double > weight_type;
-	private:
-		grid_type     m_grid;
-		weight_matrix m_wmatrix;
-		symmetric_matrix< double > m_cmatrix;
-		matrix< double > m_pmatrix;
-		double        m_secz;
-		double        m_lambda_eff;
+	const mixture_column_type& mixture() const { return m_mixture_column; }
 
-	public:
-		static vector< double > flat_reorder( const matrix<double>& indices );
+	void push_back( const mixture_type& mixture ) {
+		table_storage< detail::sample >::push_back( detail::sample( mixture ) );
+	}
+};
+}
+namespace profile {
+	using namespace boost::numeric::ublas;
 
-		restoration_functor( const vector<double>& grid, const symmetric_matrix< double >& corr_matrix, double secz, const weight_type& basewf );
+	typedef double                 value_type;
+	typedef vector<value_type>     vector_type;
+	typedef matrix<value_type>     matrix_type;
+	typedef triangular_matrix<value_type> triangular_matrix_type;
+	typedef vector_type::size_type size_type;
 
-		result_type operator() ( const vector<double>& indices ) const;
+	typedef vector_type cn2_type;
+	typedef vector_type residual_type;
+	typedef value_type  chi2_type;
 
-		double seeing( const vector<double>& cn2 ) const;
-		vector<double> indices_residual( const vector<double>& residuals ) const;
-};
+namespace detail {
+	struct sample {
+		typedef atmos::profile::cn2_type      cn2_type;
+		typedef atmos::profile::residual_type residual_type;
+		typedef atmos::profile::chi2_type     chi2_type;
+	
+		cn2_type       cn2;
+		residual_type  residual;
+		chi2_type      chi2;
+
+		sample( const cn2_type& cn2, const residual_type& residual, const chi2_type& chi2 ):
+			cn2(cn2), residual(residual), chi2(chi2) {
+		}
+		sample() {
+		}
+	};
+}
 
-class profile_storage:
-	public table_storage< profile::profile_sample > {
+class storage:
+	public table_storage< detail::sample > {
 public:
-	typedef table_storage_column< profile::profile_sample, boost::function1< profile::cn2_type,      profile::profile_sample > > cn2_column_type;
-	typedef table_storage_column< profile::profile_sample, boost::function1< profile::residual_type, profile::profile_sample > > residual_column_type;
-	typedef table_storage_column< profile::profile_sample, boost::function1< profile::chi2_type,     profile::profile_sample > > chi2_column_type;
+	typedef storage self_type;
+	typedef table_storage_column_getter< self_type, detail::sample::cn2_type >      cn2_column_type;
+	typedef table_storage_column_getter< self_type, detail::sample::residual_type > residual_column_type;
+	typedef table_storage_column_getter< self_type, detail::sample::chi2_type >     chi2_column_type;
 private:
 	cn2_column_type      m_cn2_column;
 	residual_column_type m_residual_column;
 	chi2_column_type     m_chi2_column;
 public:
-	profile_storage( size_type reserve = 60 );
+	storage( size_type reserve = 60 );
 
 	const cn2_column_type&      cn2()      const { return m_cn2_column; }
 	const residual_column_type& residual() const { return m_residual_column; }
 	const chi2_column_type&     chi2()     const { return m_chi2_column; }
 
-	void push_back( const profile::profile_sample& profile ) {
-		table_storage< profile::profile_sample >::push_back( profile );
+	void push_back( const detail::sample& profile ) {
+		table_storage< detail::sample >::push_back( profile );
 	}
 };
 
-class chi2_gt_than_filter {
+class functor {
+	public:
+		typedef detail::sample result_type;
+		typedef atmos::mixture::mixture_type mixture_type;
+		typedef atmos::mixture::mixture_covariance_type mixture_covariance_type;
+
+		typedef vector_type grid_type;
+		typedef boost::function1< vector<double>, double > weight_type;
 	private:
-		double chi2_;
+		grid_type     m_grid;
+		triangular_matrix_type m_weight_matrix;
+		triangular_matrix_type m_invweight_matrix;
+		matrix_type   m_problem_matrix;
+		double        m_secz;
+
 	public:
-		chi2_gt_than_filter( double chi2 );
-		bool operator()( const profile::profile_sample& x ) const;
+		functor( const grid_type& grid, const mixture_covariance_type& cov_matrix, double secz, const weight_type& basewf );
+
+		result_type operator() ( const mixture_type& indices ) const;
+
+		double seeing( const cn2_type& cn2 ) const;
+		vector_type indices_residual( const residual_type& residual ) const;
 };
 
+namespace filter {
+	class chi2_gt_than {
+	public:
+		typedef atmos::profile::detail::sample sample_type;
+		typedef atmos::profile::value_type value_type;
+	private:
+		value_type chi2_;
+	public:
+		chi2_gt_than( const value_type& chi2 );
+		bool operator()( const sample_type& x ) const;
+	};
+}
+}}
+
 #endif // _PROFILE_H

atmos-2.96.9.tar.gz/src/spectral_factory.cpp -> atmos-2.97.1.tar.gz/src/spectral_factory.cpp Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: spectral_factory.cpp 158 2010-09-19 08:38:25Z matwey $
+    $Id: spectral_factory.cpp 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -22,16 +22,15 @@
 #include <boost/algorithm/string.hpp>
 #include <boost/filesystem/operations.hpp>
 
-#include <boost/numeric/ublas/io.hpp>
-
+#include "io.h"
+#include "spectral_factory.h"
 #include "weif.h"
 #include "utils.h"
 
-#include "spectral_factory.h"
-
-using namespace boost::numeric::ublas;
+namespace atmos {
+namespace weight {
 
-std::pair< spectral::vector_type, spectral::vector_type > spectral::load_reaction( const boost::filesystem::path& file ) {
+std::pair< vector_type, vector_type > load_reaction( const boost::filesystem::path& file ) {
 	std::ifstream fstm( file.string().c_str() );
 	if( fstm.fail() )
 		throw std::runtime_error( "Can't open file " + file.string() );
@@ -79,7 +78,7 @@
 	return ret;
 }
 
-std::pair< spectral::vector_type, spectral::vector_type > spectral::photon_distribution( const vector_type& spectral_lambda, const vector_type& spectral_response, const vector_type& device_lambda, const vector_type& device_response ) {
+std::pair< vector_type, vector_type > photon_distribution( const vector_type& spectral_lambda, const vector_type& spectral_response, const vector_type& device_lambda, const vector_type& device_response ) {
 	std::pair< vector_type, vector_type > ret;
 
 	assert( spectral_lambda.size() == spectral_response.size() );
@@ -112,7 +111,7 @@
 	return ret;
 }
 
-void spectral::dump_weight( const boost::filesystem::path& file, const boost::function1< vector_type, double >& function, const vector_type& grid ) {
+void dump_weight( const boost::filesystem::path& file, const boost::function1< vector_type, double >& function, const vector_type& grid ) {
 	static const char* name = { "A", "B", "C", "D", "AB", "AC", "AD", "BC", "BD", "CD", "L", "T" };
 	static const unsigned int name_size = sizeof(name) / sizeof(name0);
 	const double wscale = 1e+11;
@@ -146,37 +145,38 @@
 	}
 }
 
-spectral_factory::spectral_factory(const boost::filesystem::path& wd, const boost::filesystem::path& mass_response, const boost::filesystem::path& ccd_response ):
+factory::factory(const boost::filesystem::path& wd, const boost::filesystem::path& mass_response, const boost::filesystem::path& ccd_response ):
 	m_mass_response(mass_response),
 	m_ccd_response(ccd_response),
 	m_wdir(wd),
 	m_basis_grid( lg_generator(0.0, 0.04, 0.25), 50 ),
-	m_mass_reaction( spectral::load_reaction(mass_response) ),
-	m_dimm_reaction( spectral::load_reaction(ccd_response) ) {
+	m_mass_reaction( load_reaction(mass_response) ),
+	m_dimm_reaction( load_reaction(ccd_response) ) {
 
 	if( ! (boost::filesystem::exists( wd ) && boost::filesystem::is_directory( wd )) )
 		throw std::runtime_error( wd.string() + " is not a directory or doesn't exists");
 }
-spectral_factory::~spectral_factory() {
-
+factory::~factory() {
 }
-const boost::filesystem::path spectral_factory::get_spectrum_filename( const std::string& spectype ) const {
+const boost::filesystem::path factory::get_spectrum_filename( const std::string& spectype ) const {
 	std::string spectral_type( spectype );
 	boost::to_lower( spectral_type );
 	boost::filesystem::path spect( m_wdir / (spectral_type + ".sp") );
 	return spect;
 }
-const spectral_factory::weight_type& spectral_factory::spectral( const std::string& spectype, const vector<double>& inner_diam, const vector<double>& outer_diam, double dimmbase, double dimmsize ) {
+const factory::weight_type& factory::weight( const std::string& spectype, const vector<double>& inner_diam, const vector<double>& outer_diam, double dimmbase, double dimmsize ) {
 	storage_type::const_iterator it = m_storage.find( spectype );
 	if( it != m_storage.end() )
 		return it->second;
 
-	std::pair< vector_type, vector_type > spectral_response = spectral::load_reaction( get_spectrum_filename( spectype ) );
+	std::pair< vector_type, vector_type > spectral_response = load_reaction( get_spectrum_filename( spectype ) );
 
-	std::pair< vector_type, vector_type > mass_response = spectral::photon_distribution( spectral_response.first, spectral_response.second, m_mass_reaction.first, m_mass_reaction.second );
-	std::pair< vector_type, vector_type > dimm_response = spectral::photon_distribution( spectral_response.first, spectral_response.second, m_dimm_reaction.first, m_dimm_reaction.second );
+	std::pair< vector_type, vector_type > mass_response = photon_distribution( spectral_response.first, spectral_response.second, m_mass_reaction.first, m_mass_reaction.second );
+	std::pair< vector_type, vector_type > dimm_response = photon_distribution( spectral_response.first, spectral_response.second, m_dimm_reaction.first, m_dimm_reaction.second );
 
-	std::pair< storage_type::iterator, bool > ret = m_storage.insert( std::make_pair( spectype, weight_type( m_basis_grid, utils::ublas::project( m_basis_grid , weif::massdimm_weight( mass_response.first, mass_response.second, dimm_response.first, dimm_response.second, inner_diam, outer_diam, dimmsize, dimmbase ) ) ) ) );
+	std::pair< storage_type::iterator, bool > ret = m_storage.insert( std::make_pair( spectype, weight_type( m_basis_grid, utils::ublas::project( m_basis_grid , atmos::weight::massdimm_weight( mass_response.first, mass_response.second, dimm_response.first, dimm_response.second, inner_diam, outer_diam, dimmsize, dimmbase ) ) ) ) );
 
 	return ret.first->second;
 }
+
+}}

atmos-2.96.9.tar.gz/src/spectral_factory.h -> atmos-2.97.1.tar.gz/src/spectral_factory.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: spectral_factory.h 158 2010-09-19 08:38:25Z matwey $
+    $Id: spectral_factory.h 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -23,29 +23,31 @@
 #include <string>
 
 #include <boost/filesystem/path.hpp>
-#include <boost/numeric/ublas/vector.hpp>
 #include <boost/function.hpp>
+#include <boost/numeric/ublas/vector.hpp>
 
+#include "spline.h"
 #include "nocase_less.h"
 #include "vector_generator.h"
+#include "utils.h"
 #include "weif.h"
 
-namespace spectral {
+namespace atmos {
+namespace weight {
 
 using namespace boost::numeric::ublas;
-typedef vector< double > vector_type;
+using namespace utils::algo;
+using namespace utils::ublas;
 
 std::pair< vector_type, vector_type > load_reaction( const boost::filesystem::path& file );
 std::pair< vector_type, vector_type > photon_distribution( const vector_type& spectral_lambda, const vector_type& spectral_response, const vector_type& device_lambda, const vector_type& device_response );
 
-void dump_weight( const boost::filesystem::path& file, const boost::function1< vector_type, double >& function, const vector_type& grid = utils::ublas::vector_generator< lg_generator >( lg_generator(0.0, 0.04, 0.25), 50 ) );
-
-}
+void dump_weight( const boost::filesystem::path& file, const boost::function1< vector_type, double >& function, const vector_type& grid = vector_generator< lg_generator >( lg_generator(0.0, 0.04, 0.25), 50 ) );
 
-class spectral_factory{
+class factory{
 public:
-	typedef boost::numeric::ublas::vector< double > vector_type;
-	typedef boost::numeric::ublas::vector< boost::numeric::ublas::vector<double> > vector_vector_type;
+	typedef atmos::weight::vector_type vector_type;
+	typedef vector< vector_type > vector_vector_type;
 
 	typedef utils::ublas::cubic_spline< vector_type, vector_vector_type, vector_type, vector_vector_type > weight_type;
 	typedef std::map< std::string, weight_type, utils::algo::nocase_less< std::string > > storage_type;
@@ -56,7 +58,7 @@
 	boost::filesystem::path m_ccd_response;
 	boost::filesystem::path m_wdir;
 
-	utils::ublas::vector_generator< lg_generator > m_basis_grid;
+	vector_generator< lg_generator > m_basis_grid;
 
 	storage_type m_storage;
 
@@ -65,14 +67,15 @@
 private:
 	const boost::filesystem::path get_spectrum_filename( const std::string& spectype ) const;
 public:
-	spectral_factory(const boost::filesystem::path& wd, const boost::filesystem::path& mass_response, const boost::filesystem::path& ccd_response );
-	~spectral_factory();
+	factory(const boost::filesystem::path& wd, const boost::filesystem::path& mass_response, const boost::filesystem::path& ccd_response );
+	~factory();
 
-	const weight_type& spectral( const std::string& spectype, const boost::numeric::ublas::vector<double>& inner_diam, const boost::numeric::ublas::vector<double>& outer_diam, double dimmbase, double dimmsize );
+	const weight_type& weight( const std::string& spectype, const vector_type& inner_diam, const vector_type& outer_diam, double dimmbase, double dimmsize );
 
 	const const_iterator begin() const { return m_storage.begin(); }
 	const const_iterator end()   const { return m_storage.end(); }
 };
 
-#endif //_SPECTRAL_FACTORY_H
+}}
 
+#endif //_SPECTRAL_FACTORY_H

atmos-2.96.9.tar.gz/src/table_storage.h -> atmos-2.97.1.tar.gz/src/table_storage.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: table_storage.h 143 2010-04-04 08:42:34Z matwey $
+    $Id: table_storage.h 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2010  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -75,29 +75,97 @@
 		m_storage.erase( unstable_remove_if( m_storage.begin(), m_storage.end(), pred ), m_storage.end() );
 	}
 	template<class Predicate> void remove_if( Predicate pred ) {
-		m_storage.erase( remove( m_storage.begin(), m_storage.end(), pred ), m_storage.end() );
+		using std::remove_if;
+		m_storage.erase( remove_if( m_storage.begin(), m_storage.end(), pred ), m_storage.end() );
 	}
 };
 
-template<class T, class F> class table_storage_column {
+namespace detail {
+template<class Result, class Iterator> struct table_storage_column_traits {
+	typedef Iterator subiterator;
+	typedef typename std::iterator_traits< subiterator >::reference arg_type;
+	typedef Result result_type;
+	typedef boost::function1< result_type, arg_type > function_type;
+	typedef boost::transform_iterator< function_type, subiterator, result_type > iterator;
+};
+}
+
+template<class Storage, class T> class table_storage_column_getter{
 public:
-	typedef T                subvalue_type;
-	typedef table_storage<T> table_storage_type;
-	typedef F                proxy_function_type;
+	typedef Storage storage_type;
+	typedef typename storage_type::const_iterator subiterator;
+	typedef const T& result_type;
+	typedef typename detail::table_storage_column_traits< result_type, subiterator >::function_type function_type;
+	typedef typename detail::table_storage_column_traits< result_type, subiterator >::iterator      iterator;
+private:
+	const storage_type& m_storage;
+	function_type m_function;
+public:
+	table_storage_column_getter( const storage_type& storage, const function_type& function ): m_storage(storage), m_function(function) {}
 
-	typedef typename table_storage_type::const_iterator const_subiterator;
+	iterator begin() const { return iterator( m_storage.begin(), m_function ); }
+	iterator end() const   { return iterator( m_storage.end(), m_function ); }
+};
 
-	typedef boost::transform_iterator< proxy_function_type, const_subiterator > const_iterator;
+template<class Storage, class T> class table_storage_column_setter{
+public:
+	typedef Storage storage_type;
+	typedef typename storage_type::iterator subiterator;
+	typedef T& result_type;
+	typedef typename detail::table_storage_column_traits< result_type, subiterator >::function_type function_type;
+	typedef typename detail::table_storage_column_traits< result_type, subiterator >::iterator      iterator;
 private:
-	const table_storage<T>& m_storage;
-	proxy_function_type m_proxy_function;
+	storage_type& m_storage;
+	function_type m_function;
 public:
-	table_storage_column( const table_storage<T>& storage, const F& proxy_function ):
-		m_storage( storage ), m_proxy_function( proxy_function ) {
-	}
+	table_storage_column_setter( storage_type& storage, const function_type& function ): m_storage(storage), m_function(function) {}
+
+	iterator begin() { return iterator( m_storage.begin(), m_function ); }
+	iterator end()   { return iterator( m_storage.end(), m_function ); }
+};
 
-	const_iterator begin() const { return const_iterator( m_storage.begin(), m_proxy_function ); }
-	const_iterator end()   const { return const_iterator( m_storage.end(), m_proxy_function ); }
+template<class Storage, class T> class table_storage_column_generator{
+public:
+	typedef Storage storage_type;
+	typedef typename storage_type::const_iterator subiterator;
+	typedef T result_type;
+	typedef typename detail::table_storage_column_traits< result_type, subiterator >::function_type function_type;
+	typedef typename detail::table_storage_column_traits< result_type, subiterator >::iterator      iterator;
+private:
+	const storage_type& m_storage;
+	function_type m_function;
+public:
+	table_storage_column_generator( const storage_type& storage, const function_type& function ): m_storage(storage), m_function(function) {}
+
+	iterator begin() const { return iterator( m_storage.begin(), m_function ); }
+	iterator end() const   { return iterator( m_storage.end(), m_function ); }
+};
+
+template<class Storage, class T> class table_storage_column_proxy{
+public:
+	typedef Storage storage_type;
+	typedef typename storage_type::iterator       subiterator;
+	typedef typename storage_type::const_iterator const_subiterator;
+	typedef T& result_type;
+	typedef const T& const_result_type;
+	typedef typename detail::table_storage_column_traits< result_type, subiterator >::function_type function_type;
+	typedef typename detail::table_storage_column_traits< result_type, subiterator >::iterator      iterator;
+	typedef typename detail::table_storage_column_traits< const_result_type, const_subiterator >::function_type const_function_type;
+	typedef typename detail::table_storage_column_traits< const_result_type, const_subiterator >::iterator      const_iterator;
+private:
+	storage_type& m_storage;
+	function_type m_function;
+	const_function_type m_const_function;
+public:
+	table_storage_column_proxy( storage_type& storage, const function_type& function, const const_function_type& const_function ):
+		m_storage(storage), m_function(function), m_const_function(const_function) {}
+	template<class F> table_storage_column_proxy( storage_type& storage, const F& function ):
+		m_storage(storage), m_function(function), m_const_function(function) {}
+
+	iterator begin() { return iterator( m_storage.begin(), m_function ); }
+	iterator end()   { return iterator( m_storage.end(), m_function ); }
+	const_iterator begin() const { return const_iterator( m_storage.begin(), m_const_function ); }
+	const_iterator end()   const { return const_iterator( m_storage.end(), m_const_function ); }
 };
 
 template<class T> void swap( table_storage<T>& t, table_storage<T>& u ) {

atmos-2.96.9.tar.gz/src/typesdef.h -> atmos-2.97.1.tar.gz/src/typesdef.h Changed

atmos-2.97.1.tar.gz/src/utils.cpp Added

@@ -0,0 +1,48 @@
+/*
+    $Id: utils.cpp 180 2011-01-30 15:56:33Z matwey $
+    Copyright (C) 2011  Sternberg Astronomical Institute, MSU
+
+    This program is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <cmath>
+
+#include "utils.h"
+
+namespace utils {
+namespace algo {
+
+linear_generator::linear_generator( result_type start, result_type step ): start_(start), step_(step) {
+}
+linear_generator::result_type linear_generator::operator()( size_type i ) const {
+	return ( start_ + step_ * i );
+}
+linear_generator::size_type linear_generator::size( result_type max ) const {
+	return static_cast<size_type>( std::floor( (max - start_) / step_ ) );
+}
+
+lg_generator::lg_generator( result_type start, result_type lstep, result_type gstep ): start_(start), t_(lstep), q_(result_type(1) + gstep / 2) {
+	max_ = static_cast< size_type >( std::floor( result_type(2) / gstep - start / lstep ) );
+}
+lg_generator::result_type lg_generator::operator()( size_type i ) const {
+	if( i <= max_ ) return ( start_ + t_ * i );
+	return ( start_ + t_ * max_ ) * std::pow( q_, static_cast<int>(i - max_) ); // FIXME: pow(double,int) for <4.2.1
+}
+lg_generator::size_type lg_generator::size( result_type max ) const {
+	result_type gmax = std::floor( std::log( max / ( start_ + t_ * max_ ) ) / std::log( q_ ) );
+	if( gmax > result_type(0) ) return static_cast< size_type >( gmax ) + max_ + 1;
+	return static_cast< size_type >(std::floor( (max - start_) / t_ ));
+}
+
+}}

atmos-2.96.9.tar.gz/src/utils.h -> atmos-2.97.1.tar.gz/src/utils.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: utils.h 158 2010-09-19 08:38:25Z matwey $
+    $Id: utils.h 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2010  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -16,6 +16,13 @@
     along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
 
+#ifndef _UTILS_H_
+#define _UTILS_H_
+
+#include <boost/numeric/ublas/matrix_proxy.hpp>
+#include <boost/numeric/ublas/triangular.hpp>
+#include <boost/numeric/ublas/vector.hpp>
+
 namespace utils {
 namespace ublas {
 
@@ -60,5 +67,31 @@
 	return ret;
 }
 
-} // utils
 } // ublas
+
+namespace algo {
+struct linear_generator {
+	typedef double result_type;
+	typedef size_t size_type;
+
+	linear_generator( result_type start, result_type step );
+	result_type operator()( size_type i ) const;
+	size_type size( result_type max ) const;
+private:
+	result_type start_,step_;
+};
+struct lg_generator {
+	typedef double result_type;
+	typedef size_t size_type;
+
+	lg_generator( result_type start, result_type lstep, result_type gstep );
+	result_type operator()( size_type i ) const;
+	size_type size( result_type max ) const;
+private:
+	result_type start_,t_,q_;
+	size_type max_;
+};
+} // algo
+} // utils
+
+#endif // _UTILS_H_

atmos-2.96.9.tar.gz/src/weif.cpp -> atmos-2.97.1.tar.gz/src/weif.cpp Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: weif.cpp 158 2010-09-19 08:38:25Z matwey $
+    $Id: weif.cpp 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2010  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -17,16 +17,18 @@
 */
 
 #include <cmath>
-#include <boost/function.hpp>
-
-#include <boost/numeric/ublas/io.hpp>
 
-#include "weif.h"
+#include <boost/function.hpp>
+#include <boost/numeric/ublas/vector_proxy.hpp>
 
-#include "gsl_wrap.h"
+#include "atmos.h"
 #include "fft.h"
+#include "gsl_wrap.h"
+#include "utils.h"
+#include "weif.h"
 
-namespace weif {
+namespace atmos {
+namespace weight {
 
 const double POWIND = -2.66666666667;
 const double POWDIFF = -0.666666666667;
@@ -45,8 +47,9 @@
 	double diameter3;
 	double diameter4;
 	double PIx2lambda0;
-	const sfun_type& sfun;
+	sfun_type sfun;
 	double altitude;
+	double powind;
 
 	mass_params( const sfun_type& sfun ): sfun(sfun) {}
 };
@@ -58,31 +61,29 @@
 	double dimmbase;
 	double PIx2lambda0;
 	double phi;
-	const sfun_type& sfun;
+	sfun_type sfun;
 	double altitude;
 
 	dimm_params( const sfun_type& sfun ): sfun(sfun) {}
 };
 
 double mass_integrand(double x, void * params) {
+	using namespace std;
+	
 	mass_params p = *( mass_params *) params;
-	double apert1 = 0.0, apert2 = 0.0, apert3 = 0.0, apert4 = 0;
+	double p1 = 0;
+	double p3 = 0;
 	double f = p.altitude*x*x/2;
 	const std::complex<double>& sfun_f = p.sfun(f);
 	double sf = imag( sfun_f ) * cos( p.PIx2lambda0*f ) + real( sfun_f ) * sin( p.PIx2lambda0*f );
-	sf = sf*sf*pow( x, POWIND );
-	if( p.diameter1 != p.diameter2 ) {                // First aperture:
-		apert1 = p.diameter1*j1( p.diameter1*x );
-		apert2 = p.diameter2*j1( p.diameter2*x );
-		apert1 = ( apert1 - apert2 )/( p.diameter1*p.diameter1 - p.diameter2*p.diameter2 );
+	sf = sf*sf*pow( x, p.powind );
+	if( p.diameter1 != p.diameter2 ) {
+		p1 = ( p.diameter1*j1( p.diameter1*x ) - p.diameter2*j1( p.diameter2*x ) ) / ( p.diameter1*p.diameter1 - p.diameter2*p.diameter2 );
 	}
-	if( p.diameter3 != p.diameter4) {                 // Second aperture: the same with indices 1->3, 2->4:
-		apert3 = p.diameter3*j1( p.diameter3*x );
-		apert4 = ( p.diameter4 == 0 )? 0 : p.diameter4*j1( p.diameter4*x );
-		apert3 = ( apert3 - apert4 )/( p.diameter3*p.diameter3 - p.diameter4*p.diameter4 );
+	if( p.diameter3 != p.diameter4 ) {
+		p3 = ( p.diameter3*j1( p.diameter3*x ) - p.diameter4*j1( p.diameter4*x ) ) / ( p.diameter3*p.diameter3 - p.diameter4*p.diameter4 );
 	}
-	apert1 = 2.0*( apert1 - apert3 )/x;
-	return sf*apert1*apert1;
+	return sf * 4.0 * (p1 * p3) / (x*x);
 }
 
 double dimm_integrand(double x, void * params) {
@@ -113,12 +114,12 @@
 	return it;
 }
 
-std::pair< vector<double>, vector< std::complex<double> > > spectral_band_fourier( const vector<double>& lambda, const vector<double>& dist, unsigned int center ) {
+std::pair< vector_type, vector< complex_type > > spectral_band_fourier( const vector<double>& lambda, const vector<double>& dist, unsigned int center ) {
 	const int GRIDFACT = 4 ;    // 0 - no stretch, 1 - stretch twice, 2 - stretch 4-times etc.
 	double dlambda = lambda(1) - lambda(0);
 	int dist_ftnn = 2*GRIDFACT*nexthigher(dist.size());
 	vector< std::complex<double> > dist_ft = zero_vector< std::complex<double> >( dist_ftnn );
-	vector<double> dist_freq( dist_ftnn/2 );
+	vector< double > dist_freq( dist_ftnn/2 );
 	for( vector<double>::size_type l = 0; l < dist_freq.size(); l++ ) dist_freq(l) = l/(dlambda*MKM2CM)/dist_ftnn;
 	for( vector<double>::size_type l = 0; l < dist.size(); l++ ) {      // Fill in the data array for FFT:
 		vector<double>::size_type j = ( l < center ) ? dist_ftnn + l - center : l - center;
@@ -128,10 +129,12 @@
 	return std::make_pair( dist_freq, subrange( dist_ft, 0, dist_ftnn/2 ) );
 }
 
-double mass_weight( double height, double eps1, double eps2, double eps3, double eps4, double lambda0, const boost::function1< std::complex<double>, double >& sfun, double max_freq ) {
+namespace detail {
+
+basic_weight::value_type basic_weight::do_calculate( value_type height, value_type eps1, value_type eps2, value_type eps3, value_type eps4 ) const {
 	gsl::integration_workspace w(WORKSPACE);
 	double error;
-	mass_params params(sfun);
+	mass_params params( mass_sfourier_ );
 	gsl_function F;
 	F.function = &mass_integrand;
 	F.params = &params;
@@ -139,31 +142,60 @@
 	params.diameter2 = M_PI*eps2;
 	params.diameter3 = M_PI*eps3;
 	params.diameter4 = M_PI*eps4;
-	params.PIx2lambda0 = 2*M_PI*lambda0*MKM2CM;
+	params.PIx2lambda0 = 2*M_PI * mass_leff_ * MKM2CM;
+	params.powind = POWIND + freq_power_;
 	double result;
 	params.altitude = height*KM2CM;
-	double hi_limit = std::sqrt( 2 * max_freq / params.altitude );
+	double hi_limit = std::sqrt( 2 * mass_freq_max_ / params.altitude );
 	if( std::isfinite( hi_limit ) )
 		gsl::integration_qags( &F, 0, hi_limit, 0, EPS, WORKSPACE, w, &result, &error);
 	else
 		gsl::integration_qagiu( &F, 0, 0, EPS, WORKSPACE, w, &result, &error);
-	return result*CALIBR;
+	return result*CALIBR*std::pow(M_PI*100,freq_power_);
+}
+
+basic_weight::basic_weight( const vector_type& mass_lambda, const vector_type& mass_response, const vector_type& ap_inn, const vector_type& ap_out, const value_type& freq_power ):
+	ap_inn_( ap_inn ),
+	ap_out_( ap_out ),
+	mass_leff_index_( lambda_effective( mass_lambda, element_div( mass_response, mass_lambda*MKM2CM ) ).index() ),
+	mass_leff_( mass_lambda( mass_leff_index_ ) ),
+	mass_sfourier_( spectral_band_fourier( mass_lambda, element_div( mass_response, mass_lambda*MKM2CM ), mass_leff_index_ ) ),
+	freq_power_(freq_power) {
+
+	mass_freq_max_ = mass_sfourier_.first_data()( mass_sfourier_.size() - 1 );
+}
+
+basic_weight::result_type basic_weight::operator() ( const arg_type& height ) const {
+	const size_type mass_nchan = atmos::channels;
+	vector_type ret( atmos::indices_size );
+	size_type cur_n = 0;
+
+	for( size_type i = 0; i < mass_nchan; ++i )
+		ret( cur_n++ ) = do_calculate( height, ap_out_(i), ap_inn_(i), ap_out_(i), ap_inn_(i) );
+
+	for( size_type i = 0; i < mass_nchan; ++i )
+		for( size_type j = i + 1; j < mass_nchan; j++ )
+			ret( cur_n++ ) = do_calculate( height, ap_out_(i), ap_inn_(i), ap_out_(j), ap_inn_(j) );
+
+	return ret;
+}
+
 }
 
-double dimm_weight( double height, double diam, double dbase, double phi, double lambda0, const boost::function1< std::complex<double>, double >& sfun, double max_freq ) {
+massdimm_weight::value_type massdimm_weight::do_dimm_calculate( value_type height, value_type diam, value_type dbase, value_type phi ) const {
 	gsl::integration_workspace w(WORKSPACE);
 	double error;
-	dimm_params params(sfun);
+	dimm_params params(dimm_sfourier_);
 	gsl_function F;
 	F.function = &dimm_integrand;
 	F.params = &params;
 	params.aperture = M_PI*diam;
 	params.dimmbase = 2*M_PI*dbase;
-	params.PIx2lambda0 = 2*M_PI*lambda0*MKM2CM;
+	params.PIx2lambda0 = 2*M_PI*dimm_leff_*MKM2CM;
 	params.phi = phi;
 	double result;
 	params.altitude = height*KM2CM;
-	double hi_limit = std::sqrt( 2 * max_freq / params.altitude );
+	double hi_limit = std::sqrt( 2 * dimm_freq_max_ / params.altitude );
 	if( std::isfinite( hi_limit ) )
 		gsl::integration_qags( &F, 0, hi_limit, 0, EPS, WORKSPACE, w, &result, &error);
 	else
@@ -171,57 +203,26 @@
 	return result*DIMMCONST;
 }
 
-}
-
-weif::massdimm_weight::massdimm_weight( const vector_type& mass_lambda, const vector_type& mass_response, const vector_type& dimm_lambda, const vector_type& dimm_response, const vector_type& ap_inn, const vector_type& ap_out, double dimmsize, double dimmbase ):
-	ap_inn_( ap_inn ), ap_out_( ap_out ), dimmsize_( dimmsize * 100 ), dimmbase_( dimmbase * 100 ),
-	mass_leff_index_( lambda_effective( mass_lambda, element_div( mass_response, mass_lambda*MKM2CM ) ).index() ),
+massdimm_weight::massdimm_weight( const vector_type& mass_lambda, const vector_type& mass_response, const vector_type& dimm_lambda, const vector_type& dimm_response, const vector_type& ap_inn, const vector_type& ap_out, const value_type& dimmsize, const value_type& dimmbase ):
+	basic_weight( mass_lambda, mass_response, ap_inn, ap_out ),
+	dimmsize_( dimmsize * 100 ), dimmbase_( dimmbase * 100 ),
 	dimm_leff_index_( lambda_effective( dimm_lambda, dimm_response ).index() ),
-	mass_leff_( mass_lambda( mass_leff_index_ ) ),
 	dimm_leff_( dimm_lambda( dimm_leff_index_ ) ),
-	mass_sfourier_( spectral_band_fourier( mass_lambda, element_div( mass_response, mass_lambda*MKM2CM ), mass_leff_index_ ) ),
 	dimm_sfourier_( spectral_band_fourier( dimm_lambda, dimm_response, dimm_leff_index_ ) ) {
 
-	mass_freq_max_ = mass_sfourier_.first_data()( mass_sfourier_.size() - 1 );
 	dimm_freq_max_ = dimm_sfourier_.first_data()( dimm_sfourier_.size() - 1 );
 }
 
-weif::massdimm_weight::vector_type weif::massdimm_weight::operator() ( const value_type& height ) const {
-	const size_type mass_nchan = 4;
-	vector_type ret( mass_nchan * (mass_nchan + 1) / 2 + 2 );
-	size_type cur_n = 0;
-
-	for( size_type i = 0; i < mass_nchan; ++i )
-		ret( cur_n++ ) = mass_weight( height, ap_out_(i), ap_inn_(i), 0, 0, mass_leff_, mass_sfourier_, mass_freq_max_ );
-
-	for( size_type i = 0; i < mass_nchan; ++i )
-		for( size_type j = i + 1; j < mass_nchan; j++ )
-			ret( cur_n++ ) = mass_weight( height, ap_out_(i), ap_inn_(i), ap_out_(j), ap_inn_(j), mass_leff_, mass_sfourier_, mass_freq_max_ );
-
-	ret( cur_n++ ) = dimm_weight( height, dimmsize_, dimmbase_, 0.0,    dimm_leff_, dimm_sfourier_, dimm_freq_max_ ) * 1e+11;
-	ret( cur_n++ ) = dimm_weight( height, dimmsize_, dimmbase_, M_PI/2, dimm_leff_, dimm_sfourier_, dimm_freq_max_ ) * 1e+11;
+massdimm_weight::result_type massdimm_weight::operator() ( const arg_type& height ) const {
+	vector_type ret( atmos::indices_size + atmos::dimm_indices_size );
+	size_type cur_n = atmos::indices_size;
+
+	subrange( ret, 0, atmos::indices_size ) = basic_weight::operator() ( height );
+	
+	ret( cur_n++ ) = do_dimm_calculate( height, dimmsize_, dimmbase_, 0.0 ) * 1e+11;
+	ret( cur_n++ ) = do_dimm_calculate( height, dimmsize_, dimmbase_, M_PI/2 ) * 1e+11;
 
 	return ret;
 }
 
-linear_generator::linear_generator( result_type start, result_type step ): start_(start), step_(step) {
-}
-linear_generator::result_type linear_generator::operator()( size_type i ) const {
-	return ( start_ + step_ * i );
-}
-linear_generator::size_type linear_generator::size( result_type max ) const {
-	return static_cast<size_type>( std::floor( (max - start_) / step_ ) );
-}
-
-lg_generator::lg_generator( result_type start, result_type lstep, result_type gstep ): start_(start), t_(lstep), q_(result_type(1) + gstep / 2) {
-	max_ = static_cast< size_type >( std::floor( result_type(2) / gstep - start / lstep ) );
-}
-lg_generator::result_type lg_generator::operator()( size_type i ) const {
-	if( i <= max_ ) return ( start_ + t_ * i );
-	return ( start_ + t_ * max_ ) * std::pow( q_, static_cast<int>(i - max_) ); // FIXME: pow(double,int) for <4.2.1
-}
-lg_generator::size_type lg_generator::size( result_type max ) const {
-	result_type gmax = std::floor( std::log( max / ( start_ + t_ * max_ ) ) / std::log( q_ ) );
-	if( gmax > result_type(0) ) return static_cast< size_type >( gmax ) + max_ + 1;
-	return static_cast< size_type >(std::floor( (max - start_) / t_ ));
-}
\ No newline at end of file
+}}

atmos-2.96.9.tar.gz/src/weif.h -> atmos-2.97.1.tar.gz/src/weif.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: weif.h 158 2010-09-19 08:38:25Z matwey $
+    $Id: weif.h 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2010  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -19,70 +19,66 @@
 #ifndef _WEIF_H
 #define _WEIF_H
 
-#include <boost/numeric/ublas/vector.hpp>
-#include <boost/numeric/ublas/vector_proxy.hpp>
+#include <complex>
 
-#include <utility>
+#include <boost/numeric/ublas/vector.hpp>
 
 #include "spline.h"
 
-namespace weif {
-
-	using namespace std;
+namespace atmos {
+namespace weight {
 	using namespace boost::numeric::ublas;
-	using boost::numeric::ublas::vector;
 
-	class massdimm_weight {
+	typedef double value_type;
+	typedef std::complex< value_type > complex_type;
+	typedef vector< value_type >   vector_type;
+	typedef vector_type::size_type size_type;
+	typedef vector_type            weight_type;
+	
+namespace detail {
+	class basic_weight {
 	public:
-		typedef double value_type;
-		typedef vector< value_type > vector_type;
-		typedef vector_type::size_type size_type;
-		typedef vector_type result_type;
-
-		typedef utils::ublas::cubic_spline< vector_type, vector< std::complex<double> >, vector_type, vector< std::complex<double> > > sfourier_type;
+		typedef atmos::weight::value_type value_type;
+		typedef value_type  arg_type;
+		typedef weight_type result_type;
+
+	protected:
+		typedef utils::ublas::cubic_spline< vector_type, vector< complex_type >, vector_type, vector< complex_type > > sfourier_spline_type;
+		typedef sfourier_spline_type sfourier_type;
+	private:
+		value_type do_calculate( value_type height, value_type eps1, value_type eps2, value_type eps3, value_type eps4 ) const;
 	public:
-		massdimm_weight( const vector_type& mass_lambda, const vector_type& mass_response, const vector_type& dimm_lambda, const vector_type& dimm_response, const vector_type& ap_inn, const vector_type& ap_out, double dimmsize, double dimmbase );
+		basic_weight( const vector_type& mass_lambda, const vector_type& mass_response, const vector_type& ap_inn, const vector_type& ap_out, const value_type& freq_power = value_type(0) );
 
-		vector_type operator() ( const value_type& height ) const;
+		virtual result_type operator() ( const arg_type& height ) const;
 	private:
 		vector_type ap_inn_;
 		vector_type ap_out_;
-		double dimmsize_;
-		double dimmbase_;
 		size_type mass_leff_index_;
-		size_type dimm_leff_index_;
-		double mass_leff_;
-		double dimm_leff_;
+		value_type mass_leff_;
 		sfourier_type mass_sfourier_;
-		sfourier_type dimm_sfourier_;
-		double mass_freq_max_;
-		double dimm_freq_max_;
+		value_type mass_freq_max_;
+		value_type freq_power_;
 	};
-
-	/* make_weight( crv, geometry ) */
 }
 
-    struct linear_generator {
-        typedef double result_type;
-        typedef size_t size_type;
-
-        linear_generator( result_type start, result_type step );
-        result_type operator()( size_type i ) const;
-        size_type size( result_type max ) const;
-    private:
-        result_type start_,step_;
-    };
-
-    struct lg_generator {
-        typedef double result_type;
-        typedef size_t size_type;
-
-        lg_generator( result_type start, result_type lstep, result_type gstep );
-        result_type operator()( size_type i ) const;
-        size_type size( result_type max ) const;
-    private:
-        result_type start_,t_,q_;
-        size_type max_;
-    };
+class massdimm_weight:
+	public detail::basic_weight {
+private:
+	value_type do_dimm_calculate( value_type height, value_type diam, value_type dbase, value_type phi ) const;
+public:
+	massdimm_weight( const vector_type& mass_lambda, const vector_type& mass_response, const vector_type& dimm_lambda, const vector_type& dimm_response, const vector_type& ap_inn, const vector_type& ap_out, const value_type& dimmsize, const value_type& dimmbase );
+	
+	result_type operator() ( const arg_type& height ) const;
+private:
+	value_type dimmsize_;
+	value_type dimmbase_;
+	size_type dimm_leff_index_;
+	value_type dimm_leff_;
+	sfourier_type dimm_sfourier_;
+	value_type dimm_freq_max_;
+};
+
+}}
 
 #endif // _WEIF_H

atmos-2.96.9.tar.gz/src/worksheet.cpp -> atmos-2.97.1.tar.gz/src/worksheet.cpp Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: worksheet.cpp 158 2010-09-19 08:38:25Z matwey $
+    $Id: worksheet.cpp 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -22,11 +22,9 @@
 #include <boost/algorithm/string.hpp>
 #include <boost/bind.hpp>
 #include <boost/filesystem/path.hpp>
-#include <boost/iterator.hpp>
-#include <boost/iterator/filter_iterator.hpp>
 #include <boost/format.hpp>
 #include <boost/lexical_cast.hpp>
-#include <boost/numeric/ublas/vector.hpp>
+#include <boost/numeric/ublas/vector_proxy.hpp>
 
 #include <yalpa/sun.h>
 #include <yalpa/sidtime.h>
@@ -38,7 +36,7 @@
 #include "weif.h" // for wf::*_generator
 #include "worksheet.h"
 
-#include "photometry.h"
+namespace atmos {
 
 struct worksheet::profile_grid {
 	boost::numeric::ublas::vector<double> vec;
@@ -50,7 +48,7 @@
 	public local_context< worksheet > {
 public:
 	typedef context_type::input::time_type time_type;
-	typedef context_type::spectral_factory_type::weight_type weight_type;
+	typedef context_type::weight_factory_type::weight_type weight_type;
 private:
 	std::string m_spectral;
 	yalpa::equatorial<double> m_coords;
@@ -58,7 +56,7 @@
 	object( const yalpa::equatorial<double> ecoords, const std::string& spectral, context_reference_type ctx ): local_context_type( ctx ), m_spectral(spectral), m_coords(ecoords) {}
 
 	double airmass( const time_type& time ) const;
-	const weight_type& spectral();
+	const weight_type& weight();
 	const std::string& sclass() const { return m_spectral; }
 };
 
@@ -68,27 +66,29 @@
 class worksheet::input::background_accumulator:
 	public measurement_accumulator {
 private:
-	std::auto_ptr<background_model> m_model;
-	const ublas::vector< double > m_deadtime;
-	background_nonpoisson_filter m_nonpoisson_filter;
+	std::auto_ptr< atmos::background::model::abstract > m_model;
+	const atmos::parameters::vector_type m_deadtime;
+	atmos::background::filter::nonpoisson m_nonpoisson_filter;
 public:
-	background_accumulator( const boost::posix_time::ptime& time, context_reference_type ctx ):
+	background_accumulator( const time_type& time, context_reference_type ctx ):
 		measurement_accumulator( time, ctx ),
-		m_model( ctx.m_bfactory.create_new_model(time) ),
-		m_deadtime( ctx.m_pproxy.nonlinearity() ),
+		m_model( ctx.m_background_factory.create_new_model(time) ),
+		m_deadtime( ctx.parameters().nonlinearity() ),
 		m_nonpoisson_filter(
 			ctx.m_va"background-threshold".as<double>(),
-			ctx.pproxy().nonpoisson()(3),
-			1000 * ctx.pproxy().basetime() / ctx.pproxy().exposition() ) {
+			ctx.parameters().nonpoisson()(3),
+			1000 * ctx.parameters().basetime() / ctx.parameters().exposition() ) {
 
 		if( !m_model.get() )
 			throw std::logic_error("There is not any background model");
 	};
 	virtual ~background_accumulator() {};
 	virtual void accumulate(  const time_type& time, const flux_type& mean, const vars_type& vars, const rho1_type& rho1, const rho2_type& rho2 ) {
+		using atmos::indices::flux_nonlinearity_correction;
+		
 		if( ! m_nonpoisson_filter( time, mean, vars ) ) {
 			if( m_model.get() ) {
-				m_model->accumulate( time, indices::flux_nonlinearity_correction( mean, m_deadtime ) );
+				m_model->accumulate( time, flux_nonlinearity_correction( mean, m_deadtime ) );
 			}
 		}
 	}
@@ -98,8 +98,8 @@
 		context().m_file_output.background( starttime() );
 
 		if( m_model.get() && m_model->size() > 0 ) {
-			context().m_bfactory.handle_new_model( m_model.release() );
-			context().m_file_output.fluxes( starttime(), context().m_bfactory.background( starttime() ) );
+			context().m_background_factory.handle_new_model( m_model.release() );
+			context().m_file_output.fluxes( starttime(), context().m_background_factory.background( starttime() ) );
 		}
 	}
 };
@@ -108,33 +108,35 @@
 public:
 	typedef context_type::input::size_type size_type;
 private:
-	boost::posix_time::ptime m_accumulation_begin;
-	boost::posix_time::ptime m_accumulation_end;
+	time_type m_accumulation_begin;
+	time_type m_accumulation_end;
 	double m_mean_mz;
 	size_type m_num;
 
-	dimm_storage m_dimm;
-
+	atmos::dimm::storage    m_dimm;
+	atmos::indices::storage m_indices;
+	
 	//FIXME: use params proxy lazy initialization
-	const ublas::vector< double > m_deadtime;
-	const ublas::vector< double > m_scatter;
-	const ublas::vector< double > m_npoisson;
+	const atmos::parameters::vector_type m_deadtime;
+	const atmos::parameters::vector_type m_scatter;
+	const atmos::parameters::vector_type m_npoisson;
 
-	indices_storage m_indices;
-	indices_flux_threshold_filter m_flux_filter;
+	atmos::indices::filter::flux_threshold m_flux_filter;
 public:
-	normal_accumulator( const boost::posix_time::ptime& time, context_reference_type ctx ):
+	normal_accumulator( const time_type& time, context_reference_type ctx ):
 		measurement_accumulator( time, ctx ),
 		m_mean_mz(0.0),
 		m_num(0),
-		m_deadtime( ctx.pproxy().nonlinearity() ),
-		m_scatter( ctx.pproxy().scatter() ),
-		m_npoisson( ctx.pproxy().nonpoisson() ),
+		m_deadtime( ctx.parameters().nonlinearity() ),
+		m_scatter( ctx.parameters().scatter() ),
+		m_npoisson( ctx.parameters().nonpoisson() ),
 		m_flux_filter( ctx.m_va"flux-threshold".as<double>() ) {
 
 	}
 	virtual ~normal_accumulator() {};
 	virtual void accumulate( const time_type& time, const flux_type& mean, const vars_type& vars, const rho1_type& rho1, const rho2_type& rho2 ) {
+		using namespace atmos::indices;
+		
 		if( m_accumulation_begin == boost::posix_time::not_a_date_time )
 			m_accumulation_begin = time;
 		m_accumulation_end = time;
@@ -142,38 +144,42 @@
 		m_num++;
 		m_mean_mz = m_mean_mz * ( double(m_num-1) / double(m_num) ) + context().curobject().airmass( time ) * ( 1.0/m_num );
 
-		indices::flux_type flux_corrected = indices::flux_nonlinearity_correction( mean, m_deadtime );
+		flux_type flux_corrected = flux_nonlinearity_correction( mean, m_deadtime );
 
-		const boost::tuple< vars_type, rho1_type, rho2_type >& vars_corrected = indices::vars_correction(
+		const boost::tuple< vars_type, rho1_type, rho2_type >& vars_corrected = vars_correction(
 			flux_corrected,
 			element_prod( vars, outer_prod( mean, mean ) ),
 			element_prod( rho1, outer_prod( mean, mean ) ),
 			element_prod( rho2, element_prod( mean, mean ) ), m_deadtime, m_npoisson );
 
-		flux_corrected = indices::flux_background_correction( flux_corrected, m_scatter, context().m_bfactory.background( time ) );
+		flux_corrected = flux_background_correction( flux_corrected, m_scatter, context().background_factory().background( time ) );
 
 		if( ! m_flux_filter( flux_corrected ) ) {
-			const boost::tuple< indices::indices_type, indices::indesis_type>& inds = indices::calculate( flux_corrected, vars_corrected.get<0>(), vars_corrected.get<1>(), vars_corrected.get<2>() );
+			const boost::tuple< indices_type, indesis_type>& inds = calculate( flux_corrected, vars_corrected.get<0>(), vars_corrected.get<1>(), vars_corrected.get<2>() );
 
 			m_indices.push_back( time, flux_corrected, inds.get<0>(), inds.get<1>() );
 		}
 	}
 	virtual void dimm_accumulate( const time_type& time, unsigned int framenum, const dimm_flux_type& left, const dimm_flux_type& right, const dimm_vars_type& xvars, const dimm_flux_type& yvars ) {
+		using atmos::dimm::correction;
+		
 		if( m_accumulation_begin == boost::posix_time::not_a_date_time )
 			m_accumulation_begin = time;
 		m_accumulation_end = time;
 
-		m_dimm.push_back( time, framenum, left, right, dimm::correction( xvars ), dimm::correction( yvars ) );
+		m_dimm.push_back( time, framenum, left, right, correction( xvars ), correction( yvars ) );
 	}
 	virtual void flush() {
 		using namespace utils::algo;
+		using namespace atmos;
+		using namespace boost::numeric::ublas;
 
 		const boost::posix_time::ptime& meantime = m_accumulation_begin + (m_accumulation_end-m_accumulation_begin)/2;
 
 		typedef context_type::integrals_factory_type integrals_factory_type;
 		typedef integrals_factory_type::functor_type integrals_functor_type;
 		typedef integrals_functor_type::result_type  integrals_result_type;
-		typedef restoration_functor                   restoration_functor_type;
+		typedef profile::functor                    restoration_functor_type;
 		typedef restoration_functor_type::result_type restoration_result_type;
 
 		context().m_file_output.normal( starttime(), m_mean_mz );
@@ -185,12 +191,17 @@
 			double dimm_turbulence_err = 0.0;
 
 			if( m_dimm.size() > 1 ) {
-				m_dimm.correct_low_vars();
+				dimm::vars_type xvars = average( m_dimm.xvars().begin(), m_dimm.xvars().end(), mean< dimm::vars_type >() );
+				dimm::vars_type yvars = average( m_dimm.yvars().begin(), m_dimm.yvars().end(), mean< dimm::vars_type >() );
 
-				dimm::dimm_vars_type xvars = average( m_dimm.xvars().begin(), m_dimm.xvars().end(), mean< dimm_vars_type >() );
-				dimm::dimm_vars_type yvars = average( m_dimm.yvars().begin(), m_dimm.yvars().end(), mean< dimm_vars_type >() );
-				dimm::dimm_vars_type xerrs = average( m_dimm.xvars().begin(), m_dimm.xvars().end(), error_of_mean< dimm_vars_type >(xvars) );
-				dimm::dimm_vars_type yerrs = average( m_dimm.yvars().begin(), m_dimm.yvars().end(), error_of_mean< dimm_vars_type >(yvars) );
+				/* low vars correction */
+				dimm::correct_low_vars( m_dimm.xvars().begin(), m_dimm.xvars().end(), xvars(0) );
+				dimm::correct_low_vars( m_dimm.yvars().begin(), m_dimm.yvars().end(), yvars(0) );
+				xvars = average( m_dimm.xvars().begin(), m_dimm.xvars().end(), mean< dimm::vars_type >() );
+				yvars = average( m_dimm.yvars().begin(), m_dimm.yvars().end(), mean< dimm::vars_type >() );
+				
+				dimm::vars_type xerrs = average( m_dimm.xvars().begin(), m_dimm.xvars().end(), error_of_mean< dimm::vars_type >(xvars) );
+				dimm::vars_type yerrs = average( m_dimm.yvars().begin(), m_dimm.yvars().end(), error_of_mean< dimm::vars_type >(yvars) );
 
 				double Jlong = dimm::turbulence( context().m_pproxy.dimm_scale(), context().m_pproxy.dimm_Klongitudinal(), xvars(1) );
 				double Jlong_err = dimm::turbulence( context().m_pproxy.dimm_scale(), context().m_pproxy.dimm_Klongitudinal(), xerrs(1) );
@@ -213,16 +224,14 @@
 			context().m_file_output.indices( meantime, indices_avg, indices_err, indesis_avg, indesis_err );
 			context().m_file_output.fluxes( meantime, flux_avg, flux_err );
 
-			integrals_functor_type ifun = context().m_ifactory.get_functor(
-				context().curobject().sclass(),
-				context().curobject().spectral() );
+			integrals_functor_type ifun = context().m_integrals_factory( context().curobject().sclass(), context().curobject().weight() );
 
-			integrals_storage integrals;
-			for( indices_storage::const_iterator it = m_indices.begin(); it != m_indices.end(); ++it ) {
+			integrals::storage integrals;
+			for( indices::storage::const_iterator it = m_indices.begin(); it != m_indices.end(); ++it ) {
 				integrals.push_back( ifun( it->indices, it->indesis, m_mean_mz ) );
 			}
 
-			integrals.unstable_remove_if( integrals_regular_filter() );
+			integrals.unstable_remove_if( integrals::filter::regular() );
 
 			double free_seeing_avg       = average( integrals.free_seeing().begin(),       integrals.free_seeing().end(),       mean< double >() );
 			double isoplanatic_angle_avg = average( integrals.isoplanatic_angle().begin(), integrals.isoplanatic_angle().end(), mean< double >() );
@@ -248,60 +257,34 @@
 				moment_2_avg, moment_2_err );
 
 			/* mixture */
-			typedef vector< double > mixture_type;
-			std::vector< mixture_type > mixture_samples;
-			mixture_samples.reserve( m_indices.size() );
+			mixture::storage mixture_samples( m_indices, m_dimm, m_dimm.size() > 1 ? context().m_pproxy.dimm_scale() : 0 );
 
 			vector< double > grid( context().m_va"profile-grid".as< profile_grid >().vec );
 
 			if( m_dimm.size() > 1 ) {
-				double scale2 = context().m_pproxy.dimm_scale(); scale2*=scale2;
-				dimm_storage::const_iterator dimm_current = m_dimm.begin();
-				dimm_storage::const_iterator dimm_next = dimm_current + 1;
-				indices_storage::const_iterator it = m_indices.begin();
-				while( dimm_current != m_dimm.end() ) {
-					while( it != m_indices.end() ) {
-						if( dimm_next != m_dimm.end() && dimm_next->timestamp <= it->timestamp )
-							break;
-
-						mixture_type mixture_sample( 12 );
-						project( mixture_sample, range(0,10) ) = restoration_functor::flat_reorder( it->indices );
-						mixture_sample(10) = scale2 * dimm_current->xvars(1) * 1e+11;
-						mixture_sample(11) = scale2 * dimm_current->yvars(1) * 1e+11;
-						++it;
-						mixture_samples.push_back( mixture_sample );
-					}
-					dimm_current = dimm_next;
-					++dimm_next;
-				}
-				vector< double > newgrid( grid.size() + 1 );
-				newgrid(0) = 0;
-				subrange( newgrid, 1, grid.size()+1 ) = grid;
-				grid = newgrid;
-			} else {
-				for( indices_storage::const_iterator it = m_indices.begin(); it != m_indices.end(); ++it ){
-					mixture_samples.push_back( restoration_functor::flat_reorder( it->indices ) );
-				}
+				grid.resize( grid.size() + 1 );
+				subrange( grid, 1, grid.size() ) = subrange( grid, 0, grid.size() - 1 );
+				grid(0) = 0;
 			}
+			
 			/* covariate mixture */
-			mixture_type mixture_avg = average( mixture_samples.begin(), mixture_samples.end(), mean< mixture_type >() );
-			symmetric_matrix< double > mixture_covariance = average( mixture_samples.begin(), mixture_samples.end(), covariance< mixture_type >(mixture_avg) );
-			triangular_adaptor< symmetric_matrix< double >, strict_lower > ( mixture_covariance ) *= 0.8;
+			mixture::mixture_type mixture_avg = average( mixture_samples.mixture().begin(), mixture_samples.mixture().end(), mean< mixture::mixture_type >() );
+			symmetric_matrix< double > mixture_covariance = average( mixture_samples.mixture().begin(), mixture_samples.mixture().end(), covariance< mixture::mixture_type >(mixture_avg) );
 			/* prepare grid */
 
 			restoration_functor_type rfun(
 				grid,
 				mixture_covariance,
 				m_mean_mz,
-				context().curobject().spectral()
+				context().curobject().weight()
 				);
 
-			profile_storage profile;
-			for( std::vector< mixture_type >::const_iterator it = mixture_samples.begin(); it != mixture_samples.end(); ++it ) {
-				profile.push_back( rfun( *it ) );
+			profile::storage profile;
+			for( mixture::storage::const_iterator it = mixture_samples.begin(); it != mixture_samples.end(); ++it ) {
+				profile.push_back( rfun( it->mixture ) );
 			}
 
-			profile.unstable_remove_if( chi2_gt_than_filter( context().m_va"chi2-threshold".as<double>() ) );
+			profile.unstable_remove_if( profile::filter::chi2_gt_than( context().m_va"chi2-threshold".as<double>() ) );
 
 			if( profile.size() ) {
 				double chi2_avg = average( profile.chi2().begin(), profile.chi2().end(), mean< double >() );
@@ -326,8 +309,8 @@
 	double secz2 = secz1 - 1.0;
 	return secz1 - secz2*( 0.0018167 + secz2*( 0.002875 + 0.0008083*secz2 ) );
 }
-const worksheet::object::weight_type& worksheet::object::spectral() {
-	return context().m_sfactory.spectral(
+const worksheet::object::weight_type& worksheet::object::weight() {
+	return context().m_weight_factory.weight(
 		m_spectral,
 		context().m_pproxy.inner_diam(),
 		context().m_pproxy.outer_diam(),
@@ -338,6 +321,7 @@
 void validate(boost::any& v, const std::vector<std::string>& values, worksheet::profile_grid* target_type, int) {
 	namespace po = boost::program_options;
 	using namespace boost::numeric;
+	using namespace utils::algo;
 	using namespace utils::ublas;
 
 	typedef std::vector< std::string > split_vector_type;
@@ -387,7 +371,7 @@
 		( "background-threshold,b", po::value<double>()->default_value(4.0), "set background threshold" )
 		( "chi2-threshold,x", po::value<double>()->default_value(25.0), "set chi2 threshold" )
 		( "flux-threshold",  po::value<double>()->default_value(2.0), "set flux threshold" )
-		( "profile-grid,Z",  po::value<profile_grid>()->default_value( profile_grid( utils::ublas::vector_generator< lg_generator >( lg_generator( 0.5, 1.0, 2.0 ),6 ) ), "0.5,1,2,4,8,16" ), "set for profile grid generator" )
+		( "profile-grid,Z",  po::value<profile_grid>()->default_value( profile_grid( utils::ublas::vector_generator< utils::algo::lg_generator >( utils::algo::lg_generator( 0.5, 1.0, 2.0 ),6 ) ), "0.5,1,2,4,8,16" ), "set for profile grid generator" )
 		( "mass-response,r", po::value<std::string>()->default_value("mass_maid.crv"), "file with mass detector response" )
 		( "ccd-response,c",  po::value<std::string>()->default_value("ccd.crv"), "file with ccd response" )
 		( "dump-weight,w",   "dump weight functions to the files")
@@ -414,15 +398,17 @@
 	}
 	return 0;
 }
-background_model* worksheet::select_background_model( const boost::posix_time::ptime& time ) const {
+atmos::background::model::abstract* worksheet::select_background_model( const boost::posix_time::ptime& time ) const {
+	using namespace atmos::background::model;
+	
 	double alt = sun_altitude( time ) * 180 / M_PI; // degrees
 
 	if( alt > -5.0 )
-		return new background_daytime_model( time );
+		return new daytime( time );
 	else if( alt > -12.0 )
-		return new background_twilight_model( time, boost::bind(std::mem_fun(&worksheet::sun_altitude), this, _1) );
+		return new twilight( time, boost::bind(std::mem_fun(&worksheet::sun_altitude), this, _1) );
 	// else
-	return new background_night_model( time );
+	return new night( time );
 }
 double worksheet::sun_altitude( const boost::posix_time::ptime& time ) const {
 	yalpa::ecliptic<double> sun = yalpa::sun_position( time );
@@ -433,8 +419,8 @@
 	m_va(va),
 	m_pcontainer(),
 	m_pproxy( m_pcontainer ),
-	m_sfactory( boost::filesystem::path( va"spectra-dir".as<std::string>() ), boost::filesystem::path( va"mass-response".as<std::string>() ), boost::filesystem::path( va"ccd-response".as<std::string>() ) ),
-	m_bfactory( boost::bind( std::mem_fun( &worksheet::select_background_model ), this, _1) ),
+	m_weight_factory( boost::filesystem::path( va"spectra-dir".as<std::string>() ), boost::filesystem::path( va"mass-response".as<std::string>() ), boost::filesystem::path( va"ccd-response".as<std::string>() ) ),
+	m_background_factory( boost::bind( std::mem_fun( &worksheet::select_background_model ), this, _1) ),
 	m_file_input( std::cin, *this ),
 	m_file_output( std::cout )
 {
@@ -456,8 +442,8 @@
 }
 void worksheet::post() {
 	if( m_va.count("dump-weight") ) {
-		for( spectral_factory_type::const_iterator it = m_sfactory.begin(); it != m_sfactory.end(); ++it ) {
-			spectral::dump_weight( ( boost::format("weight-%1%.wf") % it->first ).str(), it->second );
+		for( weight_factory_type::const_iterator it = weight_factory().begin(); it != weight_factory().end(); ++it ) {
+			atmos::weight::dump_weight( ( boost::format("weight-%1%.wf") % it->first ).str(), it->second );
 		}
 	}
 }
@@ -511,4 +497,4 @@
 	m_maccum->dimm_accumulate( time, framenum, left, right, xvars, yvars );
 }
 
-
+}

atmos-2.96.9.tar.gz/src/worksheet.h -> atmos-2.97.1.tar.gz/src/worksheet.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: worksheet.h 147 2010-04-17 09:25:25Z matwey $
+    $Id: worksheet.h 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -19,58 +19,45 @@
 #ifndef _WORKSHEET_H
 #define _WORKSHEET_H
 
+#include <boost/program_options.hpp>
+
 #include "background_factory.h" // for background_factory
 #include "file_input.h"       // for file_input
 #include "file_output.h"      // for file_output
-#include "params.h"           // for parameters_*
-#include "spectral_factory.h" // for spectral_factory
-#include "typesdef.h"         // for nocase_less
-
 #include "indices.h"
 #include "integrals.h"
+#include "local_context.h"
+#include "params.h"           // for parameters_*
 #include "profile.h"
+#include "spectral_factory.h" // for spectral_factory
 
-#include <boost/date_time/posix_time/posix_time.hpp>
-#include <boost/program_options.hpp>
+namespace atmos {
 
-template<class Context> class local_context {
-	public:
-		typedef Context                context_type;
-		typedef context_type&          context_reference_type;
-		typedef const context_type&    context_const_reference_type;
-		typedef local_context<Context> local_context_type;
-	private:
-		context_reference_type m_ctx;
-	protected:
-		inline context_reference_type       context()       throw() { return m_ctx; }
-		inline context_const_reference_type context() const throw() { return m_ctx; }
-	public:
-		local_context( context_reference_type ctx ) throw() : m_ctx(ctx) {}
-};
+using utils::local_context;
 
 class worksheet {
 	friend class local_context<worksheet>;
 	public:
 		typedef local_context<worksheet> context_type;
-		typedef parameters_container< std::map< std::string, std::string, utils::algo::nocase_less< std::string > > > parameters_container_type;
-		typedef parameters_proxy< parameters_container_type > parameters_proxy_type;
-		typedef spectral_factory spectral_factory_type;
-		typedef background_factory background_factory_type;
-		typedef integrals_factory integrals_factory_type;
+		typedef atmos::parameters::container      parameters_container_type;
+		typedef atmos::parameters::proxy          parameters_proxy_type;
+		typedef atmos::weight::factory            weight_factory_type;
+		typedef atmos::background::factory        background_factory_type;
+		typedef atmos::integrals::functor_factory integrals_factory_type;
 	public:
 		static int parse_program_options( int argc, char** argv, boost::program_options::variables_map& va );
 	public:
 		struct profile_grid;
 	private:
 		class input:
-			public file_input,
+			public atmos::input::file,
 			public local_context< worksheet > {
 		private:
 			class measurement_accumulator: public local_context< worksheet > {
 			private:
-				boost::posix_time::ptime m_starttime;
+				time_type m_starttime;
 			public:
-				measurement_accumulator( const boost::posix_time::ptime& time, context_reference_type ctx ) throw() : local_context_type( ctx ), m_starttime(time) {};
+				measurement_accumulator( const time_type& time, context_reference_type ctx ) throw() : local_context_type( ctx ), m_starttime(time) {};
 				virtual ~measurement_accumulator() = 0;
 				virtual void accumulate( const time_type& time, const flux_type& mean, const vars_type& vars, const rho1_type& rho1, const rho2_type& rho2 ) = 0;
 				virtual void dimm_accumulate( const time_type& time, unsigned int framenum, const dimm_flux_type& left, const dimm_flux_type& right, const dimm_vars_type& xvars, const dimm_flux_type& yvars ) = 0;
@@ -80,12 +67,12 @@
 					return dynamic_cast<const T*>(this);
 				}
 
-				inline const boost::posix_time::ptime& starttime() const { return m_starttime; }
+				inline const time_type& starttime() const { return m_starttime; }
 			};
 
 			class dummy_accumulator: public measurement_accumulator {
 			public:
-				dummy_accumulator( const boost::posix_time::ptime& time, context_reference_type ctx ) throw(): measurement_accumulator( time, ctx ) {};
+				dummy_accumulator( const time_type& time, context_reference_type ctx ) throw(): measurement_accumulator( time, ctx ) {};
 				virtual ~dummy_accumulator() {};
 				virtual void accumulate(  const time_type& time, const flux_type& mean, const vars_type& vars, const rho1_type& rho1, const rho2_type& rho2 ) {}
 				virtual void dimm_accumulate( const time_type& time, unsigned int framenum, const dimm_flux_type& left, const dimm_flux_type& right, const dimm_vars_type& xvars, const dimm_flux_type& yvars ) {}
@@ -97,7 +84,7 @@
 		private:
 			std::auto_ptr<measurement_accumulator> m_maccum;
 
-			template<class T> void switch_accumulator( const boost::posix_time::ptime& time ) {
+			template<class T> void switch_accumulator( const time_type& time ) {
 				std::auto_ptr< measurement_accumulator > new_accumulator( new T( time,context() ) );
 				std::swap( m_maccum, new_accumulator );
 				new_accumulator->flush();
@@ -116,7 +103,7 @@
 			void on_measurements( const time_type& time, const flux_type& mean, const vars_type& vars, const rho1_type& rho1, const rho2_type& rho2 );
 			void on_dimm( const time_type& time, unsigned int framenum, const dimm_flux_type& left, const dimm_flux_type& right, const dimm_vars_type& xvars, const dimm_flux_type& yvars );
 		public:
-			input( std::istream& stm, context_reference_type ctx ): file_input( stm ), local_context_type( ctx ), m_maccum( new dummy_accumulator( time_type() ,ctx) ) {}
+			input( std::istream& stm, context_reference_type ctx ): atmos::input::file( stm ), local_context_type( ctx ), m_maccum( new dummy_accumulator( time_type() ,ctx) ) {}
 		};
 		
 		class object;
@@ -125,29 +112,26 @@
 		
 		parameters_container_type m_pcontainer;
 		parameters_proxy_type     m_pproxy;
-
-		spectral_factory_type     m_sfactory;
-
-		background_factory_type   m_bfactory;
-
-		integrals_factory_type    m_ifactory;
+		weight_factory_type       m_weight_factory;
+		background_factory_type   m_background_factory;
+		integrals_factory_type    m_integrals_factory;
 		
-		input       m_file_input;
-		file_output m_file_output;
+		input        m_file_input;
+		output::file m_file_output;
 
 		std::auto_ptr<object> m_curobject;
 	private:
-		background_model* select_background_model( const boost::posix_time::ptime& time ) const;
+		background::model::abstract* select_background_model( const background::time_type& time ) const;
 	public:
 		worksheet( const boost::program_options::variables_map& va );
 		~worksheet( );
 
-		inline const parameters_proxy_type& pproxy()     const { return m_pproxy;   }
-		inline const spectral_factory_type& sfactory()   const { return m_sfactory; }
-		inline const background_factory_type& bfactory() const { return m_bfactory; }
-		inline const integrals_factory_type& ifactory()  const { return m_ifactory; }
+		inline const parameters_proxy_type&   parameters()         const { return m_pproxy;   }
+		inline const weight_factory_type&     weight_factory()     const { return m_weight_factory; }
+		inline const background_factory_type& background_factory() const { return m_background_factory; }
+		inline const integrals_factory_type&  integrals_factory()  const { return m_integrals_factory; }
 
-		double sun_altitude( const boost::posix_time::ptime& time ) const;
+		double sun_altitude( const background::time_type& time ) const;
 
 		inline object& curobject() {
 			if( ! m_curobject.get() )
@@ -160,4 +144,6 @@
 		void post();
 };
 
+}
+
 #endif // _WORKSHEET_H

atmos-2.96.9.tar.gz/third-party/lsp/include/lsp/Makefile.am -> atmos-2.97.1.tar.gz/third-party/lsp/include/lsp/Makefile.am Changed

atmos-2.96.9.tar.gz/third-party/lsp/include/lsp/Makefile.in -> atmos-2.97.1.tar.gz/third-party/lsp/include/lsp/Makefile.in Changed

atmos-2.96.9.tar.gz/third-party/lsp/include/lsp/bidiagonal_transform.h -> atmos-2.97.1.tar.gz/third-party/lsp/include/lsp/bidiagonal_transform.h Changed

atmos-2.96.9.tar.gz/third-party/lsp/include/lsp/cholesky_decomposition.h -> atmos-2.97.1.tar.gz/third-party/lsp/include/lsp/cholesky_decomposition.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: cholesky_decomposition.h 68 2009-04-20 18:46:20Z matwey $
+    $Id: cholesky_decomposition.h 80 2010-12-10 17:11:54Z matwey $
     Copyright (C) 2008  Matwey V. Kornilov <matwey.kornilov@gmail.com>
 
     This program is free software: you can redistribute it and/or modify
@@ -19,24 +19,23 @@
 #ifndef _CHOLESKY_DECOMPOSITION_H
 #define _CHOLESKY_DECOMPOSITION_H
 
-#include <limits>
-
 #include <boost/numeric/ublas/matrix.hpp>
 #include <boost/numeric/ublas/matrix_proxy.hpp>
 #include <boost/numeric/ublas/banded.hpp>
 #include <boost/numeric/ublas/triangular.hpp>
 #include <boost/numeric/ublas/vector_proxy.hpp>
 
-using namespace boost::numeric::ublas;
+#include <limits>
 
 namespace lsp{
-	
-	
+
+using namespace boost::numeric::ublas;
+
 /**
  *  @brief Upper Cholesky decomposition
  *  @paramin,out m The matrix to be decomposited. Result is stored in this object.
  *
- *  Cholesky decomposition is defined only for symmetric positive semidefined matrixes.
+ *  Cholesky decomposition is defined only for symmetric positive semidefinite matrixes.
  *  This conditions aren't checked.
  *  
  *  The decomposition defined as \f$ C = R^T R \f$ where \f$ R \f$ is upper triangular.

atmos-2.96.9.tar.gz/third-party/lsp/include/lsp/givens_rotation.h -> atmos-2.97.1.tar.gz/third-party/lsp/include/lsp/givens_rotation.h Changed

@@ -1,6 +1,6 @@
 /*
-    $Id: givens_rotation.h 69 2009-08-28 18:32:24Z matwey $
-    Copyright (C) 2008  Matwey V. Kornilov <matwey.kornilov@gmail.com>
+    $Id: givens_rotation.h 80 2010-12-10 17:11:54Z matwey $
+    Copyright (C) 2008,2010  Matwey V. Kornilov <matwey.kornilov@gmail.com>
 
     This program is free software: you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
@@ -19,14 +19,14 @@
 #ifndef _GIVENS_ROTATION_H
 #define _GIVENS_ROTATION_H
 
-#include <boost/numeric/ublas/vector.hpp>
-#include <boost/numeric/ublas/matrix.hpp>
 #include <boost/numeric/ublas/matrix_proxy.hpp>
 
-using namespace boost::numeric::ublas;
+#include <lsp/utils.h>
 
 namespace lsp{
 
+using namespace boost::numeric::ublas;
+
 /**
  *  @class givens_rotation
  *  @brief A functor for the Givens rotation transformation
@@ -133,7 +133,7 @@
  *
  */
 	template<class U> void apply ( U& x, U& y ) const {
-		U w ( x * m_c + y * m_s );
+		typename lsp::temporary_type_traits< U >::type w ( x * m_c + y * m_s );
 		y = x * ( -m_s ) + y * m_c;
 		x = w;
 	}

atmos-2.96.9.tar.gz/third-party/lsp/include/lsp/householder_transform.h -> atmos-2.97.1.tar.gz/third-party/lsp/include/lsp/householder_transform.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: householder_transform.h 60 2009-03-01 18:26:41Z matwey $
+    $Id: householder_transform.h 80 2010-12-10 17:11:54Z matwey $
     Copyright (C) 2008  Matwey V. Kornilov <matwey.kornilov@gmail.com>
 
     This program is free software: you can redistribute it and/or modify
@@ -19,18 +19,19 @@
 #ifndef _HOUSEHOLDER_TRANSFORM_H
 #define _HOUSEHOLDER_TRANSFORM_H
 
-#include <lsp/utils.h>
+#include <boost/numeric/ublas/vector.hpp>
+#include <boost/numeric/ublas/matrix.hpp>
+#include <boost/numeric/ublas/symmetric.hpp>
 
 #include <algorithm>
 #include <limits>
 
-#include <boost/numeric/ublas/vector.hpp>
-#include <boost/numeric/ublas/matrix.hpp>
-
-using namespace boost::numeric::ublas;
+#include <lsp/utils.h>
 
 namespace lsp{
 
+using namespace boost::numeric::ublas;
+
 /**
  *  @class householder_transform
  *  @brief A functor for the Householder transformation
@@ -100,7 +101,7 @@
 			m_s += ( v(p)/w )*( v(p)/w );
 			for( size_type i = l; i < m; ++i )
 				m_s += ( v(i)/w )*( v(i)/w );
-			m_s = ( v(p) < 0 ? 1 : -1 ) * w * std::pow( m_s, 0.5 );
+			m_s = ( v(p) < 0 ? 1 : -1 ) * w * std::sqrt( m_s );
 		} else {
 			m_s = 0;
 		}
@@ -160,6 +161,19 @@
  *  @return \f$ h = v_p - s \f$
  */
 	inline const value_type h() const { return m_h; }
+/**
+ *  @return \f$ u \f$
+ */
+	inline vector_type u() const {
+		vector_type ret = zero_vector< value_type >( m_v.size() );
+		ret(m_p) = m_h;
+		project(ret,range(m_l,ret.size())) = project(m_v,range(m_l,m_v.size()));
+		return ret;
+	}
+/**
+ *  @return \f$ b \equiv \f$
+ */
+	inline const value_type b() const { return m_s * m_h; }
 };
 
 };

atmos-2.96.9.tar.gz/third-party/lsp/include/lsp/least_squares.h -> atmos-2.97.1.tar.gz/third-party/lsp/include/lsp/least_squares.h Changed

atmos-2.96.9.tar.gz/third-party/lsp/include/lsp/nnls.h -> atmos-2.97.1.tar.gz/third-party/lsp/include/lsp/nnls.h Changed

atmos-2.96.9.tar.gz/third-party/lsp/include/lsp/qr_decomposition.h -> atmos-2.97.1.tar.gz/third-party/lsp/include/lsp/qr_decomposition.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: qr_decomposition.h 66 2009-04-12 14:33:06Z matwey $
+    $Id: qr_decomposition.h 80 2010-12-10 17:11:54Z matwey $
     Copyright (C) 2008  Matwey V. Kornilov <matwey.kornilov@gmail.com>
 
     This program is free software: you can redistribute it and/or modify
@@ -19,31 +19,49 @@
 #ifndef _QR_DECOMPOSITION_H
 #define _QR_DECOMPOSITION_H
 
-#include <lsp/givens_rotation.h>
-
+#include <cmath>
 #include <limits>
 
 #include <boost/numeric/ublas/matrix.hpp>
 #include <boost/numeric/ublas/matrix_proxy.hpp>
 #include <boost/numeric/ublas/storage.hpp>
 
-using namespace boost::numeric::ublas;
+#include <lsp/givens_rotation.h>
+#include <lsp/householder_transform.h>
 
 namespace lsp{
 
-namespace {
-	template<class T> static T shift( T q2, T q1, T e2, T e1 ) {
+using namespace boost::numeric::ublas;
+
+struct general_tag {};
+struct bidiag_tag {};
+struct tridiag_symmetric_tag {};
+
+template< class T, class Tag = general_tag > class qr_decomposition;
+
+namespace detail {
+	template<class T> T hypot( const T& x, const T& y ) {
 		typedef T value_type;
-		const value_type f = ( q2*q2 - q1*q1 +  e2*e2 - e1*e1 ) / ( 2*e2*q1 );
-		const value_type t = ( f < 0 ?
-              		- f + std::pow(( value_type(1)+f*f ),0.5) :
-              		- f - std::pow(( value_type(1)+f*f ),0.5) );
-			return ( q2*q2 + e2*(e2 - q1/t) );
-	};
+
+		if( std::abs(x) > std::abs(y) ) {
+			const value_type t = y/x;
+			return std::abs(x) * std::sqrt( value_type(1) + t*t );
+		}
+		if( y == value_type(0) ) {
+			return 0;
+		}
+		const value_type t = x/y;
+		return std::abs(y) * std::sqrt( value_type(1) + t*t );
+	}
+	template<class T> T wilkinson_shift( const T& a1, const T& b1, const T& a2 ) {
+		typedef T value_type;
+		const value_type d = (a1 - a2) / value_type(2);
+		const value_type q = hypot(d,b1);
+		return a2 - b1*b1 / ( d + ( d > 0 ? 1 : -1 ) * q );
+	}
 };
 
 /**
- *  @class qr_decomposition
  *  @brief A functor for the modified QR decomposition
  *
  *  Modified QR decomposition is iterative alghoritm that transformates
@@ -51,7 +69,7 @@
  *  Givens transformations. The decomposition is one of SVD parts.
  *
  */
-template<class T> class qr_decomposition {
+template<class T> class qr_decomposition<T, bidiag_tag> {
 public:
 	typedef T                                  matrix_type;   //!<
 	typedef typename matrix_type::value_type   value_type;    //!<
@@ -104,13 +122,18 @@
 	template<class M1, class M2> void apply( M1& left, M2& right, const range& cell, const regular_tag& ) const {
 		typedef givens_rotation< value_type > givens_rotation_type;
 
-		const value_type lim = std::numeric_limits< value_type >::epsilon() * norm_frobenius( m_matrix ) / m_matrix.size2();
+		const value_type lim = std::numeric_limits< value_type >::epsilon() * norm_inf( m_matrix );
 
 		for( range::const_reverse_iterator it = cell.rbegin() + 1; it != cell.rend(); ++it ) {
-			while( std::abs( m_super( *it ) ) > lim ) {
-
-				value_type en1 = ( *it != cell(0) ? m_super(*it - 1) : value_type(0) );
-				value_type e0 = m_leading( cell(0) ) - shift( m_leading(*it + 1), m_leading(*it), m_super(*it), en1 ) / m_leading( cell(0) );
+			value_type last_super; /* use differential convergation control */
+			do {
+				last_super = m_super( *it );
+
+				const value_type en1 = ( *it != cell(0) ? m_super(*it - 1) : value_type(0) );
+				const value_type a1 = en1*en1 + m_leading(*it)*m_leading(*it);
+				const value_type a2 = m_super(*it)*m_super(*it) + m_leading(*it+1)*m_leading(*it+1);
+				const value_type b1 = m_super(*it)*m_leading(*it);
+				value_type e0 = m_leading( cell(0) ) - detail::wilkinson_shift( a1, b1, a2 ) / m_leading( cell(0) );
 				value_type z = m_super( cell(0) );
 
 				givens_rotation_type gr_left( e0, z );
@@ -132,7 +155,7 @@
 
 					gr_left = givens_rotation_type( m_super(*it2-1), z );
 				}
-			}
+			} while( std::abs( last_super - m_super( *it ) ) > lim );
 			m_super( *it ) = 0;
 		}
 	}
@@ -142,7 +165,7 @@
 		if( cell.size() == 1 ) /* Scalar is in diagonal form */
 			return ;
 
-		const value_type lim = 6 * std::numeric_limits< value_type >::epsilon() * norm_frobenius( m_matrix );
+		const value_type lim = 6 * std::numeric_limits< value_type >::epsilon() * norm_inf( m_matrix );
 
 		/* Looking for the zero diagonal element */
 		for( range::const_reverse_iterator it = cell.rbegin() + 1; it != cell.rend() ; ++it ) {
@@ -220,6 +243,240 @@
 
 };
 
+/**
+ *  @brief A functor for the modified QR decomposition
+ *
+ *  Modified QR decomposition is iterative alghoritm that transformates
+ *  tridiagonal banded symmetric matrix into diagonal form by the number of
+ *  Givens transformations. The decomposition is one of Schur decomposition parts.
+ *
+ */
+template<class T> class qr_decomposition<T, tridiag_symmetric_tag> {
+public:
+	typedef T                                  matrix_type;   //!<
+	typedef typename matrix_type::value_type   value_type;    //!<
+	typedef typename matrix_type::size_type    size_type;     //!<
+	typedef matrix_vector_slice< matrix_type > diagonal_type; //!<
+private:
+	struct regular_tag {};
+	struct left_tag {};
+	struct right_tag {};
+private:
+	matrix_type& m_matrix;
+	mutable diagonal_type m_super;
+	mutable diagonal_type m_leading;
+private:
+
+	template<class M1> void apply( M1& left, const range& cell, const regular_tag& ) const {
+		typedef givens_rotation< value_type > givens_rotation_type;
+
+		const value_type lim = std::numeric_limits< value_type >::epsilon() * norm_inf( m_matrix );
+
+		for( range::const_reverse_iterator it = cell.rbegin() + 1; it != cell.rend(); ++it ) {
+			value_type last_super; /* use differential convergation control */
+			do {
+				last_super = m_super( *it );
+
+				value_type e0 = m_leading( cell(0) ) - detail::wilkinson_shift( m_leading(*it), m_super(*it), m_leading(*it+1) );
+				value_type z = m_super( cell(0) );
+
+				givens_rotation_type gr_left( e0, z );
+				for( range::const_iterator it2 = cell.begin() + 1; *it2 != *it + 2; ++it2 ) {
+
+					value_type t = m_super(*it2-1);
+					gr_left.apply( m_leading(*it2-1), t );
+					gr_left.apply( m_super(*it2-1), m_leading(*it2) );
+					gr_left.apply( m_leading(*it2-1), m_super(*it2-1) );
+					gr_left.apply( t, m_leading(*it2) );
+					gr_left.apply( row(left,*it2-1), row(left,*it2) );
+
+					if( *it2 == *it + 1 ) break;
+					z = gr_left.s() * m_super(*it2);
+					m_super(*it2) = gr_left.c() * m_super(*it2);
+
+					gr_left = givens_rotation_type( m_super(*it2-1), z );
+				}
+			} while( std::abs( last_super - m_super( *it ) ) > lim );
+			m_super( *it ) = 0;
+		}
+	}
+
+	template<class M1> void apply( M1& left, const range& cell ) const {
+
+		if( cell.size() == 1 ) /* Scalar is in diagonal form */
+			return ;
+
+		/* Looking for the zero upper diagonal element */
+		for( range::const_reverse_iterator it = cell.rbegin() + 1; it != cell.rend(); ++it ) {
+			if( m_super( *it ) == 0 ) {
+				apply( left, range( cell.start(), *it + 1 ) );
+				apply( left, range( *it + 1, cell.start() + cell.size() ) );
+				return ;
+			}
+		}
+
+		apply( left, cell, regular_tag() );
+
+	}
+
+public:
+/**
+ *  @brief An object constructor
+ *  @paramin,out matrix
+ *
+ */
+	qr_decomposition( matrix_type& matrix ):
+		m_matrix( matrix ),
+		m_super(   matrix, slice(0, 1, matrix.size2() - 1), slice(1, 1, matrix.size2() - 1) ),
+		m_leading( matrix, slice(0, 1, matrix.size2()),     slice(0, 1, matrix.size2())     ) {
+
+		//assert( matrix.upper() == 1 && matrix.lower() == 0 );
+	}
+/**
+ *  @brief Transformation operaton
+ *  @paramout left  The left matrix
+ *  @paramout right The right matrix
+ *
+ *  The routine calculates and makes transformation.
+ *  Intrinsic assumption is that the all matrix are size-suitable.
+ *
+ *  \f$ M_{left} := Q M_{left} \f$
+ *
+ *  \f$ M_{matrix} := Q M_{matrix} Q^{T} \equiv T \f$
+ *
+ */
+	template<class M1> void apply( M1& left ) const {
+		apply( left, range(0, std::min( m_matrix.size1(), m_matrix.size2() ) ) );
+	}
+
+};
+
+/**
+ *  @brief A functor for the transformation matrix into the upper triangular form using Householder transformations.
+ *
+ *  Any matrix can be transformed into the upper triangular form
+ *  by the //\f$ 2 \cdot n - 1 \f$ Householder transformations.
+ *
+ *  \f
+ *  B = Q A \quad \mbox{where} \quad A \quad \mbox{is initial matrix} \quad Q \quad \mbox{are unitary matrix}
+ *  \f
+ *
+ */
+template<class T> class qr_decomposition<T, general_tag> {
+public:
+	typedef T                                matrix_type; //!< The type of the matrix object to be trasformed
+	typedef typename matrix_type::value_type value_type;  //!< The type of the elements stored in the matrix_type
+	typedef typename matrix_type::size_type  size_type;   //!< The type for seeking in the matrix object
+	typedef vector< value_type >             vector_type; //!< The type of the vector object of matrix elements
+
+	typedef matrix_range< matrix_type > result_type;
+
+private:
+	matrix_type& m_matrix;
+
+public:
+/**
+ *  @brief An object constructor
+ *  @paramin,out matrix The reference to matrix object to be transformed
+ *
+ *  Actual transformation will be performed as soon as apply(M1& left) will be called.
+ *
+ */
+	qr_decomposition( matrix_type& matrix ):
+		m_matrix( matrix ) {
+	}
+
+/**
+ *  @brief Transformation operaton
+ *  @paramout left  The left matrix
+ *
+ *  The routine calculates and makes transformation.
+ *  Intrinsic assumption is that the all matrix are size-suitable.
+ *
+ *  \f$ M_{left} := Q M_{left} \f$
+ *
+ *  \f$ M_{matrix} := Q M_{matrix} \equiv B \f$
+ *
+ */
+        template<class M1> bool apply( M1& left ) const {
+		typedef vector< value_type >                 vector_type;
+		typedef householder_transform< vector_type > householder_transform_type;
+
+		size_type min_size = std::min( m_matrix.size1(), m_matrix.size2() );
+
+		assert( left.size1() == m_matrix.size1() );
+
+		size_type i = 0;
+
+		for( i = 0; i < min_size; ++i ) {
+			householder_transform_type hleft( i+1, i, column(m_matrix,i) );
+			hleft.apply( left, row_major_tag() );
+			hleft.apply( m_matrix, row_major_tag() );
+		}
+
+		value_type lim = matrix_error() * norm_frobenius( m_matrix );
+		size_type  rank = min_size;
+		for( i = 0; i < min_size; ++i ){
+			if( std::abs( m_matrix(i,i) ) < lim ) {
+				m_matrix(i,i) = value_type(0); /* overwhelming rounding errors */
+				rank--;
+			}
+		}
+		return (rank != min_size);
+	}
+
+	template<class M1> void remove( M1& left ) const {
+		typedef givens_rotation< value_type > givens_rotation_type;
+
+		size_type min_size = std::min( m_matrix.size1() - 1, m_matrix.size2() );
+
+		assert( left.size1() == m_matrix.size1() );
+		for( size_type i = 0; i < min_size; ++i ) {
+			if( m_matrix(i+1,i) == value_type(0) ) continue;
+
+			value_type x = m_matrix(i,i);
+			value_type y = m_matrix(i+1,i);
+			givens_rotation_type gr( x, y );
+
+			gr.apply( row(m_matrix,i), row(m_matrix,i+1) );
+			gr.apply( row(left,i), row(left,i+1) );
+
+			m_matrix(i+1,i) = value_type(0); /* overwhelming rounding errors */
+		}
+	}
+
+	template<class M1> bool insert( M1& left ) const {
+		typedef householder_transform< vector_type > householder_transform_type;
+
+		size_type min_size = std::min( m_matrix.size1(), m_matrix.size2() );
+
+		assert( left.size1()  == m_matrix.size1() );
+		assert( min_size == m_matrix.size2() );
+
+		householder_transform_type hleft( min_size, min_size - 1, column( m_matrix, min_size - 1 ) );
+
+		if( hleft.s() == value_type(0) ) /* singularity checking */
+			return true;
+
+		hleft.apply( left, row_major_tag() );
+		hleft.apply( m_matrix, row_major_tag() );
+
+		for( size_type i = min_size; i < m_matrix.size1(); ++i )
+			m_matrix(i, min_size - 1) = value_type(0); /* overwhelming rounding errors */
+
+		return false;
+	}
+
+	value_type left_error() const {
+		size_type min_size = std::min( m_matrix.size1(), m_matrix.size2() );
+		return std::abs( ( 4 * m_matrix.size1() + 32 ) * ( 2 * min_size - 1 ) * std::numeric_limits< value_type >::epsilon() );
+	}
+	value_type matrix_error() const {
+		size_type min_size = std::min( m_matrix.size1(), m_matrix.size2() );
+		return std::abs( ( 6 * m_matrix.size1() - 3 * min_size + 40 ) * ( 2 * min_size - 1 ) * std::numeric_limits< value_type >::epsilon() );
+	}
+};
+
 };
 
 #endif // _QR_DECOMPOSITION_H

atmos-2.97.1.tar.gz/third-party/lsp/include/lsp/schur_decomposition.h Added

@@ -0,0 +1,118 @@
+/*
+    $Id: schur_decomposition.h 80 2010-12-10 17:11:54Z matwey $
+    Copyright (C) 2010  Matwey V. Kornilov <matwey.kornilov@gmail.com>
+
+    This program is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef _SCHUR_DECOMPOSITION_H
+#define _SCHUR_DECOMPOSITION_H
+
+#include <boost/numeric/ublas/banded.hpp>
+#include <boost/numeric/ublas/matrix.hpp>
+#include <boost/numeric/ublas/matrix_proxy.hpp>
+#include <boost/numeric/ublas/vector.hpp>
+
+#include <lsp/tridiagonal_transform.h>
+#include <lsp/qr_decomposition.h>
+#include <lsp/utils.h>
+
+namespace lsp{
+
+using namespace boost::numeric::ublas;
+
+/**
+ *  @class schur_decomposition
+ *  @brief A functor for the Schur decomposition for symmetric matrices
+ *
+ *  Schur is a factorization of symmetric matrix that
+ *  \f
+ *  Q A Q^{T} = T \quad \mbox{where} \quad T = \left| \begin{array}{ccccc}
+ *     \lambda_1   &       &       &       &       \\
+ *           & \lambda_2   &       &       &       \\
+ *           &       &\ddots &       &       \\
+ *           &       &       &\lambda_{n-1}&       \\
+ *           &       &       &       & \lambda_n   \\
+ *  \end{array} \right|,\quad A \quad \mbox{is initial matrix,} \quad U \quad \mbox{is unitary matrix,} \quad \lambda_1 \ge \lambda_2 \ge \dots \ge \lambda_{n-1} \ge \lambda_{n}
+ *  \f
+ *
+ */
+template<class T> class schur_decomposition {
+public:
+	typedef T                                       matrix_type;               //!<
+	typedef typename matrix_type::value_type        value_type;                //!<
+	typedef typename matrix_type::size_type         size_type;                 //!<
+	typedef tridiagonal_transform< matrix_type >    tridiagonal_transform_type; //!<
+	typedef qr_decomposition< matrix_type, tridiag_symmetric_tag > qr_decomposition_type; //!<
+private:
+	matrix_type& m_matrix;
+	mutable tridiagonal_transform_type m_td_trans;
+	mutable qr_decomposition_type m_qr_decomp;
+public:
+/**
+ *  @brief An object constructor
+ *  @paramin,out matrix The reference to matrix object to be decomposited
+ *
+ *  Actual decomposition will be performed as soon as apply(M1& left) will be called.
+ *
+ */
+	schur_decomposition( matrix_type& matrix ):
+		m_matrix( matrix ),
+		m_td_trans( matrix ),
+		m_qr_decomp( matrix ) {
+	}
+
+/**
+ *  @brief Decomposition operaton
+ *  @paramout left  The left matrix
+ *
+ *  The routine decomposites the matrix.
+ *  Intrinsic assumption is that the all matrix are size-suitable.
+ *
+ *  \f$ M_{left} := Q M_{left} \f$
+ *
+ *  \f$ M_{matrix} := T \f$
+ *
+ */
+	template<class M1> void apply( M1& left ) const {
+		typedef matrix_vector_slice< matrix_type > diagonal_type;
+		typedef vector< size_type > permutation_type;
+
+		m_td_trans.apply( left );
+
+		m_qr_decomp.apply( left );
+
+		diagonal_type lambda( m_matrix,
+			slice(0, 1, std::min( m_matrix.size1(), m_matrix.size2() )),
+			slice(0, 1, std::min( m_matrix.size1(), m_matrix.size2() )) );
+		permutation_type pm( lambda.size() );
+		for( typename permutation_type::iterator it = pm.begin(); it != pm.end(); ++it ){
+			*it = it.index();
+		}
+
+		std::sort( pm.begin(), pm.end(), vector_less< diagonal_type, std::greater< typename diagonal_type::value_type > >( lambda ) );
+
+		for( typename permutation_type::size_type it = 0; it != pm.size(); ++it ){
+			if( it < pm(it) ) {
+				row(m_matrix, pm(it)).swap( row(m_matrix, it) );
+				row(left, pm(it)).swap( row(left, it) );
+			}
+		}
+	}
+};
+
+
+};
+
+#endif // _SCHUR_DECOMPOSITION_H

atmos-2.96.9.tar.gz/third-party/lsp/include/lsp/singular_decomposition.h -> atmos-2.97.1.tar.gz/third-party/lsp/include/lsp/singular_decomposition.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: singular_decomposition.h 55 2009-02-05 11:01:48Z matwey $
+    $Id: singular_decomposition.h 80 2010-12-10 17:11:54Z matwey $
     Copyright (C) 2008  Matwey V. Kornilov <matwey.kornilov@gmail.com>
 
     This program is free software: you can redistribute it and/or modify
@@ -28,10 +28,10 @@
 #include <boost/numeric/ublas/matrix_proxy.hpp>
 #include <boost/numeric/ublas/banded.hpp>
 
-using namespace boost::numeric::ublas;
-
 namespace lsp{
 
+using namespace boost::numeric::ublas;
+
 /**
  *  @class singular_decomposition
  *  @brief A functor for the singular value decomposition (SVD)
@@ -57,7 +57,7 @@
 	typedef typename matrix_type::size_type         size_type;                 //!< The type for seeking in the matrix object
 	typedef banded_adaptor< matrix_type >           banded_adaptor_type;       //!< The type of banded adaptor
 	typedef bidiagonal_transform< matrix_type >     bidiagonal_transform_type; //!< The type of a functor for bidiagonal transformation
-	typedef qr_decomposition< banded_adaptor_type > qr_decomposition_type;     //!< The type of a functor for QR decomposition
+	typedef qr_decomposition< banded_adaptor_type, bidiag_tag > qr_decomposition_type;     //!< The type of a functor for QR decomposition
 private:
 	matrix_type& m_matrix;
 	mutable bidiagonal_transform_type m_bd_trans;

atmos-2.97.1.tar.gz/third-party/lsp/include/lsp/tridiagonal_transform.h Added

@@ -0,0 +1,115 @@
+/*
+    $Id: tridiagonal_transform.h 80 2010-12-10 17:11:54Z matwey $
+    Copyright (C) 2010  Matwey V. Kornilov <matwey.kornilov@gmail.com>
+
+    This program is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef _TRIDIAGONAL_TRANSFORM_H
+#define _TRIDIAGONAL_TRANSFORM_H
+
+#include <limits>
+
+#include <boost/numeric/ublas/matrix.hpp>
+#include <boost/numeric/ublas/matrix_proxy.hpp>
+#include <boost/numeric/ublas/vector.hpp>
+#include <boost/numeric/ublas/vector_proxy.hpp>
+
+#include <lsp/householder_transform.h>
+
+namespace lsp{
+
+using namespace boost::numeric::ublas;
+
+namespace detail {
+	template<class M> symmetric_adaptor<M> make_symmetric_adaptor( M& e) {
+		return symmetric_adaptor<M>(e);
+	}
+	template<class M> symmetric_adaptor<const M> make_symmetric_adaptor( const M& e) {
+		return symmetric_adaptor<const M>(e);
+	}
+}
+	
+/**
+ *  @class trdiagonal_transform
+ *  @brief A functor for the transformation symmetric matrix into the tridiagonal form using Householder transformations.
+ *
+ */
+template<class T> class tridiagonal_transform {
+public:
+	typedef T                                matrix_type; //!< The type of the matrix object to be trasformed
+	typedef typename matrix_type::value_type value_type;  //!< The type of the elements stored in the matrix_type
+	typedef typename matrix_type::size_type  size_type;   //!< The type for seeking in the matrix object
+
+private:
+	matrix_type& m_matrix;
+	size_type    m_size;
+
+public:
+/**
+ *  @brief An object constructor
+ *  @paramin,out matrix The reference to matrix object to be transformed
+ *
+ *  Actual transformation will be performed as soon as apply(M1& left, M2& right) will be called.
+ *
+ */
+	tridiagonal_transform( matrix_type& matrix ):
+		m_matrix( matrix ),
+		m_size( matrix.size1() ) {
+
+		assert( matrix.size1() == matrix.size2() );
+	}
+
+/**
+ *  @brief Transformation operaton
+ *  @paramout left  The left matrix
+ *
+ *  The routine calculates and makes transformation.
+ *  Intrinsic assumption is that the all matrix are size-suitable.
+ *
+ *  \f$ M_{left} := Q M_{left} \f$
+  *
+ *  \f$ M_{matrix} := Q M_{matrix} Q^{T} \equiv B \f$
+ *
+ */
+	template<class M1> void apply( M1& left ) const {
+		typedef vector< value_type >                 vector_type;
+		typedef householder_transform< vector_type > householder_transform_type;
+
+		assert( left.size1()  == m_matrix.size1() );
+
+		if( m_size < 3 )
+			return ;
+
+		for( size_type i = 0; i < m_size - 1; ++i ) {
+
+			householder_transform_type hleft( i+2, i+1, column(m_matrix,i) );
+			hleft.apply( column(m_matrix,i) );
+			hleft.apply( left, row_major_tag() );
+			
+			matrix_range< matrix_type > sub_submatrix( m_matrix, range(i+1, m_size), range(i+1, m_size) );
+			value_type  c = value_type(1) / hleft.b();
+			vector_type u = project( hleft.u(), range(i+1,m_size) );
+			vector_type p = prod( sub_submatrix, u );
+			vector_type w = p + c / value_type(2) * inner_prod(p,u) * u;
+
+			sub_submatrix += detail::make_symmetric_adaptor( c * (outer_prod(u,w) + outer_prod(w,u)) );
+		}
+	}
+
+};
+
+};
+
+#endif // _TRIDIAGONAL_TRANSFORM_H

atmos-2.96.9.tar.gz/third-party/lsp/include/lsp/utils.h -> atmos-2.97.1.tar.gz/third-party/lsp/include/lsp/utils.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: utils.h 59 2009-03-01 18:25:07Z matwey $
+    $Id: utils.h 80 2010-12-10 17:11:54Z matwey $
     Copyright (C) 2008  Matwey V. Kornilov <matwey.kornilov@gmail.com>
 
     This program is free software: you can redistribute it and/or modify
@@ -24,10 +24,17 @@
 #include <cmath>
 #include <functional>
 
+#include <boost/numeric/ublas/expression_types.hpp>
+#include <boost/type_traits.hpp>
+#include <boost/utility/enable_if.hpp>
+
 /**
  *  @brief Basic namespace
  */
 namespace lsp {
+
+using namespace boost::numeric::ublas;
+
 /**
  *  @class less_abs
  *  @brief Comparsion of absoulte values
@@ -170,6 +177,19 @@
 
 };
 
+template<class T, class Enable = void> struct temporary_type_traits {
+	typedef T type;
+};
+template<class T> struct temporary_type_traits<T, typename boost::enable_if< boost::is_same< typename T::type_category, scalar_tag > >::type > {
+	typedef typename T::expression_type::value_type type;
+};
+template<class T> struct temporary_type_traits<T, typename boost::enable_if< boost::is_same< typename T::type_category, vector_tag > >::type > {
+	typedef typename T::vector_temporary_type type;
+};
+template<class T> struct temporary_type_traits<T, typename boost::enable_if< boost::is_same< typename T::type_category, matrix_tag > >::type > {
+	typedef typename T::matrix_temporary_type type;
+};
+
 };
 
 #endif // _UTILS_H