~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

Linux Cross Reference
Tina6/tina-libs/tina/math/mathQPR_analysis.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

  1 /**********
  2  *
  3  * This file is part of the TINA Open Source Image Analysis Environment
  4  * henceforth known as TINA
  5  *
  6  * TINA is free software; you can redistribute it and/or modify
  7  * it under the terms of the GNU General Public License as
  8  * published by the Free Software Foundation.
  9  *
 10  * TINA is distributed in the hope that it will be useful,
 11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
 12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 13  * GNU General Public License for more details.
 14  *
 15  * You should have received a copy of the GNU General Public License
 16  * along with TINA; if not, write to the Free Software Foundation, Inc.,
 17  * 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 18  *
 19  * ANY users of TINA who require exemption from the existing licence must
 20  * negotiate a new licence with Dr. Neil.A.Thacker, the sole agent for
 21  * the University of Manchester.
 22  *
 23  **********
 24  *
 25  * Program :    TINA
 26  * File    :  $Source: Paul Tar PhD work  $
 27  * Date    :  $Date: 4/6/14 $
 28  * Version :  $Revision: $
 29  *
 30  * Author  : Paul D Tar
 31  *
 32  * Notes : 
 33  *
 34  *
 35  *
 36  *********
 37 */
 38 
 39 #if HAVE_CONFIG_H
 40 #include <config.h>
 41 #endif
 42 
 43 #include <stdlib.h>
 44 #include <math.h>
 45 #include <time.h>
 46 #include <limits.h>
 47 
 48 #include <tina/sys/sysDef.h>
 49 #include <tina/sys/sysPro.h>
 50 #include <tina/math/mathDef.h>
 51 #include <tina/math/mathPro.h>
 52 #include <tina/math/math_QPRDef.h>
 53 #include <tina/math/mathQPR_alloc.h>
 54 #include <tina/math/mathQPR_utils.h>
 55 #include <tina/math/mathQPR_sample.h>
 56 #include <tina/math/mathQPR_io.h>
 57 
 58 #include "mathQPR_analysis.h"
 59 
 60 static void set_uniform_quantities(QPRDataBlock *block, QPRModel *model)
 61 {
 62    int j, k;
 63    double quantity;
 64 
 65    quantity = block->total_quantity / model->component_count;
 66    block->unknown = 0;
 67    block->poorly_modelled = 0;
 68    
 69    for (k=0; k<model->component_count; k++)
 70    {
 71       block->component_quantities[k] = quantity;
 72    }
 73 
 74    for (k=0; k<model->class_count; k++)
 75    {
 76       block->class_quantities[k] = 0;
 77       for (j=0; j<model->component_count; j++)
 78       {
 79          if (model->class_map->el[k][j] > 0)
 80          {
 81             block->class_quantities[k] += block->component_quantities[j];
 82          }
 83       }
 84    }
 85 }
 86 
 87 static void compute_pkX(QPRDataBlock *block, QPRModel *model)
 88 {
 89    int X, k;
 90    double mXk, mX;
 91    double mnXk, mnX;
 92 
 93    for (X=0; X<block->bin_count; X++)
 94    {
 95       if (model->components[X] != NULL)
 96       {
 97          /* Allocate data bin if required */
 98          if (block->components[X] == NULL)
 99          {
100             block->components[X] =
101                ralloc(sizeof(QPRDataBin)*model->component_count);
102             if (block->components[X] == NULL)
103             {
104                error("compute_pkX: unable to allocate components\n", non_fatal);
105                return;
106             }
107             for (k=0; k<model->component_count; k++)
108             {
109                block->components[X][k].pkX = 0;
110                block->components[X][k].pknX = 0;
111             }
112          }
113 
114          /* Bayes Theorem */
115          mX = mnX = 0;
116          for (k=0; k<model->component_count; k++)
117          {
118             mXk = model->components[X][k].pXk*block->component_quantities[k];
119             mnXk = (1-model->components[X][k].pXk)*block->component_quantities[k];
120             mX += mXk;
121             mnX += mnXk;
122          }
123          if (mX > 0)
124          {
125             for (k=0; k<model->component_count; k++)
126             {
127                mXk = model->components[X][k].pXk*
128                      block->component_quantities[k];
129                block->components[X][k].pkX = mXk / mX;
130             }
131          }
132          if (mnX > 0)
133          {
134             for (k=0; k<model->component_count; k++)
135             {
136                mnXk = (1-model->components[X][k].pXk)*
137                       block->component_quantities[k];
138                block->components[X][k].pknX = mnXk / mnX;
139             }
140          }
141       }
142       else
143       {
144          /* Pattern bin exists in data but not model */
145          if (block->components[X] != NULL)
146          {
147             for (k=0; k<model->component_count; k++)
148             {
149                block->components[X][k].pkX = 0;
150                block->components[X][k].pknX = 0;
151             }
152          }
153       }
154    }
155 }
156 
157 static void compute_quantities(QPRDataBlock *block, QPRModel *model)
158 {
159    int X, k, l;
160    double mX;
161 
162    for (k=0; k<model->component_count; k++)
163    {
164       block->component_quantities[k] = 0.0;
165    }
166    block->unknown = 0;
167    block->poorly_modelled = 0;
168 
169    /* Component quantities */
170    for (X=0; X<block->bin_count; X++)
171    {
172       if (block->components[X] != NULL)
173       {
174          for (k=0; k<model->component_count; k++)
175          {
176             block->component_quantities[k] +=
177                block->components[X][k].pkX * block->hist[X];
178          }
179       }
180    }
181 
182    /* Class quantities */
183    for (k=0; k<model->class_count; k++)
184    {
185       block->class_quantities[k] = 0.0;
186       for (l=0; l<model->component_count; l++)
187       {
188          if (model->class_map->el[k][l] > 0)
189          {
190             block->class_quantities[k] += block->component_quantities[l];
191          }
192       }
193    }
194 
195    /* Unknown quantities */
196    for (X=0; X<model->bin_count; X++)
197    {
198       mX = 0;
199       
200       if (model->components[X] != NULL)
201       {
202          for (k=0; k<model->component_count; k++)
203          {
204             mX += model->components[X][k].pXk * block->component_quantities[k];
205          }
206       }
207 
208       if (mX == 0)
209       {
210          block->unknown += block->hist[X];
211       }
212       else if (mX < 5)
213       {
214          block->poorly_modelled += block->hist[X];
215       }
216    }
217 }
218 
219 /************************** OLD COVAR CODE, PRE AREA ***************************/
220 
221 #ifdef _QPRDEBUG
222 static void compute_cov_OLD(QPRRegion *region, QPRModel *model, Bool crb, Bool ide,
223                             Bool tde, Bool amp, Bool scale)
224 {
225    int i, j, k, X;
226    double hX, hnX;
227    double pkX, pnkX, piX, pjX;
228    double pknX, pnknX, pinX, pjnX;
229    double pXi, pniX, pXj, pXk, pnXk;
230    double var_pXk, sum;
231    double QiQj, sum_data, sum_model;
232    double deltaQik, deltaQjk;
233 
234    Mat *component_delta_q, *component_amp, *component_amp_trans;
235    Mat *component_cov_times_amp;
236    Mat *class_times_component_cov, *class_trans;
237    
238    mat_zero(region->component_cov_data);
239    mat_zero(region->component_cov_model);
240    mat_zero(region->component_cov_total);
241    mat_zero(region->class_cov);
242    /* Cramer Rao lower variance bound method, i.e. var due to incoming data H(X) */
243    if (crb)
244    {
245       for (i=0; i<model->component_count; i++)
246       {
247          for (j=0; j<=i; j++)
248          {
249             sum_data = 0;
250             QiQj = region->component_quantities[i] * 
251                    region->component_quantities[j];
252             
253             for (X=0; X<region->bin_count; X++)
254             {
255                if (region->hist[X] > 0 && region->components[X] != NULL)
256                {
257                   piX = region->components[X][i].pkX;
258                   pjX = region->components[X][j].pkX;
259                   sum_data += piX*pjX*region->hist[X];
260                }
261             }
262             if (QiQj > 0)
263             {
264                region->component_cov_data->el[i][j] = sum_data / QiQj;
265                if (i != j) region->component_cov_data->el[j][i] = sum_data / QiQj;
266             }
267          }
268       }
269       mat_invert(region->component_cov_data);
270    }
271    /* Error propagation method */
272    if (ide || tde)
273    {
274       for (i=0; i<model->component_count; i++)
275       {
276          for (j=0; j<=i; j++)
277          {
278             sum_data = 0;
279             sum_model = 0;
280             for (X=0; X<region->bin_count; X++)
281             {
282                if (region->hist[X] > 0 &&
283                    region->components[X] != NULL &&
284                    model->components[X] != NULL)
285                {
286                   piX = region->components[X][i].pkX;
287                   pinX = region->components[X][i].pknX;
288                   pjX = region->components[X][j].pkX;
289                   pjnX = region->components[X][j].pknX;
290                   hX = region->hist[X]; hnX = region->total_quantity-hX;
291                   /* Single EM step variance from incoming data, i.e. H(X) */
292                   if (ide)
293                   {
294                      /* Assume Poisson var on H(X), i.e. var(H(X)) = H(X) */
295                      sum_data += piX*pjX*hX;
296                   }
297                   /* Single EM step variance from training data, i.e. P(X|k) */
298                   if (tde) 
299                   {
300                      for (k=0; k<model->component_count; k++)
301                      {
302                         deltaQik = deltaQjk = 0;
303                         pkX = region->components[X][k].pkX;
304                         pnkX = 1.0 - pkX;
305                         pknX = region->components[X][k].pknX;
306                         pnknX = 1.0 - pknX;
307                         pXk = model->components[X][k].pXk;
308                         pnXk = 1.0 - pXk;
309 
310                         if (i==k)
311                         {  
312                            if (pXk > 0) deltaQik = (pkX*pnkX*pnXk*hX)/pXk;
313                            deltaQik = deltaQik-pnknX*pknX*hnX;
314                         }
315                         else
316                         {
317                            deltaQik = pinX*pknX*hnX;
318                            if (pXk > 0) deltaQik = deltaQik-(piX*pkX*pnXk*hX)/pXk;
319                         }
320 
321                         if (j==k)
322                         {
323                            if (pXk > 0) deltaQjk = (pkX*pnkX*pnXk*hX)/pXk;
324                            deltaQjk = deltaQjk-pnknX*pknX*hnX;
325                         }
326                         else
327                         {
328                            deltaQjk = pjnX*pknX*hnX;
329                            if (pXk > 0) deltaQjk = deltaQjk-(pjX*pkX*pnXk*hX)/pXk;
330                         }
331                      
332                         var_pXk = model->components[X][k].var_pXk;
333                         
334                         sum_model += deltaQik*deltaQjk*var_pXk;
335                      }
336                   }
337                }
338             }
339 
340             region->component_cov_data->el[i][j] += sum_data;
341             region->component_cov_model->el[i][j] += sum_model;
342 
343             if (i != j)
344             {
345                region->component_cov_data->el[j][i] += sum_data;
346                region->component_cov_model->el[j][i] += sum_model;
347             }
348          }
349       }
350       /* Scaling of error to model fit */
351       if (scale)
352       {
353          mat_times(region->model_fit, component_cov_de);
354       }
355 
356       /* Amplification of error */
357       if (amp)
358       {
359          component_delta_q = mat_make(model->component_count,
360                                       model->component_count);
361          component_amp = mat_make(model->component_count,
362                                   model->component_count);
363          component_amp_trans = mat_make(model->component_count,
364                                         model->component_count);
365 
366          mat_zero(component_delta_q);
367          mat_zero(component_amp);
368          mat_zero(component_amp_trans);
369 
370          for (i=0; i<model->component_count; i++)
371          {
372             for (j=0; j<model->component_count; j++)
373             {
374                sum = 0;
375                for (X=0; X<region->bin_count; X++)
376                {
377                   if (region->hist[X] > 0 &&
378                       region->components[X] != NULL &&
379                       model->components[X] != NULL)
380                   {
381                      if (i==j)
382                      {
383                         pXi = model->components[X][i].pXk;
384                         pniX = 1-region->components[X][i].pkX;
385                         sum += pXi*pniX;
386                      }
387                      else
388                      {
389                         pXj = model->components[X][j].pXk;
390                         piX = region->components[X][i].pkX;
391                         sum += -(pXj*piX);
392                      }
393                   }
394                }
395                component_delta_q->el[i][j] = sum;
396                if (i==j) component_amp->el[i][j] = 1.0;
397             }
398          }
399 
400          component_cov_times_amp = mat_make(model->component_count,
401                                             model->component_count);
402          mat_zero(component_cov_times_amp);
403 
404          mat_diff(component_amp, component_delta_q);
405          mat_invert(component_amp);
406          mat_transp(component_amp_trans, component_amp);
407 
408          /* data error */
409          mat_prod(component_cov_times_amp, component_amp, region->component_cov_data);
410          mat_prod(region->component_cov_data, component_cov_times_amp, component_amp_trans);
411 
412          /* model error */
413          mat_prod(component_cov_times_amp, component_amp, region->component_cov_model);
414          mat_prod(region->component_cov_model, component_cov_times_amp, component_amp_trans);
415 
416          mat_free(component_delta_q);
417          mat_free(component_amp);
418          mat_free(component_amp_trans);
419          mat_free(component_cov_times_amp);
420       }
421    }
422 
423    /* Combine component cov to determine overall var on class quantities */ 
424    class_times_component_cov = mat_make(model->class_count, model->component_count);
425    class_trans = mat_make(model->component_count, model->class_count);
426 
427    mat_accum(region->component_cov_total, 1, region->component_cov_data);
428    mat_accum(region->component_cov_total, 1, region->component_cov_model);
429    
430    mat_prod(class_times_component_cov, model->class_map, region->component_cov_total);
431    mat_transp(class_trans, model->class_map);
432    mat_prod(region->class_cov, class_times_component_cov, class_trans);
433 
434    mat_free(class_times_component_cov);
435    mat_free(class_trans);
436 }
437 #endif
438 /***************************** END OF OLD COVAR CODE *****************************/
439 
440 static void compute_amp_matrix(QPRDataBlock *block, QPRModel *model,
441                                Mat *amp, Mat *amp_trans)
442 {
443    int i, j, X;
444    double sum, pXi, pXj, piX, pniX;
445    Mat *grad_q;
446 
447    grad_q = mat_make(model->component_count, model->component_count);
448 
449    for (i=0; i<model->component_count; i++)
450    {
451       for (j=0; j<model->component_count; j++)
452       {
453          sum = 0;
454          for (X=0; X<block->bin_count; X++)
455          {
456             if (block->hist[X] > 0 &&
457                 block->components[X] != NULL &&
458                 model->components[X] != NULL)
459             {
460                if (i==j)
461                {
462                   pXi = model->components[X][i].pXk;
463                   pniX = 1-block->components[X][i].pkX;
464                   sum += pXi*pniX;
465                }
466                else
467                {
468                   pXj = model->components[X][j].pXk;
469                   piX = block->components[X][i].pkX;
470                   sum += -(pXj*piX);
471                }
472             }
473          }
474          grad_q->el[i][j] = sum;
475          if (i==j) amp->el[i][j] = 1.0; else amp->el[i][j] = 0;
476       }
477    }
478 
479    mat_diff(amp, grad_q);
480    mat_invert(amp);
481    mat_transp(amp_trans, amp);
482 
483    mat_free(grad_q);
484 }
485 
486 static void compute_quantity_cov(QPRDataBlock *block, QPRModel *model, Mat *amp,
487                                  Mat *amp_trans, Bool scale)
488 {
489    int i, j, k, X;
490    double hX, hnX, pkX, pnkX, piX, pjX, pknX, pnknX, pinX, pjnX;
491    double pXk, pnXk, var_pXk;
492    double covar, deltaQik, deltaQjk;
493    Mat *data_X_cov, *model_X_cov, *cov_times_amp;
494    Mat *class_times_component_cov, *class_map_trans;
495 
496    data_X_cov = mat_make(model->component_count, model->component_count);
497    model_X_cov = mat_make(model->component_count, model->component_count);
498    cov_times_amp = mat_make(model->component_count, model->component_count);
499    class_times_component_cov = mat_make(model->class_count, model->component_count);
500    class_map_trans = mat_make(model->component_count, model->class_count);
501 
502    mat_zero(block->component_cov_data);
503    mat_zero(block->component_cov_model);
504    mat_zero(block->component_cov_total);
505    mat_zero(block->class_cov);
506 
507    for (X=0; X<block->bin_count; X++)
508    {
509       if (block->hist[X] > 0 && block->components[X] != NULL &&
510           model->components[X] != NULL)
511       {
512          mat_zero(data_X_cov);
513          mat_zero(model_X_cov);
514          hX = block->hist[X];
515          hnX = block->total_quantity-hX;
516 
517          for (i=0; i<model->component_count; i++)
518          {
519             piX = block->components[X][i].pkX;
520             pinX = block->components[X][i].pknX;
521                
522             for (j=0; j<=i; j++)
523             {
524                pjX = block->components[X][j].pkX;
525                pjnX = block->components[X][j].pknX;
526 
527                /* Single EM step covariance from incoming data, i.e. H(X) */
528                covar = piX*pjX*hX;
529                data_X_cov->el[i][j] = covar;
530                if (i != j) data_X_cov->el[j][i] = covar;
531 
532                /* Single EM step covariance from training data, i.e. P(X|k) */
533                covar = 0;
534                for (k=0; k<model->component_count; k++)
535                {
536                   deltaQik = deltaQjk = 0;
537                   pkX = block->components[X][k].pkX;
538                   pnkX = 1.0 - pkX;
539                   pknX = block->components[X][k].pknX;
540                   pnknX = 1.0 - pknX;
541                   pXk = model->components[X][k].pXk;
542                   pnXk = 1.0 - pXk;
543 
544                   if (i==k)
545                   {  
546                      if (pXk > 0) deltaQik = (pkX*pnkX*pnXk*hX)/pXk;
547                      deltaQik = deltaQik-pnknX*pknX*hnX;
548                   }
549                   else
550                   {
551                      deltaQik = pinX*pknX*hnX;
552                      if (pXk > 0) deltaQik = deltaQik-(piX*pkX*pnXk*hX)/pXk;
553                   }
554 
555                   if (j==k)
556                   {
557                      if (pXk > 0) deltaQjk = (pkX*pnkX*pnXk*hX)/pXk;
558                      deltaQjk = deltaQjk-pnknX*pknX*hnX;
559                   }
560                   else
561                   {
562                      deltaQjk = pjnX*pknX*hnX;
563                      if (pXk > 0) deltaQjk = deltaQjk-(pjX*pkX*pnXk*hX)/pXk;
564                   }
565                      
566                   var_pXk = model->components[X][k].var_pXk;
567                   covar += deltaQik*deltaQjk*var_pXk;
568                }
569 
570                model_X_cov->el[i][j] = covar;
571                if (i != j) model_X_cov->el[j][i] = covar;
572             }
573          }
574          mat_accum(block->component_cov_data, 1, data_X_cov);
575          mat_accum(block->component_cov_model, 1, model_X_cov);
576       }
577    }
578 
579    mat_prod(cov_times_amp, amp, block->component_cov_data);
580    mat_prod(block->component_cov_data, cov_times_amp, amp_trans);
581    mat_prod(cov_times_amp, amp, block->component_cov_model);
582    mat_prod(block->component_cov_model, cov_times_amp, amp_trans);
583 
584    if (scale)
585    {
586       mat_times(block->model_fit, block->component_cov_data);
587       mat_times(block->model_fit, block->component_cov_model);
588    }
589 
590    /* Combine data and model covariances across all components to determine overall class cov */
591    mat_accum(block->component_cov_total, 1, block->component_cov_data);
592    mat_accum(block->component_cov_total, 1, block->component_cov_model);
593    mat_transp(class_map_trans, model->class_map);
594    mat_prod(class_times_component_cov, model->class_map, block->component_cov_total);
595    mat_prod(block->class_cov, class_times_component_cov, class_map_trans);
596 
597    mat_free(data_X_cov);
598    mat_free(model_X_cov);
599    mat_free(cov_times_amp);
600    mat_free(class_times_component_cov);
601    mat_free(class_map_trans);
602 }
603 
604 void qpr_compute_block_cov(QPRDataBlock *block, QPRModel *model, Bool scale)
605 {
606    Mat *amp, *amp_trans;
607 
608    amp = mat_make(model->component_count, model->component_count);
609    amp_trans = mat_make(model->component_count, model->component_count);
610 
611    compute_amp_matrix(block, model, amp, amp_trans);
612    compute_quantity_cov(block, model, amp, amp_trans, scale);
613 
614    mat_free(amp);
615    mat_free(amp_trans);
616 }
617 
618 void qpr_compute_data_cov(QPRData *data, QPRModel *model, Bool scale)
619 {
620    int r;
621 
622    mat_zero(data->class_cov);
623 
624    for (r=0; r<data->block_count; r++)
625    {
626       qpr_compute_block_cov(data->blocks[r], model, scale);
627       mat_accum(data->class_cov, 1, data->blocks[r]->class_cov);
628       /********** WRONG! region covs not independent ****************/
629    }
630 }
631 
632 static double model_fit(QPRDataBlock *block, QPRModel *model, Bool include_model_var,
633                         int min_data, int max_data)
634 {
635    double mX, hX, qk, fit;
636    double diff, sum, sqrt_mX_var, sqrt_hX_var;
637    int k, X, n;
638 
639    n = 0; sum = 0;
640 
641    for (X=0; X<model->bin_count; X++)
642    {
643       hX = block->hist[X];
644       if (hX >= min_data && hX <= max_data)
645       {
646          sqrt_hX_var = 0.25; /* Assume Poisson H(X) */
647          mX = 0;
648          sqrt_mX_var = 0;
649          
650          if (model->components[X] != NULL)
651          {
652             for (k=0; k<model->component_count; k++)
653             {
654                qk = block->component_quantities[k];
655                mX += model->components[X][k].pXk*qk;
656             }
657             if (include_model_var)
658             {
659                for (k=0; k<model->component_count; k++)
660                {
661                   if (mX > 0)
662                   {
663                      qk = block->component_quantities[k];
664                      sqrt_mX_var += (qk*qk*model->components[X][k].var_pXk)/(4*mX);
665                   }
666                }
667             }
668          }
669 
670          diff = sqrt(mX) - sqrt(hX);
671          sum += (diff*diff) / (sqrt_mX_var + sqrt_hX_var);
672          n++;
673       }
674    }
675 
676    fit = sum/((double)n-model->component_count);
677 
678    if (isnan(fit) || isinf(fit)) fit = INT_MAX;
679    return fit;
680 }
681 
682 double qpr_block_model_fit(QPRDataBlock *block, QPRModel *model,
683                            Bool include_model_var)
684 {
685    return model_fit(block, model, include_model_var, 1, INT_MAX);
686 }
687 
688 double qpr_data_model_fit(QPRData *data, QPRModel *model,
689                           Bool include_model_var)
690 {
691    int r, n;
692    double fit, fit_sum;
693 
694    fit_sum = 0; n = 0;
695 
696    for (r=0; r<data->block_count; r++)
697    {
698       fit = qpr_block_model_fit(data->blocks[r], model, include_model_var);
699       if (fit != INT_MAX)
700       {
701          fit_sum += fit;
702          n++;
703       }
704    }
705 
706    fit = fit_sum/(double)n;
707    return fit;
708 }
709 
710 double qpr_analyse_block(QPRDataBlock *block, QPRModel *model, double em_limit, Bool reset)
711 {
712    double fit, old_fit, fit_change;
713 
714    if (reset) set_uniform_quantities(block, model);
715 
716    old_fit = 0;
717    do
718    {
719       compute_pkX(block, model);
720       compute_quantities(block, model);
721       fit = qpr_block_model_fit(block, model, true);
722       fit_change = fabs(old_fit-fit);
723       old_fit = fit;
724    }
725    while (fit_change > em_limit);
726 
727    block->model_fit = fit;
728    
729    return fit;
730 }
731 
732 double qpr_analyse_data(QPRData *data, QPRModel *model, double em_limit, Bool reset)
733 {
734    int r, n, k;
735    double fit, fit_sum;
736 
737    fit_sum = 0; n = 0;
738 
739    for (k=0; k<data->class_count; k++)
740    {
741       data->class_quantities[k] = 0;
742    }
743    data->unknown = 0;
744    data->poorly_modelled = 0;
745    data->total_quantity = 0;
746 
747    for (r=0; r<data->block_count; r++)
748    {
749       fit = qpr_analyse_block(data->blocks[r], model, em_limit, reset);
750       
751      /* if (fit != INT_MAX)
752       {*/
753          fit_sum += fit;
754          n++;
755      /* }*/
756 
757       /* Should check here to exclude bad regions 
758        if (....)
759       ********************************************/
760 
761       for (k=0; k<data->class_count; k++)
762       {
763          data->class_quantities[k] += data->blocks[r]->class_quantities[k];
764       }
765       data->unknown += data->blocks[r]->unknown;
766       data->poorly_modelled += data->blocks[r]->poorly_modelled;
767       data->total_quantity += data->blocks[r]->total_quantity;      
768    }
769 
770    fit = fit_sum/(double)n;
771    data->model_fit = fit;
772 
773    return fit;
774 }
775 
776 static void compute_mean_pXk(QPRData *data, QPRModel *model)
777 {
778    int r, k, X;
779    double pkX, pknX;
780    double rXk, rnXk;
781    double hXk, hnXk;
782    double hX, hnX, qk;
783    QPRDataBlock *block;
784 
785    for (X=0; X<model->bin_count; X++)
786    {
787       if (model->components[X] != NULL)
788       {
789          for (k=0; k<model->component_count; k++)
790          {
791             if (model->component_types[k] == QPR_NORMAL_COMPONENT)
792             {
793                hXk = 0;
794                hnXk = 0;
795                qk = 0;
796 
797                for (r=0; r<data->block_count; r++)
798                {
799                   block = data->blocks[r];
800 
801                   if (block->components[X] != NULL)
802                   {
803                      pkX = block->components[X][k].pkX;
804                      pknX = block->components[X][k].pknX;
805                      hX = block->hist[X];
806                      hnX = block->total_quantity - hX;
807 
808                      rXk = pkX*hX;
809                      hXk += rXk;
810 
811                      rnXk = pknX*hnX;
812                      hnXk += rnXk;
813 
814                     qk += block->component_quantities[k];
815                   }
816                }
817                model->components[X][k].pXk = (hXk == 0) ? 0 : hXk/qk;
818                model->components[X][k].var_pXk = hXk/(qk*qk);
819             }
820          }
821       }
822    }
823 }
824 
825 double qpr_extract_components(QPRData *data, QPRModel *model, double em_limit)
826 {
827    int r;
828    double mean_fit, old_mean_fit, fit_change;
829    
830    old_mean_fit = INT_MAX;
831 
832    for (r=0; r<data->block_count; r++)
833    {
834       set_uniform_quantities(data->blocks[r], model);
835    }
836 
837    do
838    {
839       qpr_analyse_data(data, model, em_limit, false);
840       compute_mean_pXk(data, model);
841       mean_fit = qpr_data_model_fit(data, model, false);
842       fit_change = fabs(old_mean_fit-mean_fit);
843       old_mean_fit = mean_fit;
844    }
845    while (fit_change > em_limit);
846 
847    return mean_fit;
848 }
849 
850 void qpr_randomise_model(QPRModel *model)
851 {
852    int X, k;
853    
854    for (X=0; X<model->bin_count; X++)
855    {
856       for (k=0; k<model->component_count; k++)
857       {
858          qpr_model_fill(model, X, k, 0, (double)rand_int(1, 100));
859       }
860    }
861    qpr_normalise_model(model, true);
862 }
863 
864 QPRModel *qpr_auto_build_model(QPRData *data, int max_loops, double em_limit,
865                                double fit_limit, void (*init_model) (QPRModel *),
866                                int build_type)
867 {
868    QPRModel *model = NULL, *best_model;
869    int loop;
870    double fit, best_fit;
871    clock_t start, end;
872    double elapsed;
873 
874    if (build_type != QPR_RAND_SEED && 
875        build_type != QPR_RAND_ICA &&
876        build_type != QPR_RAND_SEED_THEN_ICA)
877    {
878       error("qpr_auto_build_model: Invalid build_type\n", non_fatal);
879       return NULL;
880    }
881 
882    start = clock();
883    
884    best_fit = INT_MAX;
885    best_model = NULL;
886    loop = 0;
887 
888    printf("Training using %i regions\n", data->block_count);
889 
890    while (best_fit >= fit_limit && loop < max_loops)
891    {
892       if (model != NULL)
893       {
894          qpr_free_model(model);
895          model = NULL;
896       }
897       model = qpr_alloc_model(data->bin_count, data->component_count, 1);
898       if (model == NULL)
899       {
900          error("qpr_auto_build_model: Unable to alloc model\n", non_fatal);
901          return NULL;
902       }
903 
904       if (build_type == QPR_RAND_SEED ||
905           build_type == QPR_RAND_SEED_THEN_ICA)
906       {
907          printf("Seed attempt %i, ", loop);
908          (*init_model)(model);
909          qpr_analyse_data(data, model, em_limit, true);
910          fit = qpr_data_model_fit(data, model, false);
911          printf("fit %f using %i components, ", fit, data->component_count);
912       }
913       else /* build_type == QPR_RAND_ICA */
914       {
915          printf("ICA attempt %i, ", loop);
916          (*init_model)(model);
917          fit = qpr_extract_components(data, model, em_limit);
918          printf("fit %f using %i components, ", fit, data->component_count);
919       }
920 
921       if (!isnan(fit) && fit < best_fit)
922       {
923          best_fit = fit;
924          best_model = model;
925          model = NULL;
926       }
927       printf("current best fit %f\n", best_fit);
928 
929       loop++;
930    }
931 
932    if (build_type == QPR_RAND_SEED_THEN_ICA && best_model != NULL)
933    {
934       format("Components %i, seed best fit: %f, ", data->component_count, best_fit);
935 
936       printf("ICA refinement, ");
937       best_fit = qpr_extract_components(data, best_model, em_limit);
938       printf("best fit %f using %i components\n", best_fit, data->component_count);
939 
940       format("ICA best fit: %f\n", best_fit);
941    }
942    
943    if (model != NULL) qpr_free_model(model);
944 
945    end = clock();
946    elapsed = ((double) (end - start)) / CLOCKS_PER_SEC;
947    printf("Took %f seconds\n", elapsed);
948 
949    return best_model;
950 }
951 
952 void qpr_optimise_initial_quantities(QPRDataBlock *block, QPRModel *model, int runs,
953                                      double total_quantity)
954 {
955    int k, r;
956    double *best_q, best_fit;
957    double sum, fit;
958 
959    best_q = dvector_alloc(0, model->component_count);
960    best_fit = INT_MAX;
961 
962    for (r=0; r<runs; r++)
963    {
964       sum = 0;
965       for (k=0; k<model->component_count; k++)
966       {
967          block->component_quantities[k] = (double)rand_int(10, 10000);
968          sum += block->component_quantities[k];
969       }
970       for (k=0; k<model->component_count; k++)
971       {
972          block->component_quantities[k] *= (total_quantity/sum);
973       }
974 
975       fit = qpr_block_model_fit(block, model, true);
976       if (fit < best_fit)
977       {
978          best_fit = fit;
979          for (k=0; k<model->component_count; k++)
980          {
981             best_q[k] = block->component_quantities[k];
982          }
983       }
984    }
985 
986    for (k=0; k<model->component_count; k++)
987    {
988       block->component_quantities[k] = best_q[k];
989    }
990 
991    dvector_free(best_q, 0);
992 }
993 
994 

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.