SHARE
TWEET

Untitled

a guest May 19th, 2019 86 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. // Copyright (c) 2016, The Regents of the University of California (Regents).
  2. // All rights reserved.
  3.  
  4. // Redistribution and use in source and binary forms, with or without
  5. // modification, are permitted provided that the following conditions are
  6. // met:
  7.  
  8. //    1. Redistributions of source code must retain the above copyright
  9. //       notice, this list of conditions and the following disclaimer.
  10.  
  11. //    2. Redistributions in binary form must reproduce the above
  12. //       copyright notice, this list of conditions and the following
  13. //       disclaimer in the documentation and/or other materials provided
  14. //       with the distribution.
  15.  
  16. //    3. Neither the name of the copyright holder nor the names of its
  17. //       contributors may be used to endorse or promote products derived
  18. //       from this software without specific prior written permission.
  19.  
  20. // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS AS IS
  21. // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  22. // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  23. // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
  24. // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  25. // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  26. // SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  27. // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  28. // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  29. // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  30. // POSSIBILITY OF SUCH DAMAGE.
  31.  
  32. // Please contact the author(s) of this library if you have any questions.
  33. // Author: Michael Kellman   ( kellman@eecs.berkeley.edu )
  34.  
  35. #include <vector>
  36. #include <iostream>
  37. #include <cstdlib>
  38. #include <Eigen/Eigen>
  39. #include <Eigen/Sparse>
  40. #include <time.h>
  41.  
  42. #include "fftw3.h"
  43. #include <tr1/memory>
  44.  
  45. using namespace std;
  46. typedef std::shared_ptr<fftw_complex> cp;
  47.  
  48. using namespace Eigen;
  49. typedef Triplet<double> T;
  50. typedef SparseMatrix<double> SpMat;
  51.  
  52. #include "als.h"
  53.  
  54. void als(cp in, cp out, int N, double smooth, double p, int maxIter)
  55. {
  56.     // Initialize weights
  57.     vector<T> weights;
  58.     VectorXd vweights(N);
  59.     VectorXd sig(N);
  60.     weights.reserve(N);
  61.     for (int ii = 0; ii < N; ++ii)
  62.     {
  63.         weights.push_back(T(ii,ii,1));
  64.         vweights[ii] = 1;
  65.        
  66.         //setup the in vector
  67.         sig[ii] = in.get()[ii][0];
  68.     }
  69.  
  70.     // Initialize Diff matrix
  71.     vector<T> difftriplet;
  72.     difftriplet.reserve(N);
  73.     for (int ii = 0; ii < N-2; ++ii)
  74.     {
  75.         difftriplet.push_back(T(ii,ii,1));
  76.         difftriplet.push_back(T(ii,ii+1,-2));
  77.         difftriplet.push_back(T(ii,ii+2,1));
  78.     }
  79.  
  80.     SpMat D(N,N);
  81.     SpMat W(N,N);
  82.     D.setFromTriplets(difftriplet.begin(),difftriplet.end());
  83.     SpMat Dset(N,N);
  84.     Dset = smooth*D.transpose()*D;
  85.     SpMat temp(N,N);
  86.     VectorXd x(N);
  87.  
  88.     for(int ii = 0; ii < maxIter; ii++)
  89.     {
  90.         // decomposition
  91.         W.setFromTriplets(weights.begin(),weights.end());
  92.         SimplicialCholesky<SpMat> chol(W+Dset);
  93.         x = chol.solve(sig.cwiseProduct(vweights));
  94.         // update
  95.         for(int jj = 0; jj < N; jj++)
  96.         {
  97.             if(sig[jj] > x[jj])
  98.             {
  99.                 vweights[jj] = p;
  100.             }
  101.             else
  102.             {
  103.                 vweights[jj] = 1-p;
  104.             }
  105.             weights[jj] =  T(jj,jj,vweights[jj]);
  106.         }
  107.     }
  108.  
  109.     // subtract the baseline
  110.     for(int ii = 0;ii < N; ii++)
  111.     {
  112.         //out.get()[ii][0] = in.get()[ii][0] - x[ii];
  113.         out.get()[ii][0] = x[ii];
  114.         out.get()[ii][1] = 0; // This should remain 0.
  115.         //cout << in.get()[ii][0] << " : " << out.get()[ii][0] << " : " << x[ii] << endl;
  116.         //cout << out.get()[ii][0] << " : " << x[ii] << endl;
  117.     }
  118. }
  119.  
  120. void findBaseline(double *input, double *output, int length,
  121.                   double smooth, double p, double maxIter)
  122. {
  123.     cp in(new fftw_complex[length]);
  124.  
  125.     for (int i = 0; i < length; i++)
  126.     {
  127.         in.get()[i][0] = input[i];
  128.         in.get()[i][1] = 0;
  129.     }
  130.  
  131.     cp out(new fftw_complex[length]);
  132.  
  133.     als(in, out, length, smooth, p, maxIter);
  134.  
  135.     for (int i = 0; i < length; i++)
  136.     {
  137.         output[i] = out.get()[i][0];
  138.     }
  139. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top