dpfoose / FreeIModPoly

A free software implementation of the "Vancouver Raman Algorithm", an improved modified polynomial regression for autofluorescent background correction, for MATLAB/Octave, R, and the Armadillo C++ library
MIT License
3 stars 4 forks source link

Only one iteration #2

Open githubuser0xFFFF opened 6 years ago

githubuser0xFFFF commented 6 years ago

Hi Daniel,

first thank you for your work and your open source C++ IModPoly implementation.

I just tested your C++ IModPoly implementation and at first I thought it would all work well. I could see the calculated basline and the corrected spectrum returned from IModPoly function. Then I wanted to know, how the parameters poly_order, max_it and threshold influence the calculation. So I placed some debug output (two std::cout lines) line into the main loop of the calculation - no other modifications:

do { 
        //always perform at least one regression on the "shrunken" spectrum
        //Polynomial Fitting
        coefs = FreeIModPoly::OrdinaryLeastSquares(X, prev_fit);
        fit = FreeIModPoly::CalcPoly(coefs, new_abscissa);
        //Residual and dev calc (residual calculted inside CalcDev)
        dev = FreeIModPoly::CalcDev(prev_fit, fit);
        std::cout << "prev_dev" << prev_dev << " dev " << dev << std::endl;
        //error criteria
        err = FreeIModPoly::CalcErr(dev, prev_dev);
        //Reconstruction of model input
        fit += dev * ones(rows);
        //if a value in the previous fit is lower than this fit, take the previous
        uvec ind = find (prev_fit < fit);
        fit(ind) = prev_fit(ind);
        prev_fit = fit;
        prev_dev = dev;
        std::cout << "IModPoly regression: " << i << " err: " << err << " threshold "
            << threshold << " max_it " << max_it << std::endl;
        ++i;
} while (err > threshold && (no_max_it || i <= max_it));

I can see, that this loop is executed only once because prev_dev and dev are always equal. This causes the line err = FreeIModPoly::CalcErr(dev, prev_dev); to return an error of 0. And this error in turn is always smaller than the threshold of 0.05 and the while condition (err > threshold) will be false.

I tested it with different spectra, but always with the same result - the loop is executed only once. So now I'm not sure if I do soemthing wrong or if there is a bug or problem in your code. Do you have some spectra data for me for testing that will cause the loop to execute multiple times?

githubuser0xFFFF commented 6 years ago

Ok,

i verified the algorithm by using the simulated data from "Automated Autofluorescence Background Subtraction Algorithm for Biomedical Raman Spectroscopy" and I can confirm that your IModPoly C++ implementation algorithm is broken. The main loops is only executed once.