Added mathsym+tcc and boost against all advice

eo/contrib/mathsym/regression/Dataset.cpp

@@ -0,0 +1,133 @@
+/*	    
+ *             Copyright (C) 2005 Maarten Keijzer
+ *
+ *          This program is free software; you can redistribute it and/or modify
+ *          it under the terms of version 2 of the GNU General Public License as 
+ *          published by the Free Software Foundation. 
+ *
+ *          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, write to the Free Software
+ *          Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#include "Dataset.h"
+#include <fstream>
+#include <sstream>
+
+#include <iostream>
+
+using namespace std;
+
+class DataSetImpl {
+    public: 
+    vector< vector<double> > inputs;
+    vector<double> targets;
+
+    void read_data(vector<string> strings) {
+	// find the number of inputs
+	
+	istringstream cnt(strings[0]);
+	unsigned n = 0;
+	for (;;) {
+	    string s;
+	    cnt >> s;
+	    if (!cnt) break;
+	    ++n;
+	}
+
+	inputs.resize(strings.size(), vector<double>(n-1));
+	targets.resize(strings.size());
+
+	for (unsigned i = 0; i < strings.size(); ++i) {
+	    istringstream is(strings[i]);
+	    for (unsigned j = 0; j < n; ++j) {
+		
+		if (!is) {
+		    cerr << "Too few targets in record " << i << endl;
+		    exit(1);
+		}
+		
+		if (j < n-1) {
+		    is >> inputs[i][j];
+		} else {
+		    is >> targets[i];
+		}
+		
+	    }
+	}
+	
+    }
+    
+};
+
+Dataset::Dataset() { pimpl = new DataSetImpl; }
+Dataset::~Dataset() { delete pimpl; }
+Dataset::Dataset(const Dataset& that) { pimpl = new DataSetImpl(*that.pimpl); }
+Dataset& Dataset::operator=(const Dataset& that) { *pimpl = *that.pimpl; return *this; }
+
+unsigned Dataset::n_records() const { return pimpl->targets.size(); }
+unsigned Dataset::n_fields()  const { return pimpl->inputs[0].size(); }
+const std::vector<double>& Dataset::get_inputs(unsigned record) const { return pimpl->inputs[record]; }
+double Dataset::get_target(unsigned record) const { return pimpl->targets[record]; }
+
+double error(string errstr);
+
+void Dataset::load_data(std::string filename) {
+    vector<string> strings; // first load it in strings
+
+    ifstream is(filename.c_str());
+
+    for(;;) {
+	string s;
+	getline(is, s);
+	if (!is) break;
+
+	if (s[0] == '#') continue; // comment, skip
+
+	strings.push_back(s);
+    }
+   
+    is.close();
+
+    if (strings.size() == 0) {
+	error("No data could be loaded");
+    }
+    
+    pimpl->read_data(strings);
+    
+}
+
+std::vector<double> Dataset::input_minima() const {
+    vector<vector<double> >& in = pimpl->inputs;
+    
+    vector<double> mn(in[0].size(), 1e+50);
+    for (unsigned i = 0; i < in.size(); ++i) {
+	for (unsigned j = 0; j < in[i].size(); ++j) {
+	    mn[j] = std::min(mn[j], in[i][j]);
+	}
+    }
+    
+    return mn;
+}
+
+vector<double> Dataset::input_maxima() const {
+    vector<vector<double> >& in = pimpl->inputs;
+    
+    vector<double> mx(in[0].size(), -1e+50);
+    for (unsigned i = 0; i < in.size(); ++i) {
+	for (unsigned j = 0; j < in[i].size(); ++j) {
+	    mx[j] = std::max(mx[j], in[i][j]);
+	}
+    }
+    
+    return mx;
+}
+
+
+
+

eo/contrib/mathsym/regression/Dataset.h

@@ -0,0 +1,51 @@
+/*	    
+ *             Copyright (C) 2005 Maarten Keijzer
+ *
+ *          This program is free software; you can redistribute it and/or modify
+ *          it under the terms of version 2 of the GNU General Public License as 
+ *          published by the Free Software Foundation. 
+ *
+ *          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, write to the Free Software
+ *          Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#ifndef DATASET_H_
+#define DATASET_H_
+
+#include <string>
+#include <vector>
+
+class DataSetImpl;
+
+class Dataset {
+    
+    DataSetImpl* pimpl;
+    
+    Dataset& operator=(const Dataset&); // cannot assign
+    public:
+
+    Dataset();
+    ~Dataset();
+    Dataset(const Dataset&);
+
+    void load_data(std::string filename);
+    
+    unsigned n_records() const;
+    unsigned n_fields() const;
+
+    const std::vector<double>& get_inputs(unsigned record) const;
+    double get_target(unsigned record) const;
+
+    std::vector<double> input_minima() const;
+    std::vector<double> input_maxima() const;
+    
+};
+
+#endif
+

eo/contrib/mathsym/regression/ErrorMeasure.cpp

@@ -0,0 +1,237 @@
+/*	    
+ *             Copyright (C) 2005 Maarten Keijzer
+ *
+ *          This program is free software; you can redistribute it and/or modify
+ *          it under the terms of version 2 of the GNU General Public License as 
+ *          published by the Free Software Foundation. 
+ *
+ *          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, write to the Free Software
+ *          Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+
+#include <vector>
+#include <valarray>
+
+#include "ErrorMeasure.h"
+#include "Dataset.h"
+#include "Sym.h"
+#include "sym_compile.h"
+#include "TargetInfo.h"
+
+using namespace std;
+
+#ifdef INTERVAL_DEBUG
+
+#include <BoundsCheck.h>
+#include <FunDef.h>
+
+vector<double> none;
+IntervalBoundsCheck bounds(none, none);
+
+#endif
+
+
+
+static double not_a_number = atof("nan");
+
+class ErrorMeasureImpl {
+    public:
+	const Dataset& data;
+	TargetInfo train_info;
+	
+	ErrorMeasure::measure measure;
+	
+	Scaling no_scaling;
+	
+	ErrorMeasureImpl(const Dataset& d, double t_p, ErrorMeasure::measure m) : data(d), measure(m) {
+	    
+#ifdef INTERVAL_DEBUG
+	    bounds = IntervalBoundsCheck(d.input_minima(), d.input_maxima());
+#endif
+
+	    unsigned nrecords = d.n_records();
+	    unsigned cases = unsigned(t_p * nrecords);
+	    
+	    valarray<double> t(cases);
+
+	    for (unsigned i = 0; i < cases; ++i) {
+		t[i] = data.get_target(i);
+	    }
+
+	    train_info = TargetInfo(t);
+	    no_scaling = Scaling(new NoScaling);
+	}
+	
+	ErrorMeasure::result eval(const valarray<double>& y) {
+	    
+	    ErrorMeasure::result result;
+	    result.scaling = no_scaling;
+	    
+	    
+	    switch(measure) {
+		case ErrorMeasure::mean_squared:
+		    result.error = pow(train_info.targets() - y, 2.0).sum() / y.size();
+		    return result;
+		case ErrorMeasure::absolute:
+		    result.error = abs(train_info.targets() - y).sum() / y.size();
+		    return result;
+		case ErrorMeasure::mean_squared_scaled:
+		    result.scaling = ols(y, train_info);
+		    result.error  = pow(train_info.targets() - result.scaling->transform(y), 2.0).sum() / y.size();
+		    return result;
+		default: 
+		    cerr << "Unknown measure encountered: " << measure << " " << __FILE__ << " " << __LINE__ << endl;
+	    }
+	    
+	    return result;
+	}   
+	
+	unsigned train_cases() const {
+	    return train_info.targets().size();
+	}
+
+	vector<ErrorMeasure::result> multi_function_eval(const vector<Sym>& pop) {
+	    
+	    multi_function all = compile(pop);
+	    std::vector<double> y(pop.size());
+	    std::vector<double> err(pop.size()); 
+	    
+	    const std::valarray<double>& t = train_info.targets();
+	    
+	    for (unsigned i = 0; i < train_cases(); ++i) {
+		// evaluate
+		all(&data.get_inputs(i)[0], &y[0]);
+
+		for (unsigned j = 0; j < y.size(); ++j) {
+		    double diff = y[j] - t[i];
+		    if (measure == ErrorMeasure::mean_squared) { // branch prediction will probably solve this inefficiency
+			err[j] += diff * diff;
+		    } else {
+			err[j] += fabs(diff);
+		    }
+		    
+		}
+		
+	    }
+	    
+	    std::vector<ErrorMeasure::result> result(pop.size());
+
+	    double n = train_cases();
+	    Scaling no = Scaling(new NoScaling);
+	    for (unsigned i = 0; i < pop.size(); ++i) {
+		result[i].error = err[i] / n;
+		result[i].scaling = no;
+	    }
+
+	    return result;
+	    
+	}
+
+	vector<ErrorMeasure::result> single_function_eval(const vector<Sym> & pop) {
+	    
+	    vector<single_function> funcs(pop.size());
+	    compile(pop, funcs); // get one function pointer for each individual
+
+	    valarray<double> y(train_cases());
+	    vector<ErrorMeasure::result> result(pop.size());
+	    for (unsigned i = 0; i < funcs.size(); ++i) {
+		for (unsigned j = 0; j < train_cases(); ++j) {
+		    y[j] = funcs[i](&data.get_inputs(j)[0]);
+		}
+	
+#ifdef INTERVAL_DEBUG
+		//cout << "eval func " << i << " " << pop[i] << endl;
+		pair<double, double> b = bounds.calc_bounds(pop[i]);
+		
+		// check if y is in bounds
+		for (unsigned j = 0; j < y.size(); ++j) {
+		    if (y[j] < b.first -1e-4 || y[j] > b.second + 1e-4 || !finite(y[j])) {
+			cout << "Error " << y[j] << " not in " << b.first << ' ' << b.second << endl;
+			cout << "Function " << pop[i] << endl;
+			exit(1);
+		    }
+		}
+#endif
+		
+		result[i] = eval(y);
+	    }
+	    
+	    return result; 
+	}
+	
+	vector<ErrorMeasure::result> calc_error(const vector<Sym>& pop) {
+	
+	    // currently we can only accumulate simple measures such as absolute and mean_squared
+	    switch(measure) {
+		case ErrorMeasure::mean_squared:
+		case ErrorMeasure::absolute:
+		    return multi_function_eval(pop);
+		case ErrorMeasure::mean_squared_scaled:
+		    return single_function_eval(pop);
+	    }
+	    
+	    return vector<ErrorMeasure::result>();
+	}
+	
+};
+
+ErrorMeasure::result::result() {
+    error = 0.0;
+    scaling = Scaling(0);
+}
+
+bool ErrorMeasure::result::valid() const {
+    return isfinite(error);
+}
+
+ErrorMeasure::ErrorMeasure(const Dataset& data, double train_perc, measure meas) {
+    pimpl = new ErrorMeasureImpl(data, train_perc, meas);
+}
+
+ErrorMeasure::~ErrorMeasure() { delete pimpl; }
+ErrorMeasure::ErrorMeasure(const ErrorMeasure& that) { pimpl = new ErrorMeasureImpl(*that.pimpl); }
+
+
+ErrorMeasure::result ErrorMeasure::calc_error(Sym sym) {
+   
+    single_function f = compile(sym);
+    
+    valarray<double> y(pimpl->train_cases());
+     
+    for (unsigned i = 0; i < y.size(); ++i) {
+
+	y[i] = f(&pimpl->data.get_inputs(i)[0]);
+	
+	if (!finite(y[i])) {
+	    result res;
+	    res.scaling = Scaling(new NoScaling);
+	    res.error = not_a_number;
+	    return res;
+	}
+	
+    }
+   
+    return pimpl->eval(y); 
+}
+
+vector<ErrorMeasure::result> ErrorMeasure::calc_error(const vector<Sym>& syms) {
+    return pimpl->calc_error(syms);
+    
+}
+
+double ErrorMeasure::worst_performance() const {
+    
+    if (pimpl->measure == mean_squared_scaled) {
+	return pimpl->train_info.tvar();
+    }
+    
+    return 1e+20; // TODO: make this general
+}
+

eo/contrib/mathsym/regression/ErrorMeasure.h

@@ -0,0 +1,61 @@
+/*	    
+ *             Copyright (C) 2005 Maarten Keijzer
+ *
+ *          This program is free software; you can redistribute it and/or modify
+ *          it under the terms of version 2 of the GNU General Public License as 
+ *          published by the Free Software Foundation. 
+ *
+ *          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, write to the Free Software
+ *          Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#ifndef ERROR_MEASURE_H
+#define ERROR_MEASURE_H
+
+#include "Scaling.h"
+
+class ErrorMeasureImpl;
+class Sym;
+class Dataset;
+
+class ErrorMeasure {
+    
+    ErrorMeasureImpl* pimpl;
+    
+    public :
+
+	enum measure {
+	    absolute,
+	    mean_squared,
+	    mean_squared_scaled,
+	};
+	
+	struct result {
+	    double  error;
+	    Scaling scaling;
+	    
+	    result();
+	    bool valid() const;
+	};
+	
+	ErrorMeasure(const Dataset& data, double train_perc, measure meas = mean_squared);
+	
+	~ErrorMeasure();
+	ErrorMeasure(const ErrorMeasure& that);
+	ErrorMeasure& operator=(const ErrorMeasure& that);
+	
+	result calc_error(Sym sym);
+
+	std::vector<result> calc_error(const std::vector<Sym>& sym);
+	
+	double worst_performance() const;
+};
+
+#endif
+

eo/contrib/mathsym/regression/Scaling.cpp

@@ -0,0 +1,417 @@
+/*	    
+ *             Copyright (C) 2005 Maarten Keijzer
+ *
+ *          This program is free software; you can redistribute it and/or modify
+ *          it under the terms of version 2 of the GNU General Public License as 
+ *          published by the Free Software Foundation. 
+ *
+ *          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, write to the Free Software
+ *          Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#include "Scaling.h"
+#include "TargetInfo.h"
+
+using namespace std;
+
+Scaling slope(const std::valarray<double>& x, const TargetInfo& targets) {
+    
+    double xx = 0.0;
+    double xy = 0.0;
+    
+    const valarray<double>& y = targets.targets();
+    
+    for (unsigned i = 0; i < x.size(); ++i) {
+	xx += x[i] * x[i];
+	xy += x[i] * y[i];
+    }
+    
+    if (xx < 1e-7) return Scaling(new LinearScaling(0.0,0.0));
+    
+    double b = xy / xx;
+   
+    return Scaling(new LinearScaling(0.0, b));
+    
+}
+
+// Still needs proper testing with non-trivial lambda
+Scaling regularized_least_squares(const std::valarray<double>& inputs, const TargetInfo& targets, double lambda) {
+    
+    double n = inputs.size();
+    
+    valarray<double> x = inputs;
+
+    double a,b,d;
+    a=b=d=0;
+
+    for (unsigned i = 0; i < n; ++i) {
+	a += 1 + lambda;
+	b += x[i];
+	d += x[i] * x[i] + lambda;
+    }
+    
+    //invert
+    
+    double ad_bc = a*d - b * b;
+    // if ad_bc equals zero there's a problem
+   
+    if (ad_bc < 1e-17) return Scaling(new LinearScaling);
+    
+    double ai = d/ad_bc;
+    double bi = -b/ad_bc;
+    double di = a/ad_bc;
+    double ci = bi;
+    
+    // Now multiply this inverted covariance matrix (C^-1) with x' * t
+    
+    std::valarray<double> ones = x;
+    
+    // calculate C^-1 * x' )
+    for (unsigned i = 0; i < n; ++i) 
+    {
+	ones[i] = (ai + bi * x[i]);
+	x[i]    = (ci + di * x[i]);
+    }
+
+    // results are in [ones, x], now multiply with y
+
+    a = 0.0; // intercept
+    b = 0.0; // slope
+    
+    const valarray<double>& t = targets.targets();
+    
+    for (unsigned i = 0; i < n; ++i)
+    {
+	a += ones[i] * t[i];
+	b += x[i] * t[i];
+    }
+    
+    return Scaling(new LinearScaling(a,b));
+}
+
+Scaling ols(const std::valarray<double>& y, const std::valarray<double>& t) {
+    double n = y.size();
+    
+    double y_mean = y.sum() / n;
+    double t_mean = t.sum() / n;
+   
+    std::valarray<double> y_var = (y - y_mean);
+    std::valarray<double> t_var = (t - t_mean);
+    std::valarray<double> cov = t_var * y_var;
+    
+    y_var *= y_var;
+    t_var *= t_var;
+    
+    double sumvar = y_var.sum();
+    
+    if (sumvar == 0. || sumvar/n < 1e-7 || sumvar/n > 1e+7) // breakout when numerical problems are likely
+	return Scaling(new LinearScaling(t_mean,0.));
+
+    
+    double b = cov.sum() / sumvar;
+    double a = t_mean - b * y_mean;
+    
+    Scaling s = Scaling(new LinearScaling(a,b));
+
+    return s;
+}
+
+Scaling ols(const std::valarray<double>& y, const TargetInfo& targets) {
+    double n = y.size();
+    
+    double y_mean = y.sum() / n;
+    
+    std::valarray<double> y_var = (y - y_mean);
+    std::valarray<double> cov = targets.tcov_part() * y_var;
+    
+    y_var *= y_var;
+
+    double sumvar = y_var.sum();
+    
+    if (sumvar == 0. || sumvar/n < 1e-7 || sumvar/n > 1e+7) // breakout when numerical problems are likely
+	return Scaling(new LinearScaling(targets.tmean(),0.));
+
+    
+    double b = cov.sum() / sumvar;
+    double a = targets.tmean() - b * y_mean;
+    
+    if (!finite(b)) {
+	
+	cout << a << ' ' << b << endl;
+	cout << sumvar << endl;
+	cout << y_mean << endl;
+	cout << cov.sum() << endl;
+	exit(1);
+    }
+	
+    Scaling s = Scaling(new LinearScaling(a,b));
+
+    return s;
+}
+
+
+Scaling wls(const std::valarray<double>& inputs, const TargetInfo& targets) {
+    
+    std::valarray<double> x = inputs;
+    const std::valarray<double>& w = targets.weights();
+    
+    unsigned n = x.size();
+    // First calculate x'*W (as W is a diagonal matrix it's simply elementwise multiplication
+    std::valarray<double> wx = targets.weights() * x;
+    
+    // Now x'*W is contained in [w,wx], calculate x' * W * x (the covariance)
+    double a,b,d;
+    a=b=d=0.0;
+    
+    for (unsigned i = 0; i < n; ++i)
+    {
+	a += w[i];
+	b += wx[i];
+	d += x[i] * wx[i];
+    }
+ 
+    //invert
+    
+    double ad_bc = a*d - b * b;
+    // if ad_bc equals zero there's a problem
+   
+    if (ad_bc < 1e-17) return Scaling(new LinearScaling);
+    
+    double ai = d/ad_bc;
+    double bi = -b/ad_bc;
+    double di = a/ad_bc;
+    double ci = bi;
+    
+    // Now multiply this inverted covariance matrix (C^-1) with x' * W * y
+    
+    // create alias to reuse the wx we do not need anymore
+    std::valarray<double>& ones = wx;
+    
+    // calculate C^-1 * x' * W (using the fact that W is diagonal)
+    for (unsigned i = 0; i < n; ++i) 
+    {
+	ones[i] = w[i]*(ai + bi * x[i]);
+	x[i]    = w[i]*(ci + di * x[i]);
+    }
+
+    // results are in [ones, x], now multiply with y
+
+    a = 0.0; // intercept
+    b = 0.0; // slope
+    
+    const valarray<double>& t = targets.targets();
+    
+    for (unsigned i = 0; i < n; ++i)
+    {
+	a += ones[i] * t[i];
+	b += x[i] * t[i];
+    }
+    
+    return Scaling(new LinearScaling(a,b));
+}
+
+
+//Scaling med(const std::valarray<double>& inputs, const TargetInfo& targets);
+
+double mse(const std::valarray<double>& y, const TargetInfo& t) {
+
+    valarray<double> residuals = t.targets()-y;
+    residuals *= residuals;
+    double sz = residuals.size();
+    if (t.has_weights()) {
+	residuals *= t.weights();
+	sz = 1.0;
+    }
+	
+    return residuals.sum() / sz;
+}
+
+double rms(const std::valarray<double>& y, const TargetInfo& t) {
+    return sqrt(mse(y,t));
+}
+    
+double mae(const std::valarray<double>& y, const TargetInfo& t) {
+    valarray<double> residuals = abs(t.targets()-y);
+    if (t.has_weights()) residuals *= t.weights();
+    return residuals.sum() / residuals.size();
+}
+
+
+/*
+    double standard_error(const std::valarray<double>& y, const std::pair<double,double>& scaling) {
+	double a = scaling.first;
+	double b = scaling.second;
+	double n = y.size();
+	double se = sqrt( pow(a+b*y-current_set->targets,2.0).sum() / (n-2));
+	
+	double mean_y = y.sum() / n;
+	double sxx = pow( y - mean_y, 2.0).sum();
+
+	return se / sqrt(sxx);
+    }
+  
+    double scaled_mse(const std::valarray<double>& y){
+	std::pair<double,double> scaling;
+	return scaled_mse(y,scaling);
+    }
+    
+    double scaled_mse(const std::valarray<double>& y, std::pair<double, double>& scaling)
+    {
+	scaling = scale(y);
+	
+	double a = scaling.first;
+	double b = scaling.second;
+	
+	std::valarray<double> tmp = current_set->targets - a - b * y;
+	tmp *= tmp;
+	
+	if (weights.size())
+	    return (weights * tmp).sum();
+
+	return tmp.sum() / tmp.size();
+    }
+   
+    double robust_mse(const std::valarray<double>& ny, std::pair<double, double>& scaling) {
+	
+	double smse = scaled_mse(ny,scaling);
+
+	std::valarray<double> y = ny;
+	// find maximum covariance case 
+	double n = y.size();
+
+	int largest = 0;
+	
+	{
+	    double y_mean = y.sum() / n;
+	
+	    std::valarray<double> y_var = (y - y_mean);
+	    std::valarray<double> cov = tcov * y_var;
+	
+	    std::valarray<bool> maxcov = cov == cov.max();
+	
+	    for (unsigned i = 0; i < maxcov.size(); ++i) {
+		if (maxcov[i]) {
+		    largest = i;
+		    break;
+		}
+	    }
+	}
+	
+	double y_mean = (y.sum() - y[largest]) / (n-1);
+	y[largest] = y_mean; // dissappears from covariance calculation
+	
+	std::valarray<double> y_var = (y - y_mean);
+	std::valarray<double> cov = tcov * y_var;
+	y_var *= y_var;
+
+	double sumvar = y_var.sum();
+	
+	if (sumvar == 0. || sumvar/n < 1e-7 || sumvar/n > 1e+7) // breakout when numerical problems are likely
+	    return worst_performance();
+	
+	double b = cov.sum() / sumvar;
+	double a = tmean - b * y_mean;
+	
+	std::valarray<double> tmp = current_set->targets - a - b * y;
+	tmp[largest] = 0.0;
+	tmp *= tmp;
+	
+	double smse2 = tmp.sum() / (tmp.size()-1);
+	
+	static std::ofstream os("smse.txt");
+	os << smse << ' ' << smse2 << '\n';
+	
+	if (smse2 > smse) {
+	    return worst_performance();
+	    //std::cerr << "overfit? " << smse << ' ' << smse2 << '\n';
+	}
+	
+	scaling.first = a;
+	scaling.second = b;
+	
+	return smse2;
+    }
+    
+    class Sorter {
+	const std::valarray<double>& scores;
+	public:
+	    Sorter(const std::valarray<double>& _scores) : scores(_scores) {}
+	    
+	    bool operator()(unsigned i, unsigned j) const {
+		return scores[i] < scores[j];
+	    }
+    };
+    
+    double coc(const std::valarray<double>& y) {
+	std::vector<unsigned> indices(y.size());
+	for (unsigned i = 0; i < y.size(); ++i) indices[i] = i;
+	std::sort(indices.begin(), indices.end(), Sorter(y));
+	
+	const std::valarray<double>& targets = current_set->targets;
+	
+	double neg = 1.0 - targets[indices[0]];
+	double pos = targets[indices[0]];
+	
+	double cumpos = 0;
+	double cumneg = 0;
+	double sum=0;
+	
+	double last_score = y[indices[0]];
+	
+	for(unsigned i = 1; i < targets.size(); ++i) {
+	        
+	    if (fabs(y[indices[i]] - last_score) < 1e-9) { // we call it tied
+		pos += targets[indices[i]];
+		neg += 1.0 - targets[indices[i]];
+		
+		if (i < targets.size()-1)
+		    continue;
+	    }
+	    sum += pos * cumneg + (pos * neg) * 0.5;
+	    cumneg += neg;
+	    cumpos += pos;
+	    pos = targets[indices[i]];
+	    neg = 1.0 - targets[indices[i]];
+	    last_score = y[indices[i]];
+	}
+	
+	return sum / (cumneg * cumpos);
+    }
+   
+    // iterative re-weighted least squares (for parameters.classification)
+    double irls(const std::valarray<double>& scores, std::pair<double,double>& scaling) {
+	const std::valarray<double>& t = current_set->targets;
+	
+	std::valarray<double> e(scores.size());
+	std::valarray<double> u(scores.size()); 
+	std::valarray<double> w(scores.size());
+	std::valarray<double> z(scores.size());
+	
+	parameters.use_irls = false; parameters.classification=false;
+	scaling = scale(scores);
+	parameters.use_irls=true;parameters.classification=true;
+	
+	if (scaling.second == 0.0) return worst_performance(); 
+	
+	for (unsigned i = 0; i < 10; ++i) {
+	    e = exp(scaling.first + scaling.second*scores);
+	    u = e / (e + exp(-(scaling.first + scaling.second * scores)));
+	    w = u*(1.-u);
+	    z = (t-u)/w;
+	    scaling = wls(scores, u, w);
+	    //double ll = (log(u)*t + (1.-log(u))*(1.-t)).sum();
+	    //std::cout << "Scale " << i << ' ' << scaling.first << " " << scaling.second << " LL " << 2*ll << std::endl;
+	}
+
+	// log-likelihood
+	u = exp(scaling.first + scaling.second*scores) / (1 + exp(scaling.first + scaling.second*scores));
+	double ll = (log(u)*t + (1.-log(u))*(1.-t)).sum();
+	return 2*ll;
+    }
+*/

