bioem.cpp 21.1 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
#include <fstream>
#include <boost/program_options.hpp>
#include <iostream>
#include <algorithm>
#include <iterator>
#include <stdio.h>
#include <stdlib.h>
#include <string>
#include <cmath>
#include <omp.h>

#include <fftw3.h>
#include <math.h>
#include "cmodules/timer.h"

#include "param.h"
#include "bioem.h"
#include "model.h"
#include "map.h"


#include "bioem_algorithm.h"


using namespace boost;
namespace po = boost::program_options;

using namespace std;

// A helper function of Boost
template<class T>
ostream& operator<<(ostream& os, const vector<T>& v)
{
    copy(v.begin(), v.end(), ostream_iterator<T>(os, " "));
    return os;
}

bioem::bioem()
{
	
}

bioem::~bioem()
{
	
}

int bioem::configure(int ac, char* av[])
{
    /**************************************************************************************/
    /**** Configuration Routine using boost for extracting parameters, models and maps ****/
    /**************************************************************************************/
    /****** And Precalculating necessary grids, map crosscorrelations and kernels  ********/
    /*************************************************************************************/

    /*** Inizialzing default variables ***/
    std::string infile,modelfile,mapfile;
    Model.readPDB=false;
    param.writeAngles=false;
	RefMap.dumpMap = false;
	RefMap.loadMap = false;

    /*************************************************************************************/
    cout << " ++++++++++++ FROM COMMAND LINE +++++++++++\n\n";
    /*************************************************************************************/

    /********************* Command line reading input with BOOST ************************/

    try {
        po::options_description desc("Command line inputs");
        desc.add_options()
        ("Inputfile", po::value<std::string>(), "Name of input parameter file")
        ("Modelfile", po::value< std::string>() , "Name of model file")
        ("Particlesfile", po::value< std::string>(), "Name of paricles file")
        ("ReadPDB", "(Optional) If reading model file in PDB format")
		("DumpMaps", "(Optional) Dump maps after they were red from maps file")
		("LoadMapDump", "(Optional) Read Maps from dump instead of maps file")
        ("help", "(Optional) Produce help message")
        ;

        po::positional_options_description p;
        p.add("Inputfile", -1);
        p.add("Modelfile", -1);
        p.add("Particlesfile", -1);
        p.add("ReadPDB", -1);
		p.add("DumpMaps", -1);
		p.add("LoadMapDump", -1);

        po::variables_map vm;
        po::store(po::command_line_parser(ac, av).
                  options(desc).positional(p).run(), vm);
        po::notify(vm);

        if((ac < 6)) {
            std::cout << desc << std::endl;
            return 0;
        }
        if (vm.count("help")) {
            cout << "Usage: options_description [options]\n";
            cout << desc;
            return 0;
        }

        if (vm.count("Inputfile"))
        {
            cout << "Input file is: ";
            cout << vm["Inputfile"].as< std::string >()<< "\n";
            infile=vm["Inputfile"].as< std::string >();
        }
        if (vm.count("Modelfile"))
        {
            cout << "Model file is: "
                 << vm["Modelfile"].as<  std::string  >() << "\n";
            modelfile=vm["Modelfile"].as<  std::string  >();
        }

        if (vm.count("ReadPDB"))
        {
            cout << "Reading model file in PDB format.\n";
            Model.readPDB=true;
        }
        
        if (vm.count("DumpMaps"))
        {
            cout << "Dumping Maps after reading from file.\n";
            RefMap.dumpMap = true;
        }
        
        if (vm.count("LoadMapDump"))
        {
            cout << "Loading Map dump.\n";
            RefMap.loadMap = true;
        }        

        if (vm.count("Particlesfile"))
        {
            cout << "Paricle file is: "
                 << vm["Particlesfile"].as< std::string >() << "\n";
            mapfile=vm["Particlesfile"].as< std::string >();
        }
    }
    catch(std::exception& e)
    {
        cout << e.what() << "\n";
        return 1;
    }

    /********************* Reading Parameter Input ***************************/
    // copying inputfile to param class
    param.fileinput = infile.c_str();
    param.readParameters();

    /********************* Reading Model Input ******************************/
    // copying modelfile to model class
    Model.filemodel = modelfile.c_str();
    Model.readModel();

    /********************* Reading Particle Maps Input **********************/
    /********* HERE: PROBLEM if maps dont fit on the memory!! ***************/
    // copying mapfile to ref map class
    RefMap.filemap = mapfile.c_str();
    RefMap.readRefMaps();

    /****************** Precalculating Necessary Stuff *********************/
    precalculate();
	
	param.nTotGridAngles = 10;
	param.nTotCTFs = 10;
	//param.param_device.maxDisplaceCenter = 0;
	
	deviceInit();

    return(0);
}

