00001
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
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
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
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 }