Skip to content

Commit 1488fb9

Browse files
andcastilloandcastillo
committed
Debuging
Updated the dependencies. Should work with the last version of optimise-lorentzian, that is pointing to my dev version. Working for non overlaped peaks. Still needing to deal with the overlap.
1 parent d2b4816 commit 1488fb9

File tree

5 files changed

+186
-24
lines changed

5 files changed

+186
-24
lines changed

‎package.json‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -32,7 +32,7 @@
3232
"ml-stat": "^1.0.1"
3333
},
3434
"dependencies": {
35-
"ml-optimize-lorentzian": "0.0.1",
35+
"ml-optimize-lorentzian": "../optimize-lorentzian",
3636
"ml-stat": "^1.0.1",
3737
"xy-parser": "^1.1.0"
3838
}

‎src/gsdLight.js‎

Lines changed: 8 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -5,13 +5,6 @@ function gsdLight(x, y, options){
55
if (options.noiseLevel===undefined) options.noiseLevel=0;
66
if (options.maxCriteria===undefined) options.maxCriteria=true;
77

8-
//TODO Find a better way to do it.
9-
if (options.functionType === undefined) options.functionType = "lorentzian";
10-
//Lorentzian by default
11-
var widthCorrection = 1.04720;// From http://mathworld.wolfram.com/LorentzianFunction.html
12-
if(options.functionType === "gaussian")
13-
widthCorrection = 1.17741; // From https://en.wikipedia.org/wiki/Gaussian_function#Properties
14-
158

169
if (options.noiseLevel>0) {
1710
y=[].concat(y);
@@ -30,25 +23,26 @@ function gsdLight(x, y, options){
3023
var Y = new Array(size);
3124
var dY = new Array(size);
3225
var ddY = new Array(size);
33-
var dX = new Array(size);
26+
//var dX = new Array(size);
3427
var dx = x[1]-x[0];
3528

3629
for (var j = 2; j < size+2; j++) {
37-
dx = 0;
30+
/*dx = 0;
3831
dX[j-2] = x[j]-x[j-1];
3932
for(k=-1;k<=2;k++)
4033
dx+=x[j-k+1]-x[j-k];
41-
dx/=4;
34+
dx/=4;*/
35+
dx = x[j]-x[j-1];
4236
Y[j-2]=(1/35.0)*(-3*y[j-2] + 12*y[j-1] + 17*y[j] + 12*y[j+1] - 3*y[j+2]);
4337
X[j-2]=x[j];
4438
dY[j-2]=(1/(12*dx))*(y[j-2] - 8*y[j-1] + 8*y[j+1] - y[j+2]);
4539
ddY[j-2]=(1/(7*dx*dx))*(2*y[j-2] - y[j-1] - 2*y[j] - y[j+1] + 2*y[j+2]);
4640

47-
if(Math.max(dx,dX[j-2])/Math.min(dx,dX[j-2])>1.00001){
41+
/*if(Math.max(dx,dX[j-2])/Math.min(dx,dX[j-2])>1.00001){
4842
4943
Y[j-2]=0;
5044
51-
}
45+
}*/
5246

5347
}
5448

@@ -71,6 +65,7 @@ function gsdLight(x, y, options){
7165
//console.log(dx);
7266
//By the intermediate value theorem We cannot find 2 consecutive maxima or minima
7367
for (var i = 1; i < Y.length -1 ; i++){
68+
//console.log(dY[i]);
7469
if ((dY[i] < dY[i-1]) && (dY[i] <= dY[i+1])||
7570
(dY[i] <= dY[i-1]) && (dY[i] < dY[i+1])) {
7671
lastMin = X[i];
@@ -128,7 +123,7 @@ function gsdLight(x, y, options){
128123
signals.push({
129124
x: frequency,
130125
y: height,
131-
width: linewidth*widthCorrection
126+
width: linewidth//*widthCorrection
132127
})
133128
}
134129
}

‎src/index.js‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -73,7 +73,7 @@ function gsd(x, y, options){
7373
var intervals = new Array();
7474
var lastMax = null;
7575
var lastMin = null;
76-
console.log("dx "+dx);
76+
//console.log("dx "+dx);
7777
for (var i = 1; i < Y.length -1 ; i++){
7878
if ((dY[i] < dY[i-1]) && (dY[i] <= dY[i+1])||
7979
(dY[i] <= dY[i-1]) && (dY[i] < dY[i+1])) {

‎src/optimize.js‎

Lines changed: 139 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,139 @@
1+
/**
2+
* Created by acastillo on 9/6/15.
3+
*/
4+
var Opt = require("ml-optimize-lorentzian");
5+
6+
function sampleFunction(from, to, x, y, lastIndex){
7+
//console.log(from+" "+to);
8+
//console.log(lastIndex+" "+x[lastIndex[0]]);
9+
var nbPoints = x.length;
10+
var sampleX = [];
11+
var sampleY = [];
12+
var direction = Math.sign(x[1]-x[0]);//Direction of the derivative
13+
var delta = (to-from)/2;
14+
var mid = (from+to)/2;
15+
var stop = false;
16+
var index = lastIndex[0];
17+
while(!stop&&index<nbPoints){
18+
if(Math.abs(x[index]-mid)<=delta){
19+
sampleX.push(x[index]);
20+
sampleY.push(y[index]);
21+
index++;
22+
23+
}
24+
//It is outside the range.
25+
else{
26+
27+
if(Math.sign(mid-x[index])==direction){
28+
//We'll reach the mid going in the current direction
29+
index++;
30+
}
31+
else{
32+
//There is not more peaks in the current range
33+
stop=true;
34+
}
35+
}
36+
//console.log(sampleX);
37+
}
38+
lastIndex[0]=index;
39+
return [sampleX, sampleY];
40+
}
41+
42+
function optimizePeaks(peakList,x,y,n){
43+
var i, j, lastIndex=[0];
44+
var groups = groupPeaks(peakList,n);
45+
var result = [];
46+
//console.log(x[0]+" "+x[1]);
47+
for(i=0;i<groups.length;i++){
48+
//console.log(peakList[i]);
49+
var peaks = groups[i].group;
50+
if(peaks.length>1){
51+
//Multiple peaks
52+
console.log("Pending group of overlaped peaks "+peaks.length);
53+
}
54+
else{
55+
//Single peak
56+
peaks = peaks[0];
57+
var sampling = sampleFunction(peaks.x-n*peaks.width,
58+
peaks.x+n*peaks.width,x,y,lastIndex);
59+
//console.log(sampling);
60+
if(sampling[0].length>5){
61+
var error = peaks.width/10000;
62+
var opts = [ 3, 100, error, error, error, error*10, error*10, 11, 9, 1 ];
63+
//var gauss = Opt.optimizeSingleGaussian(sampling[0], sampling[1], opts, peaks);
64+
var gauss = Opt.optimizeSingleGaussian2(sampling[0], sampling[1], opts, peaks);
65+
//console.log(gauss.X2);
66+
gauss = gauss.p;
67+
//console.log("Before");
68+
//console.log(peakList[i]);
69+
result.push({x:gauss[0],y:gauss[2],width:gauss[1]*2.35482}); // From https://en.wikipedia.org/wiki/Gaussian_function#Properties}
70+
//console.log("After");
71+
//console.log(result);
72+
}
73+
}
74+
75+
}
76+
return result;
77+
}
78+
79+
function groupPeaks(peakList,nL){
80+
var group = [];
81+
var groups = [];
82+
var i, j;
83+
var limits = [peakList[0].x,nL*peakList[0].width];
84+
var upperLimit, lowerLimit;
85+
//Merge forward
86+
for(i=0;i<peakList.length;i++){
87+
//If the 2 things overlaps
88+
if(Math.abs(peakList[i].x-limits[0])<(nL*peakList[i].width+limits[1])){
89+
//Add the peak to the group
90+
group.push(peakList[i]);
91+
//Update the group limits
92+
upperLimit = limits[0]+limits[1];
93+
if(peakList[i].x+nL*peakList[i].width>upperLimit){
94+
upperLimit = peakList[i].x+nL*peakList[i].width;
95+
}
96+
lowerLimit = limits[0]-limits[1];
97+
if(peakList[i].x-nL*peakList[i].width<lowerLimit){
98+
lowerLimit = peakList[i].x-nL*peakList[i].width;
99+
}
100+
limits = [(upperLimit+lowerLimit)/2,Math.abs(upperLimit-lowerLimit)/2];
101+
102+
}
103+
else{
104+
groups.push({limits:limits,group:group});
105+
//var optmimalPeak = fitSpectrum(group,limits,spectrum);
106+
group=[peakList[i]];
107+
limits = [peakList[i].x,nL*peakList[i].width];
108+
}
109+
}
110+
groups.push({limits:limits,group:group});
111+
//Merge backward
112+
for(i =groups.length-2;i>=0;i--){
113+
//The groups overlaps
114+
if(Math.abs(groups[i].limits[0]-groups[i+1].limits[0])<
115+
(groups[i].limits[1]+groups[i+1].limits[1])/2){
116+
for(j=0;j<groups[i+1].group.length;j++){
117+
groups[i].group.push(groups[i+1].group[j]);
118+
}
119+
upperLimit = groups[i].limits[0]+groups[i].limits[1];
120+
if(groups[i+1].limits[0]+groups[i+1].limits[1]>upperLimit){
121+
upperLimit = groups[i+1].limits[0]+groups[i+1].limits[1];
122+
}
123+
lowerLimit = groups[i].limits[0]-groups[i].limits[1];
124+
if(groups[i+1].limits[0]-groups[i+1].limits[1]<lowerLimit){
125+
lowerLimit = groups[i+1].limits[0]-groups[i+1].limits[1];
126+
}
127+
//console.log(limits);
128+
groups[i].limits = [(upperLimit+lowerLimit)/2,Math.abs(upperLimit-lowerLimit)/2];
129+
130+
groups.splice(i+1,1);
131+
}
132+
}
133+
return groups;
134+
}
135+
136+
137+
138+
module.exports=optimizePeaks;
139+

‎test/ubiquitin.js‎

Lines changed: 37 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -3,33 +3,61 @@
33
var fs = require('fs');
44

55
var gsd = require("../src/gsdLight");
6+
var optimizePeaks = require("../src/optimize");
67
var parser = require('xy-parser');
78
var Stat = require('ml-stat');
89

910

1011

1112

12-
describe('Global spectra deconvolution2', function () {
13+
//describe('Global spectra deconvolution2', function () {
1314

1415

1516
//var spectrum=parser.parse(fs.readFileSync('./test/ubiquitin.txt', 'utf-8'), {arrayType: 'xxyy'});
16-
var spectrum=parser.parse(fs.readFileSync('./test/ubiquitin.txt', 'utf-8'), {arrayType: 'xxyy'});
17-
17+
var d = new Date();
18+
var n = d.getTime();
19+
var spectrum=parser.parse(fs.readFileSync('./ubiquitin.txt', 'utf-8'), {arrayType: 'xxyy'});
20+
d = new Date();
21+
console.log("Parsing time: "+(d.getTime()-n));
1822
var noiseLevel=Stat.array.median(spectrum[1].filter(function(a) {return (a>0)}));
1923
//console.log("noiseLevel "+noiseLevel);
20-
var result = gsd(spectrum[0],spectrum[1], {noiseLevel: noiseLevel, minMaxRatio:0, broadRatio:0, functionType:"gaussian"});
21-
24+
/*var x = [];
25+
var y = [];
26+
for(var i=2900;i<2923;i++){
27+
x.push(spectrum[0][i]);
28+
y.push(spectrum[1][i]);
29+
}
30+
spectrum[0]=x;
31+
spectrum[1]=y;
32+
console.log(spectrum);*/
33+
var result = gsd(spectrum[0],spectrum[1], {noiseLevel: 0, minMaxRatio:0, broadRatio:0});
34+
//console.log(result);
35+
d = new Date();
36+
console.log("Parsing + gsd time: "+(d.getTime()-n));
37+
/*var resultX = [];
2238
for(var i=0;i<result.length;i++){
23-
console.log(result[i].x+" "+result[i].width+" "+result[i].y);
39+
if(Math.abs(result[i].x-239.15)<0.2)
40+
resultX.push(result[i]);
41+
}
42+
result = resultX;
43+
console.log(resultX);*/
44+
result = optimizePeaks(result,spectrum[0],spectrum[1],2);
45+
d = new Date();
46+
console.log("Parsing + gsd + optimization time: "+(d.getTime()-n));
47+
//console.log(result);
48+
49+
for(var i=0;i<1500;i++){
50+
//if(result[i].width<0.05)
51+
console.log(result[i].x+" "+result[i].width);
2452
}
2553
console.log(result.length);
2654

2755
// Test case obtained from Pag 443, Chap 8.
28-
it('Should provide the right result ...', function () {
56+
/* it('Should provide the right result ...', function () {
2957
//var result = gsd(spectrum[0],spectrum[1], 0.001, 0.1);
3058
//console.log(result);
3159
32-
});
33-
});
60+
});*/
61+
//});
3462

3563

0 commit comments

Comments
 (0)