source: orange/source/orange/rulelearner.cpp @ 10531:8bb17eeb5895

Revision 10531:8bb17eeb5895, 79.6 KB checked in by anze <anze.staric@…>, 2 years ago (diff)

Fixed line endings.

Line 
1/*
2    This file is part of Orange.
3   
4    Copyright 1996-2010 Faculty of Computer and Information Science, University of Ljubljana
5    Contact: janez.demsar@fri.uni-lj.si
6
7    Orange is free software: you can redistribute it and/or modify
8    it under the terms of the GNU General Public License as published by
9    the Free Software Foundation, either version 3 of the License, or
10    (at your option) any later version.
11
12    Orange is distributed in the hope that it will be useful,
13    but WITHOUT ANY WARRANTY; without even the implied warranty of
14    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15    GNU General Public License for more details.
16
17    You should have received a copy of the GNU General Public License
18    along with Orange.  If not, see <http://www.gnu.org/licenses/>.
19*/
20
21#include "filter.hpp"
22#include "table.hpp"
23#include "stat.hpp"
24#include "measures.hpp"
25#include "discretize.hpp"
26#include "distvars.hpp"
27#include "classfromvar.hpp"
28#include "progress.hpp"
29
30
31#include "rulelearner.ppp"
32
33DEFINE_TOrangeVector_classDescription(PRule, "TRuleList", true, ORANGE_API)
34DEFINE_TOrangeVector_classDescription(PEVDist, "TEVDistList", true, ORANGE_API)
35
36#ifdef _MSC_VER
37#if _MSC_VER < 1300
38template<class T>
39inline T &min(const T&x, const T&y)
40{ return x<y ? x : y; }
41#endif
42#endif
43
44TRule::TRule()
45: weightID(0),
46  quality(ILLEGAL_FLOAT),
47  complexity(0),
48  coveredExamples(NULL),
49  coveredExamplesLength(-1),
50  parentRule(NULL),
51  chi(0.0),
52  requiredConditions(0),
53  baseDist(NULL)
54{}
55
56
57TRule::TRule(PFilter af, PClassifier cl, PLearner lr, PDistribution dist, PExampleTable ce, const int &w, const float &qu)
58: filter(af),
59  classifier(cl),
60  learner(lr),
61  classDistribution(dist),
62  examples(ce),
63  weightID(w),
64  quality(qu),
65  coveredExamples(NULL),
66  coveredExamplesLength(-1),
67  parentRule(NULL),
68  chi(0.0),
69  valuesFilter(NULL),
70  requiredConditions(0),
71  baseDist(dist)
72{}
73
74
75TRule::TRule(const TRule &other, bool copyData)
76: filter(other.filter? other.filter->deepCopy() : PFilter()), //
77  valuesFilter(other.valuesFilter? other.valuesFilter->deepCopy() : PFilter()),
78  classifier(other.classifier),
79  learner(other.learner),
80  complexity(other.complexity),
81  classDistribution(copyData ? other.classDistribution: PDistribution()),
82  examples(copyData ? other.examples : PExampleTable()),
83  weightID(copyData ? other.weightID : 0),
84  quality(copyData ? other.quality : ILLEGAL_FLOAT),
85  coveredExamples(copyData && other.coveredExamples && (other.coveredExamplesLength >= 0) ? (int *)memcpy(new int[other.coveredExamplesLength], other.coveredExamples, other.coveredExamplesLength) : NULL),
86  coveredExamplesLength(copyData ? other.coveredExamplesLength : -1),
87  parentRule(other.parentRule),
88  chi(other.chi),
89  baseDist(other.baseDist),
90  requiredConditions(other.requiredConditions)
91{}
92
93
94TRule::~TRule()
95{   delete coveredExamples; }
96
97bool TRule::operator ()(const TExample &ex)
98{
99  checkProperty(filter);
100  return filter->call(ex);
101}
102
103
104#define HIGHBIT 0x80000000
105
106PExampleTable TRule::operator ()(PExampleTable gen, const bool ref, const bool negate)
107{
108  checkProperty(filter);
109
110  TExampleTable *table = ref ? mlnew TExampleTable(gen, 1) : mlnew TExampleTable(PExampleGenerator(gen));
111  PExampleGenerator wtable = table;
112
113  PEITERATE(ei, gen)
114    if (filter->call(*ei) != negate)
115      table->addExample(*ei);
116
117  return wtable;
118}
119
120
121void TRule::filterAndStore(PExampleTable gen, const int &wei, const int &targetClass, const int *prevCovered, const int anExamples)
122{
123  checkProperty(filter);
124  examples=this->call(gen);
125  weightID = wei;
126  classDistribution = getClassDistribution(examples, wei);
127  if (classDistribution->abs==0)
128    return;
129
130  if (targetClass>=0)
131    classifier = mlnew TDefaultClassifier(gen->domain->classVar, TValue(targetClass), classDistribution);
132  else if (learner) {
133    classifier = learner->call(examples,wei);
134  }
135  else
136    classifier = mlnew TDefaultClassifier(gen->domain->classVar, classDistribution);
137/*  if (anExamples > 0) {
138    const int bitsInInt = sizeof(int)*8;
139    coveredExamplesLength = anExamples/bitsInInt + 1;
140    coveredExamples = (int *)malloc(coveredExamplesLength);
141    if (prevCovered) {
142      memcpy(coveredExamples, prevCovered, coveredExamplesLength);
143
144      int *cei = coveredExamples-1;
145      int mask = 0;
146      int inBit = 0;
147
148      PEITERATE(ei, gen) {
149        if (!(*cei & mask)) {
150          if (inBit)
151            *cei = *cei << inBit;
152          while(!*++cei);
153          mask = -1;
154          inBit = bitsInInt;
155        }
156
157        while( (*cei & HIGHBIT) == 0) {
158          *cei = *cei << 1;
159          *cei = mask << 1;
160          inBit--;
161        }
162
163        if (filter->call(*ei)) {
164          *cei = (*cei << 1) | 1;
165          table->addExample(*ei);
166        }
167        else
168          *cei = *cei << 1;
169
170        mask = mask << 1;
171        inBit--;
172      }
173    }
174
175    else {
176      int *cei = coveredExamples;
177      int inBit = bitsInInt;
178
179      PEITERATE(ei, gen) {
180        if (filter->call(*ei)) {
181          *cei = (*cei << 1) | 1;
182          table->addExample(*ei);
183        }
184        else
185          *cei = *cei << 1;
186
187        if (!--inBit) {
188          inBit = bitsInInt;
189          cei++;
190        }
191      }
192      *cei = *cei << inBit;
193    }
194  } */
195}
196
197
198
199bool haveEqualValues(const TRule &r1, const TRule &r2)
200{
201  const TDefaultClassifier *clsf1 = r1.classifier.AS(TDefaultClassifier);
202  const TDefaultClassifier *clsf2 = r2.classifier.AS(TDefaultClassifier);
203  if (!clsf1 || !clsf2)
204    return false;
205
206  const TDiscDistribution *dist1 = dynamic_cast<const TDiscDistribution *>(clsf1->defaultDistribution.getUnwrappedPtr());
207  const TDiscDistribution *dist2 = dynamic_cast<const TDiscDistribution *>(clsf2->defaultDistribution.getUnwrappedPtr());
208
209  float high1 = dist1->highestProb();
210  float high2 = dist2->highestProb();
211
212  for(TDiscDistribution::const_iterator d1i(dist1->begin()), d1e(dist1->end()), d2i(dist2->begin()), d2e(dist2->end());
213      (d1i!=d1e) && (d2i!=d2e);
214      d1i++, d2i++)
215    if ((*d1i == high1) && (*d2i == high2))
216      return true;
217
218  return false;
219}
220
221
222bool TRule::operator <(const TRule &other) const
223{
224  if (!haveEqualValues(*this, other))
225    return false;
226
227  bool different = false;
228
229  if (coveredExamples && other.coveredExamples) {
230    int *c1i = coveredExamples;
231    int *c2i = other.coveredExamples;
232    for(int i = coveredExamplesLength; i--; c1i++, c2i++) {
233      if (*c1i & ~*c2i)
234        return false;
235      if (*c1i != *c2i)
236        different = true;
237    }
238  }
239  else {
240    raiseError("operator not implemented yet");
241  }
242
243  return different;
244}
245
246
247bool TRule::operator <=(const TRule &other) const
248{
249  if (!haveEqualValues(*this, other))
250    return false;
251
252  if (coveredExamples && other.coveredExamples) {
253    int *c1i = coveredExamples;
254    int *c2i = other.coveredExamples;
255    for(int i = coveredExamplesLength; i--; c1i++, c2i++) {
256      if (*c1i & ~*c2i)
257        return false;
258    }
259  }
260
261  else {
262    raiseError("operator not implemented yet");
263  }
264
265  return true;
266}
267
268
269bool TRule::operator >(const TRule &other) const
270{
271  if (!haveEqualValues(*this, other))
272    return false;
273
274  bool different = false;
275  if (coveredExamples && other.coveredExamples) {
276    int *c1i = coveredExamples;
277    int *c2i = other.coveredExamples;
278    for(int i = coveredExamplesLength; i--; c1i++, c2i++) {
279      if (~*c1i & *c2i)
280        return false;
281      if (*c1i != *c2i)
282        different = true;
283    }
284  }
285
286  else {
287    raiseError("operator not implemented yet");
288  }
289
290  return different;
291}
292
293
294bool TRule::operator >=(const TRule &other) const
295{
296  if (!haveEqualValues(*this, other))
297    return false;
298
299  if (coveredExamples && other.coveredExamples) {
300    int *c1i = coveredExamples;
301    int *c2i = other.coveredExamples;
302    for(int i = coveredExamplesLength; i--; c1i++, c2i++) {
303      if (~*c1i & *c2i)
304        return false;
305    }
306  }
307
308  else {
309    raiseError("operator not implemented yet");
310  }
311
312  return true;
313}
314
315
316bool TRule::operator ==(const TRule &other) const
317{
318  if (!haveEqualValues(*this, other))
319    return false;
320
321  if (coveredExamples && other.coveredExamples) {
322    return !memcmp(coveredExamples, other.coveredExamples, coveredExamplesLength);
323  }
324
325  else {
326    raiseError("operator not implemented yet");
327  }
328
329  return false;
330}
331
332
333
334TRuleValidator_LRS::TRuleValidator_LRS(const float &a, const float &min_coverage, const int &max_rule_complexity, const float &min_quality)
335: alpha(a),
336  min_coverage(min_coverage),
337  max_rule_complexity(max_rule_complexity),
338  min_quality(min_quality)
339{}
340
341bool TRuleValidator_LRS::operator()(PRule rule, PExampleTable, const int &, const int &targetClass, PDistribution apriori) const
342{
343  const TDiscDistribution &obs_dist = dynamic_cast<const TDiscDistribution &>(rule->classDistribution.getReference());
344  if (!obs_dist.cases || obs_dist.cases < min_coverage)
345    return false;
346  if (max_rule_complexity > -1.0 && rule->complexity > max_rule_complexity)
347    return false;
348  if (min_quality>rule->quality)
349    return false;
350
351  const TDiscDistribution &exp_dist = dynamic_cast<const TDiscDistribution &>(apriori.getReference());
352  if (obs_dist.abs == exp_dist.abs) //it turns out that this happens quite often
353    return false;
354
355  if (alpha >= 1.0)
356    return true;
357
358  if (targetClass == -1) {
359    float lrs = 0.0;
360    for(TDiscDistribution::const_iterator odi(obs_dist.begin()), ode(obs_dist.end()), edi(exp_dist.begin()), ede(exp_dist.end());
361        (odi!=ode); odi++, edi++) {
362      if ((edi!=ede) && (*edi) && (*odi))
363        lrs += *odi * log(*odi / ((edi != ede) & (*edi > 0.0) ? *edi : 1e-5));
364    }
365    lrs = 2 * (lrs - obs_dist.abs * log(obs_dist.abs / exp_dist.abs));
366    return (lrs > 0.0) && (chisqprob(lrs, float(obs_dist.size()-1)) <= alpha);
367  }
368
369  float p = (targetClass < obs_dist.size()) ? obs_dist[targetClass] : 1e-5;
370  const float P = (targetClass < exp_dist.size()) && (exp_dist[targetClass] > 0.0) ? exp_dist[targetClass] : 1e-5;
371  float n = obs_dist.abs - p;
372  float N = exp_dist.abs - P;
373  if (n>=N)
374    return false;
375  if (N<=0.0)
376    N = 1e-6f;
377  if (p<=0.0)
378    p = 1e-6f;
379  if (n<=0.0)
380    n = 1e-6f;
381
382  float lrs = 2 * (p*log(p/P) + n*log(n/N) - obs_dist.abs * log(obs_dist.abs/exp_dist.abs));
383  return (lrs > 0.0) && (chisqprob(lrs, 1.0f) <= alpha);
384}
385
386
387float TRuleEvaluator_Entropy::operator()(PRule rule, PExampleTable, const int &, const int &targetClass, PDistribution apriori)
388{
389  const TDiscDistribution &obs_dist = dynamic_cast<const TDiscDistribution &>(rule->classDistribution.getReference());
390  if (!obs_dist.cases)
391    return -numeric_limits<float>::max();
392
393  if (targetClass == -1)
394    return -getEntropy(dynamic_cast<TDiscDistribution &>(rule->classDistribution.getReference()));
395
396  const TDiscDistribution &exp_dist = dynamic_cast<const TDiscDistribution &>(apriori.getReference());
397
398  float p = (targetClass < obs_dist.size()) ? obs_dist[targetClass] : 0.0;
399  const float P = (targetClass < exp_dist.size()) && (exp_dist[targetClass] > 0.0) ? exp_dist[targetClass] : 1e-5;
400
401  float n = obs_dist.abs - p;
402  float N = exp_dist.abs - P;
403  if (N<=0.0)
404    N = 1e-6f;
405  if (p<=0.0)
406    p = 1e-6f;
407  if (n<=0.0)
408    n = 1e-6f;
409
410  return ((p*log(p) + n*log(n) - obs_dist.abs * log(obs_dist.abs)) / obs_dist.abs);
411}
412
413float TRuleEvaluator_Laplace::operator()(PRule rule, PExampleTable, const int &, const int &targetClass, PDistribution apriori)
414{
415  const TDiscDistribution &obs_dist = dynamic_cast<const TDiscDistribution &>(rule->classDistribution.getReference());
416  if (!obs_dist.cases)
417    return 0;
418
419  float p;
420  if (targetClass == -1) {
421    p = float(obs_dist.highestProb());
422    return (p+1)/(obs_dist.abs+obs_dist.size());
423  }
424  p = float(obs_dist[targetClass]);
425  return (p+1)/(obs_dist.abs+2);
426}
427
428TRuleEvaluator_LRS::TRuleEvaluator_LRS(const bool &sr)
429: storeRules(sr)
430{
431  TRuleList *ruleList = mlnew TRuleList;
432  rules = ruleList;
433}
434
435float TRuleEvaluator_LRS::operator()(PRule rule, PExampleTable, const int &, const int &targetClass, PDistribution apriori)
436{
437  const TDiscDistribution &obs_dist = dynamic_cast<const TDiscDistribution &>(rule->classDistribution.getReference());
438  if (!obs_dist.cases)
439    return 0.0;
440
441  const TDiscDistribution &exp_dist = dynamic_cast<const TDiscDistribution &>(apriori.getReference());
442
443  if (obs_dist.abs >= exp_dist.abs) //it turns out that this happens quite often
444    return 0.0;
445
446  if (targetClass == -1) {
447    float lrs = 0.0;
448    for(TDiscDistribution::const_iterator odi(obs_dist.begin()), ode(obs_dist.end()), edi(exp_dist.begin()), ede(exp_dist.end());
449        (odi!=ode); odi++, edi++) {
450      if ((edi!=ede) && (*edi) && (*odi))
451        lrs += *odi * log(*odi / ((edi != ede) & (*edi > 0.0) ? *edi : 1e-5));
452    }
453    lrs = 2 * (lrs - obs_dist.abs * log(obs_dist.abs / exp_dist.abs));
454    return lrs;
455  }
456
457  float p = (targetClass < obs_dist.size()) ? obs_dist[targetClass]-0.5 : 1e-5f;
458  float P = (targetClass < exp_dist.size()) && (exp_dist[targetClass] > 0.0) ? exp_dist[targetClass] : 1e-5f;
459
460  float n = obs_dist.abs - p;
461  float N = exp_dist.abs - P;
462
463  // border cases
464  if (N<=1e-5f) return 0.0;
465  if (p<=1e-5f) return 0.0;
466  if (n<=1e-5f) n = 1e-5f;
467
468  if (p<=(p+n)*P/(P+N))
469    return 0.0;
470
471  float ep = (p+n)*P/(P+N);
472  float lrs = 2 * (p*log(p/ep) + n*log(n/(p+n))+(P-p)*log((P-p)/(P+N-p-n))+(N-n)*log((N-n)/(P+N-p-n))-(P-p)*log(P/(P+N))-N*log(N/(P+N)));
473  if (storeRules) {
474    TFilter_values *filter;
475    filter = rule->filter.AS(TFilter_values);
476    int ncond = filter->conditions->size();
477    TRuleList &rlist = rules.getReference();
478    if (!rlist.at(ncond) || rlist.at(ncond)->quality < lrs)
479        rlist.at(ncond) = rule;
480  }
481  return lrs;
482}
483
484
485TEVDist::TEVDist(const float & mu, const float & beta, PFloatList & percentiles)
486: mu(mu),
487  beta(beta),
488  percentiles(percentiles)
489{
490  maxPercentile = (float)0.95;
491  step = (float)0.1;
492}
493
494TEVDist::TEVDist()
495{
496  maxPercentile = (float)0.95;
497  step = (float)0.1;
498}
499
500// compute cumulative probability with extreme value distribution
501double TEVDist::getProb(const float & chi)
502{ 
503  // percentiles are not computed
504  if (!percentiles || percentiles->size()==0 || percentiles->at(percentiles->size()-1)<chi)
505    return 1.0-exp(-exp((double)(mu-chi)/beta));
506  // low chi
507  if (chi < percentiles->at(0)-1e-6)
508    return 1.0 - 0.05 * chi / percentiles->at(0);
509  TFloatList::const_iterator pi(percentiles->begin()), pe(percentiles->end());
510  for (int i=0; (pi+1)!=pe; pi++,i++) {
511    float a = *pi;
512    float b = *(pi+1);
513    if (chi>=(a-1e-6) && chi <= (b+1e-6))
514      if (b - a <= 0.0)
515        return maxPercentile-i*step;
516      else
517        return (maxPercentile-i*step)-step*(chi-a)/(b-a);
518  }
519  return 1.0;
520}
521
522float TEVDist::median()
523{
524  if (!percentiles || percentiles->size()==0)
525    return mu + beta*0.36651292; // log(log(2))
526  if (percentiles->size()%2 == 0)
527    return (percentiles->at(percentiles->size()/2-1)+percentiles->at(percentiles->size()/2))/2;
528  else
529    return (percentiles->at(percentiles->size()/2));
530}
531
532TEVDistGetter_Standard::TEVDistGetter_Standard(PEVDistList dists)
533: dists(dists)
534{}
535
536TEVDistGetter_Standard::TEVDistGetter_Standard()
537{}
538
539PEVDist TEVDistGetter_Standard::operator()(const PRule, const int & parentLength, const int & length) const
540{
541  // first element (correction for inter - rule optimism
542  if (!length)
543    return dists->at(0);
544  // optimism between rule length of "parentLength" and true length "length"
545  int indx = length*(length-1)/2 + parentLength + 1;
546  if (dists->size() > indx)
547    return dists->at(indx);
548  return NULL;
549}
550
551float getChi(float p, float n, float P, float N)
552{
553  p = p - 0.5;
554  n = n + 0.5;
555  if (p<=(p+n)*P/(P+N))
556      return 0.0;
557  float ep = (p+n)*P/(P+N);
558  return 2*(p*log(p/ep)+n*log(n/(p+n))+(P-p)*log((P-p)/(P+N-p-n))+(N-n)*log((N-n)/(P+N-p-n))-(P-p)*log(P/(P+N))-N*log(N/(P+N)));
559}
560
561LNLNChiSq::LNLNChiSq(PEVDist evd, const float & chi, const float & priorProb)
562: evd(evd),
563  chi(chi)
564{
565  // extreme alpha = probability of obtaining chi with such extreme value distribution
566  extremeAlpha = evd->getProb(chi);
567  if (extremeAlpha < 1.0-evd->maxPercentile)
568    extremeAlpha = -1.0;
569  // exponent used in FT cumulative function (min because it is negative)
570  exponent = min(float(log(log(1/evd->maxPercentile))),(evd->mu-chi)/evd->beta); //(evd->mu-chi)/evd->beta; //min(float(log(log(1/evd->maxPercentile))),(evd->mu-chi)/evd->beta);
571  pp = priorProb;
572}
573
574double LNLNChiSq::operator()(float chix) {
575  if (chix > 1400)
576    return -1000.0;
577
578  double chip;
579  if (chix<=0.0)
580    chip = 0.5;
581  else
582    chip = chisqprob((double)chix,1.0)*0.5;//pp; // in statc
583
584  //chip = chisqprob((double)chix,1.0); // in statc
585/*  if (chi > 5.2 && chi < 5.3)
586  printf("chip = %f, chip-e = %f\n", chip, chip-extremeAlpha); */
587
588  if (extremeAlpha > 0.0)
589    return chip-extremeAlpha;
590
591  if (chip<=0.0)
592    return -1000.0;
593
594  if (chip < 1e-6)
595    return log(chip)-exponent;
596  return log(-log(1-chip))-exponent;
597}
598
599LRInv::LRInv(float & n, float & P, float & N, float chiCorrected)
600: n(n),
601  P(P),
602  chiCorrected(chiCorrected),
603  N(N)
604{}
605
606double LRInv::operator()(float p){
607  // check how it is done in call
608  return getChi(p,n-p,P,N-P) - chiCorrected;
609}
610
611LRInvMean::LRInvMean(float correctedP, PRule rule, PRule groundRule, const int & targetClass)
612: n(rule->classDistribution->abs),
613  p(correctedP),
614  P(groundRule->classDistribution->atint(targetClass)),
615  N(groundRule->classDistribution->abs)
616{}
617
618double LRInvMean::operator()(float pc){
619  return - getChi(p,n-p,pc*N/n,N-pc*N/n) + 0.30;
620}
621
622
623LRInvE::LRInvE(PRule rule, PRule groundRule, const int & targetClass, float chiCorrected)
624: n(rule->classDistribution->abs),
625  p(rule->classDistribution->atint(targetClass)),
626  N(groundRule->classDistribution->abs),
627  chiCorrected(chiCorrected)
628{}
629
630
631double LRInvE::operator()(float P){
632  // check how it is done in call
633  P *= N/n;
634  return - getChi(p,n-p,P,N-P) + chiCorrected;
635}
636
637// Implementation of Brent's root finding method.
638float brent(const float & minv, const float & maxv, const int & maxsteps, DiffFunc * func, float threshold)
639{
640  float a = minv;
641  float b = maxv;
642  float fa = func->call(a);
643  float fb = func->call(b);
644
645
646 // float threshold = 0.01 * (maxv - minv);
647  if (fb>0 && fa>0 && fb>fa || fb<0 && fa<0 && fb<fa)
648    return a;
649    if (fb>0 && fa>0 && fb<fa || fb<0 && fa<0 && fb>fa)
650    return b;
651
652  float c = a; // c is previous value of b
653  float fe, fc = fa;
654  float m = 0.0, e = 0.0, d = 0.0;
655  int counter = 0;
656  while (1) {
657    counter += 1;
658    if (fb == fa)
659      return b;
660    else if (fb!=fc && fa!=fc)
661      d = a*fb*fc/(fa-fb)/(fa-fc)+b*fa*fc/(fb-fa)/(fb-fc)+c*fa*fb/(fc-fa)/(fc-fb);
662    else
663      d = b-fb*(b-a)/(fb-fa);
664
665    m = (a+b)/2;
666    if (d<=m && d>=b || d>=m && d<=b)
667      e = d;
668    else
669      e = m;
670    fe = func->call(e);
671    if (fe*fb<0) {
672      a = b;
673      fa = fb;
674    }
675    c = b;
676    fc = fb;
677    b = e;
678    fb = fe;
679
680/*    if (counter<maxsteps)
681        printf("a = %f, b=%f, fa=%f, fb=%f, maxv=%f, minv=%f\n",a,b,fa,fb,maxv,minv); */
682    //  return 0.0;
683    if ((b>threshold && fb*func->call(b-threshold)<=0) || fb*func->call(b+threshold)<=0)
684    {
685      return b;
686    }
687    if (abs(a-b)<threshold && fa*fb<0)
688      return (a+b)/2.;
689    if (fb*fa>0 || b>maxv || b<minv)
690      return 0.0;
691  }
692}
693
694TRuleEvaluator_mEVC::TRuleEvaluator_mEVC(const int & m, PEVDistGetter evDistGetter, PVariable probVar, PRuleValidator validator, const int & min_improved, const float & min_improved_perc, const int & optimismReduction)
695: m(m),
696  evDistGetter(evDistGetter),
697  probVar(probVar),
698  validator(validator),
699  min_improved(min_improved),
700  min_improved_perc(min_improved_perc),
701  bestRule(NULL),
702  ruleAlpha(1.0),
703  attributeAlpha(1.0),
704  optimismReduction(optimismReduction)
705{}
706
707TRuleEvaluator_mEVC::TRuleEvaluator_mEVC()
708: m(0),
709  evDistGetter(NULL),
710  probVar(NULL),
711  validator(NULL),
712  min_improved(1),
713  min_improved_perc(0),
714  bestRule(NULL),
715  ruleAlpha(1.0),
716  attributeAlpha(1.0),
717  optimismReduction(0)
718{}
719
720void TRuleEvaluator_mEVC::reset()
721{
722  bestRule = NULL;
723}
724
725
726/* The method validates rule's attributes' significance with respect to its extreme value corrected distribution. */
727bool TRuleEvaluator_mEVC::ruleAttSignificant(PRule rule, PExampleTable examples, const int & weightID, const int &targetClass, PDistribution apriori, float & aprioriProb)
728{
729  PEVDist evd;
730
731  // Should classical LRS be used, or EVC corrected?
732  bool useClassicLRS = false;
733  if (!optimismReduction)
734    useClassicLRS = true;
735  if (!useClassicLRS)
736  {
737    evd = evDistGetter->call(rule, 0, 0);
738    if (evd->mu < 1.0)
739      useClassicLRS = true;
740  }
741
742  TFilter_values *filter;
743  if (rule->valuesFilter)
744    filter = rule->valuesFilter.AS(TFilter_values);
745  else
746    filter = rule->filter.AS(TFilter_values);
747
748  // Loop through all attributes - remove each and check significance
749  bool rasig = true;
750  int i,j;
751  float chi;
752  TFilter_values *newfilter;
753  for (i=0; i<filter->conditions->size(); i++)
754  {
755      if (i<rule->requiredConditions)
756        continue;
757      // create a rule without one condition
758      TRule *newRule = new TRule();
759      PRule wnewRule = newRule;
760      wnewRule->filter = new TFilter_values();
761      wnewRule->filter->domain = examples->domain;
762      wnewRule->complexity = rule->complexity - 1;
763      newfilter = newRule->filter.AS(TFilter_values);
764      for (j=0; j<filter->conditions->size(); j++)
765        if (j!=i)
766          newfilter->conditions->push_back(filter->conditions->at(j));
767      wnewRule->filterAndStore(examples, weightID, targetClass);
768
769      // compute lrs of rule vs new rule (without one condtion)
770      if (!rule->classDistribution->abs || wnewRule->classDistribution->abs == rule->classDistribution->abs)
771        return false;
772      chi = getChi(rule->classDistribution->atint(targetClass),
773                   rule->classDistribution->abs - rule->classDistribution->atint(targetClass),
774                   wnewRule->classDistribution->atint(targetClass),
775                   wnewRule->classDistribution->abs - wnewRule->classDistribution->atint(targetClass));
776      // correct lrs with evc
777      if (!useClassicLRS) {
778        LNLNChiSq *diffFunc = new LNLNChiSq(evd,chi,aprioriProb);
779            chi = brent(0.0,chi,100, diffFunc, 0.1f);
780        delete diffFunc;
781      }
782      // check significance
783      rasig = rasig & ((chi > 0.0) && (chisqprob(chi, 1.0f) <= attributeAlpha));
784      if (!rasig)
785        return false;
786  }
787  return true;
788}
789
790float combineEPositives(float N, float P, float oldQ, float n, float q)
791{
792  if (oldQ >= P/N)
793      return q*n;
794  if (P <= 0.1)
795      return 0.0;
796  return N*oldQ/P*q*n;
797}
798
799float TRuleEvaluator_mEVC::evaluateRulePessimistic(PRule rule, PExampleTable examples, const int & weightID, const int &targetClass, PDistribution apriori, const int & rLength, const float & aprioriProb)
800{
801  // extreme value distribution; optimism from ground Rule to rule
802  PEVDist evd = evDistGetter->call(rule, 0, 0);
803
804  if (!evd || evd->mu < 0.0)
805    return -10e+6;
806  //printf("mu=%f, beta=%f\n",evd->mu,evd->beta);
807  if (evd->mu == 0.0 || !rule->parentRule)
808  {
809    // return as if rule distribution is not optimistic
810    rule->chi = getChi(rule->classDistribution->atint(targetClass), rule->classDistribution->abs - rule->classDistribution->atint(targetClass),
811                apriori->atint(targetClass), apriori->abs - apriori->atint(targetClass));
812      rule->estRF = rule->classDistribution->atint(targetClass)/rule->classDistribution->abs;
813    return (rule->classDistribution->atint(targetClass)+m*aprioriProb)/(rule->classDistribution->abs+m);
814  }
815
816  // rule's improvement chi
817  float chi = getChi(rule->classDistribution->atint(targetClass), rule->classDistribution->abs - rule->classDistribution->atint(targetClass),
818                   rule->parentRule->classDistribution->atint(targetClass), rule->parentRule->classDistribution->abs - rule->parentRule->classDistribution->atint(targetClass));
819
820  float ePos = 0.0;
821  float median = evd->median();
822  float rule_acc = rule->classDistribution->atint(targetClass)/rule->classDistribution->abs;
823  float parent_acc = rule->parentRule->classDistribution->atint(targetClass)/rule->parentRule->classDistribution->abs;
824
825   // need correcting? if rule very good, correcting will not change anything
826  if ((evd->mu-chi)/evd->beta < -100)
827      // return as if rule distribution is not optimistic
828    ePos = rule->classDistribution->atint(targetClass);
829  else if (rule_acc < parent_acc)
830    ePos = rule->classDistribution->atint(targetClass);
831  else if (chi<=median)
832    ePos = rule->classDistribution->abs * parent_acc;
833  else {
834      // compute ePos
835      LRInvE *diffFunc = new LRInvE(rule, rule->parentRule, targetClass, median);
836      ePos = brent(rule->parentRule->classDistribution->atint(targetClass)/rule->parentRule->classDistribution->abs*rule->classDistribution->atint(targetClass), rule->classDistribution->atint(targetClass), 100, diffFunc, 0.1f);
837      delete diffFunc;
838
839  }
840
841  float ePosOA; // expected positive examples considering base rule also
842  ePosOA = combineEPositives(rule->parentRule->classDistribution->abs, rule->parentRule->classDistribution->atint(targetClass), rule->parentRule->estRF, rule->classDistribution->abs, ePos/rule->classDistribution->abs);
843  rule->chi = getChi(ePosOA, rule->classDistribution->abs - ePosOA,
844                apriori->atint(targetClass), apriori->abs - apriori->atint(targetClass));
845
846  rule->estRF = ePosOA/rule->classDistribution->abs;
847  float quality = (ePosOA + m*aprioriProb)/(rule->classDistribution->abs+m);
848
849  if (quality > aprioriProb)
850    return quality;
851  return aprioriProb-0.01+0.01*rule_acc;
852}
853
854float TRuleEvaluator_mEVC::evaluateRuleM(PRule rule, PExampleTable examples, const int & weightID, const int &targetClass, PDistribution apriori, const int & rLength, const float & aprioriProb)
855{
856  if (!m & !rule->classDistribution->abs)
857    return 0.0;
858  rule->chi = getChi(rule->classDistribution->atint(targetClass), rule->classDistribution->abs - rule->classDistribution->atint(targetClass),
859                     apriori->atint(targetClass), apriori->abs - apriori->atint(targetClass));
860  float p = rule->classDistribution->atint(targetClass)+m*apriori->atint(targetClass)/apriori->abs;
861  p = p/(rule->classDistribution->abs+m);
862  return p;
863}
864
865// evaluates a rule, which can have its base conditions, base rule
866// can be evaluted in any way.
867float TRuleEvaluator_mEVC::evaluateRuleEVC(PRule rule, PExampleTable examples, const int & weightID, const int &targetClass, PDistribution apriori, const int & rLength, const float & aprioriProb)
868{
869  // extreme value distribution; optimism from parent Rule to rule
870  PEVDist evd = evDistGetter->call(rule, 0, rLength - rule->requiredConditions);
871
872  // get base distribution (normally apriori, otherwise distribution of argument)
873  PDistribution base = rule->baseDist;
874  float baseProb = base->atint(targetClass)/base->abs;
875
876  if (!evd || evd->mu < 0.0)
877    return -10e+6;
878  //printf("mu=%f, beta=%f\n",evd->mu,evd->beta);
879  if (evd->mu < 1.0001)
880  {
881    // return as if rule distribution is not optimistic
882    rule->chi = getChi(rule->classDistribution->atint(targetClass), rule->classDistribution->abs - rule->classDistribution->atint(targetClass),
883                apriori->atint(targetClass), apriori->abs - apriori->atint(targetClass));
884    rule->estRF = rule->classDistribution->atint(targetClass)/rule->classDistribution->abs;
885    rule->distP = rule->classDistribution->atint(targetClass);
886    return (rule->classDistribution->atint(targetClass)+m*aprioriProb)/(rule->classDistribution->abs+m);
887  }
888
889  // rule's chi
890  float chi = getChi(rule->classDistribution->atint(targetClass), rule->classDistribution->abs - rule->classDistribution->atint(targetClass),
891                   base->atint(targetClass), base->abs - base->atint(targetClass));
892  float ePos = 0.0;
893  float median = evd->median();
894  float rule_acc = rule->classDistribution->atint(targetClass)/rule->classDistribution->abs;
895
896  if ((evd->mu-chi)/evd->beta < -500)
897    ePos = rule->classDistribution->atint(targetClass);
898  if (rule_acc < baseProb)
899    ePos = rule->classDistribution->atint(targetClass);
900  else if (chi <= median+1e-6)
901      ePos = rule->classDistribution->abs * baseProb;
902  else {
903      // correct chi
904      LNLNChiSq *diffFunc = new LNLNChiSq(evd,chi,aprioriProb);
905      rule->chi = brent(0.0,chi,100, diffFunc, 0.1f); // this is only the correction of one step chi
906      delete diffFunc;
907
908      // compute expected number of positive examples relatively to base rule
909      if (rule->chi > 0.0)
910      {
911        // correct optimism
912
913        float baseTarget = base->atint(targetClass);
914        LRInv *diffFunc = new LRInv(rule->classDistribution->abs, baseTarget, base->abs, rule->chi); //-0.45);
915        ePos = brent(base->atint(targetClass)/base->abs*rule->classDistribution->abs, rule->classDistribution->atint(targetClass), 100, diffFunc, 0.1f);
916        //printf("epos = %4.3f\n",ePos);
917        delete diffFunc;
918        if (rule->classDistribution->abs == rule->classDistribution->atint(targetClass))
919        {
920            if (rule->parentRule)
921            {
922                float parentOpt = rule->parentRule->estRF / (rule->parentRule->classDistribution->atint(targetClass) / rule->parentRule->classDistribution->abs);
923                float ruleOpt = ePos / rule->classDistribution->atint(targetClass);
924                if (ruleOpt < parentOpt)
925                    ePos = (ruleOpt + (parentOpt - ruleOpt) * 0.366) * rule->classDistribution->atint(targetClass);
926            }
927            else
928                ePos = (0.634 * ePos/rule->classDistribution->atint(targetClass) + 0.366) * rule->classDistribution->atint(targetClass);
929        }
930      }
931      else
932        ePos = rule->classDistribution->abs * baseProb;
933
934  }
935
936  //
937  rule->chi = getChi(ePos, rule->classDistribution->abs - ePos,
938                apriori->atint(targetClass), apriori->abs - apriori->atint(targetClass));
939
940  rule->estRF = ePos/rule->classDistribution->abs;
941
942  float quality = (ePos + m*aprioriProb)/(rule->classDistribution->abs+m);
943
944  if (quality > aprioriProb)
945    return quality;
946  if (rule_acc < aprioriProb)
947    return rule_acc-0.01;
948  return aprioriProb-0.01+0.01*chi/median; 
949}
950
951float TRuleEvaluator_mEVC::operator()(PRule rule, PExampleTable examples, const int & weightID, const int &targetClass, PDistribution apriori)
952{
953  rule->chi = 0.0;
954  if (!rule->classDistribution->abs || !rule->classDistribution->atint(targetClass))
955    return 0;
956
957  // evaluate rule
958  int rLength = rule->complexity;
959  float aprioriProb = apriori->atint(targetClass)/apriori->abs;
960  float quality;
961
962  if (optimismReduction == 0)
963    quality = evaluateRuleM(rule,examples,weightID,targetClass,apriori,rLength,aprioriProb);
964  else if (optimismReduction == 1)
965    quality = evaluateRulePessimistic(rule,examples,weightID,targetClass,apriori,rLength,aprioriProb);
966  else if (optimismReduction == 2)
967    quality = evaluateRuleEVC(rule,examples,weightID,targetClass,apriori,rLength,aprioriProb);
968
969  if (quality < 0.0)
970    return quality;
971  if (!probVar || !returnExpectedProb)
972    return quality;
973
974  // get rule's probability coverage
975  int improved = 0;
976  PEITERATE(ei, rule->examples)
977    if ((*ei).getClass().intV == targetClass && quality > (*ei)[probVar].floatV)
978      improved ++;
979
980  // compute future quality = expected quality when rule is finalised
981  float bestQuality;
982  float futureQuality = 0.0;
983  if (rule->classDistribution->atint(targetClass) == rule->classDistribution->abs)
984    futureQuality = -1.0;
985  else if (improved >= min_improved &&
986           improved/rule->classDistribution->atint(targetClass) > min_improved_perc*0.01 &&
987           quality > (aprioriProb + 1e-3))
988//    futureQuality = quality;
989    futureQuality = 1.0 + quality;
990  else {
991    PDistribution oldRuleDist = rule->classDistribution;
992    float rulesTrueChi = rule->chi;
993    float rulesDistP = rule->distP;
994    rule->classDistribution = mlnew TDiscDistribution(examples->domain->classVar);
995    rule->classDistribution->setint(targetClass, oldRuleDist->atint(targetClass));
996    rule->classDistribution->abs = rule->classDistribution->atint(targetClass);
997    rule->complexity += 1;
998
999    float estRF = rule->estRF;
1000    if (optimismReduction == 0)
1001      bestQuality = evaluateRuleM(rule,examples,weightID,targetClass,apriori,rLength+1,aprioriProb);
1002    else if (optimismReduction == 1)
1003      bestQuality = evaluateRulePessimistic(rule,examples,weightID,targetClass,apriori,rLength+1,aprioriProb);
1004    else if (optimismReduction == 2)
1005      bestQuality = evaluateRuleEVC(rule,examples,weightID,targetClass,apriori,rLength+1,aprioriProb);
1006
1007    rule->estRF = estRF;
1008    rule->classDistribution = oldRuleDist;
1009    rule->chi = rulesTrueChi;
1010    rule->complexity -= 1;
1011    rule->distP = rulesDistP;
1012
1013    if (bestQuality <= quality)
1014      futureQuality = -1.0;
1015    else if (bestRule && bestQuality <= bestRule->quality)
1016      futureQuality = -1.0;
1017    else {
1018      futureQuality = 0.0; 
1019      PEITERATE(ei, rule->examples) {
1020        if ((*ei).getClass().intV != targetClass)
1021          continue;
1022        if (bestQuality <= (*ei)[probVar].floatV) {
1023          continue;
1024        }
1025      /*  float x = ((*ei)[probVar].floatV-quality); //*rule->classDistribution->abs;
1026        if ((*ei)[probVar].floatV > quality)
1027          x *= (1.0-quality)/(bestQuality-quality);
1028        x /= sqrt(quality*(1.0-quality)); // rule->classDistribution->abs*
1029        futureQuality += log(1.0-max(1e-12,1.0-2*zprob(x))); */
1030        float x;
1031        if ((*ei)[probVar].floatV > quality)
1032        {
1033            x = 1.0 - ((*ei)[probVar].floatV-quality) / (bestQuality-quality);
1034        }
1035        else
1036        {
1037            x = 1.0 + quality - (*ei)[probVar].floatV;
1038        }
1039        futureQuality += x;
1040      }
1041//      futureQuality = 1.0 - exp(futureQuality);
1042      futureQuality /= rule->classDistribution->atint(targetClass);
1043    }
1044  }
1045
1046  // store best rule as best rule and return expected quality of this rule
1047  rule->quality = quality;
1048  if (improved >= min_improved &&
1049      improved/rule->classDistribution->atint(targetClass) > min_improved_perc*0.01 &&
1050      quality > (aprioriProb + 1e-3) &&
1051      (!bestRule || (quality>bestRule->quality+1e-3)) &&
1052      (!validator || validator->call(rule, examples, weightID, targetClass, apriori))) {
1053
1054      TRule *pbestRule = new TRule(rule.getReference(), true);
1055      PRule wpbestRule = pbestRule;
1056      // check if rule is significant enough
1057      bool ruleGoodEnough = true;
1058
1059      if (ruleAlpha < 1.0)
1060        ruleGoodEnough = ruleGoodEnough & ((rule->chi > 0.0) && (chisqprob(rule->chi, 1.0f) <= ruleAlpha));
1061      if (ruleGoodEnough && attributeAlpha < 1.0)
1062        ruleGoodEnough = ruleGoodEnough & ruleAttSignificant(rule, examples, weightID, targetClass, apriori, aprioriProb);
1063      if (ruleGoodEnough)
1064      {
1065        bestRule = wpbestRule;
1066        bestRule->quality = quality;
1067        futureQuality += 1.0;
1068      }
1069  }
1070  return futureQuality;
1071}
1072
1073bool worstRule(const PRule &r1, const PRule &r2)
1074{ return    (r1->quality > r2->quality)
1075          || (r1->quality==r2->quality
1076          && r1->complexity < r2->complexity);
1077}
1078/*         || (r1->quality==r2->quality)
1079            && (   (r1->complexity < r2->complexity)
1080                || (r1->complexity == r2->complexity)
1081                   && ((int(r1.getUnwrappedPtr()) ^ int(r2.getUnwrappedPtr())) & 16) != 0
1082               ); }  */
1083
1084bool inRules(PRuleList rules, PRule rule)
1085{
1086  TRuleList::const_iterator ri(rules->begin()), re(rules->end());
1087  PExampleGenerator rulegen = rule->examples;
1088  for (; ri!=re; ri++) {
1089    PExampleGenerator rigen = (*ri)->examples;
1090    if (rigen->numberOfExamples() == rulegen->numberOfExamples()) {
1091      TExampleIterator rei(rulegen->begin()), ree(rulegen->end());
1092      TExampleIterator riei(rigen->begin()), riee(rigen->end());
1093      for (; rei != ree && !(*rei).compare(*riei); ++rei, ++riei) {
1094      }
1095        if (rei == ree)
1096          return true;
1097    }
1098  }
1099  return false;
1100}
1101
1102TRuleBeamFilter_Width::TRuleBeamFilter_Width(const int &w)
1103: width(w)
1104{}
1105
1106
1107void TRuleBeamFilter_Width::operator()(PRuleList &rules, PExampleTable, const int &)
1108{
1109  if (rules->size() > width) {
1110    sort(rules->begin(), rules->end(), worstRule);
1111
1112    TRuleList *filteredRules = mlnew TRuleList;
1113    PRuleList wFilteredRules = filteredRules;
1114
1115    int nRules = 0;
1116    TRuleList::const_iterator ri(rules->begin()), re(rules->end());
1117    while (nRules < width && ri != re) {
1118      if (!inRules(wFilteredRules,*ri)) {
1119        wFilteredRules->push_back(*ri);
1120        nRules++;
1121      }
1122      ri++;
1123    }
1124    rules =  wFilteredRules;
1125  }
1126}
1127
1128
1129inline void _selectBestRule(PRule &rule, PRule &bestRule, int &wins, TRandomGenerator &rgen)
1130{
1131  if ((rule->quality > bestRule->quality) || (rule->complexity < bestRule->complexity)) {
1132    bestRule = rule;
1133    wins = 1;
1134  }
1135  else if ((rule->complexity == bestRule->complexity) && rgen.randbool(++wins))
1136    bestRule = rule;
1137}
1138
1139
1140
1141PRuleList TRuleBeamInitializer_Default::operator()(PExampleTable data, const int &weightID, const int &targetClass, PRuleList baseRules, PRuleEvaluator evaluator, PDistribution apriori, PRule &bestRule)
1142{
1143  checkProperty(evaluator);
1144
1145  TRuleList *ruleList = mlnew TRuleList();
1146  PRuleList wruleList = ruleList;
1147
1148  TRandomGenerator rgen(data->numberOfExamples());
1149  int wins;
1150
1151  if (baseRules && baseRules->size())
1152    PITERATE(TRuleList, ri, baseRules) {
1153      TRule *newRule = mlnew TRule((*ri).getReference(), true);
1154      PRule wNewRule = newRule;
1155      ruleList->push_back(wNewRule);
1156      if (!newRule->examples)
1157      {
1158        newRule->filterAndStore(data,weightID,targetClass);
1159        newRule->baseDist = newRule->classDistribution;
1160      }
1161      newRule->quality = evaluator->call(wNewRule, data, weightID, targetClass, apriori);
1162      if (!bestRule || (newRule->quality > bestRule->quality)) {
1163        bestRule = wNewRule;
1164        wins = 1;
1165      }
1166      else
1167        if (newRule->quality == bestRule->quality)
1168          _selectBestRule(wNewRule, bestRule, wins, rgen);
1169    }
1170
1171  else {
1172     TRule *ubestRule = mlnew TRule();
1173     bestRule = ubestRule;
1174     ruleList->push_back(bestRule);
1175     ubestRule->filter = new TFilter_values();
1176     ubestRule->filter->domain = data->domain;
1177     ubestRule->filterAndStore(data, weightID,targetClass);
1178     ubestRule->baseDist = ubestRule->classDistribution;
1179     ubestRule->complexity = 0;
1180  }
1181
1182  return wruleList;
1183}
1184
1185
1186PRuleList TRuleBeamRefiner_Selector::operator()(PRule wrule, PExampleTable data, const int &weightID, const int &targetClass)
1187{
1188  if (!discretization) {
1189    discretization = mlnew TEntropyDiscretization();
1190    dynamic_cast<TEntropyDiscretization *>(discretization.getUnwrappedPtr())->forceAttribute = true;
1191    dynamic_cast<TEntropyDiscretization *>(discretization.getUnwrappedPtr())->maxNumberOfIntervals = 5;
1192  }
1193
1194  TRule &rule = wrule.getReference();
1195  TFilter_values *filter = wrule->filter.AS(TFilter_values);
1196  if (!filter)
1197    raiseError("a filter of type 'Filter_values' expected");
1198
1199  TRuleList *ruleList = mlnew TRuleList;
1200  PRuleList wRuleList = ruleList;
1201
1202  TDomainDistributions ddist(wrule->examples, wrule->weightID);
1203
1204  const TVarList &attributes = rule.examples->domain->attributes.getReference();
1205
1206  vector<bool> used(attributes.size(), false);
1207  PITERATE(TValueFilterList, vfi, filter->conditions)
1208    used[(*vfi)->position] = true;
1209
1210  vector<bool>::const_iterator ui(used.begin());
1211  TDomainDistributions::const_iterator di(ddist.begin());
1212  TVarList::const_iterator vi(attributes.begin()), ve(attributes.end());
1213  int pos = 0;
1214  for(; vi != ve; vi++, ui++, pos++, di++) {
1215    if ((*vi)->varType == TValue::INTVAR) {
1216      if (!*ui) {
1217        vector<float>::const_iterator idi((*di).AS(TDiscDistribution)->begin());
1218        for(int v = 0, e = (*vi)->noOfValues(); v != e; v++)
1219          if (*idi>0) {
1220            TRule *newRule = mlnew TRule(rule, false);
1221            ruleList->push_back(newRule);
1222            newRule->complexity++;
1223
1224            filter = newRule->filter.AS(TFilter_values);
1225
1226            TValueFilter_discrete *newCondition = mlnew TValueFilter_discrete(pos, *vi, 0);
1227            filter->conditions->push_back(newCondition);
1228
1229            TValue value = TValue(v);
1230            newCondition->values->push_back(value);
1231            newRule->filterAndStore(rule.examples, rule.weightID,targetClass);
1232            newRule->parentRule = wrule;
1233          }
1234      }
1235    }
1236
1237    else if (((*vi)->varType == TValue::FLOATVAR)) {
1238
1239        if (discretization) {
1240        PVariable discretized;
1241        try {
1242          discretized = discretization->call(rule.examples, *vi, weightID);
1243        } catch(...) {
1244          continue;
1245        }
1246        TClassifierFromVar *cfv = discretized->getValueFrom.AS(TClassifierFromVar);
1247        TDiscretizer *discretizer = cfv ? cfv->transformer.AS(TDiscretizer) : NULL;
1248        if (!discretizer)
1249          raiseError("invalid or unrecognized discretizer");
1250
1251        vector<float> cutoffs;
1252        discretizer->getCutoffs(cutoffs);
1253        if (cutoffs.size()) {
1254          TRule *newRule;
1255          newRule = mlnew TRule(rule, false);
1256          PRule wnewRule = newRule;
1257          newRule->complexity++;
1258          newRule->parentRule = wrule;
1259
1260          newRule->filter.AS(TFilter_values)->conditions->push_back(mlnew TValueFilter_continuous(pos,  TValueFilter_continuous::LessEqual,     cutoffs.front(), 0, 0));
1261          newRule->filterAndStore(rule.examples, rule.weightID,targetClass);
1262          if (wrule->classDistribution->cases > wnewRule->classDistribution->cases)
1263            ruleList->push_back(newRule);
1264
1265          for(vector<float>::const_iterator ci(cutoffs.begin()), ce(cutoffs.end()-1); ci != ce; ci++) {
1266            newRule = mlnew TRule(rule, false);
1267            wnewRule = newRule;
1268            newRule->complexity++;
1269            newRule->parentRule = wrule;
1270            filter = newRule->filter.AS(TFilter_values);
1271            filter->conditions->push_back(mlnew TValueFilter_continuous(pos,  TValueFilter_continuous::Greater, *ci, 0, 0));
1272            newRule->filterAndStore(rule.examples, rule.weightID,targetClass);
1273            if (wrule->classDistribution->cases > wnewRule->classDistribution->cases)
1274              ruleList->push_back(newRule);
1275
1276            newRule = mlnew TRule(rule, false);
1277            wnewRule = newRule;
1278            newRule->complexity++;
1279            newRule->parentRule = wrule;
1280            filter = newRule->filter.AS(TFilter_values);
1281            filter->conditions->push_back(mlnew TValueFilter_continuous(pos,  TValueFilter_continuous::LessEqual, *(ci+1), 0, 0));
1282            newRule->filterAndStore(rule.examples, rule.weightID,targetClass);
1283            if (wrule->classDistribution->cases > wnewRule->classDistribution->cases)
1284              ruleList->push_back(newRule);
1285          }
1286
1287          newRule = mlnew TRule(rule, false);
1288          ruleList->push_back(newRule);
1289          newRule->complexity++;
1290
1291          newRule->filter.AS(TFilter_values)->conditions->push_back(mlnew TValueFilter_continuous(pos,  TValueFilter_continuous::Greater, cutoffs.back(), 0, 0));
1292          newRule->filterAndStore(rule.examples, rule.weightID,targetClass);
1293          newRule->parentRule = wrule;
1294        }
1295      }
1296      else
1297        raiseWarning("discretizer not given, continuous attributes will be skipped");
1298    }
1299  }
1300
1301  if (!discretization)
1302    discretization = PDiscretization();
1303  return wRuleList;
1304}
1305
1306
1307PRuleList TRuleBeamCandidateSelector_TakeAll::operator()(PRuleList &existingRules, PExampleTable, const int &)
1308{
1309  PRuleList candidates = mlnew TRuleList(existingRules.getReference());
1310//  existingRules->clear();
1311  existingRules->erase(existingRules->begin(), existingRules->end());
1312  return candidates;
1313}
1314
1315
1316PRule TRuleBeamFinder::operator()(PExampleTable data, const int &weightID, const int &targetClass, PRuleList baseRules)
1317{
1318  // set default values if value not set
1319  bool tempInitializer = !initializer;
1320  if (tempInitializer)
1321    initializer = mlnew TRuleBeamInitializer_Default;
1322  bool tempCandidateSelector = !candidateSelector;
1323  if (tempCandidateSelector)
1324    candidateSelector = mlnew TRuleBeamCandidateSelector_TakeAll;
1325  bool tempRefiner = !refiner;
1326  if (tempRefiner)
1327    refiner = mlnew TRuleBeamRefiner_Selector;
1328/*  bool tempValidator = !validator;
1329  if (tempValidator)
1330    validator = mlnew TRuleValidator_LRS((float)0.01);
1331  bool tempRuleStoppingValidator = !ruleStoppingValidator;
1332  if (tempRuleStoppingValidator)
1333    ruleStoppingValidator = mlnew TRuleValidator_LRS((float)0.05); */
1334  bool tempEvaluator = !evaluator;
1335  if (tempEvaluator)
1336    evaluator = mlnew TRuleEvaluator_Entropy;
1337  bool tempRuleFilter = !ruleFilter;
1338  if (tempRuleFilter)
1339    ruleFilter = mlnew TRuleBeamFilter_Width;
1340
1341  // create an entropy evaluator to detect rules with pure distributions
1342  PRuleEvaluator entropy_evaluator = mlnew TRuleEvaluator_Entropy;
1343
1344  checkProperty(initializer);
1345  checkProperty(candidateSelector);
1346  checkProperty(refiner);
1347  checkProperty(evaluator);
1348  checkProperty(ruleFilter);
1349  PDistribution apriori = getClassDistribution(data, weightID);
1350
1351  TRandomGenerator rgen(data->numberOfExamples());
1352  int wins = 1;
1353
1354  PRule bestRule;
1355  PRuleList ruleList = initializer->call(data, weightID, targetClass, baseRules, evaluator, apriori, bestRule);
1356
1357  PITERATE(TRuleList, ri, ruleList) {
1358    if (!(*ri)->examples)
1359      (*ri)->filterAndStore(data, weightID,targetClass);
1360    if ((*ri)->quality == ILLEGAL_FLOAT)
1361      (*ri)->quality = evaluator->call(*ri, data, weightID, targetClass, apriori);
1362  }
1363
1364  if (!bestRule->examples)
1365    bestRule->filterAndStore(data, weightID,targetClass);
1366  if (bestRule->quality == ILLEGAL_FLOAT)
1367    bestRule->quality = evaluator->call(bestRule, data, weightID, targetClass, apriori);
1368
1369  int bestRuleLength = 0;
1370  while(ruleList->size()) {
1371    PRuleList candidateRules = candidateSelector->call(ruleList, data, weightID);
1372    PITERATE(TRuleList, ri, candidateRules) {
1373      PRuleList newRules = refiner->call(*ri, data, weightID, targetClass);
1374      PITERATE(TRuleList, ni, newRules) {
1375        // evaluate rule
1376        (*ni)->quality = evaluator->call(*ni, data, weightID, targetClass, apriori);
1377        if ((*ni)->quality > bestRule->quality && (!validator || validator->call(*ni, data, weightID, targetClass, apriori)))
1378          _selectBestRule(*ni, bestRule, wins, rgen);
1379        // if rule covers only one class or has less than 2 examples, do not add it for further specialization
1380        if (entropy_evaluator->call(*ni, data, weightID, targetClass, apriori) > -0.01 ||
1381            (*ni)->classDistribution->abs < 2)
1382            continue;
1383        // is the rule a good candidate for further specialization?
1384        if (!ruleStoppingValidator || ruleStoppingValidator->call(*ni, (*ri)->examples, weightID, targetClass, (*ri)->classDistribution)) {
1385          ruleList->push_back(*ni);
1386        }
1387      }
1388    }
1389    ruleFilter->call(ruleList,data,weightID);
1390  }
1391
1392  // set empty values if value was not set (used default)
1393  if (tempInitializer)
1394    initializer = PRuleBeamInitializer();
1395  if (tempCandidateSelector)
1396    candidateSelector = PRuleBeamCandidateSelector();
1397  if (tempRefiner)
1398    refiner = PRuleBeamRefiner();
1399/*  if (tempValidator)
1400    validator = PRuleValidator();
1401  if (tempRuleStoppingValidator)
1402    ruleStoppingValidator = PRuleValidator();  */
1403  if (tempEvaluator)
1404    evaluator = PRuleEvaluator();
1405  if (tempRuleFilter)
1406    ruleFilter = PRuleBeamFilter();
1407
1408  return bestRule;
1409}
1410
1411
1412TRuleLearner::TRuleLearner(bool se, int tc, PRuleList rl)
1413: storeExamples(se),
1414  targetClass(tc),
1415  baseRules(rl)
1416{}
1417
1418
1419PClassifier TRuleLearner::operator()(PExampleGenerator gen, const int &weightID)
1420{
1421  return this->call(gen,weightID,targetClass,baseRules);
1422}
1423
1424PClassifier TRuleLearner::operator()(PExampleGenerator gen, const int &weightID, const int &targetClass, PRuleList baseRules)
1425{
1426  // Initialize default values if values not set
1427  bool tempDataStopping = !dataStopping && !ruleStopping;
1428  if (tempDataStopping)
1429    dataStopping = mlnew TRuleDataStoppingCriteria_NoPositives;
1430
1431  bool tempRuleFinder = !ruleFinder;
1432  if (tempRuleFinder)
1433    ruleFinder = mlnew TRuleBeamFinder;
1434
1435  bool tempCoverAndRemove = !coverAndRemove;
1436  if (tempCoverAndRemove)
1437    coverAndRemove = mlnew TRuleCovererAndRemover_Default;
1438
1439  checkProperty(ruleFinder);
1440  checkProperty(coverAndRemove);
1441
1442  TExampleTable *data = mlnew TExampleTable(gen);
1443  PExampleTable wdata = data;
1444
1445  if (!dataStopping && !ruleStopping)
1446    raiseError("no stopping criteria; set 'dataStopping' and/or 'ruleStopping'");
1447
1448  TRuleList *ruleList = mlnew TRuleList;
1449  PRuleList wruleList = ruleList;
1450
1451  int currWeightID = weightID;
1452
1453  float beginwe=0.0, currentwe;
1454  if (progressCallback) {
1455    if (targetClass==-1)
1456      beginwe = wdata->weightOfExamples(weightID);
1457    else {
1458      PDistribution classDist = getClassDistribution(wdata, weightID);
1459      TDiscDistribution *ddist = classDist.AS(TDiscDistribution);
1460      beginwe = ddist->atint(targetClass);
1461    }
1462    progressCallback->call(0.0);
1463  }
1464
1465  while (!dataStopping || !dataStopping->call(wdata, currWeightID, targetClass)) {
1466    PRule rule = ruleFinder->call(wdata, currWeightID, targetClass, baseRules);
1467    if (!rule)
1468      raiseError("'ruleFinder' didn't return a rule");
1469
1470    if (ruleStopping && ruleStopping->call(ruleList, rule, wdata, currWeightID))
1471      break;
1472
1473    wdata = coverAndRemove->call(rule, wdata, currWeightID, currWeightID, targetClass);
1474    ruleList->push_back(rule);
1475
1476    if (progressCallback) {
1477      if (targetClass==-1)
1478        currentwe = wdata->weightOfExamples(weightID);
1479      else {
1480        PDistribution classDist = getClassDistribution(wdata, currWeightID);
1481        TDiscDistribution *ddist = classDist.AS(TDiscDistribution);
1482        currentwe = ddist->atint(targetClass);
1483      }
1484      progressCallback->call(1-currentwe/beginwe);
1485    }
1486  }
1487  if (progressCallback)
1488    progressCallback->call(1.0);
1489
1490
1491  // Restore values
1492  if (tempDataStopping)
1493    dataStopping = PRuleDataStoppingCriteria();
1494  if (tempRuleFinder)
1495    ruleFinder = PRuleFinder();
1496  if (tempCoverAndRemove)
1497    coverAndRemove = PRuleCovererAndRemover();
1498
1499  PRuleClassifierConstructor clConstructor =
1500    classifierConstructor ? classifierConstructor :
1501    PRuleClassifierConstructor(mlnew TRuleClassifierConstructor_firstRule());
1502  return clConstructor->call(ruleList, gen, weightID);
1503};
1504
1505
1506bool TRuleDataStoppingCriteria_NoPositives::operator()(PExampleTable data, const int &weightID, const int &targetClass) const
1507{
1508  PDistribution classDist = getClassDistribution(data, weightID);
1509  TDiscDistribution *ddist = classDist.AS(TDiscDistribution);
1510
1511  return (targetClass >= 0 ? ddist->atint(targetClass) : ddist->abs) == 0.0;
1512}
1513
1514bool TRuleStoppingCriteria_NegativeDistribution::operator()(PRuleList ruleList, PRule rule, PExampleTable data, const int &weightID) const
1515{
1516  if (rule && rule->classifier)
1517  {
1518    PDistribution aprioriDist = getClassDistribution(data, weightID);
1519    TDiscDistribution *apriori = aprioriDist.AS(TDiscDistribution);
1520
1521    const TDefaultClassifier *clsf = rule->classifier.AS(TDefaultClassifier);
1522    if (!clsf)
1523      return false;
1524    const TDiscDistribution *dist = dynamic_cast<const TDiscDistribution *>(clsf->defaultDistribution.getUnwrappedPtr());
1525    const int classVal = clsf->defaultVal.intV;
1526    if (classVal<0 || classVal>=dist->size())
1527      return false;
1528    float acc = dist->atint(clsf->defaultVal.intV)/dist->abs;
1529    float accApriori = apriori->atint(clsf->defaultVal.intV)/apriori->abs;
1530    if (accApriori>acc)
1531      return true;
1532  }
1533  return false;
1534}
1535
1536
1537PExampleTable TRuleCovererAndRemover_Default::operator()(PRule rule, PExampleTable data, const int &weightID, int &newWeight, const int &targetClass) const
1538{
1539  TExampleTable *table = mlnew TExampleTable(data, 1);
1540  PExampleGenerator wtable = table;
1541
1542  TFilter &filter = rule->filter.getReference();
1543
1544  if (targetClass < 0)
1545  {
1546    PEITERATE(ei, data)
1547      if (!filter(*ei))
1548        table->addExample(*ei);
1549  }
1550  else
1551    PEITERATE(ei, data)
1552      if (!filter(*ei) || (*ei).getClass().intV!=targetClass)
1553        table->addExample(*ei);
1554
1555
1556  newWeight = weightID;
1557  return wtable;
1558}
1559
1560// classifiers
1561PRuleClassifier TRuleClassifierConstructor_firstRule::operator ()(PRuleList rules, PExampleTable table, const int &weightID)
1562{
1563  return mlnew TRuleClassifier_firstRule(rules, table, weightID);
1564}
1565
1566
1567TRuleClassifier::TRuleClassifier(PRuleList arules, PExampleTable anexamples, const int &aweightID)
1568: rules(arules),
1569  examples(anexamples),
1570  weightID(aweightID),
1571  TClassifier(anexamples->domain->classVar,true)
1572{}
1573
1574TRuleClassifier::TRuleClassifier()
1575: TClassifier(true)
1576{}
1577
1578
1579TRuleClassifier_firstRule::TRuleClassifier_firstRule(PRuleList arules, PExampleTable anexamples, const int &aweightID)
1580: TRuleClassifier(arules, anexamples, aweightID)
1581{
1582  prior = getClassDistribution(examples, weightID);
1583}
1584
1585TRuleClassifier_firstRule::TRuleClassifier_firstRule()
1586: TRuleClassifier()
1587{}
1588
1589PDistribution TRuleClassifier_firstRule::classDistribution(const TExample &ex)
1590{
1591  checkProperty(rules);
1592  checkProperty(prior);
1593
1594  PITERATE(TRuleList, ri, rules) {
1595    if ((*ri)->call(ex))
1596      return (*ri)->classDistribution;
1597  }
1598  return prior;
1599}
1600
1601void copyTable(float **dest, float **source, int nx, int ny) {
1602  for (int i=0; i<nx; i++)
1603    memcpy(dest[i],source[i],ny*sizeof(float));
1604}
1605
1606//==============================================================================
1607// return 1 if system not solving
1608// nDim - system dimension
1609// pfMatr - matrix with coefficients
1610// pfVect - vector with free members
1611// pfSolution - vector with system solution
1612// pfMatr becames trianglular after function call
1613// pfVect changes after function call
1614//
1615// Developer: Henry Guennadi Levkin
1616//
1617//==============================================================================
1618int LinearEquationsSolving(int nDim, double* pfMatr, double* pfVect, double* pfSolution)
1619{
1620  double fMaxElem;
1621  double fAcc;
1622
1623  int i, j, k, m;
1624
1625
1626  for(k=0; k<(nDim-1); k++) // base row of matrix
1627  {
1628    // search of line with max element
1629    fMaxElem = fabs( pfMatr[k*nDim + k] );
1630    m = k;
1631    for(i=k+1; i<nDim; i++)
1632    {
1633      if(fMaxElem < fabs(pfMatr[i*nDim + k]) )
1634      {
1635        fMaxElem = pfMatr[i*nDim + k];
1636        m = i;
1637      }
1638    }
1639
1640    // permutation of base line (index k) and max element line(index m)
1641    if(m != k)
1642    {
1643      for(i=k; i<nDim; i++)
1644      {
1645        fAcc               = pfMatr[k*nDim + i];
1646        pfMatr[k*nDim + i] = pfMatr[m*nDim + i];
1647        pfMatr[m*nDim + i] = fAcc;
1648      }
1649      fAcc = pfVect[k];
1650      pfVect[k] = pfVect[m];
1651      pfVect[m] = fAcc;
1652    }
1653
1654    if( pfMatr[k*nDim + k] == 0.) return 1; // needs improvement !!!
1655
1656    // triangulation of matrix with coefficients
1657    for(j=(k+1); j<nDim; j++) // current row of matrix
1658    {
1659      fAcc = - pfMatr[j*nDim + k] / pfMatr[k*nDim + k];
1660      for(i=k; i<nDim; i++)
1661      {
1662        pfMatr[j*nDim + i] = pfMatr[j*nDim + i] + fAcc*pfMatr[k*nDim + i];
1663      }
1664      pfVect[j] = pfVect[j] + fAcc*pfVect[k]; // free member recalculation
1665    }
1666  }
1667
1668  for(k=(nDim-1); k>=0; k--)
1669  {
1670    pfSolution[k] = pfVect[k];
1671    for(i=(k+1); i<nDim; i++)
1672    {
1673      pfSolution[k] -= (pfMatr[k*nDim + i]*pfSolution[i]);
1674    }
1675    pfSolution[k] = pfSolution[k] / pfMatr[k*nDim + k];
1676  }
1677
1678  return 0;
1679}
1680
1681// extracts class value index from target in rule
1682int getClassIndex(PRule r) {
1683  const TDefaultClassifier &cl = dynamic_cast<const TDefaultClassifier &>(r->classifier.getReference());
1684  return cl.defaultVal.intV;
1685}
1686
1687// constructor 1
1688TLogitClassifierState::TLogitClassifierState(PRuleList arules, PExampleTable anexamples, const int &aweightID)
1689: rules(arules),
1690  examples(anexamples),
1691  weightID(aweightID)
1692{
1693  // initialize f, p
1694  f = new float *[examples->domain->classVar->noOfValues()-1];
1695  p = new float *[examples->domain->classVar->noOfValues()];
1696  int i;
1697  for (i=0; i<examples->domain->classVar->noOfValues()-1; i++) {
1698      f[i] = new float[examples->numberOfExamples()];
1699      p[i] = new float[examples->numberOfExamples()];
1700  }
1701  p[examples->domain->classVar->noOfValues()-1] = new float[examples->numberOfExamples()];
1702
1703  betas = new float[rules->size()];
1704  priorBetas = new float[examples->domain->classVar->noOfValues()];
1705  isExampleFixed = new bool[examples->numberOfExamples()];
1706}
1707
1708// constructor 2
1709TLogitClassifierState::TLogitClassifierState(PRuleList arules, const PDistributionList &probList, PExampleTable anexamples, const int &aweightID)
1710: rules(arules),
1711  examples(anexamples),
1712  weightID(aweightID)
1713{
1714  // initialize f, p
1715  f = new float *[examples->domain->classVar->noOfValues()-1];
1716  p = new float *[examples->domain->classVar->noOfValues()];
1717  int i, j;
1718  for (i=0; i<examples->domain->classVar->noOfValues()-1; i++) {
1719    f[i] = new float[examples->numberOfExamples()];
1720    p[i] = new float[examples->numberOfExamples()];
1721    for (j=0; j<examples->numberOfExamples(); j++) {
1722        f[i][j] = 0.0;
1723        p[i][j] = 1.0/examples->domain->classVar->noOfValues();
1724    }
1725  }
1726  p[examples->domain->classVar->noOfValues()-1] = new float[examples->numberOfExamples()];
1727  for (j=0; j<examples->numberOfExamples(); j++)
1728      p[examples->domain->classVar->noOfValues()-1][j] = 1.0/examples->domain->classVar->noOfValues();
1729
1730   // if initial example probability is given, update F and P
1731  if (probList) {
1732    double *matrix = new double [sqr(examples->domain->classVar->noOfValues()-1)];
1733    double *fVals = new double [examples->domain->classVar->noOfValues()-1];
1734    double *results = new double [examples->domain->classVar->noOfValues()-1];
1735    for (i=0; i<probList->size(); i++) {
1736      int k1, k2;
1737      TDistribution *dist = mlnew TDiscDistribution(probList->at(i)->variable);
1738      PDistribution wdist = dist;
1739      // Prepare and compute expected f - values (a linear equation)
1740      for (k1=0; k1<examples->domain->classVar->noOfValues(); k1++) {
1741        if (probList->at(i)->atint(k1) >= 1.0-1e-4)
1742          wdist->setint(k1,(float)(1.0-1e-4));
1743        else if (probList->at(i)->atint(k1) <= 1e-4)
1744          wdist->setint(k1,(float)(1e-4));
1745        else
1746          wdist->setint(k1,probList->at(i)->atint(k1));
1747      }
1748      wdist->normalize();
1749      for (k1=0; k1<examples->domain->classVar->noOfValues()-1; k1++) {
1750        fVals[k1] = -wdist->atint(k1);
1751        for (k2=0; k2<examples->domain->classVar->noOfValues()-1; k2++) {
1752          if (k1==k2)
1753            matrix[k1*(examples->domain->classVar->noOfValues()-1)+k2] = wdist->atint(k1)-1;
1754          else
1755            matrix[k1*(examples->domain->classVar->noOfValues()-1)+k2] = wdist->atint(k1);
1756        }
1757      }
1758      LinearEquationsSolving(examples->domain->classVar->noOfValues()-1, matrix, fVals, results);
1759      // store values
1760      for (k1=0; k1<examples->domain->classVar->noOfValues()-1; k1++)
1761        f[k1][i] = results[k1]>0.0 ? log(results[k1]) : -10.0;
1762      for (k1=0; k1<examples->domain->classVar->noOfValues(); k1++)
1763        p[k1][i] = wdist->atint(k1);
1764    }
1765    delete [] matrix;
1766    delete [] fVals;
1767    delete [] results;
1768  }
1769
1770  // compute rule indices
1771  i=0;
1772  ruleIndices = mlnew PIntList[rules->size()];
1773  {
1774    PITERATE(TRuleList, ri, rules) {
1775      TIntList *ruleIndicesnw = mlnew TIntList();
1776      ruleIndices[i] = ruleIndicesnw;
1777      j=0;
1778      PEITERATE(ei, examples) {
1779        if ((*ri)->call(*ei))
1780          ruleIndices[i]->push_back(j);
1781        j++;
1782      }
1783      i++;
1784    }
1785  }
1786
1787  // set initial values of betas
1788  betas = new float[rules->size()];
1789  for (i=0; i<rules->size(); i++)
1790    betas[i] = (float)0.0;
1791
1792  // Add default rules
1793  priorBetas = new float[examples->domain->classVar->noOfValues()];
1794  for (i=0; i<examples->domain->classVar->noOfValues(); i++)
1795    priorBetas[i] = 0.0;
1796
1797  // computer best rule covering
1798  PDistribution apriori = getClassDistribution(examples, weightID);
1799  isExampleFixed = new bool[examples->numberOfExamples()];
1800  for (j=0; j<examples->numberOfExamples(); j++)
1801    isExampleFixed[j] = false;
1802
1803  // priorProb and avgProb
1804  TFloatList *npriorProb = mlnew TFloatList();
1805  avgPriorProb = npriorProb;
1806  TFloatList *navgProb = mlnew TFloatList();
1807  avgProb = navgProb;
1808  TIntList *pprefixRules = mlnew TIntList();
1809  prefixRules = pprefixRules;
1810  computeAvgProbs();
1811  computePriorProbs();
1812}
1813
1814TLogitClassifierState::~TLogitClassifierState()
1815{
1816  int i;
1817  for (i=0; i<examples->domain->classVar->noOfValues()-1; i++)
1818    delete [] f[i];
1819  delete [] f;
1820
1821  for (i=0; i<examples->domain->classVar->noOfValues(); i++)
1822    delete [] p[i];
1823  delete [] p;
1824  delete [] betas;
1825  delete [] priorBetas;
1826  delete [] ruleIndices;
1827  delete [] isExampleFixed;
1828}
1829
1830void TLogitClassifierState::computeAvgProbs()
1831{
1832  // compute new rule avgProbs
1833  avgProb->clear();
1834  int classInd = 0;
1835
1836  float newAvgProb;
1837  for (int ri = 0; ri<rules->size(); ri++) {
1838    newAvgProb = 0.0;
1839    classInd = getClassIndex(rules->at(ri));
1840    PITERATE(TIntList, ind, ruleIndices[ri])
1841      newAvgProb += p[classInd][*ind];
1842    avgProb->push_back(newAvgProb/ruleIndices[ri]->size());
1843  }
1844}
1845
1846// compute new prior probs
1847void TLogitClassifierState::computePriorProbs()
1848{
1849  avgPriorProb->clear();
1850  for (int pi=0; pi<examples->domain->classVar->noOfValues(); pi++) {
1851    float newPriorProb = 0.0;
1852    for (int ei=0; ei<examples->numberOfExamples(); ei++) {
1853      newPriorProb += p[pi][ei];
1854    }
1855    avgPriorProb->push_back(newPriorProb/examples->numberOfExamples());
1856  }
1857}
1858
1859void TLogitClassifierState::copyTo(PLogitClassifierState & wstate)
1860{
1861  if (!wstate) {
1862    TLogitClassifierState *state = mlnew TLogitClassifierState(rules, examples, weightID);
1863    wstate = state;
1864    wstate->ruleIndices = mlnew PIntList[rules->size()];
1865    int i;
1866    for (i=0; i<rules->size(); i++) {
1867      TIntList * tIndices = mlnew TIntList(ruleIndices[i].getReference());
1868      wstate->ruleIndices[i] = tIndices;
1869    }
1870  }
1871
1872  wstate->eval = eval;
1873  copyTable(wstate->f, f, examples->domain->classVar->noOfValues()-1, examples->numberOfExamples());
1874  copyTable(wstate->p, p, examples->domain->classVar->noOfValues(), examples->numberOfExamples());
1875  memcpy(wstate->betas,betas,sizeof(float)*rules->size());
1876  memcpy(wstate->priorBetas,priorBetas,sizeof(float)*(examples->domain->classVar->noOfValues()-1));
1877  memcpy(wstate->isExampleFixed, isExampleFixed, sizeof(bool)*examples->numberOfExamples());
1878
1879  TFloatList *pavgProb = mlnew TFloatList(avgProb.getReference());
1880  TFloatList *pavgPriorProb = mlnew TFloatList(avgPriorProb.getReference());
1881  TIntList *pprefixRules = mlnew TIntList(prefixRules.getReference());
1882  wstate->avgProb = pavgProb;
1883  wstate->avgPriorProb = pavgPriorProb;
1884  wstate->prefixRules = pprefixRules;
1885}
1886
1887void TLogitClassifierState::newBeta(int i, float b)
1888{
1889  // set new beta
1890  float diff = b-betas[i];
1891  betas[i] = b;
1892
1893
1894  // add differences to f
1895  int classIndex = getClassIndex(rules->at(i));
1896  PITERATE(TIntList, ind, ruleIndices[i])
1897    for (int fi=0; fi<examples->domain->classVar->noOfValues()-1; fi++)
1898      if (fi == classIndex)
1899        f[fi][*ind] += diff;
1900      else if (classIndex == examples->domain->classVar->noOfValues()-1)
1901        f[fi][*ind] -= diff;
1902
1903  // compute p
1904  computePs(i);
1905  computeAvgProbs();
1906  computePriorProbs();
1907}
1908
1909void TLogitClassifierState::newPriorBeta(int i, float b)
1910{
1911  // set new beta
1912  float diff = b-priorBetas[i];
1913  priorBetas[i] = b;
1914
1915  // add differences to f
1916  for (int ei=0; ei<examples->numberOfExamples(); ei++)
1917    for (int fi=0; fi<examples->domain->classVar->noOfValues()-1; fi++)
1918      if (fi == i)
1919        f[fi][ei] += diff;
1920  computePs(-1);
1921  computeAvgProbs();
1922  computePriorProbs();
1923}
1924
1925void TLogitClassifierState::updateExampleP(int ei)
1926{
1927/*  PITERATE(TIntList, ind, frontRules)
1928      if (rules->at(*ind)->call(examples->at(ei))) {
1929        p[getClassIndex(rules->at(*ind))][ei] = rules->at(*ind)->quality;
1930        for (int ci=0; ci<examples->domain->classVar->noOfValues(); ci++)
1931            if (ci!=getClassIndex(rules->at(*ind)))
1932                p[ci][ei] = (1.0-rules->at(*ind)->quality)/(examples->domain->classVar->noOfValues()-1);
1933        return;
1934      } */
1935  if (isExampleFixed[ei])
1936      return;
1937
1938  float sum = 1.0;
1939  int pi;
1940  for (pi=0; pi<examples->domain->classVar->noOfValues()-1; pi++) {
1941    if (f[pi][ei] > 10.0)
1942        p[pi][ei] = 22000.0;
1943    else
1944        p[pi][ei] = exp(f[pi][ei]);
1945    sum += p[pi][ei];
1946  }
1947  p[examples->domain->classVar->noOfValues()-1][ei] = 1.0;
1948  for (pi=0; pi<examples->domain->classVar->noOfValues(); pi+=1) {
1949    p[pi][ei] /= sum;
1950  }
1951}
1952
1953void TLogitClassifierState::computePs(int beta_i)
1954{
1955  if (beta_i<0)
1956    for (int ei=0; ei<examples->numberOfExamples(); ei++)
1957      updateExampleP(ei);
1958  else
1959    PITERATE(TIntList, ind, ruleIndices[beta_i])
1960      updateExampleP(*ind);
1961}
1962
1963void TLogitClassifierState::setFixed(int rule_i)
1964{
1965    PITERATE(TIntList, ind, ruleIndices[rule_i])
1966      isExampleFixed[*ind] = true;
1967}
1968
1969void TLogitClassifierState::updateFixedPs(int rule_i)
1970{
1971    PITERATE(TIntList, ind, ruleIndices[rule_i])
1972    {
1973        float bestQuality = 0.0;
1974        PITERATE(TIntList, fr, prefixRules) {
1975            if (rules->at(*fr)->call(examples->at(*ind)) && rules->at(*fr)->quality > bestQuality) {
1976                bestQuality = rules->at(*fr)->quality;
1977                p[getClassIndex(rules->at(*fr))][*ind] = rules->at(*fr)->quality;
1978                for (int ci=0; ci<examples->domain->classVar->noOfValues(); ci++)
1979                if (ci!=getClassIndex(rules->at(*fr)))
1980                    p[ci][*ind] = (1.0-rules->at(*fr)->quality)/(examples->domain->classVar->noOfValues()-1);
1981                break;
1982            }
1983        }
1984    }
1985}
1986
1987float TLogitClassifierState::getBrierScore()
1988{
1989  float brier = 0.0;
1990  for (int j=0; j<examples->numberOfExamples(); j++)
1991    brier += pow(1.0-p[examples->at(j).getClass().intV][j],2.0);
1992  return brier;
1993}
1994
1995float TLogitClassifierState::getAUC()
1996{
1997  float auc = 0.0, sum_ranks;
1998  int n1, n2;
1999  vector<double> probs, ranks;
2000//  for each class
2001  for (int cl=0; cl < examples->domain->classVar->noOfValues() - 1; cl++)
2002  {
2003//    probs = prediction probabilities of class
2004    probs.clear();
2005    for (int j=0; j<examples->numberOfExamples(); j++)
2006      probs.push_back(p[cl][j]);
2007
2008//    get ranks of class examples
2009//    n1 = elements of examples from class
2010//    n2 = elements of examples not from class
2011    rankdata(probs, ranks);
2012    n1=n2=0;
2013    sum_ranks=0.0;
2014    for (int j=0; j<examples->numberOfExamples(); j++)
2015    {
2016      if (examples->at(j).getClass().intV == cl)
2017      {
2018        n1++;
2019        sum_ranks+=ranks.at(j);
2020      }
2021      else n2++;
2022    }
2023//    auc += (sum(ranks) - n1*(n1+1)/2) / (n1*n2)
2024    auc += (sum_ranks - n1*(n1+1)/2) / (n1*n2);
2025  }
2026  return auc / (examples->domain->classVar->noOfValues() - 1);
2027}
2028
2029void TLogitClassifierState::setPrefixRule(int rule_i) //, int position)
2030{
2031    prefixRules->push_back(rule_i);
2032    setFixed(rule_i);
2033    updateFixedPs(rule_i);
2034    betas[rule_i] = 0.0;
2035    computeAvgProbs();
2036    computePriorProbs();
2037}
2038
2039TRuleClassifier_logit::TRuleClassifier_logit()
2040: TRuleClassifier()
2041{}
2042
2043TRuleClassifier_logit::TRuleClassifier_logit(PRuleList arules, const float &minSignificance, const float &minBeta, const float &penalty, PExampleTable anexamples, const int &aweightID, const PClassifier &classifier, const PDistributionList &probList, bool setPrefixRules, bool optimizeBetasFlag)
2044: TRuleClassifier(arules, anexamples, aweightID),
2045  minSignificance(minSignificance),
2046  priorClassifier(classifier),
2047  setPrefixRules(setPrefixRules),
2048  optimizeBetasFlag(optimizeBetasFlag),
2049  minBeta(minBeta),
2050  penalty(penalty)
2051{
2052  initialize(probList);
2053  float step = 2.0;
2054
2055  // initialize prior betas
2056
2057  // optimize betas
2058  if (optimizeBetasFlag)
2059      optimizeBetas();
2060
2061  // find front rules
2062  if (setPrefixRules)
2063  {
2064    bool changed = setBestPrefixRule();
2065    while (changed) {
2066        if (optimizeBetasFlag)
2067            optimizeBetas();
2068        changed = setBestPrefixRule();
2069    }
2070  }
2071
2072  // prepare results in Orange-like format
2073  TFloatList *aruleBetas = mlnew TFloatList();
2074  ruleBetas = aruleBetas;
2075  TRuleList *aprefixRules = mlnew TRuleList();
2076  prefixRules = aprefixRules;
2077  int i;
2078  for (i=0; i<rules->size(); i++)
2079    ruleBetas->push_back(currentState->betas[i]);
2080  for (i=0; i<currentState->prefixRules->size(); i++)
2081    prefixRules->push_back(rules->at(currentState->prefixRules->at(i)));
2082  delete [] skipRule;
2083}
2084
2085bool TRuleClassifier_logit::setBestPrefixRule()
2086{
2087// adding a prefix rule should improve brier score and it should not decrease AUC
2088
2089  PLogitClassifierState tempState;
2090  currentState->copyTo(tempState);
2091  int bestRuleI = -1;
2092  float best_brier_improvement = 0.0;
2093  float bestNewQuality = 0.0;
2094
2095  PDistribution apriori = getClassDistribution(examples, weightID);
2096  float apriori_prob, new_positive, new_covered, new_rf_rule, rf_rule;
2097  float oldQuality, newQuality, brier_improvement;
2098  float current_auc = currentState->getAUC();
2099  float current_brier = currentState->getBrierScore();
2100
2101  for (int i=0; i<rules->size(); i++) {
2102    // compute new coverage
2103    new_positive = 0.0; new_covered = 0.0;
2104    for (int j=0; j<examples->numberOfExamples(); j++)
2105        if (rules->at(i)->call(examples->at(j)) && !tempState->isExampleFixed[j])
2106        {
2107            if (getClassIndex(rules->at(i)) == examples->at(j).getClass().intV)
2108                new_positive ++;
2109            new_covered ++;
2110        }
2111    if (new_covered < 1)
2112        continue;
2113
2114/*  // rule should cover at least 50% of its own examples (so that conditions still make some sense)
2115    if (new_covered < rules->at(i)->classDistribution->abs / 2)
2116      continue; */
2117
2118    apriori_prob = apriori->atint(getClassIndex(rules->at(i))) / apriori->abs;
2119    new_rf_rule = new_positive / new_covered;
2120    if (new_rf_rule <= apriori_prob)
2121        continue;
2122
2123    // corrected quality of the rule (due to rules in the prefix set already)
2124    oldQuality = rules->at(i)->quality;
2125    rf_rule = rules->at(i)->classDistribution->atint(getClassIndex(rules->at(i)));
2126    rf_rule /= rules->at(i) -> classDistribution->abs;
2127    newQuality = oldQuality * (new_positive + 2*rf_rule) / (new_covered + 2) / rf_rule;
2128
2129    // check if improvement is good
2130    rules->at(i)->quality = newQuality;
2131    currentState->setPrefixRule(i);
2132    rules->at(i)->quality = oldQuality;
2133
2134    brier_improvement = (current_brier - currentState->getBrierScore()) / new_covered;
2135
2136    if (currentState->getAUC() >= current_auc && 
2137        brier_improvement > best_brier_improvement)
2138    {
2139        best_brier_improvement = brier_improvement;
2140        bestRuleI = i;
2141        bestNewQuality = newQuality;
2142    }
2143    tempState->copyTo(currentState);
2144  }
2145  if (bestRuleI > -1)
2146  {
2147      rules->at(bestRuleI)->quality = bestNewQuality;
2148      currentState->setPrefixRule(bestRuleI);
2149      skipRule[bestRuleI] = true;
2150      // compute new class distribution for this rule
2151      TExampleTable * newexamples = mlnew TExampleTable(examples->domain);
2152      PExampleGenerator pnewexamples = PExampleGenerator(newexamples);
2153      for (int ei=0; ei<examples->numberOfExamples(); ei++)
2154        if (!tempState->isExampleFixed[ei])
2155          newexamples->addExample(examples->at(ei));
2156      rules->at(bestRuleI)->filterAndStore(pnewexamples, rules->at(bestRuleI)->weightID, getClassIndex(rules->at(bestRuleI)));
2157      return true;
2158  }
2159  return false;
2160}
2161
2162void TRuleClassifier_logit::optimizeBetas()
2163{
2164  bool minSigChange = true;
2165  bool minBetaChange = true;
2166
2167  while (minSigChange || minBetaChange)
2168  {
2169      // learn initial model
2170      updateRuleBetas(2.0);
2171
2172      minSigChange = false;
2173      minBetaChange = false;
2174      for (int i=0; i<rules->size(); i++)
2175      {
2176         if (skipRule[i] || rules->at(i)->classDistribution->abs == prior->abs)
2177             continue;
2178         if (currentState->betas[i] < minBeta)
2179         {
2180             skipRule[i] = true;
2181             minBetaChange = true;
2182             currentState->newBeta(i,0.0);
2183         }
2184      }
2185
2186      // min significance check, the same story as minBeta
2187      // loop through rules, check if significant, if not set to 0
2188      for (int i=0; i<rules->size(); i++)
2189      {
2190         if (skipRule[i] || rules->at(i)->classDistribution->abs == prior->abs)
2191             continue;
2192         float oldBeta = currentState->betas[i];
2193         currentState->newBeta(i,0.0);
2194         if (currentState->avgProb->at(i) < (wSatQ->at(i) - 0.01))
2195             currentState->newBeta(i,oldBeta);
2196         else
2197         {
2198             skipRule[i] = true;
2199             minSigChange = true;
2200         }
2201      }
2202  }
2203}
2204
2205// Init current state
2206void TRuleClassifier_logit::initialize(const PDistributionList &probList)
2207{
2208  // compute prior distribution of learning examples
2209  prior = getClassDistribution(examples, weightID);
2210  domain = examples->domain;
2211  // set initial state
2212  TLogitClassifierState *ncurrentState = new TLogitClassifierState(rules, probList, examples, weightID);
2213  currentState = ncurrentState;
2214  // compute standard deviations of rules
2215  TFloatList *sd = new TFloatList();
2216  wsd = sd;
2217  TFloatList *sig = new TFloatList();
2218  wsig = sig;
2219  PITERATE(TRuleList, ri, rules) {
2220    float maxDiff = (*ri)->classDistribution->atint(getClassIndex(*ri))/(*ri)->classDistribution->abs;
2221    maxDiff -= (*ri)->quality;
2222    wsig->push_back(maxDiff);
2223    float n = (*ri)->examples->numberOfExamples();
2224    float a = n*(*ri)->quality;
2225    float b = n*(1.0-(*ri)->quality);
2226    float expab = log(a)+log(b)-2*log(a+b)-log(a+b+1);
2227    //wsd->push_back(exp(0.5*expab));
2228    // temprorarily, wsd will be expected beta for a rule
2229    float rf = prior->atint(getClassIndex(*ri))/prior->abs;
2230    if ((*ri)->classDistribution->abs < prior->abs)
2231        wsd->push_back(log((*ri)->quality/(1-(*ri)->quality))-log(rf/(1-rf)));
2232    else if (getClassIndex(*ri) < (*ri)->examples->domain->classVar->noOfValues() - 1)
2233        wsd->push_back(log(rf/(1-rf)));
2234    else
2235        wsd->push_back(0.0);
2236   // printf("%f %f %f %f\n", n, a, b, exp(0.5*expab));
2237  }
2238
2239  // compute satisfiable qualities
2240  TFloatList *satQ = new TFloatList();
2241  wSatQ = satQ;
2242  if (minSignificance >= 0.5 || minSignificance <= 0.0)
2243    PITERATE(TRuleList, ri, rules)
2244      //wSatQ->push_back((*ri)->quality);
2245      wSatQ->push_back(1);
2246  else
2247  {
2248    float n,a,b,error,high,low,av,avf;
2249    PITERATE(TRuleList, ri, rules) {
2250      n = (*ri)->examples->numberOfExamples();
2251      a = n*(*ri)->quality;
2252      b = n*(1.0-(*ri)->quality);
2253      error = 1.0;
2254      high = 1.0;
2255      low = (float) 0.0;
2256      while (error > 0.001)
2257      {
2258        av = (high+low)/2.0;
2259        avf = (float)betai(double(a),double(b),double(av));
2260        if (avf < minSignificance)
2261          low = av;
2262        else
2263          high = av;
2264        error = abs(avf-minSignificance);
2265      }
2266      wSatQ->push_back(av);
2267    }
2268  }
2269
2270  // Compute average example coverage and set index of examples covered by rule
2271  float **coverages = new float* [examples->domain->classVar->noOfValues()];
2272  int j=0,k=0;
2273  for (j=0; j<examples->domain->classVar->noOfValues(); j++)  {
2274    coverages[j] = new float[examples->numberOfExamples()];
2275    for (k=0; k<examples->numberOfExamples(); k++)
2276      coverages[j][k] = 0.0;
2277  }
2278  int i=0;
2279  {
2280    PITERATE(TRuleList, ri, rules) {
2281      j=0;
2282      PEITERATE(ei, examples) {
2283        if ((*ri)->call(*ei)) {
2284          //int vv = (*ei).getClass().intV;
2285          //if ((*ei).getClass().intV == getClassIndex(*ri))
2286          coverages[getClassIndex(*ri)][j] += 1.0;
2287        }
2288        j++;
2289      }
2290      i++;
2291    }
2292  }
2293
2294  // compute coverages of rules
2295  TFloatList *avgCov = new TFloatList();
2296  wavgCov = avgCov;
2297  for (i=0; i<rules->size(); i++) {
2298    float newCov = 0.0;
2299    float counter = 0.0;
2300    PITERATE(TIntList, ind, currentState->ruleIndices[i])
2301    {
2302      //if (getClassIndex(rules->at(i)) == examples->at(*ind).getClass().intV) {
2303      newCov += coverages[getClassIndex(rules->at(i))][*ind];
2304      counter++;
2305      //}
2306    }
2307    //printf(" ... %f ... %f \n", newCov, counter);
2308    if (counter) {
2309      wavgCov->push_back(newCov/counter);
2310    }
2311    else
2312      wavgCov->push_back(0.0);
2313  }
2314
2315  skipRule = new bool [rules->size()];
2316  for (i=0; i<rules->size(); i++)
2317      skipRule[i] = false;
2318}
2319
2320void TRuleClassifier_logit::updateRuleBetas(float step_)
2321{
2322    PLogitClassifierState finalState;
2323    currentState->copyTo(finalState);
2324
2325    float step = 0.1f;
2326    float error = 1e+20f;
2327    float old_error = 1e+21f;
2328
2329    int nsteps = 0;
2330    while (old_error > error || step > 0.00001)
2331    {
2332        // reduce step if improvement failed
2333        nsteps++;
2334        if (old_error < error && nsteps > 20 || nsteps > 1000)
2335        {
2336            nsteps = 0;
2337            step /= 10;
2338            finalState->copyTo(currentState);
2339        }
2340        else
2341            currentState->copyTo(finalState);
2342
2343        old_error = error;
2344        // update betas
2345        for (int i=0; i<rules->size(); i++) {
2346            if (skipRule[i])
2347                continue;
2348
2349            float der = 0.0;
2350            if (currentState->avgProb->at(i) > rules->at(i)->quality)
2351            {
2352                der += pow(rules->at(i)->quality - currentState->avgProb->at(i),2);
2353                der -= 2 * (rules->at(i)->quality - currentState->avgProb->at(i)) * currentState->betas[i];
2354            }
2355            else
2356                der -= 2 * (rules->at(i)->quality - currentState->avgProb->at(i));
2357            der += 2 * penalty * currentState->betas[i];
2358//            if (currentState->avgProb->at(i) > rules->at(i)->quality)
2359//                der = max(0.01f / step, der);
2360            currentState->newBeta(i,max(0.0f, currentState->betas[i]-step*der));
2361        }
2362        // estimate new error
2363        error = 0;
2364        for (int i=0; i<rules->size(); i++) {
2365//            if (currentState->avgProb->at(i) < rules->at(i)->quality)
2366            if (currentState->avgProb->at(i) > rules->at(i)->quality)
2367                error += pow(rules->at(i)->quality - currentState->avgProb->at(i),2) * currentState->betas[i];
2368            else
2369                error += pow(rules->at(i)->quality - currentState->avgProb->at(i),2);
2370            error += penalty*pow(currentState->betas[i],2);
2371        }
2372        //printf("error = %4.4f\n", error);
2373        // print betas
2374        //for (int i=0; i<rules->size(); i++)
2375        //    printf("%4.2f,",currentState->betas[i]);
2376        //printf("\n");
2377    }
2378
2379    finalState->copyTo(currentState);
2380}
2381
2382void TRuleClassifier_logit::addPriorClassifier(const TExample &ex, double * priorFs) {
2383  // initialize variables
2384  double *matrix = new double [sqr(examples->domain->classVar->noOfValues()-1)];
2385  double *fVals = new double [examples->domain->classVar->noOfValues()-1];
2386  double *results = new double [examples->domain->classVar->noOfValues()-1];
2387  int k1, k2;
2388  TDistribution *dist = mlnew TDiscDistribution(domain->classVar);
2389  PDistribution wdist = dist;
2390
2391  PDistribution classifierDist = priorClassifier->classDistribution(ex);
2392  // correct probablity if equals 1.0
2393  for (k1=0; k1<examples->domain->classVar->noOfValues(); k1++) {
2394    if (classifierDist->atint(k1) >= 1.0-1e-4)
2395      wdist->setint(k1,(float)(1.0-1e-4));
2396    else if (classifierDist->atint(k1) <= 1e-4)
2397      wdist->setint(k1,(float)(1e-4));
2398    else
2399      wdist->setint(k1,classifierDist->atint(k1));
2400  }
2401  wdist->normalize();
2402
2403
2404  // create matrix
2405  for (k1=0; k1<examples->domain->classVar->noOfValues()-1; k1++) {
2406    fVals[k1] = -wdist->atint(k1);
2407    for (k2=0; k2<examples->domain->classVar->noOfValues()-1; k2++) {
2408      if (k1==k2)
2409        matrix[k1*(examples->domain->classVar->noOfValues()-1)+k2] = (wdist->atint(k1)-1);
2410      else
2411        matrix[k1*(examples->domain->classVar->noOfValues()-1)+k2] = wdist->atint(k1);
2412    }
2413  }
2414  // solve equation
2415  LinearEquationsSolving(examples->domain->classVar->noOfValues()-1, matrix, fVals, results);
2416  for (k1=0; k1<examples->domain->classVar->noOfValues()-1; k1++)
2417    priorFs[k1] = results[k1]>0.0 ? log(results[k1]) : -10.0;
2418  // clean up
2419  delete [] matrix;
2420  delete [] fVals;
2421  delete [] results;
2422}
2423
2424PDistribution TRuleClassifier_logit::classDistribution(const TExample &ex)
2425{
2426  checkProperty(rules);
2427  checkProperty(prior);
2428  checkProperty(domain);
2429  TExample cexample(domain, ex);
2430
2431  TDiscDistribution *dist = mlnew TDiscDistribution(domain->classVar);
2432  PDistribution res = dist;
2433
2434  // if front rule triggers, use it first
2435  bool foundPrefixRule = false;
2436  float bestQuality = 0.0;
2437  PITERATE(TRuleList, rs, prefixRules) {   
2438    if ((*rs)->call(ex) && (*rs)->quality > bestQuality) {
2439        bestQuality = (*rs)->quality;
2440        dist->setint(getClassIndex(*rs),(*rs)->quality);
2441            for (int ci=0; ci<examples->domain->classVar->noOfValues(); ci++)
2442            if (ci!=getClassIndex(*rs))
2443                dist->setint(ci,(1.0-(*rs)->quality)/(examples->domain->classVar->noOfValues()-1));
2444        foundPrefixRule = true;
2445        break;
2446      }
2447  }
2448  if (foundPrefixRule)
2449    return dist;
2450
2451  // if correcting a classifier, use that one first then
2452  double * priorFs = new double [examples->domain->classVar->noOfValues()-1];
2453  if (priorClassifier)
2454    addPriorClassifier(ex, priorFs);
2455  else
2456    for (int k=0; k<examples->domain->classVar->noOfValues()-1; k++)
2457      priorFs[k] = 0.0;
2458  // compute return probabilities
2459  for (int i=0; i<res->noOfElements()-1; i++) {
2460    float f = priorFs[i];
2461    TFloatList::const_iterator b(ruleBetas->begin()), be(ruleBetas->end());
2462    TRuleList::iterator r(rules->begin()), re(rules->end());
2463    for (; r!=re; r++, b++)
2464      if ((*r)->call(cexample)) {
2465        if (getClassIndex(*r) == i)
2466            f += (*b);
2467        else if (getClassIndex(*r) == res->noOfElements()-1)
2468            f -= (*b);
2469      }
2470    dist->addint(i,exp(f));
2471  }
2472  dist->addint(res->noOfElements()-1,1.0);
2473  dist->normalize();
2474
2475  delete [] priorFs;
2476  return res;
2477}
Note: See TracBrowser for help on using the repository browser.