xydata.cpp

Go to the documentation of this file.
00001 
00005 /*
00006  * Author: Steven Ludtke, 04/10/2003 (sludtke@bcm.edu)
00007  * Copyright (c) 2000-2006 Baylor College of Medicine
00008  *
00009  * This software is issued under a joint BSD/GNU license. You may use the
00010  * source code in this file under either license. However, note that the
00011  * complete EMAN2 and SPARX software packages have some GPL dependencies,
00012  * so you are responsible for compliance with the licenses of these packages
00013  * if you opt to use BSD licensing. The warranty disclaimer below holds
00014  * in either instance.
00015  *
00016  * This complete copyright notice must be included in any revised version of the
00017  * source code. Additional authorship citations may be added, but existing
00018  * author citations must be preserved.
00019  *
00020  * This program is free software; you can redistribute it and/or modify
00021  * it under the terms of the GNU General Public License as published by
00022  * the Free Software Foundation; either version 2 of the License, or
00023  * (at your option) any later version.
00024  *
00025  * This program is distributed in the hope that it will be useful,
00026  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00027  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00028  * GNU General Public License for more details.
00029  *
00030  * You should have received a copy of the GNU General Public License
00031  * along with this program; if not, write to the Free Software
00032  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
00033  *
00034  * */
00035 
00036 #include "xydata.h"
00037 #include <algorithm>
00038 #ifndef WIN32
00039 #include <sys/param.h>
00040 #else
00041 #include <windows.h>
00042 #define MAXPATHLEN (MAX_PATH*4)
00043 #endif
00044 #include <cfloat>
00045 #include <cmath>
00046 #include <cstdio>
00047 #include "log.h"
00048 #include "exception.h"
00049 
00050 using namespace EMAN;
00051 
00052 XYData::XYData()
00053         : ymin(FLT_MAX), ymax(-FLT_MAX), mean_x_spacing(1.0)
00054 {
00055 }
00056 
00057 void XYData::update()
00058 {
00059         std::sort(data.begin(), data.end());
00060 
00061         ymin = FLT_MAX;
00062         ymax = -FLT_MAX;
00063 
00064         typedef vector < Pair >::const_iterator Ptype;
00065         for (Ptype p = data.begin(); p != data.end(); p++) {
00066                 if (p->y > ymax) {
00067                         ymax = p->y;
00068                 }
00069                 if (p->y < ymin) {
00070                         ymin = p->y;
00071                 }
00072         }
00073 
00074         size_t n = data.size();
00075         mean_x_spacing = (data[n - 1].x - data[0].x) / (float) n;
00076 }
00077 
00078 
00079 int XYData::read_file(const string & filename)
00080 {
00081         FILE *in = fopen(filename.c_str(), "rb");
00082         if (!in) {
00083                 LOGERR("cannot open xydata file '%s'", filename.c_str());
00084                 return 1;
00085         }
00086 
00087         char buf[MAXPATHLEN];
00088         char tmp_str[MAXPATHLEN];
00089 
00090         while (fgets(buf, MAXPATHLEN, in)) {
00091                 if (buf[0] != '#') {
00092                         float x = 0;
00093                         float y = 0;
00094 
00095                         if (sscanf(buf, " %f%[^.0-9-]%f", &x, tmp_str, &y) != 3) {
00096                                 break;
00097                         }
00098                         data.push_back(Pair(x, y));
00099                 }
00100         }
00101 
00102         fclose(in);
00103         in = 0;
00104 
00105         update();
00106 
00107         return 0;
00108 }
00109 
00110 int XYData::write_file(const string & filename) const
00111 {
00112         FILE *out = fopen(filename.c_str(), "wb");
00113         if (!out) {
00114                 LOGERR("cannot open xydata file '%s' to write", filename.c_str());
00115                 return 1;
00116         }
00117 
00118         typedef vector < Pair >::const_iterator Ptype;
00119         for (Ptype p = data.begin(); p != data.end(); p++) {
00120                 fprintf(out, "%1.6g\t%1.6g\n", p->x, p->y);
00121         }
00122 
00123         fclose(out);
00124         out = 0;
00125 
00126         return 0;
00127 }
00128 
00129 
00130 float XYData::calc_correlation(XYData * xy, float minx, float maxx) const
00131 {
00132         size_t n = data.size();
00133         float x0 = data[0].x;
00134         float xn = data[n - 1].x;
00135 
00136         if (maxx <= minx || minx >= xn || maxx <= x0) {
00137                 LOGERR("incorrect minx, maxx=%f,%f for this XYData range [%f,%f]",
00138                            minx, maxx, x0, xn);
00139                 return 0;
00140         }
00141 
00142         float scc = 0;
00143         float norm1 = 0;
00144         float norm2 = 0;
00145 
00146         xy->update();
00147         for (size_t i = 0; i < n; i++) {
00148                 float x = data[i].x;
00149                 if (x >= minx && x <= maxx && xy->is_validx(x)) {
00150                         float selfy = data[i].y;
00151                         float xyy = xy->get_yatx(x);
00152 
00153                         scc += selfy * xyy;
00154                         norm1 += selfy * selfy;
00155                         norm2 += xyy * xyy;
00156                 }
00157         }
00158 
00159         float result = scc / sqrt(norm1 * norm2);
00160         return result;
00161 }
00162 
00163 
00164 float XYData::get_yatx(float x,bool outzero)
00165 {
00166         if (data.size()==0 || mean_x_spacing==0) return 0.0;
00167 
00168         int nx = (int) data.size();
00169         // Do the range checking up front before we get into trouble
00170         if (x<data[0].x) return outzero?0.0f:data[0].y;
00171         if (x>data[nx-1].x) return outzero?0.0f:data[nx-1].y;
00172         
00173         int s = (int) floor((x - data[0].x) / mean_x_spacing);
00174         if (s>nx-1) s=nx-1;
00175 
00176         // These deal with possibly nonuniform x values. A btree would be better, but this is simple
00177         while (s>0 && data[s].x > x) s--;
00178         while (s<(nx-1) && data[s + 1].x < x ) s++;
00179         if (s>nx-2) return outzero?0.0f:data[nx-1].y;
00180         
00181         float f = (x - data[s].x) / (data[s + 1].x - data[s].x);
00182         float y = data[s].y * (1 - f) + data[s + 1].y * f;
00183         return y;
00184 }
00185 
00186 void XYData::set_xy_list(const vector<float>& xlist, const vector<float>& ylist)
00187 {
00188         if(xlist.size() != ylist.size()) {
00189                 throw InvalidParameterException("xlist and ylist size does not match!");
00190         }
00191 
00192         for(unsigned int i=0; i<xlist.size(); ++i) {
00193                 data.push_back(Pair(xlist[i], ylist[i]));
00194         }
00195 }
00196 
00197 void XYData::set_size(size_t n)
00198 {
00199         data.resize(n, Pair(0.0f, 0.0f));
00200 }
00201 
00202 vector<float> XYData::get_xlist() const
00203 {
00204         vector<float> xlist;
00205         vector<Pair>::const_iterator cit;
00206         for(cit=data.begin(); cit!=data.end(); ++cit) {
00207                 xlist.push_back( (*cit).x);
00208         }
00209 
00210         return xlist;
00211 }
00212 
00213 vector<float> XYData::get_ylist() const
00214 {
00215         vector<float> ylist;
00216         vector<Pair>::const_iterator cit;
00217         for(cit=data.begin(); cit!=data.end(); ++cit) {
00218                 ylist.push_back( (*cit).y);
00219         }
00220 
00221         return ylist;
00222 }

Generated on Tue Jun 11 12:40:27 2013 for EMAN2 by  doxygen 1.4.7