int bioem::precalculate()
{
    /**************************************************************************************/
    /* Precalculating Routine of Orientation grids, Map crosscorrelations and CTF Kernels */
    /**************************************************************************************/

    // Generating Grids of orientations
    param.CalculateGridsParam();

    //Inizialzing crosscorrelations of Maps
    memset(RefMap.sum_RefMap, 0, BIOEM_MAX_MAPS * sizeof(*RefMap.sum_RefMap));
    memset(RefMap.sumsquare_RefMap, 0, BIOEM_MAX_MAPS * sizeof(*RefMap.sum_RefMap));

    myfloat_t sum,sumsquare;

    //Precalculating cross-correlations of maps
    for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap ; iRefMap++)
    {
        calcross_cor(RefMap.Ref[iRefMap],sum,sumsquare);
        //Storing Crosscorrelations in Map class
        RefMap.sum_RefMap[iRefMap]=sum;
        RefMap.sumsquare_RefMap[iRefMap]=sumsquare;
    }

    // Precalculating CTF Kernels stored in class Param
    param.CalculateRefCTF();

    return(0);
}


int bioem::run()
{
    /**************************************************************************************/
    /**** Main BioEM routine, projects, convolutes and compares with Map using OpenMP ****/
    /**************************************************************************************/

    /**** If we want to control the number of threads -> omp_set_num_threads(XX); ******/
    /****************** Declarying class of Probability Pointer  *************************/
    pProb = new bioem_Probability[RefMap.ntotRefMap];

    printf("\tInitializing\n");
    // Inizialzing Probabilites to zero and constant to -Infinity
    for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
    {
        pProb[iRefMap].Total=0.0;
        pProb[iRefMap].Constoadd=-9999999;
        pProb[iRefMap].max_prob=-9999999;
        for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient ++)
		{
            pProb[iRefMap].forAngles[iOrient]=0.0;
            pProb[iRefMap].ConstAngle[iOrient]=-99999999;
        }
    }
    /**************************************************************************************/
	deviceStartRun();

    /******************************** MAIN CYCLE ******************************************/
	
    /*** Declaring Private variables for each thread *****/
    mycomplex_t* proj_mapFFT;
    bioem_map conv_map;

	//allocating fftw_complex vector
    proj_mapFFT= (mycomplex_t *) fftw_malloc(sizeof(mycomplex_t) *4*param.param_device.NumberPixels*param.param_device.NumberPixels);

	HighResTimer timer;

	printf("\tMain Loop (GridAngles %d, CTFs %d, RefMaps %d, Shifts (%d/%d)²), Pixels %d²\n", param.nTotGridAngles, param.nTotCTFs, RefMap.ntotRefMap, 2 * param.param_device.maxDisplaceCenter + param.param_device.GridSpaceCenter, param.param_device.GridSpaceCenter, param.param_device.NumberPixels);
	printf("\tInner Loop Count (%d %d %d) %lld\n", param.param_device.maxDisplaceCenter, param.param_device.GridSpaceCenter, param.param_device.NumberPixels, (long long int) (param.param_device.NumberPixels * param.param_device.NumberPixels * (2 * param.param_device.maxDisplaceCenter / param.param_device.GridSpaceCenter + 1) * (2 * param.param_device.maxDisplaceCenter / param.param_device.GridSpaceCenter + 1)));
    //#pragma omp parallel for
    for (int iProjectionOut = 0; iProjectionOut < param.nTotGridAngles; iProjectionOut++)
    {
        /***************************************************************************************/
        /***** Creating Projection for given orientation and transforming to Fourier space *****/
		timer.ResetStart();
        createProjection(iProjectionOut, proj_mapFFT);
		printf("Time Projection %d: %f\n", iProjectionOut, timer.GetCurrentElapsedTime());

        /***************************************************************************************/
        /***** **** Internal Loop over convolutions **** *****/
        for (int iConv = 0; iConv < param.nTotCTFs; iConv++)
        {
			printf("\t\tConvolution %d %d\n", iProjectionOut, iConv);
            /*** Calculating convolutions of projection map and crosscorrelations ***/
			timer.ResetStart();
            createConvolutedProjectionMap(iProjectionOut,iConv,proj_mapFFT,conv_map);
			printf("Time Convolution %d %d: %f\n", iProjectionOut, iConv, timer.GetCurrentElapsedTime());

            /***************************************************************************************/
            /*** Comparing each calculated convoluted map with all experimental maps ***/
			timer.ResetStart();
			compareRefMaps(iProjectionOut, iConv, conv_map);
			const double compTime = timer.GetCurrentElapsedTime();
			const int nShifts = 2 * param.param_device.maxDisplaceCenter / param.param_device.GridSpaceCenter + 1;
			const double nFlops = (double) RefMap.ntotRefMap * (double) nShifts * (double) nShifts *
				(((double) param.param_device.NumberPixels - (double) param.param_device.maxDisplaceCenter / 2.) * ((double) param.param_device.NumberPixels - (double) param.param_device.maxDisplaceCenter / 2.) * 5. + 25.) / compTime;
			const double nGBs = (double) RefMap.ntotRefMap * (double) nShifts * (double) nShifts *
				(((double) param.param_device.NumberPixels - (double) param.param_device.maxDisplaceCenter / 2.) * ((double) param.param_device.NumberPixels - (double) param.param_device.maxDisplaceCenter / 2.) * 2. + 8.) * (double) sizeof(myfloat_t) / compTime;
			const double nGBs2 = (double) RefMap.ntotRefMap * ((double) param.param_device.NumberPixels * (double) param.param_device.NumberPixels + 8.) * (double) sizeof(myfloat_t) / compTime;

			printf("Time Comparison %d %d: %f sec (%f GFlops, %f GB/s (cached), %f GB/s)\n", iProjectionOut, iConv, compTime, nFlops / 1000000000., nGBs / 1000000000., nGBs2 / 1000000000.);
        }

    }
    //deallocating fftw_complex vector
    fftw_free(proj_mapFFT);
	
	deviceFinishRun();

    /************* Writing Out Probabilities ***************/

    /*** Angular Probability ***/

    // if(param.writeAngles){
    ofstream angProbfile;
    angProbfile.open ("ANG_PROB_iRefMap");
    // }

    ofstream outputProbFile;
    outputProbFile.open ("Output_Probabilities");

    for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
    {
        /**** Total Probability ***/
        outputProbFile << "RefMap " << iRefMap << " Probability  "  << log(pProb[iRefMap].Total)+pProb[iRefMap].Constoadd+0.5*log(M_PI)+(1-param.param_device.Ntotpi*0.5)*(log(2*M_PI)+1)+log(param.param_device.volu) << " Constant " << pProb[iRefMap].Constoadd  << "\n";

        outputProbFile << "RefMap " << iRefMap << " Maximizing Param: ";

        /*** Param that maximize probability****/
        outputProbFile << (pProb[iRefMap].max_prob + 0.5 * log(M_PI) + (1 - param.param_device.Ntotpi * 0.5) * (log(2 * M_PI) + 1) + log(param.param_device.volu)) << " ";
        outputProbFile << param.angles[pProb[iRefMap].max_prob_orient].pos[0] << " ";
        outputProbFile << param.angles[pProb[iRefMap].max_prob_orient].pos[1] << " ";
        outputProbFile << param.angles[pProb[iRefMap].max_prob_orient].pos[2] << " ";
        outputProbFile << param.CtfParam[pProb[iRefMap].max_prob_conv].pos[0] << " ";
        outputProbFile << param.CtfParam[pProb[iRefMap].max_prob_conv].pos[1] << " ";
        outputProbFile << param.CtfParam[pProb[iRefMap].max_prob_conv].pos[2] << " ";
		outputProbFile << pProb[iRefMap].max_prob_cent_x << " ";
		outputProbFile << pProb[iRefMap].max_prob_cent_y;
        outputProbFile << "\n";

        /*** For individual files***/ //angProbfile.open ("ANG_PROB_"iRefMap);

        if(param.writeAngles)
		{
            for (int iProjectionOut = 0; iProjectionOut < param.nTotGridAngles; iProjectionOut++)
            {
                angProbfile << " " << iRefMap << " " << param.angles[iProjectionOut].pos[0] << " " << param.angles[iProjectionOut].pos[1] << " " << param.angles[iProjectionOut].pos[2] << " " << log(pProb[iRefMap].forAngles[iProjectionOut])+pProb[iRefMap].ConstAngle[iProjectionOut]+0.5*log(M_PI)+(1-param.param_device.Ntotpi*0.5)*(log(2*M_PI)+1)+log(param.param_device.volu) << " " << log(param.param_device.volu) << "\n";

            }
        }

    }

    angProbfile.close();
    outputProbFile.close();

    //Deleting allocated pointers

    if (pProb)
    {
        delete[] pProb;
        pProb = NULL;
    }

    if (param.refCTF)
    {
        delete[] param.refCTF;
        param.refCTF =NULL;
    }

    return(0);
}

