source: orange/source/orange/distvars.cpp @ 8103:84dd307db8a3

Revision 8103:84dd307db8a3, 33.4 KB checked in by ales_erjavec <ales.erjavec@…>, 3 years ago (diff)

Fixed randomFloat() methods (fixes #861).

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
22#include <math.h>
23#include "crc.h"
24
25#include "stat.hpp"
26#include "random.hpp"
27#include "values.hpp"
28#include "vars.hpp"
29#include "stladdon.hpp"
30#include "domain.hpp"
31#include "examplegen.hpp"
32#include "examples.hpp"
33
34
35#include "distvars.ppp"
36
37DEFINE_TOrangeVector_classDescription(PDistribution, "TDistributionList", true, ORANGE_API)
38
39
40#define CHECKVALTYPE(valType) \
41 if (! (   (valType==TValue::INTVAR) && supportsDiscrete \
42        || (valType==TValue::FLOATVAR) && supportsContinuous)) \
43   raiseError("invalid value type");
44
45#define NOT_IMPLEMENTED(x) { raiseError("'%s' is not implemented", x); throw 0; /*just to avoid warnings*/ }
46
47
48TDistribution::TDistribution()
49: unknowns(0.0),
50  abs(0.0),
51  cases(0.0),
52  normalized(false),
53  supportsDiscrete(false),
54  supportsContinuous(false)
55{}
56
57
58TDistribution::TDistribution(PVariable var)
59: variable(var),
60  unknowns(0.0),
61  abs(0.0),
62  cases(0.0),
63  normalized(false),
64  supportsDiscrete(false),
65  supportsContinuous(false)
66{}
67
68
69TDistribution *TDistribution::create(PVariable var)
70{ if (!var)
71    return NULL;
72  if (var->varType==TValue::INTVAR)
73    return mlnew TDiscDistribution(var);
74  if (var->varType==TValue::FLOATVAR)
75    return mlnew TContDistribution(var);
76
77  ::raiseErrorWho("Distribution", "unknown value type");
78  return NULL; // to make compiler happy
79}
80
81
82TDistribution *TDistribution::fromGenerator(PExampleGenerator gen, const int &position, const int &weightID)
83{
84  if (position >= gen->domain->variables->size())
85    ::raiseErrorWho("Distribution", "index %i out of range", position);
86
87  PVariable var = gen->domain->variables->at(position);
88
89  if (var->varType == TValue::INTVAR)
90    return mlnew TDiscDistribution(gen, position, weightID);
91
92  if (var->varType == TValue::FLOATVAR)
93    return mlnew TContDistribution(gen, position, weightID);
94
95  ::raiseErrorWho("Distribution", "unknown value type");
96  return NULL; // to make compiler happy
97}
98
99
100TDistribution *TDistribution::fromGenerator(PExampleGenerator gen, PVariable var, const int &weightID)
101{
102  if (var->varType == TValue::INTVAR)
103    return mlnew TDiscDistribution(gen, var, weightID);
104
105  if (var->varType == TValue::FLOATVAR)
106    return mlnew TContDistribution(gen, var, weightID);
107
108  ::raiseErrorWho("Distribution", "unknown value type");
109  return NULL; // to make compiler happy
110}
111
112
113
114// General
115
116float TDistribution::compatibility(const TSomeValue &) const
117NOT_IMPLEMENTED("compatibility")
118
119bool TDistribution::compatible(const TSomeValue &) const
120NOT_IMPLEMENTED("compatible")
121
122int TDistribution::compare(const TSomeValue &) const
123NOT_IMPLEMENTED("compare")
124
125TDistribution &TDistribution::operator += (const TDistribution &)
126NOT_IMPLEMENTED("+=")
127
128TDistribution &TDistribution::operator -= (const TDistribution &)
129NOT_IMPLEMENTED("-=")
130
131TDistribution &TDistribution::operator *= (const TDistribution &)
132NOT_IMPLEMENTED("*=")
133
134TDistribution &TDistribution::operator *= (const float &)
135NOT_IMPLEMENTED("*=")
136
137
138// Discrete
139
140const float &TDistribution::atint(const int &)
141NOT_IMPLEMENTED("atint(int)")
142
143const float &TDistribution::atint(const int &) const
144NOT_IMPLEMENTED("atint(int)")
145
146void TDistribution::addint(const int &, const float &)
147NOT_IMPLEMENTED("addint(int, float)")
148
149void TDistribution::setint(const int &, const float &)
150NOT_IMPLEMENTED("add(int, float)")
151
152int TDistribution::randomInt()
153NOT_IMPLEMENTED("randomInt()")
154
155int TDistribution::randomInt(const long &)
156NOT_IMPLEMENTED("randomInt(long)")
157
158int TDistribution::highestProbIntIndex() const
159NOT_IMPLEMENTED("highestProbIntIndex()")
160
161int TDistribution::highestProbIntIndex(const long &) const
162NOT_IMPLEMENTED("highestProbIntIndex(int)")
163
164int TDistribution::highestProbIntIndex(const TExample &) const
165NOT_IMPLEMENTED("highestProbIntIndex(TExample)")
166
167float TDistribution::p(const int &) const
168NOT_IMPLEMENTED("p(int)")
169
170int TDistribution::noOfElements() const
171NOT_IMPLEMENTED("noOfElements()")
172
173// Continuous
174
175const float &TDistribution::atfloat(const float &)
176NOT_IMPLEMENTED("atfloat(float)")
177
178const float &TDistribution::atfloat(const float &) const
179NOT_IMPLEMENTED("atfloat(float)")
180
181void TDistribution::addfloat(const float &, const float &)
182NOT_IMPLEMENTED("addfloat(float, float)")
183
184void TDistribution::setfloat(const float &, const float &)
185NOT_IMPLEMENTED("add(float, float)")
186
187float TDistribution::randomFloat()
188NOT_IMPLEMENTED("randomFloat()")
189
190float TDistribution::randomFloat(const long &)
191NOT_IMPLEMENTED("randomFloat(long)")
192
193float TDistribution::highestProbFloatIndex() const
194NOT_IMPLEMENTED("highestProbFloatIndex()")
195
196float TDistribution::average() const
197NOT_IMPLEMENTED("average()")
198
199float TDistribution::dev() const
200NOT_IMPLEMENTED("dev()")
201
202float TDistribution::var() const
203NOT_IMPLEMENTED("dev()")
204
205float TDistribution::error() const
206NOT_IMPLEMENTED("error()")
207
208float TDistribution::percentile(const float &) const
209NOT_IMPLEMENTED("percentile(float)")
210
211float TDistribution::p(const float &) const
212NOT_IMPLEMENTED("p(float)")
213
214
215
216
217TDistribution &TDistribution::operator +=(PDistribution other)
218{ return operator += (other.getReference()); }
219
220
221TDistribution &TDistribution::operator -=(PDistribution other)
222{ return operator -= (other.getReference()); }
223
224TDistribution &TDistribution::operator *=(PDistribution other)
225{ return operator *= (other.getReference()); }
226
227
228
229float TDistribution::operator -  (const TSomeValue &v) const 
230{ return 1-compatibility(v); }
231
232
233float TDistribution::operator || (const TSomeValue &v) const
234{ return 1-compatibility(v); }
235
236
237const float &TDistribution::operator[](const TValue &val) const 
238{ if (val.isSpecial()) {
239    if (variable)
240      raiseError("undefined value of attribute '%s'", variable->get_name().c_str());
241    else
242      raiseError("undefined attribute value");
243  }
244  CHECKVALTYPE(val.varType);
245  return (val.varType==TValue::INTVAR) ? atint(int(val)) : atfloat(float(val));
246}
247
248
249const float &TDistribution::operator[](const TValue &val)
250{ if (val.isSpecial()) {
251    if (variable)
252      raiseError("undefined value of attribute '%s'", variable->get_name().c_str());
253    else
254      raiseError("undefined attribute value");
255  }
256  CHECKVALTYPE(val.varType);
257  return (val.varType==TValue::INTVAR) ? atint(int(val)) : atfloat(float(val));
258}
259
260
261void TDistribution::add(const TValue &val, const float &p)
262{ 
263  if (!val.svalV || !variable || !variable->distributed) {
264    if (val.isSpecial()) {
265      unknowns += p;
266      if (!val.svalV || !val.svalV.is_derived_from(TDistribution))
267        return;
268    }
269    else {
270      CHECKVALTYPE(val.varType);
271      if (val.varType==TValue::INTVAR)
272        addint(val.intV, p);
273      else
274        addfloat(val.floatV, p);
275      return;
276    }
277  }
278
279  if (!val.svalV)
280    unknowns += p;
281
282  const TDiscDistribution *ddist = val.svalV.AS(TDiscDistribution);
283  if (ddist) {
284    if (!supportsDiscrete || variable && ddist->variable && (variable!=ddist->variable))
285      raiseError("invalid value type");
286    int i = 0;
287    const_PITERATE(TDiscDistribution, ddi, ddist)
288      addint(i++, *ddi*p);
289    return;
290  }
291
292  const TContDistribution *cdist = val.svalV.AS(TContDistribution);
293  if (cdist) {
294    if (!supportsContinuous || variable && ddist->variable && (variable!=ddist->variable))
295      raiseError("invalid value type");
296    const_PITERATE(TContDistribution, cdi, cdist)
297      addfloat((*cdi).first, (*cdi).second*p);
298    return;
299  }
300
301  raiseError("invalid value type");
302}
303
304
305void TDistribution::set(const TValue &val, const float &p)
306{ if (!val.isSpecial()) {
307    CHECKVALTYPE(val.varType);
308    if (val.varType==TValue::INTVAR)
309      setint(val.intV, p);
310    else
311      setfloat(val.floatV, p);
312  }
313}
314
315
316TValue TDistribution::highestProbValue() const 
317{ if (supportsDiscrete)
318    return TValue(highestProbIntIndex());
319  else if (supportsContinuous)
320    return TValue(highestProbFloatIndex());
321  else
322    return TValue();
323}
324
325
326TValue TDistribution::highestProbValue(const long &random) const 
327{ if (supportsDiscrete)
328    return TValue(highestProbIntIndex(random));
329  else if (supportsContinuous)
330    return TValue(highestProbFloatIndex());
331  else
332    return TValue();
333}
334
335
336TValue TDistribution::highestProbValue(const TExample &exam) const 
337{ if (supportsDiscrete)
338    return TValue(highestProbIntIndex(exam));
339  else if (supportsContinuous)
340    return TValue(highestProbFloatIndex());
341  else
342    return TValue();
343}
344
345
346TValue TDistribution::randomValue()
347{ if (supportsDiscrete)
348    return TValue(randomInt());
349  else if (supportsContinuous)
350    return TValue(randomFloat());
351  else 
352    return TValue();
353}
354
355
356TValue TDistribution::randomValue(const long &random)
357{ if (supportsDiscrete)
358    return TValue(randomInt(random));
359  else if (supportsContinuous)
360    return TValue(randomFloat(random));
361  else 
362    return TValue();
363}
364
365
366float TDistribution::p(const TValue &val) const
367{ if (val.isSpecial()) {
368    if (variable)
369      raiseError("undefined value of attribute '%s'", variable->get_name().c_str());
370    else
371      raiseError("undefined attribute value");
372  }
373  CHECKVALTYPE(val.varType);
374  return (val.varType==TValue::INTVAR) ? p(int(val)) : p(float(val));
375}
376
377
378
379
380TDiscDistribution::TDiscDistribution() 
381{ supportsDiscrete = true; };
382
383
384TDiscDistribution::TDiscDistribution(PVariable var) 
385: TDistribution(var)
386{ if (var->varType!=TValue::INTVAR)
387     raiseError("attribute '%s' is not discrete", var->get_name().c_str());
388  distribution = vector<float>(var->noOfValues(), 0.0);
389  supportsDiscrete = true;
390}
391
392
393TDiscDistribution::TDiscDistribution(int values, float value) 
394: distribution(vector<float>(values, value))
395{ cases = abs = value*values;
396  supportsDiscrete = true;
397}
398
399
400TDiscDistribution::TDiscDistribution(const vector<float> &f) 
401: distribution(f)
402{ abs = 0.0;
403  for (const_iterator fi(begin()), fe(end()); fi!=fe; abs += *(fi++));
404  cases = abs;
405  supportsDiscrete = true;
406}
407
408
409TDiscDistribution::TDiscDistribution(const float *f, const int &len)
410: distribution(f, f+len)
411{ abs = 0.0;
412  for (const_iterator fi(begin()), fe(end()); fi!=fe; abs += *(fi++));
413  cases = abs;
414  supportsDiscrete = true;
415}
416
417
418TDiscDistribution::TDiscDistribution(PDistribution other) 
419: TDistribution(other.getReference())
420{ supportsDiscrete = true; }
421
422
423TDiscDistribution::TDiscDistribution(PDiscDistribution other) 
424: TDistribution(other.getReference())
425{ supportsDiscrete = true; }
426
427
428TDiscDistribution::TDiscDistribution(PExampleGenerator gen, const int &position, const int &weightID)
429{
430  supportsDiscrete = true;
431
432  if (position >= gen->domain->variables->size())
433    raiseError("index %i out of range", position);
434
435  variable = gen->domain->variables->at(position);
436  if (variable->varType != TValue::INTVAR)
437    raiseError("attribute '%s' is not discrete", variable->get_name().c_str());
438
439  distribution = vector<float>(variable->noOfValues(), 0.0);
440
441  PEITERATE(ei, gen)
442    add((*ei)[position], WEIGHT(*ei));
443}
444
445
446TDiscDistribution::TDiscDistribution(PExampleGenerator gen, PVariable var, const int &weightID)
447: TDistribution(var)
448{
449  supportsDiscrete = true;
450
451  if (variable->varType != TValue::INTVAR)
452    raiseError("attribute '%s' is not discrete", variable->get_name().c_str());
453
454  distribution = vector<float>(var->noOfValues(), 0.0);
455
456  int position = gen->domain->getVarNum(variable, false);
457  if (position != ILLEGAL_INT)
458    PEITERATE(ei, gen)
459      add((*ei)[position], WEIGHT(*ei));
460  else
461    if (variable->getValueFrom)
462      PEITERATE(ei, gen)
463        add(variable->computeValue(*ei), WEIGHT(*ei));
464    else
465      raiseError("attribute '%s' not in domain and cannot be computed", variable->get_name().c_str());
466}
467
468
469const float &TDiscDistribution::atint(const int &v)
470{ int ms = v + 1 - size();
471  if (ms>0) {
472    reserve(v+1);
473    while (ms--)
474      push_back(0.0);
475  }
476  return distribution[v]; 
477}
478
479
480const float &TDiscDistribution::atint(const int &v) const
481{ if (!size())
482    raiseError("empty distribution");
483  if ((v < 0) || (v >= int(size()))) 
484    raiseError("value %i out of range 0-%i", v, size()-1);
485  return at(v); 
486}
487
488
489void TDiscDistribution::addint(const int &v, const float &w)
490{ if ((v<0) || (v>1e6))
491    raiseError("invalid value");
492
493  int ms = v+1 - size();
494  if (ms>0) {
495    reserve(v+1);
496    while (ms--)
497      push_back(0.0);
498  }
499
500  float &val = distribution[v];
501  val += w;
502  abs += w;
503  cases += w;
504  normalized = false;
505}
506
507
508void TDiscDistribution::setint(const int &v, const float &w)
509{ if ((v<0) || (v>1e6))
510    raiseError("invalid value");
511
512  int ms = v+1 - size();
513  if (ms>0) {
514    reserve(v+1);
515    while (ms--)
516      push_back(0.0);
517  }
518
519  float &val=distribution[v];
520  abs += w-val;
521  cases += w-val;
522  val = w;
523  normalized = false;
524}
525
526
527TDistribution &TDiscDistribution::adddist(const TDistribution &other, const float &factor)
528{
529  const TDiscDistribution *mother=dynamic_cast<const TDiscDistribution *>(&other);
530  if (!mother)
531    raiseError("wrong type of distribution for +=");
532
533  int ms = mother->size() - size();
534  if (ms>0) {
535    reserve(mother->size());
536    while (ms--)
537      push_back(0.0);
538  }
539 
540  iterator ti = begin();
541  const_iterator oi = mother->begin(), oe = mother->end();
542  while(oi!=oe)
543    *(ti++) += *(oi++) * factor;
544  abs += mother->abs * factor;
545  cases += mother->cases;
546  unknowns += mother->unknowns;
547  normalized = false;
548  return *this;
549}
550
551
552TDistribution &TDiscDistribution::operator -=(const TDistribution &other)
553{
554  const TDiscDistribution *mother=dynamic_cast<const TDiscDistribution *>(&other);
555  if (!mother)
556    raiseError("wrong type of distribution for -=");
557
558  int ms = mother->size() - size();
559  if (ms>0) {
560    reserve(mother->size());
561    while (ms--)
562      push_back(0.0);
563  }
564 
565  iterator ti = begin();
566  const_iterator oi = mother->begin(), oe = mother->end();
567  while(oi!=oe)
568    *(ti++) -= *(oi++);
569  abs -= mother->abs;
570  cases -= mother->cases;
571  unknowns -= mother->unknowns;
572  normalized = false;
573  return *this;
574}
575
576
577TDistribution &TDiscDistribution::adddist(PDistribution other, const float &factor)
578{ return adddist(other.getReference(), 1.0); }
579
580
581TDistribution &TDiscDistribution::operator +=(const TDistribution &other)
582{ return adddist(other, 1.0); }
583
584
585TDistribution &TDiscDistribution::operator +=(PDistribution other)
586{ return adddist(other.getReference(), 1.0); }
587
588
589TDistribution &TDiscDistribution::operator -=(PDistribution other)
590{ return operator -= (other.getReference()); }
591
592
593
594TDistribution &TDiscDistribution::operator *=(const float &weight)
595{ for(iterator di(begin()); di!=end(); (*(di++)) *= weight);
596  abs *= weight;
597  normalized = false;
598  return *this;
599}
600
601
602TDistribution &TDiscDistribution::operator *=(const TDistribution &other)
603{ 
604  const TDiscDistribution *mother=dynamic_cast<const TDiscDistribution *>(&other);
605  if (!mother)
606    raiseError("wrong type of distribution for *=");
607
608  abs = 0.0;
609  iterator di = begin(), de = end();
610  const_iterator di2 = mother->begin(), de2 = mother->end();
611  while ((di!=de) && (di2!=de2))
612    abs += (*(di++) *= *(di2++));
613
614  if (di!=de)
615    erase(di, de);
616
617  normalized = false;
618  return *this;
619}
620
621
622TDistribution &TDiscDistribution::operator *= (PDistribution other)
623{ return operator *= (other.getReference()); }
624
625
626TDistribution &TDiscDistribution::operator /=(const TDistribution &other)
627{ const TDiscDistribution *mother=dynamic_cast<const TDiscDistribution *>(&other);
628  if (!mother)
629    raiseError("wrong type of distribution for /=");
630
631  abs = 0.0;
632  iterator di = begin(), de = end();
633  const_iterator di2 = mother->begin(), de2 = mother->end();
634  for (; (di!=de) && (di2!=de2); di++, di2++) {
635    if ((-1e-20 < *di2) && (*di2 < 1e-20)) {
636      if ((*di<-1e-20) || (*di>1e-20))
637        raiseError("division by zero in /=");
638    }
639    else
640      abs += (*di /= *di2);
641  }
642
643  if (di!=de)
644    erase(di, de);
645
646  normalized = false;
647  return *this;
648}
649
650
651TDistribution &TDiscDistribution::operator /= (PDistribution other)
652{ return operator /= (other.getReference()); }
653
654
655TDistribution &TDiscDistribution::mul(const TDistribution &other, const float &weight)
656{ const TDiscDistribution *mother=dynamic_cast<const TDiscDistribution *>(&other);
657  if (!mother)
658    raiseError("wrong type of distribution for -=");
659
660  abs = 0.0;
661  iterator di = begin(), de = end();
662  const_iterator di2 = mother->begin(), de2 = mother->end();
663  while ((di!=de) && (di2!=de2))
664    abs += (*(di++) *= weight * *(di2++));
665
666  if (di!=de)
667    erase(di, de);
668
669  normalized = false;
670  return *this;
671}
672
673
674TDistribution &TDiscDistribution::mul(PDistribution other, const float &weight)
675{ return mul(other.getReference(), weight); }
676
677
678/*  Returns normalized scalar products of distributions of 'other' and 'this'.
679    The result corresponds to a probability that two random values chosen according
680    to the given distributions are same. */
681float TDiscDistribution::compatibility(const TSomeValue &ot) const
682{ const TDiscDistribution *dv=dynamic_cast<const TDiscDistribution *>(&ot);
683  if (dv) {
684    float sum=0;
685    for(const_iterator i1=begin(), i2=dv->begin();
686        (i1!=end());
687        sum += *(i1++) * *(i2++))
688    return sum/abs/dv->abs;
689  }
690
691  const TValue *vv=dynamic_cast<const TValue *>(&ot);
692  if (   (vv) 
693      || (vv->varType==TValue::INTVAR))
694    return (vv->intV>int(size())) ? 0.0 : operator[](vv->intV)/abs;
695     
696  raiseError("can't compare values of different types");
697  return 0.0; // to make compilers happy
698}
699
700
701/*  Declared only since it is abstract in TSomeValue.
702    Definition is somewhat artificial: compare does a lexicographical comparison of probabilities. */
703int  TDiscDistribution::compare(const TSomeValue &ot) const
704{ const TDiscDistribution *dv=dynamic_cast<const TDiscDistribution *>(&ot);
705  if (!dv)
706    raiseError("can't compare values of different types");
707
708  const_iterator i1=begin(), i2=dv->begin();
709  for( ; (i1!=end()) && (*i1==*i2); i1++, i2++);
710  if (i1==end())
711    return 0;
712  else 
713    if (*i1<*i2)
714      return -1;
715  return 1;
716}
717
718
719/*  Declared only since it is abstract in TSomeValue.
720    Definitions is somewhat artificial: compatible returns true if compatibility>0 (i.e. if there
721    is a non-xero probability that a random values with given distributions are same). */
722bool  TDiscDistribution::compatible (const TSomeValue &ot) const
723{ return (compatibility(ot)>0); }
724
725
726void TDiscDistribution::normalize()
727{ if (!normalized) {
728    if (abs) {
729      this_ITERATE(dvi)
730        *dvi /= abs;
731      abs=1.0;
732    }
733    else 
734      if (size()) {
735        float p = 1.0/float(size());
736        this_ITERATE(dvi)
737          *dvi = p;
738        abs = 1.0;
739      }
740   normalized = true;
741  }
742}
743
744
745int TDiscDistribution::highestProbIntIndex() const
746{
747  if (!size())
748    return 0;
749
750  int wins = 1;
751  int best = 0;
752  float bestP = operator[](0);
753  int i, e;
754
755  unsigned long crc;
756  INIT_CRC(crc);
757
758  for(i = 1, e = int(size()); --e; i++) {
759    const float &P = operator[](i);
760    add_CRC(P, crc);
761
762    if (P > bestP) {
763      best = i;
764      bestP = P;
765      wins = 1;
766    }
767    else if (P==bestP)
768      wins++;
769  }
770
771  if (wins==1)
772    return best;
773
774  FINISH_CRC(crc);
775  crc &= 0x7fffffff;
776
777  for(i = 0, wins = 1 + crc % wins; wins; i++)
778    if (operator[](i)==bestP)
779      wins--;
780
781  return i-1;
782}
783
784
785int TDiscDistribution::highestProbIntIndex(const long &random) const
786{
787  if (!size())
788    return 0;
789
790  int wins = 1;
791  int best = 0;
792  float bestP = operator[](0);
793  int i, e;
794
795  for(i = 1, e = int(size()); --e; i++)
796    if (operator[](i) > bestP) {
797      best = i;
798      bestP = operator[](i);
799      wins = 1;
800    }
801    else if (operator[](i)==bestP)
802      wins++;
803
804  if (wins==1)
805    return best;
806
807  for(i = 0, wins = 1 + random % wins; wins; i++)
808    if (operator[](i)==bestP)
809      wins--;
810
811  return i-1;
812}
813
814
815int TDiscDistribution::highestProbIntIndex(const TExample &exam) const
816{
817  if (!size())
818    return 0;
819
820  int wins = 1;
821  int best = 0;
822  float bestP = operator[](0);
823  int i, e;
824
825  for(i = 1, e = int(size()); --e; i++)
826    if (operator[](i) > bestP) {
827      best = i;
828      bestP = operator[](i);
829      wins = 1;
830    }
831    else if (operator[](i)==bestP)
832      wins++;
833
834  if (wins==1)
835    return best;
836
837  wins = 1 + exam.sumValues() % wins;
838
839  i = 0;   
840  while (wins)
841    if (operator[](i++)==bestP)
842      wins--;
843
844  return i-1;
845}
846
847
848float TDiscDistribution::highestProb() const
849{
850  float best=-1;
851  for(int i=0, isize = size(); i<isize; i++)
852    if (operator[](i) > best)
853      best=i;
854  if (best>=0)
855    return operator[](best);
856  else
857    return size() ? 1.0/size() : 0.0;
858}
859
860
861bool TDiscDistribution::noDeviation() const
862{ const_this_ITERATE(dvi)
863    if (*dvi)
864      return *dvi == abs;
865  return size()==1;
866}
867 
868
869int TDiscDistribution::randomInt(const long &random)
870{ 
871  float ri = (random & 0x7fffffff) / float(0x7fffffff);
872  if (!abs || !size())
873    raiseError("cannot return a random element of an empty distribution");
874  ri = fmod(ri, abs);
875  const_iterator di(begin());
876  while (ri > *di)
877    ri -= *(di++);
878  return int(di-begin());
879}
880
881
882int TDiscDistribution::randomInt()
883{ 
884  if (!abs || !size())
885    raiseError("cannot return a random element of an empty distribution");
886
887  if (!randomGenerator)
888    randomGenerator = mlnew TRandomGenerator;
889
890  float ri = randomGenerator->randfloat(abs);
891  const_iterator di(begin());
892  while (ri > *di)
893    ri -= *(di++);
894  return int(di-begin());
895}
896
897
898float TDiscDistribution::p(const int &x) const
899{ if (!abs) 
900    return size() ? 1.0/size() : 0.0;
901  if (x>=size())
902    return 0.0;
903  return atint(x)/abs; 
904}
905
906int TDiscDistribution::noOfElements() const
907{ return size(); }
908
909
910int TDiscDistribution::sumValues() const
911{ unsigned long crc;
912  INIT_CRC(crc);
913
914  const_this_ITERATE(dvi)
915      add_CRC(*dvi, crc);
916
917  FINISH_CRC(crc);
918  return int(crc & 0x7fffffff);
919}
920
921
922TContDistribution::TContDistribution()
923: sum(0.0),
924  sum2(0.0)
925{ supportsContinuous = true; }
926
927
928TContDistribution::TContDistribution(const map<float, float> &dist)
929: distribution(dist), 
930  sum(0.0),
931  sum2(0.0)
932{ abs = 0.0;
933  this_ITERATE(di) {
934    abs+=(*di).second;
935    sum+=(*di).second*(*di).first;
936    sum2+=(*di).second*(*di).first*(*di).first;
937  }
938  cases = abs;
939  supportsContinuous = true; 
940}
941
942
943TContDistribution::TContDistribution(PVariable var) 
944: TDistribution(var),
945  sum(0.0),
946  sum2(0.0)
947{ if (var->varType!=TValue::FLOATVAR)
948     raiseError("attribute '%s' is not continuous", var->get_name().c_str());
949  supportsContinuous = true; 
950}
951
952
953TContDistribution::TContDistribution(PExampleGenerator gen, const int &position, const int &weightID)
954: sum(0.0),
955  sum2(0.0)
956{
957  supportsContinuous = true;
958
959  if (position >= gen->domain->variables->size())
960    raiseError("index %i out of range", position);
961
962  variable = gen->domain->variables->at(position);
963  if (variable->varType != TValue::FLOATVAR)
964    raiseError("attribute '%s' is not continuous", variable->get_name().c_str());
965
966  PEITERATE(ei, gen)
967    add((*ei)[position], WEIGHT(*ei));
968}
969
970
971TContDistribution::TContDistribution(PExampleGenerator gen, PVariable var, const int &weightID)
972: TDistribution(var),
973  sum(0.0),
974  sum2(0.0)
975{
976  supportsContinuous = true;
977
978  if (variable->varType != TValue::FLOATVAR)
979    raiseError("attribute '%s' is not continuous", variable->get_name().c_str());
980
981  int position = gen->domain->getVarNum(variable, false);
982  if (position != ILLEGAL_INT)
983    PEITERATE(ei, gen)
984      add((*ei)[position], WEIGHT(*ei));
985  else
986    if (variable->getValueFrom)
987      PEITERATE(ei, gen)
988        add(variable->computeValue(*ei), WEIGHT(*ei));
989    else
990      raiseError("attribute '%s' not in domain and cannot be computed", variable->get_name().c_str());
991}
992
993
994const float &TContDistribution::atfloat(const float &v)
995{ if (find(v)!=end())
996    distribution[v]=0;
997  return distribution[v]; 
998}
999
1000
1001const float &TContDistribution::atfloat(const float &v) const
1002{ const_iterator vi=find(v);
1003  if (vi==end())
1004    raiseError("value %5.3f does not exist", v);
1005  return (*vi).second;
1006}
1007
1008
1009void TContDistribution::addfloat(const float &v, const float &w)
1010{ 
1011  iterator vi=find(v);
1012  if (vi==end())
1013    distribution[v]=w;
1014  else
1015    (*vi).second+=w;
1016
1017  abs += w;
1018  cases += w;
1019  sum += w * v;
1020  sum2 += w * v*v;
1021  normalized = false;
1022}
1023
1024
1025void TContDistribution::setfloat(const float &v, const float &w)
1026{ 
1027  iterator vi=find(v);
1028  if (vi==end()) {
1029    distribution[v]=w;
1030    abs += w;
1031    cases += w;
1032    sum += w * v;
1033    sum += w * v*v;
1034  }
1035  else {
1036    float dif = w - (*vi).second;
1037    abs += dif;
1038    cases += w;
1039    sum += dif * v;
1040    sum2 += dif * v*v;
1041    (*vi).second += w;
1042  }
1043 
1044  normalized = false;
1045}
1046
1047
1048TDistribution &TContDistribution::operator +=(const TDistribution &other)
1049{
1050  const TContDistribution *mother = dynamic_cast<const TContDistribution *>(&other);
1051  if (!mother)
1052    raiseError("wrong distribution type for +=");
1053
1054  const_PITERATE(TContDistribution, oi, mother) 
1055    addfloat((*oi).first, (*oi).second);
1056
1057  unknowns += mother->unknowns;
1058
1059  return *this;
1060}
1061
1062
1063TDistribution &TContDistribution::operator -=(const TDistribution &other)
1064{
1065  const TContDistribution *mother = dynamic_cast<const TContDistribution *>(&other);
1066  if (!mother)
1067    raiseError("wrong distribution type for -=");
1068
1069  const_PITERATE(TContDistribution, oi, mother) 
1070    addfloat((*oi).first, -(*oi).second);
1071
1072  unknowns -= mother->unknowns;
1073
1074  return *this;
1075}
1076
1077
1078TDistribution &TContDistribution::operator +=(PDistribution other)
1079{ return operator += (other.getReference()); }
1080
1081
1082TDistribution &TContDistribution::operator -=(PDistribution other)
1083{ return operator -= (other.getReference()); }
1084
1085
1086
1087TDistribution &TContDistribution::operator *=(const float &weight)
1088{ for(iterator i(begin()), e(end()); i!=e; (*(i++)).second*=weight);
1089  abs *= weight;
1090  sum *= weight;
1091  sum2 *= weight;
1092  normalized = false;
1093  return *this;
1094}
1095
1096
1097float TContDistribution::highestProbFloatIndex() const
1098{
1099  // Could use sumValues here, but it's too expensive; this should work for distributions that are distributed enough
1100  long sum = 0;
1101  { const_this_ITERATE(i)
1102      sum += *(long *)(&(*i).first) + *(long *)(&(*i).second);
1103  }
1104
1105  TSimpleRandomGenerator rg(sum);
1106
1107  int wins=0;
1108  const_iterator best;
1109  const_this_ITERATE(i)
1110    if (   (wins==0) && ((wins=1)==1)
1111        || ((*i).second >  (*best).second) && ((wins=1)==1)
1112        || ((*i).second == (*best).second) && rg.randbool(++wins))
1113      best = i;
1114
1115  if (!wins)
1116    raiseError("cannot compute the modus of an empty distribution");
1117
1118  return (*best).first;
1119}
1120
1121
1122float TContDistribution::highestProb() const
1123{
1124  long sum = 0;
1125  { const_this_ITERATE(i)
1126      sum += *(long *)(&(*i).first) + *(long *)(&(*i).second);
1127   }
1128
1129  TSimpleRandomGenerator rg(sum);
1130
1131  int wins=0;
1132  const_iterator best;
1133  const_this_ITERATE(i)
1134    if (   (wins==0) && ((wins=1)==1)
1135        || ((*i).second >  (*best).second) && ((wins=1)==1)
1136        || ((*i).second == (*best).second) && rg.randbool(++wins))
1137      best = i;
1138
1139  if (wins)
1140    return (*best).second;
1141  else
1142    return size() ? 1.0/size() : 0.0;
1143}
1144
1145
1146bool TContDistribution::noDeviation() const
1147{ return size()==1;
1148}
1149
1150
1151float TContDistribution::average() const
1152{ if (!abs)
1153    if (variable)
1154      raiseError("cannot compute average ('%s' has no defined values)", variable->get_name().c_str());
1155    else
1156      raiseError("cannot compute average (attribute has no defined values)");
1157
1158  return sum/abs ; 
1159}
1160
1161
1162float TContDistribution::dev() const
1163{ 
1164  if (abs<=1e-7)
1165    if (variable)
1166      raiseError("cannot compute standard deviation ('%s' has no defined values)", variable->get_name().c_str());
1167    else
1168      raiseError("cannot compute standard deviation (attribute has no defined values)");
1169
1170  const float res = sqrt((sum2-sum*sum/abs)/abs);
1171  return res > 0 ? res : 0.0;
1172}
1173 
1174float TContDistribution::var() const
1175{
1176  if (!abs)
1177    if (variable)
1178      raiseError("cannot compute variance ('%s' has no defined values)", variable->get_name().c_str());
1179    else
1180      raiseError("cannot compute variance (attribute has no defined values)");
1181
1182  const float res = (sum2-sum*sum/abs)/abs;
1183  return res > 0 ? res : 0.0;
1184}
1185 
1186float TContDistribution::error() const
1187{ if (abs <= 1.0)
1188    return 0.0;
1189  const float res = sqrt((sum2-sum*sum/abs)/(abs-1) / abs);
1190  return res > 0 ? res : 0.0;
1191}
1192
1193
1194float TContDistribution::percentile(const float &perc) const
1195{ if ((perc<0) || (perc>100))
1196    raiseError("invalid percentile");
1197
1198  if (!size())
1199    raiseError("empty distribution");
1200
1201  if (perc==0.0)
1202    return (*begin()).first;
1203 
1204  if (perc==100.0) {
1205    const_iterator li(end());
1206    return (*--li).first;
1207  }
1208
1209  float togo = abs*perc/100.0;
1210  const_iterator ths(begin()), prev, ee(end());
1211
1212  if (ths == ee)
1213    raiseError("empty distribution");
1214
1215  while ((ths != ee) && (togo > 0)) {
1216    togo -= (*ths).second;
1217    prev = ths;
1218    ths++;
1219  }
1220
1221  if ((togo < 0) || (ths == ee))
1222    return (*prev).first;
1223
1224  // togo==0.0 && ths!=ee
1225  return ((*prev).first + (*ths).first) / 2.0;
1226}
1227
1228
1229void TContDistribution::normalize()
1230{ if (!normalized) {
1231    if (abs) {
1232      this_ITERATE(dvi)
1233        (*dvi).second /= abs;
1234      sum /= abs;
1235      sum2 /= abs;
1236      abs = 1.0;
1237    }
1238    else if (size()) {
1239      float p = 1.0/float(size());
1240      sum = 0.0;
1241      sum2 = 0.0;
1242      this_ITERATE(dvi) {
1243        (*dvi).second = p;
1244        sum += (*dvi).first;
1245        sum2 += sqr((*dvi).first);
1246      }
1247      sum /= abs;
1248      sum2 /= abs;
1249      abs = 1.0;
1250    }
1251
1252    normalized = true;
1253  }
1254}
1255
1256
1257float TContDistribution::randomFloat()
1258{
1259  if (!randomGenerator)
1260    randomGenerator = mlnew TRandomGenerator;
1261
1262  float ri = randomGenerator->randfloat(abs);
1263  const_iterator di(begin());
1264  while (ri > (*di).second)
1265    ri -= (*(di++)).second;
1266  return (*di).first;
1267}
1268
1269
1270float TContDistribution::randomFloat(const long &random)
1271{ 
1272  float ri = (random & 0x7fffffff) / float(0x7fffffff);
1273  const_iterator di(begin());
1274  while (ri > (*di).second)
1275    ri -= (*(di++)).second;
1276  return (*di).first;
1277}
1278
1279
1280float TContDistribution::p(const float &x) const
1281{ const_iterator rb = upper_bound(x);
1282  if (rb==end())
1283    return 0.0;
1284  if ((*rb).first==x)
1285    return (*rb).second;
1286  if (rb==begin())
1287    return 0.0;
1288  const_iterator lb = rb;
1289  lb--;
1290
1291  return (*lb).second + (x - (*lb).first) * ((*rb).second - (*lb).second) / ((*rb).first - (*lb).first);
1292}
1293
1294
1295int TContDistribution::sumValues() const
1296{ unsigned long crc;
1297  INIT_CRC(crc);
1298
1299  const_this_ITERATE(dvi) {
1300    add_CRC((*dvi).first, crc);
1301    add_CRC((*dvi).second, crc);
1302  }
1303
1304  FINISH_CRC(crc);
1305  return int(crc & 0x7fffffff);
1306}
1307
1308
1309TGaussianDistribution::TGaussianDistribution(const float &amean, const float &asigma, const float &anabs)
1310: mean(amean),
1311  sigma(asigma)
1312{
1313  abs = anabs;
1314  normalized = true;
1315  supportsContinuous = true; 
1316}
1317
1318
1319TGaussianDistribution::TGaussianDistribution(PDistribution dist)
1320: mean(dist->average()),
1321  sigma(sqrt(dist->dev()))
1322{
1323 abs = dist->abs;
1324 normalized = true; 
1325 supportsContinuous = true; 
1326}
1327
1328
1329
1330float TGaussianDistribution::average() const
1331{ return mean; }
1332
1333
1334float TGaussianDistribution::var() const
1335{ return sigma*sigma; }
1336 
1337
1338float TGaussianDistribution::dev() const
1339{ return sigma; }
1340 
1341
1342float TGaussianDistribution::error() const
1343{ return sigma; }
1344 
1345
1346void TGaussianDistribution::normalize()
1347{ abs = 1.0; }
1348
1349
1350float TGaussianDistribution::highestProbFloatIndex() const
1351{ return mean; }
1352
1353
1354#define pi 3.1415926535897931
1355
1356float TGaussianDistribution::highestProb() const
1357{ return abs * 1/(sigma * sqrt(2*pi)); }
1358
1359
1360float TGaussianDistribution::randomFloat()
1361{ 
1362  if (!randomGenerator)
1363    randomGenerator = mlnew TRandomGenerator;
1364
1365  return (float)gasdev((double)mean, (double)sigma, randomGenerator.getReference());
1366}
1367
1368
1369float TGaussianDistribution::randomFloat(const long &random)
1370{ 
1371  TRandomGenerator rg(random);
1372  return (float)gasdev((double)mean, (double)sigma, rg);
1373}
1374
1375
1376float TGaussianDistribution::p(const float &x) const
1377{ return abs * exp(-sqr((x-mean)/2/sigma)) / (sigma*sqrt(2*pi)); }
1378
1379
1380bool TGaussianDistribution::noDeviation() const
1381{ return sigma==0.0; }
1382
1383
1384int TGaussianDistribution::sumValues() const
1385{ unsigned long crc;
1386  INIT_CRC(crc);
1387  add_CRC(mean, crc);
1388  add_CRC(sigma, crc);
1389  FINISH_CRC(crc);
1390  return int(crc & 0x7fffffff);
1391}
1392
1393
1394TDomainDistributions::TDomainDistributions()
1395{}
1396
1397
1398TDomainDistributions::TDomainDistributions(PExampleGenerator gen, const long weightID, bool skipDiscrete, bool skipContinuous)
1399{
1400  reserve(gen->domain->variables->size());
1401  PITERATE(TVarList, vi, gen->domain->variables) {
1402    bool compute;
1403    if ((*vi)->varType == TValue::INTVAR)
1404      compute = !skipDiscrete;
1405    else if ((*vi)->varType == TValue::FLOATVAR)
1406      compute = !skipContinuous;
1407    else
1408      compute = false;
1409   
1410    push_back(compute ? TDistribution::create(*vi) : PDistribution());
1411  }
1412
1413  for(TExampleIterator fi(gen->begin()); fi; ++fi) {
1414    TExample::iterator ei=(*fi).begin();
1415    float weight=WEIGHT(*fi);
1416    for(iterator di=begin(); di!=end(); di++, ei++)
1417      if (*di)
1418        (*di)->add(*ei, weight);
1419  }
1420}
1421
1422
1423void TDomainDistributions::normalize()
1424{ this_ITERATE(di)
1425    (*di)->normalize(); 
1426}
1427
1428
1429PDistribution getClassDistribution(PExampleGenerator gen, const long &weightID)
1430{ if (!gen)
1431    raiseErrorWho("getClassDistribution", "no examples");
1432
1433  if (!gen->domain || !gen->domain->classVar)
1434    raiseErrorWho("getClassDistribution", "invalid example generator or class-less domain");
1435
1436  PDistribution classDist = TDistribution::create(gen->domain->classVar);
1437  TDistribution *uclassdist = const_cast<TDistribution *>(classDist.getUnwrappedPtr());
1438  PEITERATE(ei, gen)
1439    uclassdist->add((*ei).getClass(), WEIGHT(*ei));
1440  return classDist;
1441}
1442
1443#undef NOT_IMPLEMENTED
1444#undef CHECKVALTYPE
Note: See TracBrowser for help on using the repository browser.