1 /* IsotopeDistributionCalculator.c */
2 /**********************************************************************************************************
3 Copyright (c) 2002-2013 Abdul-Rahman Allouche. All rights reserved
4
5 Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated
6 documentation files (the Gabedit), to deal in the Software without restriction, including without limitation
7 the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software,
8 and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
9
10 The above copyright notice and this permission notice shall be included in all copies or substantial portions
11 of the Software.
12
13 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED
14 TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
15 THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
16 CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
17 DEALINGS IN THE SOFTWARE.
18 ************************************************************************************************************/
19
20 #include "../../Config.h"
21 #include <stdlib.h>
22 #include <math.h>
23 #include <ctype.h>
24 #include "../Common/Global.h"
25 #include "../Utils/Constants.h"
26 #include "../Utils/AtomsProp.h"
27 #include "../Utils/Utils.h"
28 #include "../Utils/Constants.h"
29 #include "../Utils/UtilsInterface.h"
30 #include "../Files/FileChooser.h"
31 #include "../Files/FolderChooser.h"
32 #include "../Files/GabeditFolderChooser.h"
33 #include "../Common/Help.h"
34 #include "../Common/StockIcons.h"
35 #include "../Utils/GabeditTextEdit.h"
36 #include "../IsotopeDistribution/IsotopeDistributionCalculator.h"
37
38 #define DEBUGFLAG 0
39
40 /************************************************************************************************************/
41 static gint cmp_2_isodata(gconstpointer a, gconstpointer b);
42 static void free_element_data(ElementData* e);
43 static ElementData get_element_data(gchar* symbol);
44 static void free_elements(ElementData* elements, gint nElements);
45 static ElementData* get_elements(gint nElements, gint* nAtoms, gchar** symbols);
46 /* static void print_elements(ElementData* elements, gint nElements);*/
47 static void cut_peaks(GList* peaks, gdouble abundPrecision);
48 static gint summarize_peaks(GList* peaks, gdouble precision);
49 static GList *add_peak(GList* oldPeaks, IsotopeData* newPeak);
50 /* static GList *add_peak_zero(GList* oldPeaks);*/
51 static IsotopeData* new_iso(gdouble mass, gdouble abundance);
52 static IsotopeData* free_iso(IsotopeData* is);
53 static GList* compute_peaks(gint nElements, ElementData* elements, gdouble massPrecision, gdouble abundancePrecision, gchar** error);
54 static gchar* parse_formula(gchar* formula, gchar*** symbolsP, gint* nElementsP, gint** nAtomsP);
55 /************************************************************************************************************/
cmp_2_isodata(gconstpointer a,gconstpointer b)56 static gint cmp_2_isodata(gconstpointer a, gconstpointer b)
57 {
58 if(a && b && ((IsotopeData* )a)->mass > ((IsotopeData* )b)->mass) return 1;
59 else if(a && b && ((IsotopeData* )a)->mass == ((IsotopeData* )b)->mass) return 0;
60 else return -1;
61 }
62 /************************************************************************************************************/
free_element_data(ElementData * e)63 static void free_element_data(ElementData* e)
64 {
65 if(e && e->isotopes) g_free(e->isotopes);
66 }
67 /************************************************************************************************************/
get_element_data(gchar * symbol)68 static ElementData get_element_data(gchar* symbol)
69 {
70 SAtomsProp prop = prop_atom_get(symbol);
71 ElementData e;
72 gint i;
73
74 e.nIsotopes = prop.nIsotopes;
75 e.isotopes = NULL;
76 if( e.nIsotopes>0) e.isotopes = g_malloc(sizeof(IsotopeData)*e.nIsotopes);
77 for(i=0;i<e.nIsotopes;i++)
78 {
79 e.isotopes[i].mass = prop.rMass[i];
80 e.isotopes[i].abundance = prop.abundances[i]/100.0;
81 }
82 prop_atom_free(&prop);
83 e.nAtoms = 1;
84 return e;
85 }
86 /************************************************************************************************************/
free_elements(ElementData * elements,gint nElements)87 static void free_elements(ElementData* elements, gint nElements)
88 {
89 if(elements)
90 {
91 gint i;
92 for(i=0;i<nElements;i++)
93 free_element_data(&elements[i]);
94 g_free(elements);
95 }
96 }
97 /************************************************************************************************************/
get_elements(gint nElements,gint * nAtoms,gchar ** symbols)98 static ElementData* get_elements(gint nElements, gint* nAtoms, gchar** symbols)
99 {
100 gint i;
101 ElementData* elements = NULL;
102 if(nElements<1) return elements;
103 if(!nAtoms) return elements;
104 if(!symbols) return elements;
105 elements = g_malloc(nElements*sizeof(ElementData));
106 for(i=0;i<nElements;i++)
107 {
108 elements[i] = get_element_data(symbols[i]);
109 elements[i].nAtoms = nAtoms[i];
110 }
111 return elements;
112 }
113 /************************************************************************************************************/
get_str_elements(ElementData * elements,gint nElements)114 static gchar* get_str_elements(ElementData* elements, gint nElements)
115 {
116 gint i, j;
117 gchar* str = NULL;
118 gchar* dum = NULL;
119 str = g_strdup_printf("List of elements of the molecule\n");
120 dum = str;
121 str = g_strdup_printf("%s================================\n",str);
122 g_free(dum);
123 for(i=0;i<nElements;i++)
124 {
125 dum = str;
126 str = g_strdup_printf("%snAtoms = %d\n", str,elements[i].nAtoms);
127 g_free(dum);
128 for(j=0;j<elements[i].nIsotopes;j++)
129 {
130 dum = str;
131 str = g_strdup_printf("%s M = %f A = %f\n", str, elements[i].isotopes[j].mass, elements[i].isotopes[j].abundance);
132 g_free(dum);
133 }
134 }
135 dum = str;
136 str = g_strdup_printf("%s\n",str);
137 g_free(dum);
138 return str;
139 }
140 /************************************************************************************************************/
141 /*
142 static void print_elements(ElementData* elements, gint nElements)
143 {
144 gint i, j;
145 printf("List of elements of the molecule\n");
146 printf("================================\n");
147 for(i=0;i<nElements;i++)
148 {
149 printf("nAtoms = %d\n", elements[i].nAtoms);
150 for(j=0;j<elements[i].nIsotopes;j++)
151 printf("\tMass = %f Abundance = %f\n", elements[i].isotopes[j].mass, elements[i].isotopes[j].abundance);
152 }
153 printf("\n");
154 }
155 */
156 /************************************************************************************************************/
cut_peaks(GList * peaks,gdouble abundPrecision)157 static void cut_peaks(GList* peaks, gdouble abundPrecision)
158 {
159 /* don't remove the first elements */
160 while(peaks && peaks->next)
161 {
162 IsotopeData* data =((IsotopeData*)(peaks->next->data));
163 if(data && data->abundance<abundPrecision/100)
164 {
165 peaks = g_list_remove(peaks, data);
166 g_free(data);
167 continue;
168 }
169 else peaks=peaks->next;
170 }
171 }
172 /************************************************************************************************************/
summarize_peaks(GList * peaks,gdouble precision)173 static gint summarize_peaks(GList* peaks, gdouble precision)
174 {
175 GList* iter;
176 GList* dum;
177 gint nRemoved = 0;
178 for(iter=peaks; iter != NULL; iter=iter->next)
179 {
180 while( iter->next && fabs( ((IsotopeData*)(iter->next->data))->mass - ((IsotopeData*)(iter->data))->mass) < precision)
181 {
182 dum=iter->next;
183 iter->next=dum->next;
184 if(iter->next)
185 iter->next->prev=iter;
186 ((IsotopeData*)(iter->data))->abundance += ((IsotopeData*)(dum->data))->abundance;
187 free_iso((IsotopeData*)(dum->data));
188 }
189 nRemoved++;
190 }
191 return nRemoved;
192 }
193 /************************************************************************************************************/
add_peak(GList * oldPeaks,IsotopeData * newPeak)194 static GList *add_peak(GList* oldPeaks, IsotopeData* newPeak)
195 {
196
197 if(!oldPeaks) return g_list_append(oldPeaks, newPeak);
198 if(oldPeaks->data && ((IsotopeData*)(oldPeaks->data))->mass<1e-12)
199 {
200 free_isotope_distribution(oldPeaks);
201 oldPeaks = NULL;
202 return g_list_append(oldPeaks, newPeak);
203
204 }
205 return g_list_insert_sorted(oldPeaks, newPeak, cmp_2_isodata);
206 }
207 /************************************************************************************************************/
208 /*
209 static GList *add_peak_zero(GList* oldPeaks)
210 {
211 IsotopeData* is = new_iso(0,0);
212 if(!is) return oldPeaks;
213 return add_peak(oldPeaks, is);
214 }
215 */
216 /************************************************************************************************************/
new_iso(gdouble mass,gdouble abundance)217 static IsotopeData* new_iso(gdouble mass, gdouble abundance)
218 {
219 IsotopeData* is = g_malloc(sizeof(IsotopeData));
220 is->mass = mass;
221 is->abundance = abundance;
222 return is;
223 }
224 /************************************************************************************************************/
free_iso(IsotopeData * is)225 static IsotopeData* free_iso(IsotopeData* is)
226 {
227 if(is) g_free(is);
228 return NULL;
229 }
230 /************************************************************************************************************/
free_isotope_distribution(GList * isotopeDistribution)231 GList* free_isotope_distribution(GList* isotopeDistribution)
232 {
233 GList* peaks = isotopeDistribution;
234 while(peaks)
235 {
236 free_iso((IsotopeData*)peaks->data);
237 peaks = peaks->next;
238 }
239 g_list_free(peaks);
240 return NULL;
241 }
242 /************************************************************************************************************/
compute_peaks(gint nElements,ElementData * elements,gdouble massPrecision,gdouble abundancePrecision,gchar ** error)243 static GList* compute_peaks(gint nElements, ElementData* elements, gdouble massPrecision, gdouble abundancePrecision, gchar** error)
244 {
245 GList* peaks = NULL;
246 GList* newPeaks = NULL;
247 GList* iIter = NULL;
248 IsotopeData* is;
249 gint i;
250 gint j;
251 gint k;
252 gint nRemoved;
253 gint updateCounter = 0;
254 gint frequenceUpdate = 0;
255 if(*error) *error = NULL;
256 if(nElements<1 || elements[0].nAtoms<1 || elements[0].nIsotopes<1) return peaks;
257 is = new_iso(0,1);
258 peaks = add_peak(NULL, is);
259 if(!peaks)
260 {
261 if(*error) *error = g_strdup(_("No enough memory"));
262 return peaks;
263 }
264 frequenceUpdate = 100*nElements;
265 cancelIsotopeDistribution = FALSE;
266 for(i = 0; i<nElements; i++)
267 {
268 for(j = 0; j<elements[i].nAtoms; j++)
269 {
270 is = new_iso(0,1);
271 newPeaks = add_peak(NULL, is);
272 if(!is)
273 {
274 free_isotope_distribution(peaks);
275 if(*error) *error = g_strdup(_("No enough memory"));
276 return NULL;
277 }
278 for(iIter=peaks; iIter != NULL; iIter = iIter->next)
279 {
280 for(k=0;k<elements[i].nIsotopes;k++)
281 {
282 is = new_iso(0,0);
283 if(!is)
284 {
285 free_isotope_distribution(peaks);
286 if(*error) *error = g_strdup(_("No enough memory"));
287 return NULL;
288 }
289 is->mass = ((IsotopeData*)(iIter->data))->mass + elements[i].isotopes[k].mass;
290 is->abundance = ((IsotopeData*)(iIter->data))->abundance * elements[i].isotopes[k].abundance;
291 newPeaks = add_peak(newPeaks, is);
292 if(!newPeaks) { free_isotope_distribution(peaks); return NULL;}
293 if( (updateCounter++) > frequenceUpdate)
294 {
295 updateCounter = 0;
296 while( gtk_events_pending() ) gtk_main_iteration();
297 }
298 if(cancelIsotopeDistribution) break;
299 }
300 if(cancelIsotopeDistribution) break;
301 }
302 if(cancelIsotopeDistribution) break;
303 free_isotope_distribution(peaks);
304 peaks = newPeaks;
305 nRemoved = summarize_peaks(peaks,massPrecision);
306 }
307 if(cancelIsotopeDistribution) break;
308 cut_peaks(peaks,abundancePrecision);
309 }
310 if(cancelIsotopeDistribution)
311 {
312 free_isotope_distribution(peaks);
313 peaks = NULL;
314 if(*error) *error = g_strdup(_("Calculation canceled"));
315 }
316 return peaks;
317 }
318 /************************************************************************************************************/
compute_isotope_distribution(gint nElements,gint * nAtoms,gchar ** symbols,gdouble massPrecision,gdouble abundancePrecision,gchar ** error,gchar ** info)319 GList* compute_isotope_distribution(gint nElements, gint* nAtoms, gchar** symbols, gdouble massPrecision, gdouble abundancePrecision, gchar** error,gchar** info)
320 {
321 ElementData* elements = NULL;
322 GList* peaks = NULL;
323 if(*error) *error = NULL;
324 if(nElements<1)
325 {
326 if(error) *error = g_strdup(_("Number of elements can not been <1 !"));
327 return peaks;
328 }
329 elements = get_elements(nElements, nAtoms, symbols);
330 /* print_elements(elements, nElements);*/
331 if(info)*info = get_str_elements(elements, nElements);
332
333
334 peaks = compute_peaks(nElements, elements, massPrecision, abundancePrecision, error);
335 free_elements(elements, nElements);
336 return peaks;
337 }
338 /************************************************************************************************************/
parse_formula(gchar * formula,gchar *** symbolsP,gint * nElementsP,gint ** nAtomsP)339 static gchar* parse_formula(gchar* formula, gchar*** symbolsP, gint* nElementsP, gint** nAtomsP)
340 {
341 gchar** symbols = NULL;
342 gint nElements;
343 gint* nAtoms = NULL;
344 gchar* c = NULL;
345 gchar* listSymbols = NULL;
346 gchar* listNumbers = NULL;
347 gint l = 0;
348 gint i,j;
349 gchar** nA;
350 gint ns=0;
351 gint nn=0;
352
353 *symbolsP = NULL;
354 *nElementsP = 0;
355 *nAtomsP = NULL;
356
357 /* printf("formula=%s\n",formula );*/
358 if(!formula) return g_strdup(_("Formula not valid"));
359 if(strlen(formula)<1) return g_strdup(_("Formula not valid"));
360 if(islower(*formula)) return g_strdup(_("The first character of your formula is a lower case !"));
361 if(isdigit(*formula)) return g_strdup(_("The first character of your formula is a digit !"));
362
363 l =strlen(formula);
364 listSymbols = g_malloc((2*l+1)*sizeof(gchar));
365 for(i=0; i<2*l;i++) listSymbols[i] = ' ';
366 listSymbols[2*l-1] = '\0';
367 j = 0;
368 listSymbols[j] = formula[j];
369 for(i=1; i<l;i++)
370 {
371 if(isdigit(formula[i]))
372 {
373 if(listSymbols[j]!=' ') listSymbols[++j] = ' ';
374 continue;
375 }
376 if(!islower(formula[i]))
377 {
378 if(listSymbols[j]!=' ') listSymbols[++j] = ' ';
379 listSymbols[++j] = formula[i];
380 continue;
381 }
382 listSymbols[++j] = formula[i];
383 }
384 /* printf("j=%d\n",j);*/
385 listSymbols[j+1] = '\0';
386 delete_last_spaces(listSymbols);
387 delete_first_spaces(listSymbols);
388 /* printf("ListSymb=%s#\n",listSymbols);*/
389 symbols = g_strsplit(listSymbols," ",-1);
390 g_free(listSymbols);
391 i = 0;
392 for(c=symbols[i];c;c=symbols[++i])
393 {
394 /* printf("c=%s\n",c);*/
395 }
396 ns = i;
397 /* printf("ns=%d\n",ns);*/
398
399 listNumbers = g_malloc((2*l+1)*sizeof(gchar));
400 for(i=0; i<2*l;i++) listNumbers[i] = ' ';
401 listNumbers[2*l] = '\0';
402
403 j = 0;
404 for(i=1; i<l;i++)
405 {
406 if(isdigit(formula[i]))
407 {
408 listNumbers[j++] = formula[i];
409 continue;
410 }
411 if(!islower(formula[i]))
412 {
413 if(!isdigit(formula[i-1])) listNumbers[j++] = '1';
414 listNumbers[j++] = ' ';
415 continue;
416 }
417 }
418 if(!isdigit(formula[l-1])) listNumbers[j++] = '1';
419
420 /* printf("j=%d\n",j);*/
421 listNumbers[j+1] = '\0';
422 delete_last_spaces(listNumbers);
423 delete_first_spaces(listNumbers);
424 /* printf("ListNumb=%s#\n",listNumbers);*/
425 nA = g_strsplit(listNumbers," ",-1);
426 g_free(listNumbers);
427 i = 0;
428 for(c=nA[i];c;c=nA[++i])
429 {
430 /* printf("c=%s\n",c);*/
431 }
432 nn = i;
433 /* printf("ns=%d\n",nn);*/
434 if(ns<1 || ns!=nn)
435 {
436 free_one_string_table(symbols, i);
437 free_one_string_table(nA, j);
438 return g_strdup(_("Your formula is not valid !"));
439 }
440 nElements = ns;
441 for(i=0;i<nElements;i++)
442 {
443 /*printf("symbol=%s#\n",symbols[i]);*/
444 if(!test_atom_define(symbols[i]))
445 {
446 gchar* mess = g_strdup_printf(_("%s is not a known atom for Gabedit !"),symbols[i]);
447 free_one_string_table(symbols, i);
448 free_one_string_table(nA, j);
449 return mess;
450 }
451 }
452 nAtoms = g_malloc(nElements*sizeof(gint));
453 for(i=0;i<nElements;i++)
454 {
455 nAtoms[i] = atoi(nA[i]);
456 if(nAtoms[i]<1) nAtoms[i] = 1;
457 }
458 nA = free_one_string_table(nA, nElements);
459 *symbolsP = symbols;
460 *nElementsP = nElements;
461 *nAtomsP = nAtoms;
462 return NULL;
463 }
464 /************************************************************************************************************/
compute_isotope_distribution_from_formula(gchar * formula,gdouble massPrecision,gdouble abundancePrecision,gchar ** error,gchar ** info)465 GList* compute_isotope_distribution_from_formula(gchar* formula, gdouble massPrecision, gdouble abundancePrecision, gchar** error, gchar** info)
466 {
467 gchar** symbols = NULL;
468 gint nElements;
469 gint* nAtoms = NULL;
470 GList* peaks = NULL;
471 gchar* message;
472
473 if(*error) *error = NULL;
474 message = parse_formula(formula, &symbols, &nElements, &nAtoms);
475 if(message)
476 {
477 if(error) *error = message;
478 else g_free(message);
479 return peaks;
480 }
481 peaks = compute_isotope_distribution(nElements, nAtoms, symbols, massPrecision, abundancePrecision, error, info);
482 return peaks;
483 }
484 /************************************************************************************************************/
get_example_isotope_distribution()485 GList* get_example_isotope_distribution()
486 {
487 GList* peaks = NULL;
488 IsotopeData* is = new_iso(1,1);
489 peaks = add_peak(peaks, is);
490 is = new_iso(2,99);
491 peaks = add_peak(peaks, is);
492 return peaks;
493 }
494 /*****************************************************************************/
get_sum_abundance_from_list(GList * peaks)495 gdouble get_sum_abundance_from_list(GList* peaks)
496 {
497 GList* p = peaks;
498 gdouble s = 0;
499
500 s = 0;
501 for(p = peaks; p != NULL; p = p->next)
502 {
503 IsotopeData* is = (IsotopeData*) (p->data);
504 if(is) s += is->abundance;
505 }
506 return s;
507 }
508 /*****************************************************************************/
get_max_abundance_from_list(GList * peaks)509 gdouble get_max_abundance_from_list(GList* peaks)
510 {
511 GList* p = peaks;
512 gdouble max = 0;
513
514 max = 0;
515 for(p = peaks; p != NULL; p = p->next)
516 {
517 IsotopeData* is = (IsotopeData*) (p->data);
518 if(is && is->abundance>max) max = is->abundance;
519 }
520 return max;
521 }
522