350
int bioem::compareRefMaps(int iProjectionOut, int iConv, const bioem_map& conv_map, const int startMap)
351
352
{
#pragma omp parallel for
353
    for (int iRefMap = startMap; iRefMap < RefMap.ntotRefMap; iRefMap ++)
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
    {
        compareRefMapShifted<-1>(iRefMap,iProjectionOut,iConv,conv_map, pProb, param.param_device, RefMap);
    }
    return(0);
}

int bioem::createProjection(int iMap,mycomplex_t* mapFFT)
{

    /**************************************************************************************/
    /****  BioEM Create Projection routine in Euler angle predefined grid****************
     ********************* and turns projection into Fourier space **********************/
    /**************************************************************************************/

    myfloat3_t RotatedPointsModel[Model.nPointsModel];
    myfloat_t rotmat[3][3];
    myfloat_t alpha, gam,beta;
    fftw_plan plan;
    mycomplex_t* localproj;
    int totnumPixFFT=2*param.param_device.NumberPixels;

    localproj= (mycomplex_t *) fftw_malloc(sizeof(mycomplex_t) *4*param.param_device.NumberPixels*param.param_device.NumberPixels);
    memset(localproj,0,4*param.param_device.NumberPixels*param.param_device.NumberPixels*sizeof(*localproj));

    alpha=param.angles[iMap].pos[0];
    beta=param.angles[iMap].pos[1];
    gam=param.angles[iMap].pos[2];

    /**** To see how things are going: cout << "Id " << omp_get_thread_num() <<  " Angs: " << alpha << " " << beta << " " << gam << "\n"; ***/

    /********** Creat Rotation with pre-defiend grid of orientations**********/

    rotmat[0][0]=cos(gam)*cos(alpha)-cos(beta)*sin(alpha)*sin(gam);
    rotmat[0][1]=cos(gam)*sin(alpha)+cos(beta)*cos(alpha)*sin(gam);
    rotmat[0][2]=sin(gam)*sin(beta);
    rotmat[1][0]=-sin(gam)*cos(alpha)-cos(beta)*sin(alpha)*cos(gam);
    rotmat[1][1]=-sin(gam)*sin(alpha)+cos(beta)*cos(alpha)*cos(gam);
    rotmat[1][2]=cos(gam)*sin(beta);
    rotmat[2][0]=sin(beta)*sin(alpha);
    rotmat[2][1]=-sin(beta)*cos(alpha);
    rotmat[2][2]=cos(beta);


    for(int n=0; n< Model.nPointsModel; n++)
	{
        RotatedPointsModel[n].pos[0]=0.0;
        RotatedPointsModel[n].pos[1]=0.0;
        RotatedPointsModel[n].pos[2]=0.0;
    }
    for(int n=0; n< Model.nPointsModel; n++)
    {
        for(int k=0; k< 3; k++)
        {
            for(int j=0; j< 3; j++)
            {
                RotatedPointsModel[n].pos[k]+=rotmat[k][j]*Model.PointsModel[n].pos[j];
            }
        }
    }

    int i, j;

    /************ Projection over the Z axis********************/
    for(int n=0; n< Model.nPointsModel; n++)
    {
        //Getting pixel that represents coordinates & shifting the start at to Numpix/2,Numpix/2 )
        i=floor(RotatedPointsModel[n].pos[0]/(*param.pixelSize)+float(param.param_device.NumberPixels)/2.0+0.5);
        j=floor(RotatedPointsModel[n].pos[1]/(*param.pixelSize)+float(param.param_device.NumberPixels)/2.0+0.5);

        localproj[i*2*param.param_device.NumberPixels+j+param.param_device.NumberPixels*param.param_device.NumberPixels+int(param.param_device.NumberPixels/2.0)][0]+=Model.densityPointsModel[n];

    }

    /**** Output Just to check****/
    if(iMap==10)
    {
        ofstream myexamplemap;
        ofstream myexampleRot;
        myexamplemap.open ("MAP_i10");
        myexampleRot.open ("Rot_i10");
        myexamplemap << "ANGLES " << alpha << " " << beta << " " << gam << "\n";
        for(int k=0; k<2*param.param_device.NumberPixels; k++)
        {
            for(int j=0; j<2*param.param_device.NumberPixels; j++) myexamplemap << "\nMAP " << k << " " << j<< " " <<localproj[k*2*param.param_device.NumberPixels+j][0];
        }
        myexamplemap << " \n";
        for(int n=0; n< Model.nPointsModel; n++)myexampleRot << "\nCOOR " << RotatedPointsModel[n].pos[0] << " " << RotatedPointsModel[n].pos[1] << " " << RotatedPointsModel[n].pos[2];
        myexamplemap.close();
        myexampleRot.close();
    }

    /***** Converting projection to Fourier Space for Convolution later with kernel****/
    /********** Omp Critical is necessary with FFTW*******/
    //#pragma omp critical
    {
        plan = fftw_plan_dft_2d(totnumPixFFT,totnumPixFFT,localproj,mapFFT,FFTW_FORWARD,FFTW_ESTIMATE);
        fftw_execute(plan);
        fftw_destroy_plan(plan);
        fftw_free(localproj);
    }

    return(0);
}


