EMAN2
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         if (data.size()==0) return;
00060         
00061         std::sort(data.begin(), data.end());
00062 
00063         ymin = FLT_MAX;
00064         ymax = -FLT_MAX;
00065 
00066         typedef vector < Pair >::const_iterator Ptype;
00067         for (Ptype p = data.begin(); p != data.end(); p++) {
00068                 if (p->y > ymax) {
00069                         ymax = p->y;
00070                 }
00071                 if (p->y < ymin) {
00072                         ymin = p->y;
00073                 }
00074         }
00075 
00076         size_t n = data.size();
00077         mean_x_spacing = (data[n - 1].x - data[0].x) / (float) n;
00078 }
00079 
00080 
00081 int XYData::read_file(const string & filename)
00082 {
00083         FILE *in = fopen(filename.c_str(), "rb");
00084         if (!in) {
00085                 throw FileAccessException(filename.c_str());
00086 //              LOGERR("cannot open xydata file '%s'", filename.c_str());
00087 //              return 1;
00088         }
00089 
00090         char buf[MAXPATHLEN];
00091         char tmp_str[MAXPATHLEN];
00092 
00093         while (fgets(buf, MAXPATHLEN, in)) {
00094                 if (buf[0] != '#') {
00095                         float x = 0;
00096                         float y = 0;
00097 
00098                         if (sscanf(buf, " %f%[^.0-9-]%f", &x, tmp_str, &y) != 3) {
00099                                 break;
00100                         }
00101                         data.push_back(Pair(x, y));
00102                 }
00103         }
00104 
00105         fclose(in);
00106         in = 0;
00107 
00108         update();
00109 
00110         return 0;
00111 }
00112 
00113 int XYData::write_file(const string & filename) const
00114 {
00115         FILE *out = fopen(filename.c_str(), "wb");
00116         if (!out) {
00117                 LOGERR("cannot open xydata file '%s' to write", filename.c_str());
00118                 return 1;
00119         }
00120 
00121         typedef vector < Pair >::const_iterator Ptype;
00122         for (Ptype p = data.begin(); p != data.end(); p++) {
00123                 fprintf(out, "%1.6g\t%1.6g\n", p->x, p->y);
00124         }
00125 
00126         fclose(out);
00127         out = 0;
00128 
00129         return 0;
00130 }
00131 
00132 
00133 float XYData::calc_correlation(XYData * xy, float minx, float maxx) const
00134 {
00135         size_t n = data.size();
00136         float x0 = data[0].x;
00137         float xn = data[n - 1].x;
00138 
00139         if (maxx <= minx || minx >= xn || maxx <= x0) {
00140                 LOGERR("incorrect minx, maxx=%f,%f for this XYData range [%f,%f]",
00141                            minx, maxx, x0, xn);
00142                 return 0;
00143         }
00144 
00145         float scc = 0;
00146         float norm1 = 0;
00147         float norm2 = 0;
00148 
00149         xy->update();
00150         for (size_t i = 0; i < n; i++) {
00151                 float x = data[i].x;
00152                 if (x >= minx && x <= maxx && xy->is_validx(x)) {
00153                         float selfy = data[i].y;
00154                         float xyy = xy->get_yatx(x);
00155 
00156                         scc += selfy * xyy;
00157                         norm1 += selfy * selfy;
00158                         norm2 += xyy * xyy;
00159                 }
00160         }
00161 
00162         float result = scc / sqrt(norm1 * norm2);
00163         return result;
00164 }
00165 
00166 
00167 float XYData::get_yatx(float x,bool outzero)
00168 {
00169         if (data.size()==0 || mean_x_spacing==0) return 0.0;
00170 
00171         int nx = (int) data.size();
00172         // Do the range checking up front before we get into trouble
00173         if (x<data[0].x) return outzero?0.0f:data[0].y;
00174         if (x>data[nx-1].x) return outzero?0.0f:data[nx-1].y;
00175         
00176         int s = (int) floor((x - data[0].x) / mean_x_spacing);
00177         if (s>nx-1) s=nx-1;
00178 
00179         // These deal with possibly nonuniform x values. A btree would be better, but this is simple
00180         while (s>0 && data[s].x > x) s--;
00181         while (s<(nx-1) && data[s + 1].x < x ) s++;
00182         if (s>nx-2) return outzero?0.0f:data[nx-1].y;
00183         
00184         float f = (x - data[s].x) / (data[s + 1].x - data[s].x);
00185         float y = data[s].y * (1 - f) + data[s + 1].y * f;
00186         return y;
00187 }
00188 
00189 float XYData::get_yatx_smooth(float x,int smoothing)
00190 {
00191         if (data.size()==0) return 0.0;
00192         if (data.size()==1) return data[0].y;
00193         if (smoothing!=1) throw InvalidParameterException("Only smoothing==1 (linear) currently supported");
00194         
00195         int nx = (int) data.size();
00196         
00197         int s = nx/2; 
00198         if (mean_x_spacing>0) s=(int) floor((x - data[0].x) / mean_x_spacing);
00199         if (s>nx-2) s=nx-2;
00200         else if (s<0) s=0;
00201         else {
00202                 // These deal with possibly nonuniform x values. A btree would be better, but this is simple, and there usually won't be THAT many points
00203                 while (s>0 && data[s].x > x) s--;
00204                 while (s<(nx-2) && data[s + 1].x < x ) s++;
00205         }
00206         
00207         float f = 0,y=0;
00208         if (data[s + 1].x != data[s].x) {
00209                 f= (x - data[s].x) / (data[s + 1].x - data[s].x);
00210                 y = data[s].y * (1 - f) + data[s + 1].y * f;
00211         }
00212         else {
00213                 int s2=s;
00214                 while (data[s2].x==data[s].x) {
00215                         if (s2<nx-1) s2++;
00216                         if (s>0) s--;
00217                         if (s==0 &&s2==nx-1) return data[nx/2].y;
00218                 }
00219                 f= (x - data[s].x) / (data[s2].x - data[s].x);
00220                 y = data[s].y * (1 - f) + data[s2].y * f;
00221                 
00222         }
00223 //      printf("%d %1.2f x %1.2f %1.2f %1.2f y %1.2f %1.2f\n",s,f,x,data[s].x,data[s+1].x,data[s].y,data[s+1].y);
00224         return y;
00225 }
00226 
00227 // FIXME: clearly a slow way to do this
00228 void XYData::insort(float x, float y) {
00229         data.push_back(Pair(x,y));
00230         update();
00231         
00232         
00233         //      int nx = (int) data.size();
00234 //      if (nx==0) { data.push_back(Pair(x,y)); return; }
00235 //      
00236 //      int s = (int) floor((x - data[0].x) / mean_x_spacing);
00237 //      if (s>nx) s=nx;
00238 //      else if (s<0) s=0;
00239 //      else {
00240 //      // These deal with possibly nonuniform x values. A btree would be better, but this is simple, and there usually won't be THAT many points
00241 //              while (s>0 && data[s-1].x > x) s--;
00242 //              while (s<nx && data[s].x <= x ) s++;
00243 //      }
00244 //      data.insert(data.begin()+s,Pair(x,y));
00245 }       
00246 
00247 // data must be sorted before this will work
00248 void XYData::dedupx() {
00249         float acnt=1.0;
00250         for (std::vector<Pair>::iterator it = data.begin() ; it+1 < data.end(); ++it) {
00251                 while (it+1<data.end() && it->x==(it+1)->x) {
00252 //                      printf("%d\t%d\t%1.0f\t%f\t%f\n",it-data.begin(),data.end()-it,acnt,it->x,(it+1)->x);
00253                         it->y+=(it+1)->y;
00254                         acnt+=1.0;
00255                         data.erase(it+1);
00256                 }
00257                 if (acnt>1.0) {
00258                         it->y/=acnt;
00259                         acnt=1.0;
00260                 }
00261         }
00262 }
00263 
00264 void XYData::set_xy_list(const vector<float>& xlist, const vector<float>& ylist)
00265 {
00266         if(xlist.size() != ylist.size()) {
00267                 throw InvalidParameterException("xlist and ylist size does not match!");
00268         }
00269 
00270         for(unsigned int i=0; i<xlist.size(); ++i) {
00271                 data.push_back(Pair(xlist[i], ylist[i]));
00272         }
00273 }
00274 
00275 void XYData::set_size(size_t n)
00276 {
00277         data.resize(n, Pair(0.0f, 0.0f));
00278 }
00279 
00280 vector<float>  XYData::get_state() const {
00281         vector<float> list;
00282         vector<Pair>::const_iterator cit;
00283         for(cit=data.begin(); cit!=data.end(); ++cit) {
00284                 list.push_back( (*cit).x);
00285                 list.push_back( (*cit).y);
00286         }
00287 
00288         return list;
00289         
00290 }
00291 
00292 void  XYData::set_state(vector<float> list) {
00293         if(list.size()%2==1) {
00294                 throw InvalidParameterException("Invalid pickled data");
00295         }
00296 
00297         data.clear();
00298         for(unsigned int i=0; i<list.size(); i+=2) {
00299                 data.push_back(Pair(list[i], list[i+1]));
00300         }
00301 
00302         update();
00303 }
00304 
00305 vector<float> XYData::get_xlist() const
00306 {
00307         vector<float> xlist;
00308         vector<Pair>::const_iterator cit;
00309         for(cit=data.begin(); cit!=data.end(); ++cit) {
00310                 xlist.push_back( (*cit).x);
00311         }
00312 
00313         return xlist;
00314 }
00315 
00316 vector<float> XYData::get_ylist() const
00317 {
00318         vector<float> ylist;
00319         vector<Pair>::const_iterator cit;
00320         for(cit=data.begin(); cit!=data.end(); ++cit) {
00321                 ylist.push_back( (*cit).y);
00322         }
00323 
00324         return ylist;
00325 }