eo/contrib/mathsym/regression/Scaling.h

@@ -0,0 +1,93 @@
+/*	    
+ *             Copyright (C) 2005 Maarten Keijzer
+ *
+ *          This program is free software; you can redistribute it and/or modify
+ *          it under the terms of version 2 of the GNU General Public License as 
+ *          published by the Free Software Foundation. 
+ *
+ *          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, write to the Free Software
+ *          Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#ifndef SCALING_H_
+#define SCALING_H_
+
+#include "shared_ptr.h"
+
+#include <valarray>
+#include <iostream>
+#include <string>
+
+class TargetInfo;
+
+class ScalingBase {
+    public:
+
+    std::valarray<double> apply(const std::valarray<double>& x) { 
+	std::valarray<double> xtmp = x;
+	transform(xtmp);
+	return xtmp;
+    }
+	
+    virtual double transform(double input) const = 0;
+    virtual void transform(std::valarray<double>& inputs) const = 0;
+    virtual std::ostream& print(std::ostream& os, std::string str) const = 0;
+    virtual std::valarray<double> transform(const std::valarray<double>& inputs) const = 0;
+};
+
+typedef shared_ptr<ScalingBase> Scaling;
+
+class LinearScaling : public ScalingBase {
+    
+    double a,b;
+    
+    public:
+    LinearScaling() : a(0.0), b(1.0) {}
+    LinearScaling(double _a, double _b) : a(_a), b(_b) {}
+
+    double transform(double input) const { input *=b; input += a; return input; }
+    void transform(std::valarray<double>& inputs) const { inputs *= b; inputs += a; }
+    std::valarray<double> transform(const std::valarray<double>& inputs) const { 
+	std::valarray<double> y = a + b * inputs;
+	return y;
+    }
+    
+    double intercept() const { return a; }
+    double slope()     const { return b; }
+    
+    std::ostream& print(std::ostream& os, std::string str) const {
+	os.precision(16);
+	os << a << " + " << b << " * " << str;
+	return os;
+    }
+};
+
+class NoScaling : public ScalingBase{
+    void transform(std::valarray<double>&) const {}
+    double transform(double input) const { return input; }
+    std::valarray<double> transform(const std::valarray<double>& inputs) const { return inputs; }
+    std::ostream& print(std::ostream& os, std::string str) const { return os << str; }
+};
+
+extern Scaling slope(const std::valarray<double>& inputs, const TargetInfo& targets); // slope only
+extern Scaling ols(const std::valarray<double>& inputs, const TargetInfo& targets);
+extern Scaling wls(const std::valarray<double>& inputs, const TargetInfo& targets);
+extern Scaling med(const std::valarray<double>& inputs, const TargetInfo& targets);
+
+extern Scaling ols(const std::valarray<double>& inputs, const std::valarray<double>& outputs);
+
+extern double mse(const std::valarray<double>& y, const TargetInfo& t);
+extern double rms(const std::valarray<double>& y, const TargetInfo& t);
+extern double mae(const std::valarray<double>& y, const TargetInfo& t);
+
+// Todo Logistic Scaling
+
+#endif
+
+

eo/contrib/mathsym/regression/TargetInfo.cpp

@@ -0,0 +1,138 @@
+/*	    
+ *             Copyright (C) 2005 Maarten Keijzer
+ *
+ *          This program is free software; you can redistribute it and/or modify
+ *          it under the terms of version 2 of the GNU General Public License as 
+ *          published by the Free Software Foundation. 
+ *
+ *          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, write to the Free Software
+ *          Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#include "TargetInfo.h"
+
+using namespace std;
+
+TargetInfo::TargetInfo(const TargetInfo& org) { operator=(org); }
+
+TargetInfo& TargetInfo::operator=(const TargetInfo& org) {
+    _targets.resize(org._targets.size());
+    _weights.resize(org._weights.size());
+    _tcov_part.resize(org._tcov_part.size());
+    
+    _targets = org._targets;
+    _weights = org._weights;
+    _tcov_part = org._tcov_part;
+
+    _tmean = org._tmean;
+    _tvar  = org._tvar;
+    _tstd  = org._tstd;
+    _tmed = org._tmed;
+    return *this;
+}
+      
+
+TargetInfo::TargetInfo(const std::valarray<double>& t) {
+    _weights.resize(0);
+    _targets.resize(t.size());
+    _targets = t;
+    
+    _tmean = _targets.sum()/_targets.size();
+    
+    _tcov_part.resize(_targets.size());
+    _tcov_part = _targets;
+    _tcov_part -= _tmean;
+	
+    std::valarray<double> tmp = _tcov_part;
+    tmp = _tcov_part;
+    tmp *= tmp;
+	
+    _tvar = tmp.sum() / (tmp.size()-1);
+    _tstd = sqrt(_tvar);
+    _tmed = 0;
+}
+
+TargetInfo::TargetInfo(const std::valarray<double>& t, const std::valarray<double>& w) {
+
+    _targets.resize(t.size());
+    _weights.resize(w.size());
+
+    _targets = t;
+    _weights = w;
+    
+    double sumw = _weights.sum();
+	// scale weights so that they'll add up to 1
+    _weights /= sumw;
+	
+    _tmean = (_targets * _weights).sum();
+    _tcov_part.resize(_targets.size());
+    _tcov_part = _targets;
+    _tcov_part -= _tmean;
+
+    _tvar = (pow(_targets - _tmean, 2.0) * _weights).sum();
+    _tstd = sqrt(_tvar);
+    _tmed = 0.;
+}
+
+// calculate the members, now in the context of a mask
+void TargetInfo::set_training_mask(const std::valarray<bool>& tmask) {
+    
+    TargetInfo tmp;
+    
+    if (has_weights() ) {
+	tmp = TargetInfo( _targets[tmask], _weights[tmask]);
+    } else {
+	tmp = TargetInfo( _targets[tmask] );
+    }
+    
+    _tcov_part.resize(tmp._tcov_part.size());
+    _tcov_part = tmp._tcov_part;
+
+    _tmean = tmp._tmean;
+    _tvar  = tmp._tvar;
+    _tstd  = tmp._tstd;
+    _tmed =  tmp._tmed;
+
+    _training_mask.resize(tmask.size());
+    _training_mask = tmask;
+}
+
+struct SortOnTargets
+{
+    const valarray<double>& t;
+    SortOnTargets(const valarray<double>& v) : t(v) {}
+
+    bool operator()(int i, int j) const {
+	return fabs(t[i]) < fabs(t[j]);
+    }
+};
+    
+vector<int> TargetInfo::sort() {
+    
+    vector<int> ind(_targets.size());
+    for (unsigned i = 0; i < ind.size(); ++i) { ind[i] = i; }
+
+    std::sort(ind.begin(), ind.end(), SortOnTargets(_targets));
+
+    valarray<double> tmptargets = _targets;
+    valarray<double> tmpweights = _weights;
+    valarray<double> tmpcov     = _tcov_part;
+    
+    for (unsigned i = 0; i < ind.size(); ++i) 
+    {
+	_targets[i] = tmptargets[ ind[i] ];
+	_tcov_part[i] = tmpcov[ ind[i] ];	
+	if (_weights.size()) _weights[i] = tmpweights[ ind[i] ];
+    }
+
+    return ind;
+}
+
+
+

eo/contrib/mathsym/regression/TargetInfo.h

@@ -0,0 +1,65 @@
+/*	    
+ *             Copyright (C) 2005 Maarten Keijzer
+ *
+ *          This program is free software; you can redistribute it and/or modify
+ *          it under the terms of version 2 of the GNU General Public License as 
+ *          published by the Free Software Foundation. 
+ *
+ *          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, write to the Free Software
+ *          Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#ifndef TARGETINFO_H_
+#define TARGETINFO_H_
+
+#include <valarray>
+#include <vector>
+
+class TargetInfo {
+    std::valarray<double> _targets;
+    std::valarray<double> _weights;
+    std::valarray<bool>   _training_mask;
+
+    // some stuff for ols
+    std::valarray<double> _tcov_part;
+    double _tmean;
+    double _tvar;
+    double _tstd;
+    double _tmed;
+
+    public:
+    TargetInfo() {}
+
+    TargetInfo(const std::valarray<double>& t);
+    TargetInfo(const std::valarray<double>& t, const std::valarray<double>& w);
+
+    TargetInfo(const TargetInfo& org);
+    TargetInfo& operator=(const TargetInfo& org);
+    ~TargetInfo() {}
+
+    const std::valarray<double>& targets() const { return _targets; }
+    const std::valarray<double>& weights() const { return _weights; }
+    const std::valarray<bool>&   mask() const { return _training_mask; }
+
+    void set_training_mask(const std::valarray<bool>& mask);
+
+    bool has_weights() const { return _weights.size(); }
+    bool has_mask()    const { return _training_mask.size(); }
+
+    std::vector<int> sort();
+
+    const std::valarray<double>& tcov_part() const { return _tcov_part; }
+    double tmean() const { return _tmean; }
+    double tvar()  const { return _tvar; }
+    double tstd()  const { return _tstd; }
+    double devmedian() const { return _tmed; }
+};
+
+#endif
+