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
00175
00176 while (s>0 && data[s].x > x) s--;
00177 while (s<(nx-1) && data[s + 1].x < x ) s++;
00178
00179 float f = (x - data[s].x) / (data[s + 1].x - data[s].x);
00180 float y = data[s].y * (1 - f) + data[s + 1].y * f;
00181 return y;
00182 }
00183
00184 void XYData::set_xy_list(const vector<float>& xlist, const vector<float>& ylist)
00185 {
00186 if(xlist.size() != ylist.size()) {
00187 throw InvalidParameterException("xlist and ylist size does not match!");
00188 }
00189
00190 for(unsigned int i=0; i<xlist.size(); ++i) {
00191 data.push_back(Pair(xlist[i], ylist[i]));
00192 }
00193 }
00194
00195 void XYData::set_size(size_t n)
00196 {
00197 data.resize(n, Pair(0.0f, 0.0f));
00198 }
00199
00200 vector<float> XYData::get_xlist() const
00201 {
00202 vector<float> xlist;
00203 vector<Pair>::const_iterator cit;
00204 for(cit=data.begin(); cit!=data.end(); ++cit) {
00205 xlist.push_back( (*cit).x);
00206 }
00207
00208 return xlist;
00209 }
00210
00211 vector<float> XYData::get_ylist() const
00212 {
00213 vector<float> ylist;
00214 vector<Pair>::const_iterator cit;
00215 for(cit=data.begin(); cit!=data.end(); ++cit) {
00216 ylist.push_back( (*cit).y);
00217 }
00218
00219 return ylist;
00220 }