int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,bioem_map& Mapconv)
{
    /**************************************************************************************/
    /****  BioEM Create Convoluted Projection Map routine, multiplies in Fourier **********
    **************** calculated Projection with convoluted precalculated Kernel**********
    *************** and Backtransforming it to real Space ******************************/
    /**************************************************************************************/

    fftw_plan plan;
    mycomplex_t* localmultFFT;
    mycomplex_t* localconvFFT;
    int totnumPixFFT=2*param.param_device.NumberPixels;
    localmultFFT= (mycomplex_t *) fftw_malloc(sizeof(mycomplex_t)*totnumPixFFT*totnumPixFFT);
    localconvFFT= (mycomplex_t *) fftw_malloc(sizeof(mycomplex_t)*totnumPixFFT*totnumPixFFT);


    /**** Multiplying FFTmap with corresponding kernel ****/

    for(int i=0; i < 2*param.param_device.NumberPixels ; i++ )
    {
        for(int j=0; j < 2*param.param_device.NumberPixels ; j++ )
        {   //Projection*CONJ(KERNEL)
            localmultFFT[i*2*param.param_device.NumberPixels+j][0]=lproj[i*2*param.param_device.NumberPixels+j][0]*param.refCTF[iConv].cpoints[i*2*param.param_device.NumberPixels+j][0]+lproj[i*2*param.param_device.NumberPixels+j][1]*param.refCTF[iConv].cpoints[i*2*param.param_device.NumberPixels+j][1];
            localmultFFT[i*2*param.param_device.NumberPixels+j][1]=lproj[i*2*param.param_device.NumberPixels+j][1]*param.refCTF[iConv].cpoints[i*2*param.param_device.NumberPixels+j][0]-lproj[i*2*param.param_device.NumberPixels+j][0]*param.refCTF[iConv].cpoints[i*2*param.param_device.NumberPixels+j][1];
            // cout << "GG " << i << " " << j << " " << param.refCTF[iConv].cpoints[i*2*param.param_device.NumberPixels+j][0] << " " <<param.refCTF[iConv].cpoints[i*2*param.param_device.NumberPixels+j][1] <<" " <<lproj[i*2*param.param_device.NumberPixels+j][0] <<" " <<lproj[i*2*param.param_device.NumberPixels+j][1] << "\n";
        }
    }

    /**** Bringing convoluted Map to real Space ****/
    //#pragma omp critical
    {
        plan = fftw_plan_dft_2d(totnumPixFFT,totnumPixFFT,localmultFFT,localconvFFT,FFTW_BACKWARD,FFTW_ESTIMATE);
        fftw_execute(plan);
    }


    /****Asigning convolution fftw_complex to bioem_map ****/
    for(int i=0; i < param.param_device.NumberPixels ; i++ )
    {
        for(int j=0; j < param.param_device.NumberPixels ; j++ )
        {
            Mapconv.points[i][j]=localconvFFT[i*2*param.param_device.NumberPixels+j+param.param_device.NumberPixels*param.param_device.NumberPixels+int(param.param_device.NumberPixels/2.0)][0];
        }
    }

    /**** Freeing fftw_complex created (dont know if omp critical is necessary) ****/
    //#pragma omp critical
    {
        fftw_destroy_plan(plan);
        fftw_free(localconvFFT);
        fftw_free(localmultFFT);
    }

    return(0);
}

int bioem::calcross_cor(bioem_map& localmap,myfloat_t& sum,myfloat_t& sumsquare)
{
    /*********************** Routine to calculate Cross correlations***********************/

    sum=0.0;
    sumsquare=0.0;
    for (int i = 0; i < param.param_device.NumberPixels; i++)
    {
        for (int j = 0; j < param.param_device.NumberPixels; j++)
        {
            // Calculate Sum of pixels
            sum += localmap.points[i][j];
            // Calculate Sum of pixels squared
            sumsquare += localmap.points[i][j]*localmap.points[i][j];
        }
    }
    return(0);
}

int bioem::deviceInit()
{
	return(0);
}

int bioem::deviceStartRun()
{
	return(0);
}

int bioem::deviceFinishRun()
{
	return(0);
}