source: orange/source/orange/rulelearner.cpp @ 8735:065a34c267f2

Revision 8735:065a34c267f2, 81.6 KB checked in by matejd <matejd@…>, 3 years ago (diff)

Moved over code from qtgraph branch, turned primitives into a module (plot.primitives); added Visualize Qt folder to setup.py packages

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(-1),
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      }
1031      futureQuality = 1.0 - exp(futureQuality);
1032    }
1033  }
1034
1035  // store best rule as best rule and return expected quality of this rule
1036  rule->quality = quality;
1037  if (improved >= min_improved &&
1038      improved/rule->classDistribution->atint(targetClass) > min_improved_perc*0.01 &&
1039      quality > (aprioriProb + 1e-3) &&
1040      (!bestRule || (quality>bestRule->quality+1e-3)) &&
1041      (!validator || validator->call(rule, examples, weightID, targetClass, apriori))) {
1042
1043      TRule *pbestRule = new TRule(rule.getReference(), true);
1044      PRule wpbestRule = pbestRule;
1045      // check if rule is significant enough
1046      bool ruleGoodEnough = true;
1047
1048      if (ruleAlpha < 1.0)
1049        ruleGoodEnough = ruleGoodEnough & ((rule->chi > 0.0) && (chisqprob(rule->chi, 1.0f) <= ruleAlpha));
1050      if (ruleGoodEnough && attributeAlpha < 1.0)
1051        ruleGoodEnough = ruleGoodEnough & ruleAttSignificant(rule, examples, weightID, targetClass, apriori, aprioriProb);
1052      if (ruleGoodEnough)
1053      {
1054        bestRule = wpbestRule;
1055        bestRule->quality = quality;
1056        futureQuality += 1.0;
1057      }
1058  }
1059  return futureQuality;
1060}
1061
1062bool worstRule(const PRule &r1, const PRule &r2)
1063{ return    (r1->quality > r2->quality)
1064          || (r1->quality==r2->quality
1065          && r1->complexity < r2->complexity);
1066}
1067/*         || (r1->quality==r2->quality)
1068            && (   (r1->complexity < r2->complexity)
1069                || (r1->complexity == r2->complexity)
1070                   && ((int(r1.getUnwrappedPtr()) ^ int(r2.getUnwrappedPtr())) & 16) != 0
1071               ); }  */
1072
1073bool inRules(PRuleList rules, PRule rule)
1074{
1075  TRuleList::const_iterator ri(rules->begin()), re(rules->end());
1076  PExampleGenerator rulegen = rule->examples;
1077  for (; ri!=re; ri++) {
1078    PExampleGenerator rigen = (*ri)->examples;
1079    if (rigen->numberOfExamples() == rulegen->numberOfExamples()) {
1080      TExampleIterator rei(rulegen->begin()), ree(rulegen->end());
1081      TExampleIterator riei(rigen->begin()), riee(rigen->end());
1082      for (; rei != ree && !(*rei).compare(*riei); ++rei, ++riei) {
1083      }
1084        if (rei == ree)
1085          return true;
1086    }
1087  }
1088  return false;
1089}
1090
1091TRuleBeamFilter_Width::TRuleBeamFilter_Width(const int &w)
1092: width(w)
1093{}
1094
1095
1096void TRuleBeamFilter_Width::operator()(PRuleList &rules, PExampleTable, const int &)
1097{
1098  if (rules->size() > width) {
1099    sort(rules->begin(), rules->end(), worstRule);
1100
1101    TRuleList *filteredRules = mlnew TRuleList;
1102    PRuleList wFilteredRules = filteredRules;
1103
1104    int nRules = 0;
1105    TRuleList::const_iterator ri(rules->begin()), re(rules->end());
1106    while (nRules < width && ri != re) {
1107      if (!inRules(wFilteredRules,*ri)) {
1108        wFilteredRules->push_back(*ri);
1109        nRules++;
1110      }
1111      ri++;
1112    }
1113    rules =  wFilteredRules;
1114  }
1115}
1116
1117
1118inline void _selectBestRule(PRule &rule, PRule &bestRule, int &wins, TRandomGenerator &rgen)
1119{
1120  if ((rule->quality > bestRule->quality) || (rule->complexity < bestRule->complexity)) {
1121    bestRule = rule;
1122    wins = 1;
1123  }
1124  else if ((rule->complexity == bestRule->complexity) && rgen.randbool(++wins))
1125    bestRule = rule;
1126}
1127
1128
1129
1130PRuleList TRuleBeamInitializer_Default::operator()(PExampleTable data, const int &weightID, const int &targetClass, PRuleList baseRules, PRuleEvaluator evaluator, PDistribution apriori, PRule &bestRule)
1131{
1132  checkProperty(evaluator);
1133
1134  TRuleList *ruleList = mlnew TRuleList();
1135  PRuleList wruleList = ruleList;
1136
1137  TRandomGenerator rgen(data->numberOfExamples());
1138  int wins;
1139
1140  if (baseRules && baseRules->size())
1141    PITERATE(TRuleList, ri, baseRules) {
1142      TRule *newRule = mlnew TRule((*ri).getReference(), true);
1143      PRule wNewRule = newRule;
1144      ruleList->push_back(wNewRule);
1145      if (!newRule->examples)
1146      {
1147        newRule->filterAndStore(data,weightID,targetClass);
1148        newRule->baseDist = newRule->classDistribution;
1149      }
1150      newRule->quality = evaluator->call(wNewRule, data, weightID, targetClass, apriori);
1151      if (!bestRule || (newRule->quality > bestRule->quality)) {
1152        bestRule = wNewRule;
1153        wins = 1;
1154      }
1155      else
1156        if (newRule->quality == bestRule->quality)
1157          _selectBestRule(wNewRule, bestRule, wins, rgen);
1158    }
1159
1160  else {
1161     TRule *ubestRule = mlnew TRule();
1162     bestRule = ubestRule;
1163     ruleList->push_back(bestRule);
1164     ubestRule->filter = new TFilter_values();
1165     ubestRule->filter->domain = data->domain;
1166     ubestRule->filterAndStore(data, weightID,targetClass);
1167     ubestRule->baseDist = ubestRule->classDistribution;
1168     ubestRule->complexity = 0;
1169  }
1170
1171  return wruleList;
1172}
1173
1174
1175PRuleList TRuleBeamRefiner_Selector::operator()(PRule wrule, PExampleTable data, const int &weightID, const int &targetClass)
1176{
1177  if (!discretization) {
1178    discretization = mlnew TEntropyDiscretization();
1179    dynamic_cast<TEntropyDiscretization *>(discretization.getUnwrappedPtr())->forceAttribute = true;
1180    dynamic_cast<TEntropyDiscretization *>(discretization.getUnwrappedPtr())->maxNumberOfIntervals = 5;
1181  }
1182
1183  TRule &rule = wrule.getReference();
1184  TFilter_values *filter = wrule->filter.AS(TFilter_values);
1185  if (!filter)
1186    raiseError("a filter of type 'Filter_values' expected");
1187
1188  TRuleList *ruleList = mlnew TRuleList;
1189  PRuleList wRuleList = ruleList;
1190
1191  TDomainDistributions ddist(wrule->examples, wrule->weightID);
1192
1193  const TVarList &attributes = rule.examples->domain->attributes.getReference();
1194
1195  vector<bool> used(attributes.size(), false);
1196  PITERATE(TValueFilterList, vfi, filter->conditions)
1197    used[(*vfi)->position] = true;
1198
1199  vector<bool>::const_iterator ui(used.begin());
1200  TDomainDistributions::const_iterator di(ddist.begin());
1201  TVarList::const_iterator vi(attributes.begin()), ve(attributes.end());
1202  int pos = 0;
1203  for(; vi != ve; vi++, ui++, pos++, di++) {
1204    if ((*vi)->varType == TValue::INTVAR) {
1205      if (!*ui) {
1206        vector<float>::const_iterator idi((*di).AS(TDiscDistribution)->begin());
1207        for(int v = 0, e = (*vi)->noOfValues(); v != e; v++)
1208          if (*idi>0) {
1209            TRule *newRule = mlnew TRule(rule, false);
1210            ruleList->push_back(newRule);
1211            newRule->complexity++;
1212
1213            filter = newRule->filter.AS(TFilter_values);
1214
1215            TValueFilter_discrete *newCondition = mlnew TValueFilter_discrete(pos, *vi, 0);
1216            filter->conditions->push_back(newCondition);
1217
1218            TValue value = TValue(v);
1219            newCondition->values->push_back(value);
1220            newRule->filterAndStore(rule.examples, rule.weightID,targetClass);
1221            newRule->parentRule = wrule;
1222          }
1223      }
1224    }
1225
1226    else if (((*vi)->varType == TValue::FLOATVAR)) {
1227
1228        if (discretization) {
1229        PVariable discretized;
1230        try {
1231          discretized = discretization->call(rule.examples, *vi, weightID);
1232        } catch(...) {
1233          continue;
1234        }
1235        TClassifierFromVar *cfv = discretized->getValueFrom.AS(TClassifierFromVar);
1236        TDiscretizer *discretizer = cfv ? cfv->transformer.AS(TDiscretizer) : NULL;
1237        if (!discretizer)
1238          raiseError("invalid or unrecognized discretizer");
1239
1240        vector<float> cutoffs;
1241        discretizer->getCutoffs(cutoffs);
1242        if (cutoffs.size()) {
1243          TRule *newRule;
1244          newRule = mlnew TRule(rule, false);
1245          PRule wnewRule = newRule;
1246          newRule->complexity++;
1247          newRule->parentRule = wrule;
1248
1249          newRule->filter.AS(TFilter_values)->conditions->push_back(mlnew TValueFilter_continuous(pos,  TValueFilter_continuous::LessEqual,     cutoffs.front(), 0, 0));
1250          newRule->filterAndStore(rule.examples, rule.weightID,targetClass);
1251          if (wrule->classDistribution->cases > wnewRule->classDistribution->cases)
1252            ruleList->push_back(newRule);
1253
1254          for(vector<float>::const_iterator ci(cutoffs.begin()), ce(cutoffs.end()-1); ci != ce; ci++) {
1255            newRule = mlnew TRule(rule, false);
1256            wnewRule = newRule;
1257            newRule->complexity++;
1258            newRule->parentRule = wrule;
1259            filter = newRule->filter.AS(TFilter_values);
1260            filter->conditions->push_back(mlnew TValueFilter_continuous(pos,  TValueFilter_continuous::Greater, *ci, 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            newRule = mlnew TRule(rule, false);
1266            wnewRule = newRule;
1267            newRule->complexity++;
1268            newRule->parentRule = wrule;
1269            filter = newRule->filter.AS(TFilter_values);
1270            filter->conditions->push_back(mlnew TValueFilter_continuous(pos,  TValueFilter_continuous::LessEqual, *(ci+1), 0, 0));
1271            newRule->filterAndStore(rule.examples, rule.weightID,targetClass);
1272            if (wrule->classDistribution->cases > wnewRule->classDistribution->cases)
1273              ruleList->push_back(newRule);
1274          }
1275
1276          newRule = mlnew TRule(rule, false);
1277          ruleList->push_back(newRule);
1278          newRule->complexity++;
1279
1280          newRule->filter.AS(TFilter_values)->conditions->push_back(mlnew TValueFilter_continuous(pos,  TValueFilter_continuous::Greater, cutoffs.back(), 0, 0));
1281          newRule->filterAndStore(rule.examples, rule.weightID,targetClass);
1282          newRule->parentRule = wrule;
1283        }
1284      }
1285      else
1286        raiseWarning("discretizer not given, continuous attributes will be skipped");
1287    }
1288  }
1289
1290  if (!discretization)
1291    discretization = PDiscretization();
1292  return wRuleList;
1293}
1294
1295
1296PRuleList TRuleBeamCandidateSelector_TakeAll::operator()(PRuleList &existingRules, PExampleTable, const int &)
1297{
1298  PRuleList candidates = mlnew TRuleList(existingRules.getReference());
1299//  existingRules->clear();
1300  existingRules->erase(existingRules->begin(), existingRules->end());
1301  return candidates;
1302}
1303
1304
1305PRule TRuleBeamFinder::operator()(PExampleTable data, const int &weightID, const int &targetClass, PRuleList baseRules)
1306{
1307  // set default values if value not set
1308  bool tempInitializer = !initializer;
1309  if (tempInitializer)
1310    initializer = mlnew TRuleBeamInitializer_Default;
1311  bool tempCandidateSelector = !candidateSelector;
1312  if (tempCandidateSelector)
1313    candidateSelector = mlnew TRuleBeamCandidateSelector_TakeAll;
1314  bool tempRefiner = !refiner;
1315  if (tempRefiner)
1316    refiner = mlnew TRuleBeamRefiner_Selector;
1317/*  bool tempValidator = !validator;
1318  if (tempValidator)
1319    validator = mlnew TRuleValidator_LRS((float)0.01);
1320  bool tempRuleStoppingValidator = !ruleStoppingValidator;
1321  if (tempRuleStoppingValidator)
1322    ruleStoppingValidator = mlnew TRuleValidator_LRS((float)0.05); */
1323  bool tempEvaluator = !evaluator;
1324  if (tempEvaluator)
1325    evaluator = mlnew TRuleEvaluator_Entropy;
1326  bool tempRuleFilter = !ruleFilter;
1327  if (tempRuleFilter)
1328    ruleFilter = mlnew TRuleBeamFilter_Width;
1329
1330  // create an entropy evaluator to detect rules with pure distributions
1331  PRuleEvaluator entropy_evaluator = mlnew TRuleEvaluator_Entropy;
1332
1333  checkProperty(initializer);
1334  checkProperty(candidateSelector);
1335  checkProperty(refiner);
1336  checkProperty(evaluator);
1337  checkProperty(ruleFilter);
1338  PDistribution apriori = getClassDistribution(data, weightID);
1339
1340  TRandomGenerator rgen(data->numberOfExamples());
1341  int wins = 1;
1342
1343  PRule bestRule;
1344  PRuleList ruleList = initializer->call(data, weightID, targetClass, baseRules, evaluator, apriori, bestRule);
1345
1346  PITERATE(TRuleList, ri, ruleList) {
1347    if (!(*ri)->examples)
1348      (*ri)->filterAndStore(data, weightID,targetClass);
1349    if ((*ri)->quality == ILLEGAL_FLOAT)
1350      (*ri)->quality = evaluator->call(*ri, data, weightID, targetClass, apriori);
1351  }
1352
1353  if (!bestRule->examples)
1354    bestRule->filterAndStore(data, weightID,targetClass);
1355  if (bestRule->quality == ILLEGAL_FLOAT)
1356    bestRule->quality = evaluator->call(bestRule, data, weightID, targetClass, apriori);
1357
1358  int bestRuleLength = 0;
1359  while(ruleList->size()) {
1360    PRuleList candidateRules = candidateSelector->call(ruleList, data, weightID);
1361    PITERATE(TRuleList, ri, candidateRules) {
1362      PRuleList newRules = refiner->call(*ri, data, weightID, targetClass);
1363      PITERATE(TRuleList, ni, newRules) {
1364        // evaluate rule
1365        (*ni)->quality = evaluator->call(*ni, data, weightID, targetClass, apriori);
1366        if ((*ni)->quality > bestRule->quality && (!validator || validator->call(*ni, data, weightID, targetClass, apriori)))
1367          _selectBestRule(*ni, bestRule, wins, rgen);
1368        // if rule covers only one class or has less than 2 examples, do not add it for further specialization
1369        if (entropy_evaluator->call(*ni, data, weightID, targetClass, apriori) > -0.01 ||
1370            (*ni)->classDistribution->abs < 2)
1371            continue;
1372        // is the rule a good candidate for further specialization?
1373        if (!ruleStoppingValidator || ruleStoppingValidator->call(*ni, (*ri)->examples, weightID, targetClass, (*ri)->classDistribution)) {
1374          ruleList->push_back(*ni);
1375        }
1376      }
1377    }
1378    ruleFilter->call(ruleList,data,weightID);
1379  }
1380
1381  // set empty values if value was not set (used default)
1382  if (tempInitializer)
1383    initializer = PRuleBeamInitializer();
1384  if (tempCandidateSelector)
1385    candidateSelector = PRuleBeamCandidateSelector();
1386  if (tempRefiner)
1387    refiner = PRuleBeamRefiner();
1388/*  if (tempValidator)
1389    validator = PRuleValidator();
1390  if (tempRuleStoppingValidator)
1391    ruleStoppingValidator = PRuleValidator();  */
1392  if (tempEvaluator)
1393    evaluator = PRuleEvaluator();
1394  if (tempRuleFilter)
1395    ruleFilter = PRuleBeamFilter();
1396
1397  return bestRule;
1398}
1399
1400
1401TRuleLearner::TRuleLearner(bool se, int tc, PRuleList rl)
1402: storeExamples(se),
1403  targetClass(tc),
1404  baseRules(rl)
1405{}
1406
1407
1408PClassifier TRuleLearner::operator()(PExampleGenerator gen, const int &weightID)
1409{
1410  return this->call(gen,weightID,targetClass,baseRules);
1411}
1412
1413PClassifier TRuleLearner::operator()(PExampleGenerator gen, const int &weightID, const int &targetClass, PRuleList baseRules)
1414{
1415  // Initialize default values if values not set
1416  bool tempDataStopping = !dataStopping && !ruleStopping;
1417  if (tempDataStopping)
1418    dataStopping = mlnew TRuleDataStoppingCriteria_NoPositives;
1419
1420  bool tempRuleFinder = !ruleFinder;
1421  if (tempRuleFinder)
1422    ruleFinder = mlnew TRuleBeamFinder;
1423
1424  bool tempCoverAndRemove = !coverAndRemove;
1425  if (tempCoverAndRemove)
1426    coverAndRemove = mlnew TRuleCovererAndRemover_Default;
1427
1428  checkProperty(ruleFinder);
1429  checkProperty(coverAndRemove);
1430
1431  TExampleTable *data = mlnew TExampleTable(gen);
1432  PExampleTable wdata = data;
1433
1434  if (!dataStopping && !ruleStopping)
1435    raiseError("no stopping criteria; set 'dataStopping' and/or 'ruleStopping'");
1436
1437  TRuleList *ruleList = mlnew TRuleList;
1438  PRuleList wruleList = ruleList;
1439
1440  int currWeightID = weightID;
1441
1442  float beginwe=0.0, currentwe;
1443  if (progressCallback) {
1444    if (targetClass==-1)
1445      beginwe = wdata->weightOfExamples(weightID);
1446    else {
1447      PDistribution classDist = getClassDistribution(wdata, weightID);
1448      TDiscDistribution *ddist = classDist.AS(TDiscDistribution);
1449      beginwe = ddist->atint(targetClass);
1450    }
1451    progressCallback->call(0.0);
1452  }
1453
1454  while (!dataStopping || !dataStopping->call(wdata, currWeightID, targetClass)) {
1455    PRule rule = ruleFinder->call(wdata, currWeightID, targetClass, baseRules);
1456    if (!rule)
1457      raiseError("'ruleFinder' didn't return a rule");
1458
1459    if (ruleStopping && ruleStopping->call(ruleList, rule, wdata, currWeightID))
1460      break;
1461
1462    wdata = coverAndRemove->call(rule, wdata, currWeightID, currWeightID, targetClass);
1463    ruleList->push_back(rule);
1464
1465    if (progressCallback) {
1466      if (targetClass==-1)
1467        currentwe = wdata->weightOfExamples(weightID);
1468      else {
1469        PDistribution classDist = getClassDistribution(wdata, currWeightID);
1470        TDiscDistribution *ddist = classDist.AS(TDiscDistribution);
1471        currentwe = ddist->atint(targetClass);
1472      }
1473      progressCallback->call(1-currentwe/beginwe);
1474    }
1475  }
1476  if (progressCallback)
1477    progressCallback->call(1.0);
1478
1479
1480  // Restore values
1481  if (tempDataStopping)
1482    dataStopping = PRuleDataStoppingCriteria();
1483  if (tempRuleFinder)
1484    ruleFinder = PRuleFinder();
1485  if (tempCoverAndRemove)
1486    coverAndRemove = PRuleCovererAndRemover();
1487
1488  PRuleClassifierConstructor clConstructor =
1489    classifierConstructor ? classifierConstructor :
1490    PRuleClassifierConstructor(mlnew TRuleClassifierConstructor_firstRule());
1491  return clConstructor->call(ruleList, gen, weightID);
1492};
1493
1494
1495bool TRuleDataStoppingCriteria_NoPositives::operator()(PExampleTable data, const int &weightID, const int &targetClass) const
1496{
1497  PDistribution classDist = getClassDistribution(data, weightID);
1498  TDiscDistribution *ddist = classDist.AS(TDiscDistribution);
1499
1500  return (targetClass >= 0 ? ddist->atint(targetClass) : ddist->abs) == 0.0;
1501}
1502
1503bool TRuleStoppingCriteria_NegativeDistribution::operator()(PRuleList ruleList, PRule rule, PExampleTable data, const int &weightID) const
1504{
1505  if (rule && rule->classifier)
1506  {
1507    PDistribution aprioriDist = getClassDistribution(data, weightID);
1508    TDiscDistribution *apriori = aprioriDist.AS(TDiscDistribution);
1509
1510    const TDefaultClassifier *clsf = rule->classifier.AS(TDefaultClassifier);
1511    if (!clsf)
1512      return false;
1513    const TDiscDistribution *dist = dynamic_cast<const TDiscDistribution *>(clsf->defaultDistribution.getUnwrappedPtr());
1514    const int classVal = clsf->defaultVal.intV;
1515    if (classVal<0 || classVal>=dist->size())
1516      return false;
1517    float acc = dist->atint(clsf->defaultVal.intV)/dist->abs;
1518    float accApriori = apriori->atint(clsf->defaultVal.intV)/apriori->abs;
1519    if (accApriori>acc)
1520      return true;
1521  }
1522  return false;
1523}
1524
1525
1526PExampleTable TRuleCovererAndRemover_Default::operator()(PRule rule, PExampleTable data, const int &weightID, int &newWeight, const int &targetClass) const
1527{
1528  TExampleTable *table = mlnew TExampleTable(data, 1);
1529  PExampleGenerator wtable = table;
1530
1531  TFilter &filter = rule->filter.getReference();
1532
1533  if (targetClass < 0)
1534  {
1535    PEITERATE(ei, data)
1536      if (!filter(*ei))
1537        table->addExample(*ei);
1538  }
1539  else
1540    PEITERATE(ei, data)
1541      if (!filter(*ei) || (*ei).getClass().intV!=targetClass)
1542        table->addExample(*ei);
1543
1544
1545  newWeight = weightID;
1546  return wtable;
1547}
1548
1549// classifiers
1550PRuleClassifier TRuleClassifierConstructor_firstRule::operator ()(PRuleList rules, PExampleTable table, const int &weightID)
1551{
1552  return mlnew TRuleClassifier_firstRule(rules, table, weightID);
1553}
1554
1555
1556TRuleClassifier::TRuleClassifier(PRuleList arules, PExampleTable anexamples, const int &aweightID)
1557: rules(arules),
1558  examples(anexamples),
1559  weightID(aweightID),
1560  TClassifier(anexamples->domain->classVar,true)
1561{}
1562
1563TRuleClassifier::TRuleClassifier()
1564: TClassifier(true)
1565{}
1566
1567
1568TRuleClassifier_firstRule::TRuleClassifier_firstRule(PRuleList arules, PExampleTable anexamples, const int &aweightID)
1569: TRuleClassifier(arules, anexamples, aweightID)
1570{
1571  prior = getClassDistribution(examples, weightID);
1572}
1573
1574TRuleClassifier_firstRule::TRuleClassifier_firstRule()
1575: TRuleClassifier()
1576{}
1577
1578PDistribution TRuleClassifier_firstRule::classDistribution(const TExample &ex)
1579{
1580  checkProperty(rules);
1581  checkProperty(prior);
1582
1583  PITERATE(TRuleList, ri, rules) {
1584    if ((*ri)->call(ex))
1585      return (*ri)->classDistribution;
1586  }
1587  return prior;
1588}
1589
1590void copyTable(float **dest, float **source, int nx, int ny) {
1591  for (int i=0; i<nx; i++)
1592    memcpy(dest[i],source[i],ny*sizeof(float));
1593}
1594
1595//==============================================================================
1596// return 1 if system not solving
1597// nDim - system dimension
1598// pfMatr - matrix with coefficients
1599// pfVect - vector with free members
1600// pfSolution - vector with system solution
1601// pfMatr becames trianglular after function call
1602// pfVect changes after function call
1603//
1604// Developer: Henry Guennadi Levkin
1605//
1606//==============================================================================
1607int LinearEquationsSolving(int nDim, double* pfMatr, double* pfVect, double* pfSolution)
1608{
1609  double fMaxElem;
1610  double fAcc;
1611
1612  int i, j, k, m;
1613
1614
1615  for(k=0; k<(nDim-1); k++) // base row of matrix
1616  {
1617    // search of line with max element
1618    fMaxElem = fabs( pfMatr[k*nDim + k] );
1619    m = k;
1620    for(i=k+1; i<nDim; i++)
1621    {
1622      if(fMaxElem < fabs(pfMatr[i*nDim + k]) )
1623      {
1624        fMaxElem = pfMatr[i*nDim + k];
1625        m = i;
1626      }
1627    }
1628
1629    // permutation of base line (index k) and max element line(index m)
1630    if(m != k)
1631    {
1632      for(i=k; i<nDim; i++)
1633      {
1634        fAcc               = pfMatr[k*nDim + i];
1635        pfMatr[k*nDim + i] = pfMatr[m*nDim + i];
1636        pfMatr[m*nDim + i] = fAcc;
1637      }
1638      fAcc = pfVect[k];
1639      pfVect[k] = pfVect[m];
1640      pfVect[m] = fAcc;
1641    }
1642
1643    if( pfMatr[k*nDim + k] == 0.) return 1; // needs improvement !!!
1644
1645    // triangulation of matrix with coefficients
1646    for(j=(k+1); j<nDim; j++) // current row of matrix
1647    {
1648      fAcc = - pfMatr[j*nDim + k] / pfMatr[k*nDim + k];
1649      for(i=k; i<nDim; i++)
1650      {
1651        pfMatr[j*nDim + i] = pfMatr[j*nDim + i] + fAcc*pfMatr[k*nDim + i];
1652      }
1653      pfVect[j] = pfVect[j] + fAcc*pfVect[k]; // free member recalculation
1654    }
1655  }
1656
1657  for(k=(nDim-1); k>=0; k--)
1658  {
1659    pfSolution[k] = pfVect[k];
1660    for(i=(k+1); i<nDim; i++)
1661    {
1662      pfSolution[k] -= (pfMatr[k*nDim + i]*pfSolution[i]);
1663    }
1664    pfSolution[k] = pfSolution[k] / pfMatr[k*nDim + k];
1665  }
1666
1667  return 0;
1668}
1669
1670// extracts class value index from target in rule
1671int getClassIndex(PRule r) {
1672  const TDefaultClassifier &cl = dynamic_cast<const TDefaultClassifier &>(r->classifier.getReference());
1673  return cl.defaultVal.intV;
1674}
1675
1676// constructor 1
1677TLogitClassifierState::TLogitClassifierState(PRuleList arules, PExampleTable anexamples, const int &aweightID)
1678: rules(arules),
1679  examples(anexamples),
1680  weightID(aweightID)
1681{
1682  // initialize f, p
1683  f = new float *[examples->domain->classVar->noOfValues()-1];
1684  p = new float *[examples->domain->classVar->noOfValues()];
1685  int i;
1686  for (i=0; i<examples->domain->classVar->noOfValues()-1; i++) {
1687      f[i] = new float[examples->numberOfExamples()];
1688      p[i] = new float[examples->numberOfExamples()];
1689  }
1690  p[examples->domain->classVar->noOfValues()-1] = new float[examples->numberOfExamples()];
1691
1692  betas = new float[rules->size()];
1693  priorBetas = new float[examples->domain->classVar->noOfValues()];
1694  isExampleFixed = new bool[examples->numberOfExamples()];
1695}
1696
1697// constructor 2
1698TLogitClassifierState::TLogitClassifierState(PRuleList arules, const PDistributionList &probList, PExampleTable anexamples, const int &aweightID)
1699: rules(arules),
1700  examples(anexamples),
1701  weightID(aweightID)
1702{
1703  // initialize f, p
1704  f = new float *[examples->domain->classVar->noOfValues()-1];
1705  p = new float *[examples->domain->classVar->noOfValues()];
1706  int i, j;
1707  for (i=0; i<examples->domain->classVar->noOfValues()-1; i++) {
1708    f[i] = new float[examples->numberOfExamples()];
1709    p[i] = new float[examples->numberOfExamples()];
1710    for (j=0; j<examples->numberOfExamples(); j++) {
1711        f[i][j] = 0.0;
1712        p[i][j] = 1.0/examples->domain->classVar->noOfValues();
1713    }
1714  }
1715  p[examples->domain->classVar->noOfValues()-1] = new float[examples->numberOfExamples()];
1716  for (j=0; j<examples->numberOfExamples(); j++)
1717      p[examples->domain->classVar->noOfValues()-1][j] = 1.0/examples->domain->classVar->noOfValues();
1718
1719   // if initial example probability is given, update F and P
1720  if (probList) {
1721    double *matrix = new double [sqr(examples->domain->classVar->noOfValues()-1)];
1722    double *fVals = new double [examples->domain->classVar->noOfValues()-1];
1723    double *results = new double [examples->domain->classVar->noOfValues()-1];
1724    for (i=0; i<probList->size(); i++) {
1725      int k1, k2;
1726      TDistribution *dist = mlnew TDiscDistribution(probList->at(i)->variable);
1727      PDistribution wdist = dist;
1728      // Prepare and compute expected f - values (a linear equation)
1729      for (k1=0; k1<examples->domain->classVar->noOfValues(); k1++) {
1730        if (probList->at(i)->atint(k1) >= 1.0-1e-4)
1731          wdist->setint(k1,(float)(1.0-1e-4));
1732        else if (probList->at(i)->atint(k1) <= 1e-4)
1733          wdist->setint(k1,(float)(1e-4));
1734        else
1735          wdist->setint(k1,probList->at(i)->atint(k1));
1736      }
1737      wdist->normalize();
1738      for (k1=0; k1<examples->domain->classVar->noOfValues()-1; k1++) {
1739        fVals[k1] = -wdist->atint(k1);
1740        for (k2=0; k2<examples->domain->classVar->noOfValues()-1; k2++) {
1741          if (k1==k2)
1742            matrix[k1*(examples->domain->classVar->noOfValues()-1)+k2] = wdist->atint(k1)-1;
1743          else
1744            matrix[k1*(examples->domain->classVar->noOfValues()-1)+k2] = wdist->atint(k1);
1745        }
1746      }
1747      LinearEquationsSolving(examples->domain->classVar->noOfValues()-1, matrix, fVals, results);
1748      // store values
1749      for (k1=0; k1<examples->domain->classVar->noOfValues()-1; k1++)
1750        f[k1][i] = results[k1]>0.0 ? log(results[k1]) : -10.0;
1751      for (k1=0; k1<examples->domain->classVar->noOfValues(); k1++)
1752        p[k1][i] = wdist->atint(k1);
1753    }
1754    delete [] matrix;
1755    delete [] fVals;
1756    delete [] results;
1757  }
1758
1759  // compute rule indices
1760  i=0;
1761  ruleIndices = mlnew PIntList[rules->size()];
1762  {
1763    PITERATE(TRuleList, ri, rules) {
1764      TIntList *ruleIndicesnw = mlnew TIntList();
1765      ruleIndices[i] = ruleIndicesnw;
1766      j=0;
1767      PEITERATE(ei, examples) {
1768        if ((*ri)->call(*ei))
1769          ruleIndices[i]->push_back(j);
1770        j++;
1771      }
1772      i++;
1773    }
1774  }
1775
1776  // set initial values of betas
1777  betas = new float[rules->size()];
1778  for (i=0; i<rules->size(); i++)
1779    betas[i] = (float)0.0;
1780
1781  // Add default rules
1782  priorBetas = new float[examples->domain->classVar->noOfValues()];
1783  for (i=0; i<examples->domain->classVar->noOfValues(); i++)
1784    priorBetas[i] = 0.0;
1785
1786  // computer best rule covering
1787  PDistribution apriori = getClassDistribution(examples, weightID);
1788  isExampleFixed = new bool[examples->numberOfExamples()];
1789  for (j=0; j<examples->numberOfExamples(); j++)
1790    isExampleFixed[j] = false;
1791
1792  // priorProb and avgProb
1793  TFloatList *npriorProb = mlnew TFloatList();
1794  avgPriorProb = npriorProb;
1795  TFloatList *navgProb = mlnew TFloatList();
1796  avgProb = navgProb;
1797  TIntList *pprefixRules = mlnew TIntList();
1798  prefixRules = pprefixRules;
1799  computeAvgProbs();
1800  computePriorProbs();
1801}
1802
1803TLogitClassifierState::~TLogitClassifierState()
1804{
1805  int i;
1806  for (i=0; i<examples->domain->classVar->noOfValues()-1; i++)
1807    delete [] f[i];
1808  delete [] f;
1809
1810  for (i=0; i<examples->domain->classVar->noOfValues(); i++)
1811    delete [] p[i];
1812  delete [] p;
1813  delete [] betas;
1814  delete [] priorBetas;
1815  delete [] ruleIndices;
1816  delete [] isExampleFixed;
1817}
1818
1819void TLogitClassifierState::computeAvgProbs()
1820{
1821  // compute new rule avgProbs
1822  avgProb->clear();
1823  int classInd = 0;
1824
1825  float newAvgProb;
1826  for (int ri = 0; ri<rules->size(); ri++) {
1827    newAvgProb = 0.0;
1828    classInd = getClassIndex(rules->at(ri));
1829    PITERATE(TIntList, ind, ruleIndices[ri])
1830      newAvgProb += p[classInd][*ind];
1831    avgProb->push_back(newAvgProb/ruleIndices[ri]->size());
1832  }
1833}
1834
1835// compute new prior probs
1836void TLogitClassifierState::computePriorProbs()
1837{
1838  avgPriorProb->clear();
1839  for (int pi=0; pi<examples->domain->classVar->noOfValues(); pi++) {
1840    float newPriorProb = 0.0;
1841    for (int ei=0; ei<examples->numberOfExamples(); ei++) {
1842      newPriorProb += p[pi][ei];
1843    }
1844    avgPriorProb->push_back(newPriorProb/examples->numberOfExamples());
1845  }
1846}
1847
1848void TLogitClassifierState::copyTo(PLogitClassifierState & wstate)
1849{
1850  if (!wstate) {
1851    TLogitClassifierState *state = mlnew TLogitClassifierState(rules, examples, weightID);
1852    wstate = state;
1853    wstate->ruleIndices = mlnew PIntList[rules->size()];
1854    int i;
1855    for (i=0; i<rules->size(); i++) {
1856      TIntList * tIndices = mlnew TIntList(ruleIndices[i].getReference());
1857      wstate->ruleIndices[i] = tIndices;
1858    }
1859  }
1860
1861  wstate->eval = eval;
1862  copyTable(wstate->f, f, examples->domain->classVar->noOfValues()-1, examples->numberOfExamples());
1863  copyTable(wstate->p, p, examples->domain->classVar->noOfValues(), examples->numberOfExamples());
1864  memcpy(wstate->betas,betas,sizeof(float)*rules->size());
1865  memcpy(wstate->priorBetas,priorBetas,sizeof(float)*(examples->domain->classVar->noOfValues()-1));
1866  memcpy(wstate->isExampleFixed, isExampleFixed, sizeof(bool)*examples->numberOfExamples());
1867
1868  TFloatList *pavgProb = mlnew TFloatList(avgProb.getReference());
1869  TFloatList *pavgPriorProb = mlnew TFloatList(avgPriorProb.getReference());
1870  TIntList *pprefixRules = mlnew TIntList(prefixRules.getReference());
1871  wstate->avgProb = pavgProb;
1872  wstate->avgPriorProb = pavgPriorProb;
1873  wstate->prefixRules = pprefixRules;
1874}
1875
1876void TLogitClassifierState::newBeta(int i, float b)
1877{
1878  // set new beta
1879  float diff = b-betas[i];
1880  betas[i] = b;
1881
1882
1883  // add differences to f
1884  int classIndex = getClassIndex(rules->at(i));
1885  PITERATE(TIntList, ind, ruleIndices[i])
1886    for (int fi=0; fi<examples->domain->classVar->noOfValues()-1; fi++)
1887      if (fi == classIndex)
1888        f[fi][*ind] += diff;
1889      else if (classIndex == examples->domain->classVar->noOfValues()-1)
1890        f[fi][*ind] -= diff;
1891
1892  // compute p
1893  computePs(i);
1894  computeAvgProbs();
1895  computePriorProbs();
1896}
1897
1898void TLogitClassifierState::newPriorBeta(int i, float b)
1899{
1900  // set new beta
1901  float diff = b-priorBetas[i];
1902  priorBetas[i] = b;
1903
1904  // add differences to f
1905  for (int ei=0; ei<examples->numberOfExamples(); ei++)
1906    for (int fi=0; fi<examples->domain->classVar->noOfValues()-1; fi++)
1907      if (fi == i)
1908        f[fi][ei] += diff;
1909  computePs(-1);
1910  computeAvgProbs();
1911  computePriorProbs();
1912}
1913
1914void TLogitClassifierState::updateExampleP(int ei)
1915{
1916/*  PITERATE(TIntList, ind, frontRules)
1917      if (rules->at(*ind)->call(examples->at(ei))) {
1918        p[getClassIndex(rules->at(*ind))][ei] = rules->at(*ind)->quality;
1919        for (int ci=0; ci<examples->domain->classVar->noOfValues(); ci++)
1920            if (ci!=getClassIndex(rules->at(*ind)))
1921                p[ci][ei] = (1.0-rules->at(*ind)->quality)/(examples->domain->classVar->noOfValues()-1);
1922        return;
1923      } */
1924  if (isExampleFixed[ei])
1925      return;
1926
1927  float sum = 1.0;
1928  int pi;
1929  for (pi=0; pi<examples->domain->classVar->noOfValues()-1; pi++) {
1930    if (f[pi][ei] > 10.0)
1931        p[pi][ei] = 22000.0;
1932    else
1933        p[pi][ei] = exp(f[pi][ei]);
1934    sum += p[pi][ei];
1935  }
1936  p[examples->domain->classVar->noOfValues()-1][ei] = 1.0;
1937  for (pi=0; pi<examples->domain->classVar->noOfValues(); pi+=1) {
1938    p[pi][ei] /= sum;
1939  }
1940}
1941
1942void TLogitClassifierState::computePs(int beta_i)
1943{
1944  if (beta_i<0)
1945    for (int ei=0; ei<examples->numberOfExamples(); ei++)
1946      updateExampleP(ei);
1947  else
1948    PITERATE(TIntList, ind, ruleIndices[beta_i])
1949      updateExampleP(*ind);
1950}
1951
1952void TLogitClassifierState::setFixed(int rule_i)
1953{
1954    PITERATE(TIntList, ind, ruleIndices[rule_i])
1955      isExampleFixed[*ind] = true;
1956}
1957
1958void TLogitClassifierState::updateFixedPs(int rule_i)
1959{
1960    PITERATE(TIntList, ind, ruleIndices[rule_i])
1961  {
1962    float bestQuality = 0.0;
1963        PITERATE(TIntList, fr, prefixRules) {
1964              if (rules->at(*fr)->call(examples->at(*ind)) && rules->at(*fr)->quality > bestQuality) {
1965          bestQuality = rules->at(*fr)->quality;
1966                  p[getClassIndex(rules->at(*fr))][*ind] = rules->at(*fr)->quality;
1967                  for (int ci=0; ci<examples->domain->classVar->noOfValues(); ci++)
1968                      if (ci!=getClassIndex(rules->at(*fr)))
1969                          p[ci][*ind] = (1.0-rules->at(*fr)->quality)/(examples->domain->classVar->noOfValues()-1);
1970                  break;
1971              }
1972      }
1973  }
1974}
1975
1976float TLogitClassifierState::getBrierScore()
1977{
1978  float brier = 0.0;
1979  for (int j=0; j<examples->numberOfExamples(); j++)
1980    brier += pow(1.0-p[examples->at(j).getClass().intV][j],2.0);
1981  return brier;
1982}
1983
1984float TLogitClassifierState::getAUC()
1985{
1986  float auc = 0.0, sum_ranks;
1987  int n1, n2;
1988  vector<double> probs, ranks;
1989//  for each class
1990  for (int cl=0; cl < examples->domain->classVar->noOfValues() - 1; cl++)
1991  {
1992//    probs = prediction probabilities of class
1993    probs.clear();
1994    for (int j=0; j<examples->numberOfExamples(); j++)
1995      probs.push_back(p[cl][j]);
1996
1997//    get ranks of class examples
1998//    n1 = elements of examples from class
1999//    n2 = elements of examples not from class
2000    rankdata(probs, ranks);
2001    n1=n2=0;
2002    sum_ranks=0.0;
2003    for (int j=0; j<examples->numberOfExamples(); j++)
2004    {
2005      if (examples->at(j).getClass().intV == cl)
2006      {
2007        n1++;
2008        sum_ranks+=ranks.at(j);
2009      }
2010      else n2++;
2011    }
2012//    auc += (sum(ranks) - n1*(n1+1)/2) / (n1*n2)
2013    auc += (sum_ranks - n1*(n1+1)/2) / (n1*n2);
2014  }
2015  return auc / (examples->domain->classVar->noOfValues() - 1);
2016}
2017
2018void TLogitClassifierState::setPrefixRule(int rule_i) //, int position)
2019{
2020    prefixRules->push_back(rule_i);
2021    setFixed(rule_i);
2022    updateFixedPs(rule_i);
2023    betas[rule_i] = 0.0;
2024  computeAvgProbs();
2025  computePriorProbs();
2026}
2027
2028TRuleClassifier_logit::TRuleClassifier_logit()
2029: TRuleClassifier()
2030{}
2031
2032TRuleClassifier_logit::TRuleClassifier_logit(PRuleList arules, const float &minSignificance, const float &minBeta, PExampleTable anexamples, const int &aweightID, const PClassifier &classifier, const PDistributionList &probList, bool setPrefixRules, bool optimizeBetasFlag)
2033: TRuleClassifier(arules, anexamples, aweightID),
2034  minSignificance(minSignificance),
2035  priorClassifier(classifier),
2036  setPrefixRules(setPrefixRules),
2037  optimizeBetasFlag(optimizeBetasFlag),
2038  minBeta(minBeta)
2039{
2040  initialize(probList);
2041  float step = 2.0;
2042  minStep = (float)0.01;
2043
2044  // initialize prior betas
2045
2046  // optimize betas
2047  if (optimizeBetasFlag)
2048      optimizeBetas();
2049
2050  // find front rules
2051  if (setPrefixRules)
2052  {
2053      bool changed = setBestPrefixRule();
2054      while (changed) {
2055        if (optimizeBetasFlag)
2056            optimizeBetas();
2057        changed = setBestPrefixRule();
2058      }
2059  }
2060
2061  // prepare results in Orange-like format
2062  TFloatList *aruleBetas = mlnew TFloatList();
2063  ruleBetas = aruleBetas;
2064  TRuleList *aprefixRules = mlnew TRuleList();
2065  prefixRules = aprefixRules;
2066  int i;
2067  for (i=0; i<rules->size(); i++)
2068    ruleBetas->push_back(currentState->betas[i]);
2069  for (i=0; i<currentState->prefixRules->size(); i++)
2070    prefixRules->push_back(rules->at(currentState->prefixRules->at(i)));
2071  delete [] skipRule;
2072}
2073
2074bool TRuleClassifier_logit::setBestPrefixRule()
2075{
2076// adding a prefix rule should improve brier score and it should not decrease AUC
2077
2078  PLogitClassifierState tempState;
2079  currentState->copyTo(tempState);
2080  int bestRuleI = -1;
2081  float best_brier_improvement = 0.0;
2082  float bestNewQuality = 0.0;
2083
2084  PDistribution apriori = getClassDistribution(examples, weightID);
2085  float apriori_prob, new_positive, new_covered, new_rf_rule, rf_rule;
2086  float oldQuality, newQuality, brier_improvement;
2087  float current_auc = currentState->getAUC();
2088  float current_brier = currentState->getBrierScore();
2089
2090  for (int i=0; i<rules->size(); i++) {
2091    // compute new coverage
2092    new_positive = 0.0; new_covered = 0.0;
2093    for (int j=0; j<examples->numberOfExamples(); j++)
2094        if (rules->at(i)->call(examples->at(j)) && !tempState->isExampleFixed[j])
2095        {
2096            if (getClassIndex(rules->at(i)) == examples->at(j).getClass().intV)
2097                new_positive ++;
2098            new_covered ++;
2099        }
2100    if (new_covered < 1)
2101        continue;
2102
2103/*  // rule should cover at least 50% of its own examples (so that conditions still make some sense)
2104    if (new_covered < rules->at(i)->classDistribution->abs / 2)
2105      continue; */
2106
2107    apriori_prob = apriori->atint(getClassIndex(rules->at(i))) / apriori->abs;
2108    new_rf_rule = new_positive / new_covered;
2109    if (new_rf_rule <= apriori_prob)
2110        continue;
2111
2112    // corrected quality of the rule (due to rules in the prefix set already)
2113    oldQuality = rules->at(i)->quality;
2114    rf_rule = rules->at(i)->classDistribution->atint(getClassIndex(rules->at(i)));
2115    rf_rule /= rules->at(i) -> classDistribution->abs;
2116    newQuality = oldQuality * (new_positive + 2*rf_rule) / (new_covered + 2) / rf_rule;
2117
2118    // check if improvement is good
2119    rules->at(i)->quality = newQuality;
2120    currentState->setPrefixRule(i);
2121    rules->at(i)->quality = oldQuality;
2122
2123    brier_improvement = (current_brier - currentState->getBrierScore()) / new_covered;
2124
2125    if (currentState->getAUC() >= current_auc && 
2126        brier_improvement > best_brier_improvement)
2127    {
2128        best_brier_improvement = brier_improvement;
2129        bestRuleI = i;
2130        bestNewQuality = newQuality;
2131    }
2132    tempState->copyTo(currentState);
2133  }
2134  if (bestRuleI > -1)
2135  {
2136      rules->at(bestRuleI)->quality = bestNewQuality;
2137      currentState->setPrefixRule(bestRuleI);
2138      skipRule[bestRuleI] = true;
2139      // compute new class distribution for this rule
2140      TExampleTable * newexamples = mlnew TExampleTable(examples->domain);
2141      PExampleGenerator pnewexamples = PExampleGenerator(newexamples);
2142      for (int ei=0; ei<examples->numberOfExamples(); ei++)
2143        if (!tempState->isExampleFixed[ei])
2144          newexamples->addExample(examples->at(ei));
2145      rules->at(bestRuleI)->filterAndStore(pnewexamples, rules->at(bestRuleI)->weightID, getClassIndex(rules->at(bestRuleI)));
2146      return true;
2147  }
2148  return false;
2149}
2150
2151void TRuleClassifier_logit::optimizeBetas()
2152{
2153  bool minSigChange = true;
2154  bool minBetaChange = true;
2155
2156  while (minSigChange || minBetaChange)
2157  {
2158      // learn initial model
2159      updateRuleBetas(2.0);
2160
2161      minSigChange = false;
2162      minBetaChange = false;
2163      for (int i=0; i<rules->size(); i++)
2164      {
2165         if (skipRule[i] || rules->at(i)->classDistribution->abs == prior->abs)
2166             continue;
2167         if (currentState->betas[i] < minBeta)
2168         {
2169             skipRule[i] = true;
2170             minBetaChange = true;
2171             currentState->newBeta(i,0.0);
2172         }
2173      }
2174
2175      // min significance check, the same story as minBeta
2176      // loop through rules, check if significant, if not set to 0
2177      for (int i=0; i<rules->size(); i++)
2178      {
2179         if (skipRule[i] || rules->at(i)->classDistribution->abs == prior->abs)
2180             continue;
2181         float oldBeta = currentState->betas[i];
2182         currentState->newBeta(i,0.0);
2183         if (currentState->avgProb->at(i) < (wSatQ->at(i) - 0.01))
2184             currentState->newBeta(i,oldBeta);
2185         else
2186         {
2187             skipRule[i] = true;
2188             minSigChange = true;
2189         }
2190      }
2191  }
2192}
2193
2194// Init current state
2195void TRuleClassifier_logit::initialize(const PDistributionList &probList)
2196{
2197  // compute prior distribution of learning examples
2198  prior = getClassDistribution(examples, weightID);
2199  domain = examples->domain;
2200  // set initial state
2201  TLogitClassifierState *ncurrentState = new TLogitClassifierState(rules, probList, examples, weightID);
2202  currentState = ncurrentState;
2203  // compute standard deviations of rules
2204  TFloatList *sd = new TFloatList();
2205  wsd = sd;
2206  TFloatList *sig = new TFloatList();
2207  wsig = sig;
2208  PITERATE(TRuleList, ri, rules) {
2209    float maxDiff = (*ri)->classDistribution->atint(getClassIndex(*ri))/(*ri)->classDistribution->abs;
2210      maxDiff -= (*ri)->quality;
2211      wsig->push_back(maxDiff);
2212   float n = (*ri)->examples->numberOfExamples();
2213    float a = n*(*ri)->quality;
2214    float b = n*(1.0-(*ri)->quality);
2215    float expab = log(a)+log(b)-2*log(a+b)-log(a+b+1);
2216    //wsd->push_back(exp(0.5*expab));
2217    // temprorarily, wsd will be expected beta for a rule
2218    float rf = prior->atint(getClassIndex(*ri))/prior->abs;
2219    if ((*ri)->classDistribution->abs < prior->abs)
2220        wsd->push_back(log((*ri)->quality/(1-(*ri)->quality))-log(rf/(1-rf)));
2221    else if (getClassIndex(*ri) < (*ri)->examples->domain->classVar->noOfValues() - 1)
2222        wsd->push_back(log(rf/(1-rf)));
2223    else
2224        wsd->push_back(0.0);
2225   // printf("%f %f %f %f\n", n, a, b, exp(0.5*expab));
2226  }
2227
2228  // compute satisfiable qualities
2229  TFloatList *satQ = new TFloatList();
2230  wSatQ = satQ;
2231  if (minSignificance >= 0.5 || minSignificance <= 0.0)
2232    PITERATE(TRuleList, ri, rules)
2233      //wSatQ->push_back((*ri)->quality);
2234      wSatQ->push_back(1);
2235  else
2236  {
2237    float n,a,b,error,high,low,av,avf;
2238    PITERATE(TRuleList, ri, rules) {
2239      n = (*ri)->examples->numberOfExamples();
2240      a = n*(*ri)->quality;
2241      b = n*(1.0-(*ri)->quality);
2242      error = 1.0;
2243      high = 1.0;
2244      low = (float) 0.0;
2245      while (error > 0.001)
2246      {
2247        av = (high+low)/2.0;
2248        avf = (float)betai(double(a),double(b),double(av));
2249        if (avf < minSignificance)
2250          low = av;
2251        else
2252          high = av;
2253        error = abs(avf-minSignificance);
2254      }
2255      wSatQ->push_back(av);
2256    }
2257  }
2258
2259  // Compute average example coverage and set index of examples covered by rule
2260  float **coverages = new float* [examples->domain->classVar->noOfValues()];
2261  int j=0,k=0;
2262  for (j=0; j<examples->domain->classVar->noOfValues(); j++)  {
2263    coverages[j] = new float[examples->numberOfExamples()];
2264    for (k=0; k<examples->numberOfExamples(); k++)
2265      coverages[j][k] = 0.0;
2266  }
2267  int i=0;
2268  {
2269    PITERATE(TRuleList, ri, rules) {
2270      j=0;
2271      PEITERATE(ei, examples) {
2272        if ((*ri)->call(*ei)) {
2273          //int vv = (*ei).getClass().intV;
2274              //if ((*ei).getClass().intV == getClassIndex(*ri))
2275                coverages[getClassIndex(*ri)][j] += 1.0;
2276        }
2277          j++;
2278      }
2279      i++;
2280    }
2281  }
2282
2283  // compute coverages of rules
2284  TFloatList *avgCov = new TFloatList();
2285  wavgCov = avgCov;
2286  for (i=0; i<rules->size(); i++) {
2287    float newCov = 0.0;
2288    float counter = 0.0;
2289    PITERATE(TIntList, ind, currentState->ruleIndices[i])
2290    {
2291      //if (getClassIndex(rules->at(i)) == examples->at(*ind).getClass().intV) {
2292      newCov += coverages[getClassIndex(rules->at(i))][*ind];
2293      counter++;
2294      //}
2295    }
2296    //printf(" ... %f ... %f \n", newCov, counter);
2297    if (counter) {
2298      wavgCov->push_back(newCov/counter);
2299    }
2300    else
2301      wavgCov->push_back(0.0);
2302  }
2303
2304  skipRule = new bool [rules->size()];
2305  for (i=0; i<rules->size(); i++)
2306      skipRule[i] = false;
2307}
2308
2309void TRuleClassifier_logit::updateRuleBetas(float step_)
2310{
2311    PLogitClassifierState finalState;
2312    currentState->copyTo(finalState);
2313
2314    float step = 0.1f;
2315    float gamma = 0.01f;
2316    float error = 1e+20f;
2317    float old_error = 1e+21f;
2318
2319    while (old_error > error || step > 0.00001)
2320    {
2321        // reduce step if improvement failed
2322        if (old_error < error)
2323        {
2324            step /= 10;
2325            finalState->copyTo(currentState);
2326        }
2327        else
2328            currentState->copyTo(finalState);
2329
2330        old_error = error;
2331        // update betas
2332        for (int i=0; i<rules->size(); i++) {
2333            if (skipRule[i])
2334                continue;
2335
2336            float der = 0.0;
2337            if (currentState->avgProb->at(i) < rules->at(i)->quality)
2338                der -= rules->at(i)->quality - currentState->avgProb->at(i);
2339            der += 2*gamma*currentState->betas[i];
2340            currentState->newBeta(i,max(0.0f, currentState->betas[i]-step*der));
2341        }
2342        // estimate new error
2343        error = 0;
2344        for (int i=0; i<rules->size(); i++) {
2345            if (currentState->avgProb->at(i) < rules->at(i)->quality)
2346                error += rules->at(i)->quality - currentState->avgProb->at(i);
2347            error += gamma*pow(currentState->betas[i],2);
2348        }
2349        //printf("error = %4.4f\n", error);
2350        // print betas
2351        //for (int i=0; i<rules->size(); i++)
2352        //    printf("%4.2f,",currentState->betas[i]);
2353        //printf("\n");
2354    }
2355    finalState->copyTo(currentState);
2356}
2357
2358/*
2359void TRuleClassifier_logit::updateRuleBetas2(float step_)
2360{
2361
2362  stabilizeAndEvaluate(step_,-1);
2363  PLogitClassifierState finalState, tempState;
2364  currentState->copyTo(finalState);
2365
2366  float step = 2.0;
2367  int changed;
2368  float worst_underestimate, underestimate;
2369  float auc = currentState->getAUC();
2370  float brier = currentState->getBrierScore();
2371  float temp_auc, temp_brier;
2372  int worst_rule_index;
2373  while (step > 0.001)
2374  {
2375      step /= 2;
2376      changed = 0;
2377      while (changed < 100)
2378      {
2379        changed = 0;
2380        worst_underestimate = (float)0.01;
2381        worst_rule_index = -1;
2382        // find rule with greatest underestimate in probability
2383        for (int i=0; i<rules->size(); i++) {
2384            if (currentState->avgProb->at(i) >= rules->at(i)->quality)
2385                continue;
2386            if (skipRule[i])
2387                continue;
2388
2389            underestimate = (rules->at(i)->quality - currentState->avgProb->at(i));//*rules->at(i)->classDistribution->abs;
2390            // if under estimate error is big enough
2391            if (underestimate > worst_underestimate)
2392            {
2393                worst_underestimate = underestimate;
2394                worst_rule_index = i;
2395            }
2396        }
2397        if (worst_rule_index > -1)
2398        {
2399            currentState->newBeta(worst_rule_index,currentState->betas[worst_rule_index]+step);
2400            if (currentState->avgProb->at(worst_rule_index) > rules->at(worst_rule_index)->quality)
2401            {
2402                finalState->copyTo(currentState);
2403                changed = 100;
2404            }
2405            else
2406            {
2407              stabilizeAndEvaluate(step,-1);
2408              temp_auc = currentState->getAUC();
2409              temp_brier = currentState->getBrierScore();
2410              if (temp_auc >= auc && temp_brier < brier)
2411              {
2412                currentState->copyTo(finalState);
2413                changed = 0;
2414                auc = temp_auc;
2415                brier = temp_brier;
2416              }
2417              else
2418                changed ++;
2419            }
2420         // }
2421        }
2422        else
2423        {
2424          changed = 100;
2425          finalState->copyTo(currentState);
2426        }
2427      }
2428  }
2429  finalState->copyTo(currentState);
2430}
2431*/
2432
2433void TRuleClassifier_logit::stabilizeAndEvaluate(float & step, int last_changed_rule_index)
2434{
2435    PLogitClassifierState tempState;
2436    currentState->copyTo(tempState);
2437    bool changed = true;
2438    while (changed)
2439    {
2440        changed = false;
2441        for (int i=0; i<rules->size(); i++)
2442        {
2443            if (currentState->avgProb->at(i) > (rules->at(i)->quality + 0.01) && currentState->betas[i] > 0.0 &&
2444                i != last_changed_rule_index)
2445            {
2446                float new_beta = currentState->betas[i]-step > 0 ? currentState->betas[i]-step : 0.0;
2447                currentState->newBeta(i,new_beta);
2448                if (currentState->avgProb->at(i) < rules->at(i)->quality + 1e-6)
2449                {
2450                    tempState->copyTo(currentState);
2451                }
2452                else
2453                {
2454                    currentState->copyTo(tempState);
2455                    changed = true;
2456                }
2457            }
2458        }
2459    }
2460}
2461
2462
2463void TRuleClassifier_logit::addPriorClassifier(const TExample &ex, double * priorFs) {
2464  // initialize variables
2465  double *matrix = new double [sqr(examples->domain->classVar->noOfValues()-1)];
2466  double *fVals = new double [examples->domain->classVar->noOfValues()-1];
2467  double *results = new double [examples->domain->classVar->noOfValues()-1];
2468  int k1, k2;
2469  TDistribution *dist = mlnew TDiscDistribution(domain->classVar);
2470  PDistribution wdist = dist;
2471
2472  PDistribution classifierDist = priorClassifier->classDistribution(ex);
2473  // correct probablity if equals 1.0
2474  for (k1=0; k1<examples->domain->classVar->noOfValues(); k1++) {
2475    if (classifierDist->atint(k1) >= 1.0-1e-4)
2476      wdist->setint(k1,(float)(1.0-1e-4));
2477    else if (classifierDist->atint(k1) <= 1e-4)
2478      wdist->setint(k1,(float)(1e-4));
2479    else
2480      wdist->setint(k1,classifierDist->atint(k1));
2481  }
2482  wdist->normalize();
2483
2484
2485  // create matrix
2486  for (k1=0; k1<examples->domain->classVar->noOfValues()-1; k1++) {
2487    fVals[k1] = -wdist->atint(k1);
2488    for (k2=0; k2<examples->domain->classVar->noOfValues()-1; k2++) {
2489      if (k1==k2)
2490        matrix[k1*(examples->domain->classVar->noOfValues()-1)+k2] = (wdist->atint(k1)-1);
2491      else
2492        matrix[k1*(examples->domain->classVar->noOfValues()-1)+k2] = wdist->atint(k1);
2493    }
2494  }
2495  // solve equation
2496  LinearEquationsSolving(examples->domain->classVar->noOfValues()-1, matrix, fVals, results);
2497  for (k1=0; k1<examples->domain->classVar->noOfValues()-1; k1++)
2498    priorFs[k1] = results[k1]>0.0 ? log(results[k1]) : -10.0;
2499  // clean up
2500  delete [] matrix;
2501  delete [] fVals;
2502  delete [] results;
2503}
2504
2505PDistribution TRuleClassifier_logit::classDistribution(const TExample &ex)
2506{
2507  checkProperty(rules);
2508  checkProperty(prior);
2509  checkProperty(domain);
2510  TExample cexample(domain, ex);
2511
2512  TDiscDistribution *dist = mlnew TDiscDistribution(domain->classVar);
2513  PDistribution res = dist;
2514
2515  // if front rule triggers, use it first
2516  bool foundPrefixRule = false;
2517  float bestQuality = 0.0;
2518  PITERATE(TRuleList, rs, prefixRules) {
2519      if ((*rs)->call(ex) && (*rs)->quality > bestQuality) {
2520      bestQuality = (*rs)->quality;
2521          dist->setint(getClassIndex(*rs),(*rs)->quality);
2522          for (int ci=0; ci<examples->domain->classVar->noOfValues(); ci++)
2523            if (ci!=getClassIndex(*rs))
2524                dist->setint(ci,(1.0-(*rs)->quality)/(examples->domain->classVar->noOfValues()-1));
2525      foundPrefixRule = true;
2526      break;
2527      }
2528  }
2529  if (foundPrefixRule)
2530    return dist;
2531
2532  // if correcting a classifier, use that one first then
2533  double * priorFs = new double [examples->domain->classVar->noOfValues()-1];
2534  if (priorClassifier)
2535    addPriorClassifier(ex, priorFs);
2536  else
2537    for (int k=0; k<examples->domain->classVar->noOfValues()-1; k++)
2538      priorFs[k] = 0.0;
2539  // compute return probabilities
2540  for (int i=0; i<res->noOfElements()-1; i++) {
2541    float f = priorFs[i];
2542    TFloatList::const_iterator b(ruleBetas->begin()), be(ruleBetas->end());
2543    TRuleList::iterator r(rules->begin()), re(rules->end());
2544    for (; r!=re; r++, b++)
2545      if ((*r)->call(cexample)) {
2546        if (getClassIndex(*r) == i)
2547            f += (*b);
2548        else if (getClassIndex(*r) == res->noOfElements()-1)
2549          f -= (*b);
2550      }
2551    dist->addint(i,exp(f));
2552  }
2553  dist->addint(res->noOfElements()-1,1.0);
2554  dist->normalize();
2555
2556  delete [] priorFs;
2557  return res;
2558}
Note: See TracBrowser for help on using the repository browser.