PLearn 0.1
RandomVar.cc
Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // PLearn (A C++ Machine Learning Library)
00004 // Copyright (C) 1998 Pascal Vincent
00005 // Copyright (C) 1999-2002 Pascal Vincent, Yoshua Bengio and University of Montreal
00006 //
00007 
00008 // Redistribution and use in source and binary forms, with or without
00009 // modification, are permitted provided that the following conditions are met:
00010 // 
00011 //  1. Redistributions of source code must retain the above copyright
00012 //     notice, this list of conditions and the following disclaimer.
00013 // 
00014 //  2. Redistributions in binary form must reproduce the above copyright
00015 //     notice, this list of conditions and the following disclaimer in the
00016 //     documentation and/or other materials provided with the distribution.
00017 // 
00018 //  3. The name of the authors may not be used to endorse or promote
00019 //     products derived from this software without specific prior written
00020 //     permission.
00021 // 
00022 // THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR
00023 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
00024 // OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN
00025 // NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
00026 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
00027 // TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00028 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00029 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00030 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00031 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00032 // 
00033 // This file is part of the PLearn library. For more information on the PLearn
00034 // library, go to the PLearn Web site at www.plearn.org
00035 
00036 
00037  
00038 
00039 /* *******************************************************      
00040  * $Id: RandomVar.cc 3994 2005-08-25 13:35:03Z chapados $
00041  * AUTHORS: Pascal Vincent & Yoshua Bengio
00042  * This file is part of the PLearn library.
00043  ******************************************************* */
00044 
00045 
00047 #include <plearn/var/AbsVariable.h>
00048 #include <plearn/var/ConcatColumnsVariable.h>
00049 #include <plearn/var/ConcatRowsVariable.h>
00050 #include <plearn/var/DeterminantVariable.h>
00051 #include <plearn/var/EqualVariable.h>
00052 #include <plearn/var/ExpVariable.h>
00053 #include <plearn/var/ExtendedVariable.h>
00054 #include <plearn/var/LeftPseudoInverseVariable.h>
00055 #include <plearn/var/LogSumVariable.h>
00056 #include <plearn/var/LogVariable.h>
00057 #include <plearn/var/MatrixSumOfVariable.h>
00058 #include <plearn/var/PowVariable.h>
00059 #include <plearn/var/ProductVariable.h>
00060 #include <plearn/var/RightPseudoInverseVariable.h>
00061 #include <plearn/var/SoftmaxVariable.h>
00062 #include <plearn/var/SumOfVariable.h>
00063 #include <plearn/var/SumVariable.h>
00064 
00065 #include <plearn/base/general.h>
00066 #include "RandomVar.h"
00067 //#include "NaryVariable.h"
00068 //#include "ConjugateGradientOptimizer.h" // Not in the PLearn CVS repository.
00069 #include <plearn/math/plapack.h>
00070 #include <plearn/var/Var_operators.h>
00071 //#include <cmath>
00072 
00073 namespace PLearn {
00074 using namespace std;
00075 
00076 RandomVar::RandomVar() 
00077     :PP<RandomVariable>(new NonRandomVariable(1)) {}
00078 RandomVar::RandomVar(RandomVariable* v) 
00079     :PP<RandomVariable>(v) {}
00080 RandomVar::RandomVar(const RandomVar& other) 
00081     :PP<RandomVariable>(other) {}
00082 RandomVar::RandomVar(int length, int width)
00083     :PP<RandomVariable>(new NonRandomVariable(length,width)) {}
00084 
00085 RandomVar::RandomVar(const Vec& vec) 
00086     :PP<RandomVariable>(new NonRandomVariable(Var(vec))) {}
00087 RandomVar::RandomVar(const Mat& mat) 
00088     :PP<RandomVariable>(new NonRandomVariable(Var(mat))) {}
00089 RandomVar::RandomVar(const Var& v) 
00090     :PP<RandomVariable>(new NonRandomVariable(v)) {}
00091 
00092 RandomVar::RandomVar(const RVArray& rvars) 
00093     :PP<RandomVariable>(new JointRandomVariable(rvars)) {}
00094 
00095 RandomVar RandomVar::operator[](RandomVar index)
00096 { return new RandomElementOfRandomVariable((*this),index); }
00097 
00098 void RandomVar::operator=(const RVArray& vars) 
00099 { *this = RandomVar(vars); }
00100 
00101 void RandomVar::operator=(real f)
00102 {
00103     if (!(*this)->isNonRandom())
00104         PLERROR("RandomVar: can't assign values to a truly random RV");
00105     (*this)->value->value.fill(f);
00106 }
00107 
00108 void RandomVar::operator=(const Vec& v)
00109 {
00110     if (!(*this)->isNonRandom())
00111         PLERROR("RandomVar: can't assign a values to a truly random RV");
00112     (*this)->value->value << v;
00113 }
00114 
00115 void RandomVar::operator=(const Mat& m)
00116 { 
00117     if (!(*this)->isNonRandom())
00118         PLERROR("RandomVar: can't assign a Var to a truly random RV");
00119     Var v(m);
00120     (*this)->value = v;
00121 }
00122 
00123 void RandomVar::operator=(const Var& v)
00124 {
00125     if (!(*this)->isNonRandom())
00126         PLERROR("RandomVar: can't assign a Var to a truly random RV");
00127     (*this)->value = v;
00128 }
00129 
00130 RVInstance RandomVar::operator==(const Var& v) const
00131 {
00132     return RVInstance(*this,v);
00133 }
00134 
00135 // make an array of RV's
00136 RVArray RandomVar::operator&(const RandomVar& v) const
00137 {
00138     return RVArray(*this,v,10);
00139 }
00140 
00141 ConditionalExpression RandomVar::operator|(RVArray a) const
00142 {
00143     RVInstanceArray rvia(a.size());
00144     for (int i=0;i<a.size();i++)
00145     {
00146         rvia[i].V = a[i];
00147         rvia[i].v = a[i]->value;
00148     }
00149     return ConditionalExpression(RVInstance(*this,Var((*this)->length())),rvia);
00150 }
00151   
00152 ConditionalExpression RandomVar::operator|(RVInstanceArray rhs) const
00153 {
00154     return ConditionalExpression(RVInstance(*this,Var((*this)->length())),rhs);
00155 }
00156 
00157 #if 0
00158 
00159 RandomVar RandomVar::operator[](RandomVar index)
00160 { return matRandomVarElement(*this,index); }
00161 
00162 Vec RandomVar::operator[](int i)
00163 { return matRandomVarElement(*this,i); }
00164 
00165 RandomVar RandomVar::operator()(RandomVar i, RandomVar j)
00166 {
00167     return 
00168         new RandomElementOfRandomVariable(*this,i*((real)((*this)->value->matValue.width()))+j);
00169 }
00170 
00171 real& RandomVar::operator()(int i, int j)
00172 {
00173     return new RandomVarElement(*this,i*((real)((*this)->value->matValue.width()))+j);
00174 }
00175 
00176 #endif
00177 
00178 
00179 
00182 int RandomVariable::rv_counter = 0;
00183 
00184 RandomVariable::RandomVariable(int thelength, int thewidth)
00185     :rv_number(rv_counter++), value(thelength,thewidth),marked(false), 
00186      EMmark(false), pmark(false), learn_the_parameters(0)
00187 {
00188 }
00189 
00190 RandomVariable::RandomVariable(const Vec& the_value)
00191     :rv_number(rv_counter++), value(the_value),marked(false), EMmark(false)
00192     , pmark(false), learn_the_parameters(0)
00193 {
00194 }
00195 
00196 RandomVariable::RandomVariable(const Mat& the_value)
00197     :rv_number(rv_counter++), value(the_value),marked(false), EMmark(false),
00198      pmark(false), learn_the_parameters(0)
00199 {
00200 }
00201 
00202 RandomVariable::RandomVariable(const Var& the_value)
00203     :rv_number(rv_counter++), value(the_value),marked(false), EMmark(false),
00204      pmark(false), learn_the_parameters(0)
00205 {
00206 }
00207 
00208 RandomVariable::RandomVariable(const RVArray& the_parents, int thelength)
00209     :rv_number(rv_counter++),
00210      parents(the_parents), value(thelength), marked(false), EMmark(false),
00211      pmark(false), learn_the_parameters(new bool[the_parents.size()])
00212 {
00213 }
00214 
00215 RandomVariable::RandomVariable(const RVArray& the_parents, int thelength, 
00216                                int thewidth)
00217     : rv_number(rv_counter++),
00218       parents(the_parents), value(thelength,thewidth),marked(false), 
00219       EMmark(false), pmark(false),
00220       learn_the_parameters(new bool[the_parents.size()])
00221 {
00222 }
00223 
00224 
00225 RandomVar RandomVariable::subVec(int start, int length) { 
00226     return new SubVecRandomVariable(this,start,length);
00227 }
00228 
00229 real RandomVariable::EM(const RVArray& parameters_to_learn,
00230                         VarArray& prop_path, VarArray& observedVars, VMat distr, int n_samples, 
00231                         int max_n_iterations, real relative_improvement_threshold,
00232                         bool accept_worsening_likelihood)
00233 {
00234     real avgnegloglik = 0;
00235     real previous_nll=FLT_MAX, nll_change;
00236     bool EMfinished= !(max_n_iterations>0);
00237     int n_epochs=0;
00238 
00239     EMTrainingInitialize(parameters_to_learn);
00240     clearEMmarks();
00241     while (!EMfinished) {
00242         avgnegloglik=epoch(prop_path, observedVars, distr,n_samples);
00243         cout << "EM epoch NLL = " << avgnegloglik << endl;
00244         nll_change = (previous_nll - avgnegloglik)/fabs(previous_nll);
00245         if (nll_change < -1e-4  &&  !accept_worsening_likelihood)
00246             printf("%s %s from %f to %f\n", "RandomVariable::EM",
00247                    "An EM epoch yielded worse negative log-likelihood,",
00248                    previous_nll, avgnegloglik);
00249         n_epochs++;
00250         EMfinished = canStopEM() &&
00251             ((n_epochs >= max_n_iterations) ||
00252              (fabs(nll_change) <= relative_improvement_threshold) ||
00253              (!accept_worsening_likelihood &&
00254               nll_change  <= relative_improvement_threshold));
00255         previous_nll=avgnegloglik;
00256     }
00257     return avgnegloglik;
00258 }
00259 
00260 real RandomVariable::epoch(VarArray& prop_path, 
00261                            VarArray& observed_vars, 
00262                            const VMat& distr, int n_samples,
00263                            bool do_EM_learning)
00264 {
00265     real avg_cost = 0;
00266     if (do_EM_learning) 
00267     {
00268         EMEpochInitialize();
00269         clearEMmarks();
00270     }
00271     for (int i=0;i<n_samples;i++)
00272     {
00273         Vec sam(distr->width());
00274         distr->getRow(i,sam);
00275         observed_vars << sam;
00276         prop_path.fprop(); // computes logP in last element of prop_path
00277         Var logp = prop_path.last();
00278 
00279 #if 0
00280         // debugging
00281         cout << "at example i=" << i << endl;
00282         VarArray sources = logp->sources();
00283         logp->unmarkAncestors();
00284         sources.printInfo();
00285         prop_path.printInfo();
00286 #endif
00287         
00288         avg_cost -= logp->value[0];
00289         if (do_EM_learning)
00290             // last = LHS observed value
00291             EMBprop(observed_vars.last()->value,1.0);
00292     }
00293 
00294     if (do_EM_learning) 
00295     {
00296         EMUpdate();
00297         clearEMmarks();
00298     }
00299     return avg_cost / n_samples;
00300 }
00301 
00302 void RandomVariable::unmarkAncestors() 
00303 { 
00304     if (pmark)
00305     {
00306         marked=false; 
00307         pmark=false;
00308         for (int i=0;i<parents.size();i++)
00309             parents[i]->unmarkAncestors();
00310     }
00311 }
00312 
00313 void RandomVariable::clearEMmarks()
00314 { 
00315     if (EMmark)
00316     {
00317         EMmark=false; 
00318         for (int i=0;i<parents.size();i++)
00319             parents[i]->clearEMmarks();
00320     }
00321 }
00322 
00323 void RandomVariable::setKnownValues()
00324 {
00325     if (!pmark && !marked)
00326     {
00327         pmark=true;
00328         bool all_parents_marked=true;
00329         for (int i=0;i<parents.size();i++)
00330         {
00331             parents[i]->setKnownValues();
00332             all_parents_marked &= parents[i]->isMarked();
00333         }
00334         setValueFromParentsValue();
00335         if (all_parents_marked)
00336             marked=true;
00337     }
00338 }
00339 
00340 void RandomVariable::EMUpdate()
00341 {
00342     if (EMmark) return;
00343     EMmark=true;
00344     for (int i=0;i<parents.size();i++)
00345         if (!parents[i]->isConstant()) parents[i]->EMUpdate();
00346 }
00347 
00348 bool RandomVariable::canStopEM()
00349 {
00350     // propagate to parents
00351     bool can=true;
00352     for (int i=0;i<parents.size() && !can;i++)
00353         can = parents[i]->canStopEM();
00354     return can;
00355 }
00356 
00357 void RandomVariable::
00358 EMTrainingInitialize(const RVArray& parameters_to_learn)
00359 {
00360     if (EMmark) return;
00361     EMmark=true;
00362     int n_can_learn=0;
00363     int n_random=0;
00364     for (int i=0;i<parents.size();i++)
00365     {
00366         if (parameters_to_learn.contains(parents[i]))
00367         {
00368             if (!parents[i]->isConstant())
00369                 PLERROR("Trying to learn a parameter that is not constant!");
00370             learn_the_parameters[i] = true;
00371             n_can_learn++;
00372         }
00373         else learn_the_parameters[i] = false;
00374         if (!parents[i]->isNonRandom())
00375             n_random++;
00376     }
00377     if (n_can_learn>0 && n_random>0)
00378         PLERROR("RandomVariable: can't learn some parameter if others are random");
00379 
00380     for (int i=0;i<parents.size();i++)
00381         parents[i]->EMTrainingInitialize(parameters_to_learn);
00382 }
00383 
00384 void RandomVariable::EMEpochInitialize()
00385 {
00386     if (EMmark) return;
00387     EMmark=true;
00388     for (int i=0;i<parents.size();i++)
00389         parents[i]->EMEpochInitialize();
00390 }
00391 
00392 Var RandomVariable::P(const Var& obs, const RVInstanceArray& RHS)
00393 { return exp(logP(obs,RHS,0)); }
00394 
00395 
00396 RandomVariable::~RandomVariable() 
00397 { if (learn_the_parameters) delete learn_the_parameters; }
00398 
00399 Var RandomVariable::ElogP(const Var& obs, RVArray& parameters_to_learn, 
00400                           const RVInstanceArray& RHS)
00401 {
00402     PLERROR("ElogP not implemented for this class (%s)",classname());
00403     return Var(0);
00404 }
00405 
00408 RandomVar operator*(RandomVar a, RandomVar b)
00409 {
00410     if (a->isScalar() || b->isScalar()); // scalar times something
00411     else if (a->isVec() && b->isVec()) // vec times vec
00412     {
00413         if (a->length()*a->width() != b->length()*b->width())
00414             PLERROR("In RandomVar operator*(RandomVar a, RandomVar b) cannot do a dot product between 2 vecs with different sizes");
00415     }
00416     else if (a->isRowVec()) // rowvec times mat
00417     {
00418         if (a->length() != b->width())
00419             PLERROR("In RandomVar operator*(RandomVar a, RandomVar b) in case rowvec times mat: a->length() != b->width()");
00420     }
00421     else if (b->isRowVec()) // mat times rowvec
00422     {
00423         if (b->length() != a->width())
00424             PLERROR("In RandomVar operator*(RandomVar a, RandomVar b) in case mat times rowvec: b->length() != a->width()");
00425     }
00426     else
00427         PLERROR("In RandomVar operator*(RandomVar a, RandomVar b) This case is not handled (but maybe it should be...)");
00428 
00429     return new ProductRandomVariable(a,b);
00430 }
00431 
00432 RandomVar operator+(RandomVar a, RandomVar b)
00433 {
00434     return new PlusRandomVariable(a,b);
00435 }
00436 
00437 RandomVar operator-(RandomVar a, RandomVar b)
00438 {
00439     return new MinusRandomVariable(a,b);
00440 }
00441 
00442 RandomVar operator/(RandomVar a, RandomVar b)
00443 {
00444     return new ElementWiseDivisionRandomVariable(a,b);
00445 }
00446 
00447 // exponential function applied element-by-element
00448 RandomVar exp(RandomVar x) { return new ExpRandomVariable(x); }
00449 
00450 // natural logarithm function applied element-by-element
00451 RandomVar log(RandomVar x) { return new LogRandomVariable(x); }
00452 
00453 RandomVar extend(RandomVar v, real extension_value, int n_extend)
00454 { return new ExtendedRandomVariable(v,extension_value,n_extend); }
00455 
00456 RandomVar hconcat(const RVArray& a)
00457 { return new ConcatColumnsRandomVariable(a); }
00458 
00459 
00460 real EM(ConditionalExpression conditional_expression,
00461         RVArray parameters_to_learn,
00462         VMat distr, int n_samples, int max_n_iterations, 
00463         real relative_improvement_threshold,
00464         bool accept_worsening_likelihood,
00465         bool compute_final_train_NLL)
00466 {
00467     // assign the value fields of the RV's to those provided by user
00468     RandomVar& LHS=conditional_expression.LHS.V;
00469     Var& lhs_observation=conditional_expression.LHS.v;
00470     VarArray prop_input_vars;
00472     // NOTE NOTE NOTE:
00473     //
00474     // THE ORDER OF THE VALUES IN THE DISTRIBUTION MUST BE:
00475     // (1) conditioning variables (RHS), (2) output variables (LHS)
00477     VarArray distr_observed_vars = 
00478         conditional_expression.RHS.instances() & (VarArray)lhs_observation;
00479     // note that we don't use LHS->value to put the LHS_observation
00480     // in case some distribution need to compare the observation
00481     // with a function of their parents that is put in their value field.
00482     Var logp = logP(conditional_expression,false);
00483     VarArray prop_path = 
00484         propagationPath(distr_observed_vars & parameters_to_learn.values(), logp);
00485     // do the actual EM training with multiple epochs
00486     real train_NLL =
00487         LHS->EM(parameters_to_learn,prop_path,distr_observed_vars, distr, 
00488                 n_samples, max_n_iterations, relative_improvement_threshold,
00489                 accept_worsening_likelihood);
00490     if (compute_final_train_NLL)
00491         train_NLL = LHS->epoch(prop_path,distr_observed_vars,distr,
00492                                n_samples,false);
00493     LHS->unmarkAncestors();
00494     return train_NLL;
00495 }
00496 
00497 real EM(ConditionalExpression conditional_expression,
00498         RVArray parameters_to_learn,
00499         VMat distr, int n_samples, int max_n_iterations, 
00500         real relative_improvement_threshold,
00501         bool compute_final_train_NLL)
00502 {
00503     // assign the value fields of the RV's to those provided by user
00504     RandomVar& LHS=conditional_expression.LHS.V;
00505     Var& lhs_observation=conditional_expression.LHS.v;
00506     VarArray prop_input_vars;
00508     // NOTE NOTE NOTE:
00509     //
00510     // THE ORDER OF THE VALUES IN THE DISTRIBUTION MUST BE:
00511     // (1) conditioning variables (RHS), (2) output variables (LHS)
00513     VarArray distr_observed_vars = 
00514         conditional_expression.RHS.instances() & (VarArray)lhs_observation;
00515     // note that we don't use LHS->value to put the LHS_observation
00516     // in case some distribution need to compare the observation
00517     // with a function of their parents that is put in their value field.
00518     RVInstanceArray params_to_learn(parameters_to_learn.size());
00519     for (int i=0;i<parameters_to_learn.size();i++)
00520     {
00521         params_to_learn[i].V=parameters_to_learn[i];
00522         // params_to_learn[i].v will hold "new" params in EM epoch
00523         params_to_learn[i].v=Var(parameters_to_learn[i]->length());
00524         // initialize "new" params with current value
00525         params_to_learn[i].v->value << params_to_learn[i].V->value->value;
00526     }
00527     Var elogp = ElogP(conditional_expression,params_to_learn,false);
00528     VarArray new_params = params_to_learn.instances();
00529     VarArray current_params = parameters_to_learn.values();
00530     VarArray prop_path = 
00531         propagationPath(distr_observed_vars & current_params & new_params, elogp);
00532     // do the actual EM training with N epochs, N=#free parameters
00533     // (this is because, assuming that maximizing the auxiliary function
00534     // is solvable analytically, i.e. it is quadratic in the parameters,
00535     // then N iterations of conjugate gradiends should suffice. With
00536     // numerical errors, we can tolerate a bit more...
00537     int n_free_params = new_params.nelems();
00538     int max_n_Q_iterations = 1 + (int)(n_free_params*1.5);
00539     Vec params(n_free_params);
00540     // again, assuming solvable max Q, a specialized but faster CG is enough
00541     Var totalElogP = meanOf(elogp,distr_observed_vars,
00542                             distr,n_samples,new_params);
00543     PLERROR("In EM (RandomVar.cc), code using ConjugateGradientOptimizer is now commented out");
00544     max_n_Q_iterations = max_n_Q_iterations; // TODO Remove this (just to make the compiler happy).
00545     /* REMOVED because not in the PLearn CVS repository.
00546      
00547     ConjugateGradientOptimizer opt(new_params, totalElogP,
00548     0.001,0.001,max_n_Q_iterations);
00549     */
00550 
00551     // the outer loop is over EM iterations
00552     real avgnegloglik = 0;
00553     real previous_nll=FLT_MAX, nll_change;
00554     bool EMfinished= !(max_n_iterations>0);
00555     int n_epochs=0;
00556 
00557     while (!EMfinished) {
00558         // the inner loop is over "new params" optimization of totalElogP 
00559         PLERROR("In EM (RandomVar.cc), code using ConjugateGradientOptimizer is now commented out");
00560 //    opt.optimize(); COMMENTED (same as above)
00561         avgnegloglik = - totalElogP->value[0];
00562         cout << "EM epoch -Q = " << avgnegloglik << endl;
00563         nll_change = (previous_nll - avgnegloglik)/fabs(previous_nll);
00564         if (nll_change < -1e-4)
00565             printf("%s %s from %f to %f\n", "RandomVariable::EM",
00566                    "An EM epoch yielded worse negative log-likelihood,",
00567                    previous_nll, avgnegloglik);
00568         n_epochs++;
00569         EMfinished = 
00570             ((n_epochs >= max_n_iterations) ||
00571              (fabs(nll_change) <= relative_improvement_threshold) ||
00572              nll_change  <= relative_improvement_threshold);
00573         previous_nll=avgnegloglik;
00574 
00575         // copy the "new params" to the "current params"
00576         new_params >> params;
00577         current_params << params;
00578     }
00579 
00580     if (compute_final_train_NLL)
00581     {
00582         Var logp = logP(conditional_expression,false);
00583         Var totalLogP = meanOf(logp,distr_observed_vars,
00584                                distr,n_samples,current_params);
00585         totalLogP->fprop();
00586         avgnegloglik = - totalLogP->value[0];
00587     }
00588 
00589     LHS->unmarkAncestors();
00590     return avgnegloglik;
00591 }
00592 
00593 Var P(ConditionalExpression conditional_expression, 
00594       bool clearMarksUponReturn) 
00595 {
00596     RandomVar& LHS = conditional_expression.LHS.V;
00597     RVInstanceArray& RHS = conditional_expression.RHS;
00598     // traverse the tree of ancestors of this node
00599     // and mark nodes which are deterministic descendents of RHS
00600     // and of non-random variables
00601     // while setting their "value" field to this Var function of them.
00602     LHS->markRHSandSetKnownValues(RHS);
00603 
00604     Var p = LHS->P(conditional_expression.LHS.v,RHS);
00605 
00606     if (clearMarksUponReturn)
00607         // put the network back in its original state
00608         LHS->unmarkAncestors();
00609 
00610     // make sure that all the paths which do not
00611     // depend on x, y, and the tunable_parameters are correctly computed
00612     p->fprop_from_all_sources();
00613     return p;
00614 }
00615 
00616 Var logP(ConditionalExpression conditional_expression, bool clearMarksUponReturn,
00617          RVInstanceArray* parameters_to_learn) 
00618 {
00619     RandomVar& LHS = conditional_expression.LHS.V;
00620     RVInstanceArray& RHS = conditional_expression.RHS;
00621     // traverse the tree of ancestors of this node
00622     // and mark nodes which are deterministic descendents of RHS
00623     // and of non-random variables
00624     // while setting their "value" field to this Var function of them.
00625     LHS->markRHSandSetKnownValues(RHS);
00626 
00627     Var logp = LHS->logP(conditional_expression.LHS.v,RHS,parameters_to_learn);
00628 
00629     if (clearMarksUponReturn)
00630     {
00631         // put the network back in its original state
00632         LHS->unmarkAncestors();
00633         for (int i=0;i<RHS.size();i++) RHS[i].V->unmark();
00634     }
00635 
00636     // make sure that all the paths which do not
00637     // depend on x, y, and the tunable_parameters are correctly computed
00638     logp->fprop_from_all_sources();
00639     return logp;
00640 }
00641 
00642 Var ElogP(ConditionalExpression conditional_expression, 
00643           RVInstanceArray& parameters_to_learn,
00644           bool clearMarksUponReturn)
00645 { 
00646     return logP(conditional_expression,clearMarksUponReturn,&parameters_to_learn);
00647 }
00648 
00649 
00650 // integrate the RV over the given hiddenRV
00651 // and return the resulting new RandomVariable. This 
00652 // may be difficult to do in general...
00653 RandomVar marginalize(const RandomVar& RV, const RandomVar& hiddenRV)
00654 { 
00655     PLERROR("marginalize not implemented yet..."); 
00656     return RandomVar();
00657 }
00658 
00659 // Sample an instance from the given conditional expression,
00660 // of the form (LHS|RHS) where LHS is a RandomVar and
00661 // RHS is a RVInstanceArray, e.g. (X==x && Z==z && W==w).
00662 // THIS IS A VERY INEFFICIENT IMPLEMENTATION IF TO
00663 // BE CALLED MANY TIMES.
00664 Vec sample(ConditionalExpression conditional_expression)
00665 {
00666     Var instance = Sample(conditional_expression);
00667     instance->fprop_from_all_sources();
00668     return instance->value;
00669 }
00670 
00671 // Sample N instances from the given conditional expression,
00672 // of the form (LHS|RHS) where LHS is a RandomVar and
00673 // RHS is a RVInstanceArray, e.g. (X==x && Z==z && W==w).
00674 // Put the N instances in the rows of the given Nxd matrix.
00675 // THIS ALSO SHOWS HOW TO REPEATEDLY SAMPLE IN AN EFFICIENT
00676 // MANNER (rather than call "Vec sample(ConditionalExpression)").
00677 void sample(ConditionalExpression conditional_expression,Mat& samples)
00678 {
00679     if (samples.length()==0) return;
00680     Var instance = Sample(conditional_expression);
00681     instance->fprop_from_all_sources();
00682     samples(0) << instance->value;
00683     if (samples.length()>0)
00684     {
00685         VarArray path;
00686         instance->random_sources().setMark(); // mark the random sources
00687         instance->markPath(); // mark successors of the random sources
00688         instance->buildPath(path); // extract path from the random sources to instance
00689         // and clear marks
00690         for (int i=1;i<samples.length();i++)
00691         {
00692             path.fprop();
00693             samples(i) << instance->value;
00694         }
00695     }
00696 }
00697 
00698 // Return a Var which depends functionally on the RHS instances
00699 // and the value of other RandomVars which are non-random and
00700 // influence the LHS.
00701 Var Sample(ConditionalExpression conditional_expression)
00702 {
00703     RVInstanceArray& RHS = conditional_expression.RHS;
00704     RandomVar& LHS = conditional_expression.LHS.V;
00705     LHS->markRHSandSetKnownValues(RHS);
00706     LHS->unmarkAncestors();
00707     return LHS->value;
00708 }
00709 
00710 // multivariate d-dimensional diagonal normal with NON-RANDOM and CONSTANT 
00711 // parameters (default means = 0, default standard deviations = 1)
00712 RandomVar normal(real mean, real standard_dev, int d,
00713                  real minimum_standard_deviation)
00714 {
00715     RandomVar means(d); 
00716     means->value->value.fill(mean);
00717     RandomVar logvar(d);
00718     real variance = standard_dev*standard_dev-
00719         minimum_standard_deviation*minimum_standard_deviation;
00720     if (variance<=0)
00721         PLERROR("normal: variance should be positive");
00722     logvar->value->value.fill((real)log((double)variance));
00723     return new DiagonalNormalRandomVariable(means,logvar,
00724                                             minimum_standard_deviation);
00725 }
00726 
00727 // diagonal normal with general parameters
00728 // given by the provided RandomVar's
00729 RandomVar normal(RandomVar mean, RandomVar log_variance,
00730                  real minimum_standard_deviation)
00731 {
00732     return new DiagonalNormalRandomVariable(mean,log_variance,
00733                                             minimum_standard_deviation);
00734 }
00735 
00736 RandomVar mixture(RVArray components, RandomVar log_weights)
00737 {
00738     return new MixtureRandomVariable(components,log_weights);
00739 }
00740 
00741 RandomVar multinomial(RandomVar log_probabilities)
00742 {
00743     return new MultinomialRandomVariable(log_probabilities);
00744 }
00745 
00748 ConditionalExpression::
00749 ConditionalExpression(RVInstance lhs, RVInstanceArray rhs) 
00750     :LHS(lhs), RHS(rhs) {}
00751 
00752 ConditionalExpression::
00753 ConditionalExpression(RVInstance lhs)
00754     :LHS(lhs), RHS() {}
00755 
00756 ConditionalExpression::
00757 ConditionalExpression(RandomVar lhs)
00758     :LHS(lhs,Var(lhs->length())), RHS() {}
00759 
00760 // build from multiple LHS RVInstances: make one RVInstance
00761 // from the joint of the RVs and the vconcat of the instances.
00762 ConditionalExpression::ConditionalExpression(RVInstanceArray lhs)
00763     :LHS(lhs.random_variables(),vconcat(lhs.instances())), RHS() {}
00764 
00767 RVInstance::RVInstance(const RandomVar& VV, const Var& vv) :V(VV), v(vv) 
00768 {
00769     if (VV->length()!=vv->length())
00770         PLERROR("Associating a RandomVar of length %d to a Var of length %d",
00771                 VV->length(),vv->length());
00772 }
00773 
00774 RVInstance::RVInstance() {}
00775 
00776 RVInstanceArray RVInstance::operator&&(RVInstance rvi)
00777 {
00778     return RVInstanceArray(*this,rvi);
00779 }
00780 
00781 ConditionalExpression RVInstance::operator|(RVInstanceArray a)
00782 {
00783     return ConditionalExpression(*this,a);
00784 }
00785 
00786 // swap the v with the V->value
00787 void RVInstance::swap_v_and_Vvalue()
00788 { Var tmp = v; v = V->value; V->value = tmp; }
00789 
00790 
00793 RVInstanceArray::RVInstanceArray()
00794     : Array<RVInstance>(0,0)
00795 {}
00796 
00797 RVInstanceArray::RVInstanceArray(int n,int n_extra)
00798     : Array<RVInstance>(n,n_extra)
00799 {}
00800 
00801 RVInstanceArray::RVInstanceArray(const Array<RVInstance>& va)
00802     : Array<RVInstance>(va) {} 
00803 
00804 RVInstanceArray::RVInstanceArray(const RVInstance& v, int n_extra)
00805     : Array<RVInstance>(1,n_extra)
00806 { (*this)[0] = v; }
00807 
00808 RVInstanceArray::RVInstanceArray(const RVInstance& v1, const RVInstance& v2, int n_extra)
00809     : Array<RVInstance>(2,n_extra)
00810 { 
00811     (*this)[0] = v1; 
00812     (*this)[1] = v2; 
00813 }
00814 
00815 RVInstanceArray::RVInstanceArray(const RVInstance& v1, const RVInstance& v2, 
00816                                  const RVInstance& v3, int n_extra)
00817     : Array<RVInstance>(3,n_extra)
00818 { 
00819     (*this)[0] = v1; 
00820     (*this)[1] = v2; 
00821     (*this)[2] = v3; 
00822 }
00823 
00824 int RVInstanceArray::length() const {
00825     int l=0;
00826     for (int i=0;i<size();i++)
00827         l += (*this)[i].V->length();
00828     return l;
00829 }
00830 
00831 VarArray RVInstanceArray::values() const {
00832     VarArray vals(size());
00833     for (int i=0;i<size();i++)
00834         vals[i]=(*this)[i].V->value;
00835     return vals;
00836 }
00837 
00838 VarArray RVInstanceArray::instances() const {
00839     VarArray vals(size());
00840     for (int i=0;i<size();i++)
00841         vals[i]=(*this)[i].v;
00842     return vals;
00843 }
00844 
00845 RVArray RVInstanceArray::random_variables() const {
00846     RVArray vars(size());
00847     for (int i=0;i<size();i++)
00848         vars[i]=(*this)[i].V;
00849     return vars;
00850 }
00851 
00852 RVInstanceArray RVInstanceArray::operator&&(RVInstance rvi)
00853 {
00854     return PLearn::operator&(*this,(RVInstanceArray)rvi);
00855     //return this->operator&((RVInstanceArray)rvi);
00856 }
00857 
00858 ConditionalExpression RVInstanceArray::operator|(RVInstanceArray RHS)
00859 { 
00860     return ConditionalExpression(RVInstance(random_variables(),
00861                                             vconcat(instances())),RHS);
00862 }
00863 
00864 int RVInstanceArray::compareRVnumbers(const RVInstance* rvi1, 
00865                                       const RVInstance* rvi2)
00866 {
00867     return rvi1->V->rv_number - rvi2->V->rv_number;
00868 }
00869 
00870 // sorts in-place the elements by V->rv_number (topological order of 
00871 // the graphical model) (in the order: ancestors -> descendants)
00872 void RVInstanceArray::sort()
00873 {
00874     RVInstance* array = data();
00875     qsort(array,size(),sizeof(RVInstance),(compare_function)compareRVnumbers);
00876 }
00877 
00880 RVArray::RVArray()
00881     : Array<RandomVar>(0,0)
00882 {}
00883 
00884 RVArray::RVArray(int n,int n_extra)
00885     : Array<RandomVar>(n,n_extra)
00886 {}
00887 
00888 RVArray::RVArray(const Array<RandomVar>& va)
00889     : Array<RandomVar>(va) {} 
00890 
00891 RVArray::RVArray(const RandomVar& v, int n_extra)
00892     : Array<RandomVar>(1,n_extra)
00893 { (*this)[0] = v; }
00894 
00895 RVArray::RVArray(const RandomVar& v1, const RandomVar& v2, int n_extra)
00896     : Array<RandomVar>(2,n_extra)
00897 { 
00898     (*this)[0] = v1; 
00899     (*this)[1] = v2; 
00900 }
00901 
00902 RVArray::RVArray(const RandomVar& v1, const RandomVar& v2, const RandomVar& v3,
00903                  int n_extra)
00904     : Array<RandomVar>(3,n_extra)
00905 { 
00906     (*this)[0] = v1; 
00907     (*this)[1] = v2; 
00908     (*this)[2] = v3; 
00909 }
00910 
00911 int RVArray::length() const
00912 {
00913     int l=0;
00914     for (int i=0;i<size();i++)
00915         l += (*this)[i]->length();
00916     return l;
00917 }
00918 
00919 VarArray RVArray::values() const
00920 {
00921     VarArray vals(size());
00922     for (int i=0;i<size();i++)
00923         vals[i]=(*this)[i]->value;
00924     return vals;
00925 }
00926 
00927 RandomVar RVArray::operator[](RandomVar index) 
00928 { return new RVArrayRandomElementRandomVariable(*this, index); }
00929 
00930 int RVArray::compareRVnumbers(const RandomVar* v1, const RandomVar* v2)
00931 {
00932     return (*v1)->rv_number - (*v2)->rv_number;
00933 }
00934 
00935 // sorts in-place the elements by rv_number (topological order of 
00936 // the graphical model) (in the order: ancestors -> descendants)
00937 void RVArray::sort()
00938 {
00939     RandomVar* array = data();
00940     qsort(array,size(),sizeof(RandomVar),(compare_function)compareRVnumbers);
00941 }
00942 
00945 StochasticRandomVariable::StochasticRandomVariable(int length)
00946     :RandomVariable(length) {}
00947 
00948 StochasticRandomVariable::StochasticRandomVariable(const RVArray& parameters,
00949                                                    int length)
00950     :RandomVariable(parameters,length)
00951 {}
00952 
00953 StochasticRandomVariable::StochasticRandomVariable(const RVArray& parameters,
00954                                                    int length, int width)
00955     :RandomVariable(parameters,length,width)
00956 {
00957 }
00958 
00959 void StochasticRandomVariable::setKnownValues()
00960 {
00961     if (!marked && !pmark)
00962     {
00963         pmark=true;
00964         // a StochasticRandomVariable cannot be non-random
00965         // unless it is a "dirac" distribution (i.e., isNonRandom()==true).
00966         for (int i=0;i<parents.size();i++)
00967             parents[i]->setKnownValues();
00968         setValueFromParentsValue();
00969         if (isNonRandom()) marked=true;
00970     }
00971 }
00972 
00975 // these are only used by NonRandomVariable
00976 FunctionalRandomVariable::FunctionalRandomVariable(int thelength)
00977     :RandomVariable(thelength) {}
00978 FunctionalRandomVariable::FunctionalRandomVariable(int thelength, int thewidth)
00979     :RandomVariable(thelength,thewidth) {}
00980 FunctionalRandomVariable::FunctionalRandomVariable(const Vec& the_value)
00981     :RandomVariable(the_value) {}
00982 FunctionalRandomVariable::FunctionalRandomVariable(const Mat& the_value)
00983     :RandomVariable(the_value) {}
00984 FunctionalRandomVariable::FunctionalRandomVariable(const Var& the_value)
00985     :RandomVariable(the_value) {}
00986 
00987 // only used by the other sub-classes
00988 FunctionalRandomVariable::FunctionalRandomVariable(const RVArray& the_parents,
00989                                                    int length)
00990     :RandomVariable(the_parents,length) {}
00991 
00992 FunctionalRandomVariable::FunctionalRandomVariable(const RVArray& the_parents,
00993                                                    int length, int width)
00994     :RandomVariable(the_parents,length,width) {}
00995 
00996 Var FunctionalRandomVariable::logP(const Var& obs, const RVInstanceArray& RHS,
00997                                    RVInstanceArray* parameters_to_learn)
00998 {
00999     // gather the unobserved parents
01000     int np=parents.size();
01001     RVInstanceArray unobserved_parents(0,np);
01002     for (int i=0;i<np;i++)
01003         if (!(parents[i]->isMarked() || parents[i]->isNonRandom()))
01004             unobserved_parents &= RVInstance(parents[i],Var(parents[i]->length()));
01005     // simplest case first
01006     int nup = unobserved_parents.size();
01007     if (nup==0)
01008     {
01009         if (isDiscrete())
01010             return isequal(value,obs);
01011         // else
01012         return isequal(value,obs)*FLT_MAX;
01013     }
01014     // else
01015     Var *JacobianCorrection=0;
01016     if (invertible(obs,unobserved_parents,&JacobianCorrection))
01017     {
01018         Var logp(1);
01019         // sort the unobserved parents in topological order of graph
01020         unobserved_parents.sort();
01021         bool first = true;
01022         RVInstanceArray RHS(0,RHS.size()+unobserved_parents.size());
01023         RVInstanceArray xRHS(RHS);
01024         for (int i=0;i<nup;i++)
01025         {
01026             // note these are still symbolic computations to build Var logp
01027             if (first)
01028                 logp = unobserved_parents[i].V->logP(unobserved_parents[i].v,xRHS,
01029                                                      parameters_to_learn);
01030             else
01031                 logp = logp + 
01032                     unobserved_parents[i].V->logP(unobserved_parents[i].v,xRHS,
01033                                                   parameters_to_learn);
01034             first = false;
01035             // add the visited parents to the RHS, e.g. to compute
01036             //   P(P1=p1,P2=p2,P3=p3) = P(P1=p1)*P(P2=p2|P1=p1)*P(P3=p3|P2=p2,P1=p1)
01037             //
01038             xRHS &= unobserved_parents[i];
01039         }
01040         if (JacobianCorrection)
01041             return logp + *JacobianCorrection;
01042         return logp;
01043     }
01044     // else
01045     return 
01046         PLearn::logP(ConditionalExpression
01047                      (RVInstance(marginalize(this, unobserved_parents.random_variables()), obs),
01048                       RHS),true,parameters_to_learn);
01049 }
01050 
01051 bool FunctionalRandomVariable::invertible(const Var& obs, 
01052                                           RVInstanceArray& unobserved_parents, Var** JacobianCorrection)
01053 {
01054     PLERROR("FunctionalRandomVariable::invertible() should not be called\n"
01055             "Either the sub-class should re-implement logP() or re-define\n"
01056             "invertible() appropriately.");
01057     return false;
01058 }
01059 
01060 // note that stochastic RV's are never non-random,
01061 // but a fonctional RV is non-random if it has
01062 // no parents (i.e., it is a NonRandomVariable) or if all its
01063 // parents are non-random. Note also that the marked field
01064 // sometimes means conditionally non-random during calls
01065 // to functions such as logP.
01066 bool FunctionalRandomVariable::isNonRandom()
01067 { 
01068     bool non_random=true;
01069     for (int i=0;i<parents.size() && non_random;i++)
01070         non_random = parents[i]->isNonRandom();
01071     return non_random;
01072 }
01073 
01074 bool FunctionalRandomVariable::isDiscrete()
01075 {
01076     bool all_discrete = true;
01077     for (int i=0;i<parents.size() && all_discrete;i++)
01078         all_discrete = parents[i]->isDiscrete();
01079     return all_discrete;
01080 }
01081 
01082 
01085 NonRandomVariable::NonRandomVariable(int thelength)
01086     :FunctionalRandomVariable(thelength) {}
01087 
01088 NonRandomVariable::NonRandomVariable(int thelength, int thewidth)
01089     :FunctionalRandomVariable(thelength,thewidth) {}
01090 
01091 NonRandomVariable::NonRandomVariable(const Var& v)
01092     :FunctionalRandomVariable(v) {}
01093 
01096 JointRandomVariable::JointRandomVariable(const RVArray& variables)
01097     :FunctionalRandomVariable(variables,variables.length()) 
01098 { 
01099     if (variables.size()==0)
01100         PLERROR("JointRandomVariables(RVArray) expects an array with >0 elements");
01101 }
01102 
01103 void JointRandomVariable::setValueFromParentsValue() 
01104 {
01105     if (marked) return;
01106     VarArray values(parents.size());
01107     for (int i=0;i<parents.size();i++)
01108         values[i]=parents[i]->value;
01109     value = vconcat(values);
01110 }
01111 
01112 bool JointRandomVariable::invertible(const Var& obs, 
01113                                      RVInstanceArray& unobserved_parents, Var** JacobianCorrection)
01114 {
01115     int p=0;
01116     int j=0;
01117     int nun=unobserved_parents.size();
01118     for (int i=0;i<parents.size();i++)
01119     {
01120         if (j==nun)
01121             PLERROR("JointRandomVariable::invertible ==> logic error");
01122         int l = parents[i]->length();
01123         if (unobserved_parents[j].V==parents[i])
01124             unobserved_parents[j++].v = obs->subVec(p,l);
01125         p+=l;
01126     }
01127     return true;
01128 }
01129 
01130 void JointRandomVariable::EMBprop(const Vec obs, real posterior)
01131 {
01132     int p=0;
01133     for (int i=0;i<parents.size();i++)
01134     {
01135         int l = parents[i]->length();
01136         // watch for redundant computation!
01137         parents[i]->EMBprop(obs.subVec(p,l),posterior);
01138         p+=l;
01139     }
01140 }
01141 
01144 RandomElementOfRandomVariable::RandomElementOfRandomVariable(const RandomVar& v, 
01145                                                              const RandomVar& index)
01146     :FunctionalRandomVariable(v&index,1)
01147 { 
01148     if (index->length()!=1)
01149         PLERROR("RandomElementOfRandomVariables expects an index RandomVar of length 1");
01150 }
01151 
01152 void RandomElementOfRandomVariable::setValueFromParentsValue() 
01153 {
01154     if (marked) return;
01155     value = v()->value[index()->value];
01156 }
01157 
01158 bool RandomElementOfRandomVariable::
01159 invertible(const Var& obs, 
01160            RVInstanceArray& unobserved_parents,
01161            Var** JacobianCorrection)
01162 {
01163     // IT IS POSSIBLE TO COMPUTE A LOGP WITHOUT INTEGRATING
01164     // IF v() IS OBSERVED
01165     if (v()==unobserved_parents[0].V && unobserved_parents.size()==1)
01166         PLERROR("RandomElementOFRandomVariable could compute logP but not implemented");
01167     return false;
01168 }
01169 
01170 void RandomElementOfRandomVariable::EMBprop(const Vec obs, real posterior)
01171 {
01172 }
01173 
01176 RVArrayRandomElementRandomVariable::
01177 RVArrayRandomElementRandomVariable(const RVArray& table, const RandomVar& index)
01178     :FunctionalRandomVariable((RVArray)(table&index),table[0]->length())
01179 { 
01180     int l=table[0]->length();
01181     for (int i=1;i<table.size();i++)
01182         if (table[i]->length()!=l)
01183             PLERROR("RVArrayRandomElementRandomVariables expect all the elements of table\n"
01184                     " to have the same length (%-th has length %d while 0-th has length %d",
01185                     i,table[i]->length(),l);
01186     if (index->length()!=1)
01187         PLERROR("RVArrayRandomElementRandomVariables expect an index RandomVar of length 1");
01188 }
01189 
01190 void RVArrayRandomElementRandomVariable::setValueFromParentsValue() 
01191 {
01192     if (marked) return;
01193     int n = parents.size()-1;
01194     VarArray parents_values(n);
01195     for (int i=0;i<n;i++)
01196         parents_values[i]=parents[i]->value;
01197     value = parents_values[index()->value];
01198 }
01199 
01200 Var RVArrayRandomElementRandomVariable::
01201 logP(const Var& obs, const RVInstanceArray& RHS,
01202      RVInstanceArray* parameters_to_learn)
01203 {
01204     int n = parents.size()-1;
01205     if (index()->isMarked() || index()->isNonRandom())
01206         // special case where index is observed, just pass to selected parent
01207     {
01208         VarArray parents_logp(n);
01209         for (int i=0;i<n;i++)
01210             parents_logp[i]=parents[i]->logP(obs,RHS,parameters_to_learn);
01211         return parents_logp[index()->value];
01212     }
01213     // otherwise, build a Mixture, with log_weights = logP(index()==i)
01214     VarArray log_weights(n);
01215     RVArray components(n);
01216     for (int i=0;i<n;i++)
01217     {
01218         Var indx(1);
01219         indx = (real)i;
01220         log_weights[i] = index()->logP(indx,RHS,parameters_to_learn);
01221         components[i] = parents[i];
01222     }
01223     RandomVar mixt = mixture(components,RandomVar(vconcat(log_weights)));
01224     return mixt->logP(obs,RHS,parameters_to_learn);
01225 }
01226 
01227 void RVArrayRandomElementRandomVariable::EMBprop(const Vec obs, real posterior)
01228 {
01229 }
01230 
01231 
01234 NegRandomVariable::NegRandomVariable(RandomVariable* input)
01235     :FunctionalRandomVariable(input->length()) {}
01236 
01237 void NegRandomVariable::setValueFromParentsValue()
01238 {
01239     if (marked) return;
01240     value = -parents[0]->value;
01241 }
01242 
01243 bool NegRandomVariable::invertible(const Var& obs, 
01244                                    RVInstanceArray& unobserved_parents,
01245                                    Var** JacobianCorrection)
01246 {
01247     unobserved_parents[0].v = -obs;
01248     return true;
01249 }
01250 
01251 void NegRandomVariable::EMBprop(const Vec obs, real posterior)
01252 {
01253     if (!parents[0]->isConstant())
01254         parents[0]->EMBprop(-obs,posterior);
01255 }
01256 
01259 ExpRandomVariable::ExpRandomVariable(RandomVar& input)
01260     :FunctionalRandomVariable(input,input->length()) {}
01261 
01262 void ExpRandomVariable::setValueFromParentsValue()
01263 {
01264     if (marked) return;
01265     value = exp(parents[0]->value);
01266 }
01267 
01268 bool ExpRandomVariable::invertible(const Var& obs, 
01269                                    RVInstanceArray& unobserved_parents,
01270                                    Var** JacobianCorrection)
01271 {
01272     unobserved_parents[0].v = log(obs);
01273     return true;
01274 }
01275 
01276 void ExpRandomVariable::EMBprop(const Vec obs, real posterior)
01277 {
01278     if (!parents[0]->isConstant())
01279         parents[0]->EMBprop(log(obs),posterior);
01280 }
01281 
01284 LogRandomVariable::LogRandomVariable(RandomVar& input)
01285     :FunctionalRandomVariable(input,input->length()) {}
01286 
01287 void LogRandomVariable::setValueFromParentsValue()
01288 {
01289     if (marked) return;
01290     value = log(parents[0]->value);
01291 }
01292 
01293 bool LogRandomVariable::invertible(const Var& obs, 
01294                                    RVInstanceArray& unobserved_parents,
01295                                    Var** JacobianCorrection)
01296 {
01297     unobserved_parents[0].v = exp(obs);
01298     return true;
01299 }
01300 
01301 void LogRandomVariable::EMBprop(const Vec obs, real posterior)
01302 {
01303     if (!parents[0]->isConstant())
01304         parents[0]->EMBprop(exp(obs),posterior);
01305 }
01306 
01309 PlusRandomVariable::PlusRandomVariable(RandomVar input1, RandomVar input2)
01310     : FunctionalRandomVariable(input1 & input2, 
01311                                MAX(input1->length(),input2->length())),
01312       parent_to_learn(parents[0]), other_parent(parents[0]),
01313       numerator(value->length()), difference(value->length())
01314 {
01315     if(input1->length() != input2->length() &&
01316        input1->length() !=1 && input2->length()!=1)
01317         PLERROR("PlusRandomVariable(RandomVariable* in1, RandomVariable* in2) in1 and"
01318                 "in2 must have the same length or one of them must be of length 1");
01319 }
01320   
01321 void PlusRandomVariable::setValueFromParentsValue()
01322 {
01323     if (marked) return;
01324     value = X0()->value + X1()->value;
01325 }
01326 
01327 bool PlusRandomVariable::invertible(const Var& obs, 
01328                                     RVInstanceArray& unobserved_parents,
01329                                     Var** JacobianCorrection)
01330 {
01331     if (unobserved_parents.size()==2)
01332         return false; // can't invert if two parents are unobserved
01333     if (unobserved_parents[0].V == X0())
01334         unobserved_parents[0].v = obs - X1()->value;
01335     else
01336         unobserved_parents[0].v = obs - X0()->value;
01337     return true;
01338     
01339 }
01340 
01341 void PlusRandomVariable::
01342 EMTrainingInitialize(const RVArray& parameters_to_learn)
01343 {
01344     RandomVariable::EMTrainingInitialize(parameters_to_learn);
01345     if (learn_X0() && learn_X1())
01346         PLERROR("PlusRandomVariable: can't learn both X0 and X1");
01347     if (learn_X0() || learn_X1())
01348     {
01349         learn_something=true;
01350         if (learn_X0())
01351         {
01352             parent_to_learn = X0();
01353             other_parent = X1();
01354         }
01355         else
01356         {
01357             parent_to_learn = X1();
01358             other_parent = X0();
01359         }
01360     }
01361 }
01362 
01363 void PlusRandomVariable::EMEpochInitialize()
01364 {
01365     if (EMmark) return;
01366     RandomVariable::EMEpochInitialize();
01367     if (learn_something)
01368     {
01369         numerator.clear();
01370         denom = 0.;
01371     }
01372 }
01373 
01374 void PlusRandomVariable::EMBprop(const Vec obs, real posterior)
01375 {
01376     if (learn_something)
01377     {
01378         // numerator += posterior * (obs - other_parent->value->value);
01379         substract(obs,other_parent->value->value,difference);
01380         multiplyAcc(numerator, difference,posterior);
01381         denom += posterior;
01382         if (!other_parent->isConstant())
01383         {
01384             // propagate to other parent
01385             substract(obs,parent_to_learn->value->value,difference);
01386             other_parent->EMBprop(difference,posterior);
01387         }
01388     }
01389     else
01390     {
01391         if (!X1()->isConstant())
01392         {
01393             substract(obs,X0()->value->value,difference);
01394             X1()->EMBprop(difference,posterior);
01395         }
01396         if (!X0()->isConstant())
01397         {
01398             substract(obs,X1()->value->value,difference);
01399             X0()->EMBprop(difference,posterior);
01400         }
01401     }
01402 }
01403 
01404 void PlusRandomVariable::EMUpdate()
01405 {
01406     if (EMmark) return;
01407     EMmark=true;
01408     if (learn_something && denom>0)
01409         // new value = numerator / denom
01410         multiply(numerator,real(1.0/denom),parent_to_learn->value->value);
01411     if (!learn_X0() && !X0()->isConstant())
01412         X0()->EMUpdate();
01413     if (!learn_X1() && !X1()->isConstant())
01414         X1()->EMUpdate();
01415 }
01416 
01419 MinusRandomVariable::MinusRandomVariable(RandomVar input1, RandomVar input2)
01420     : FunctionalRandomVariable(input1 & input2, 
01421                                MAX(input1->length(),input2->length())),
01422       parent_to_learn(parents[0]), other_parent(parents[0]),
01423       numerator(value->length()), difference(value->length())
01424 {
01425     if(input1->length() != input2->length() &&
01426        input1->length() !=1 && input2->length()!=1)
01427         PLERROR("MinusRandomVariable(RandomVariable* in1, RandomVariable* in2) in1 and"
01428                 "in2 must have the same length or one of them must be of length 1");
01429 }
01430   
01431 void MinusRandomVariable::setValueFromParentsValue()
01432 {
01433     if (marked) return;
01434     value = X0()->value - X1()->value;
01435 }
01436 
01437 bool MinusRandomVariable::invertible(const Var& obs, 
01438                                      RVInstanceArray& unobserved_parents,
01439                                      Var** JacobianCorrection)
01440 {
01441     if (unobserved_parents.size()==2)
01442         return false; // can't invert if two parents are unobserved
01443     if (unobserved_parents[0].V == X0())
01444         unobserved_parents[0].v = obs + X1()->value;
01445     else
01446         unobserved_parents[0].v = X0()->value - obs;
01447     return true;
01448     
01449 }
01450 
01451 void MinusRandomVariable::
01452 EMTrainingInitialize(const RVArray& parameters_to_learn)
01453 {
01454     RandomVariable::EMTrainingInitialize(parameters_to_learn);
01455     if (learn_X0() && learn_X1())
01456         PLERROR("MinusRandomVariable: can't learn both X0 and X1");
01457     if (learn_X0() || learn_X1())
01458     {
01459         learn_something=true;
01460         if (learn_X0())
01461         {
01462             parent_to_learn = X0();
01463             other_parent = X1();
01464         }
01465         else
01466         {
01467             parent_to_learn = X1();
01468             other_parent = X0();
01469         }
01470     }
01471 }
01472 
01473 void MinusRandomVariable::EMEpochInitialize()
01474 {
01475     if (EMmark) return;
01476     RandomVariable::EMEpochInitialize();
01477     if (learn_something)
01478     {
01479         numerator.clear();
01480         denom = 0.0;
01481     }
01482 }
01483 
01484 void MinusRandomVariable::EMBprop(const Vec obs, real posterior)
01485 {
01486     if (learn_something)
01487     {
01488         if (learn_X0())
01489             // numerator += posterior * (obs + other_parent->value->value);
01490             add(obs,other_parent->value->value,difference);
01491         else
01492             // numerator += posterior * (other_parent->value->value - obs);
01493             substract(other_parent->value->value,obs,difference);
01494 
01495         multiplyAcc(numerator, difference,posterior);
01496         denom += posterior;
01497         if (!other_parent->isConstant())
01498         {
01499             // propagate to other parent
01500             if (learn_X0())
01501                 add(obs,parent_to_learn->value->value,difference);
01502             else
01503                 substract(parent_to_learn->value->value,obs,difference);
01504             other_parent->EMBprop(difference,posterior);
01505         }
01506     }
01507     else
01508     {
01509         if (!X1()->isConstant())
01510         {
01511             substract(X0()->value->value,obs,difference);
01512             X1()->EMBprop(difference,posterior);
01513         }
01514         if (!X0()->isConstant())
01515         {
01516             add(obs,X1()->value->value,difference);
01517             X0()->EMBprop(difference,posterior);
01518         }
01519     }
01520 }
01521 
01522 void MinusRandomVariable::EMUpdate()
01523 {
01524     if (EMmark) return;
01525     EMmark=true;
01526     if (learn_something && denom>0)
01527         // new value = numerator / denom
01528         multiply(numerator,real(1.0/denom),parent_to_learn->value->value);
01529     if (!learn_X0() && !X0()->isConstant())
01530         X0()->EMUpdate();
01531     if (!learn_X1() && !X1()->isConstant())
01532         X1()->EMUpdate();
01533 }
01534 
01537 ElementWiseDivisionRandomVariable::
01538 ElementWiseDivisionRandomVariable(RandomVar input1, RandomVar input2)
01539     : FunctionalRandomVariable(input1 & input2, 
01540                                MAX(input1->length(),input2->length()))
01541 {
01542     if(input1->length() != input2->length() &&
01543        input1->length() !=1 && input2->length()!=1)
01544         PLERROR("ElementWiseDivisionRandomVariable(RandomVariable* in1, RandomVariable* in2) in1 and"
01545                 "in2 must have the same length or one of them must be of length 1");
01546 }
01547   
01548 void ElementWiseDivisionRandomVariable::setValueFromParentsValue()
01549 {
01550     if (marked) return;
01551     value = X0()->value / X1()->value;
01552 }
01553 
01554 bool ElementWiseDivisionRandomVariable::invertible(const Var& obs, 
01555                                                    RVInstanceArray& unobserved_parents,
01556                                                    Var** JacobianCorrection)
01557 {
01558     if (unobserved_parents.size()==2)
01559         return false; // can't invert if two parents are unobserved
01560     if (unobserved_parents[0].V == X0())
01561         unobserved_parents[0].v = obs * X1()->value;
01562     else
01563         unobserved_parents[0].v = X0()->value / obs;
01564     return true;
01565 }
01566 
01567 void ElementWiseDivisionRandomVariable::
01568 EMTrainingInitialize(const RVArray& parameters_to_learn)
01569 {
01570 }
01571 
01572 void ElementWiseDivisionRandomVariable::EMEpochInitialize()
01573 {
01574 }
01575 
01576 void ElementWiseDivisionRandomVariable::EMBprop(const Vec obs, real posterior)
01577 {
01578 }
01579 
01580 void ElementWiseDivisionRandomVariable::EMUpdate()
01581 {
01582     PLERROR("ElementWiseDivisionRandomVariable::EMUpdate() not implemented");
01583 }
01584 
01587 ProductRandomVariable::ProductRandomVariable(RandomVar input1, 
01588                                              RandomVar input2)
01589     : FunctionalRandomVariable(input1 & input2, input1->value->matValue.length(),
01590                                input2->value->matValue.width()),
01591       m(input1->value->matValue.length()), n(input1->value->matValue.width()),
01592       l(input2->value->matValue.width()), learn_something(false)
01593 {
01594     if (n != input2->value->matValue.length())
01595         PLERROR("ProductRandomVariable(X0,X1): X0(%d,%d)'s width (%d) must match"
01596                 "X1(%d,%d)'s length (%d)", input1->value->matValue.length(),
01597                 input1->value->matValue.width(), input1->value->matValue.width(),
01598                 input2->value->matValue.length(), input2->value->matValue.width(),
01599                 input2->value->matValue.length());
01600     scalars = (m==1 && n==1 && l==1);
01601 }
01602 
01603 void ProductRandomVariable::setValueFromParentsValue()
01604 {
01605     if (marked) return;
01606     value = Var(new ProductVariable(X0()->value,
01607                                     X1()->value));
01608 }
01609 
01610 bool ProductRandomVariable::invertible(const Var& obs, 
01611                                        RVInstanceArray& unobserved_parents,
01612                                        Var** JacobianCorrection)
01613 {
01614     if (unobserved_parents.size()==2)
01615         return false; // can't invert if two parents are unobserved
01616     if (unobserved_parents[0].V == X0())
01617     {
01618         unobserved_parents[0].v = 
01619             product(obs,rightPseudoInverse(X1()->value));
01620         **JacobianCorrection = log(abs(det(X1()->value)));
01621     }
01622     else
01623     {
01624         unobserved_parents[0].v = 
01625             product( leftPseudoInverse(X0()->value), obs);
01626         **JacobianCorrection = log(abs(det(X0()->value)));
01627     }
01628     return true;
01629     
01630 }
01631 
01632 void ProductRandomVariable::
01633 EMTrainingInitialize(const RVArray& parameters_to_learn)
01634 {
01635     RandomVariable::EMTrainingInitialize(parameters_to_learn);
01636     if (learn_X0() && learn_X1())
01637         PLERROR("ProductRandomVariable: can't learn both X0 and X1");
01638     if (learn_X0() || learn_X1())
01639     {
01640         denom.resize(n,n);
01641         tmp2.resize(n,n);
01642         tmp4.resize(n);
01643         learn_something=true;
01644         if (learn_X0())
01645         {
01646             X0numerator.resize(m,n);
01647             tmp1.resize(m,n);
01648             if (!X1()->isNonRandom())
01649             {
01650                 tmp3.resize(m,l);
01651                 vtmp3 = tmp3.toVec();
01652             }
01653         }
01654         else
01655         {
01656             X1numerator.resize(n,l);
01657             tmp1.resize(n,l);
01658             if (!X0()->isNonRandom())
01659             {
01660                 tmp3.resize(n,m);
01661                 vtmp3 = tmp3.toVec();
01662             }
01663         }
01664     }
01665     else
01666         learn_something=false;
01667 }
01668 
01669 void ProductRandomVariable::EMEpochInitialize()
01670 {
01671     if (EMmark) return;
01672     RandomVariable::EMEpochInitialize();
01673     if (learn_something)
01674     {
01675         denom.clear();
01676         if (learn_X0())
01677             X0numerator.clear();
01678         else
01679             X1numerator.clear();
01680     }
01681 }
01682 
01683 void ProductRandomVariable::EMBprop(const Vec obs, real posterior)
01684 {
01685     if (learn_something)
01686     {
01687         if (learn_X0())
01688         {
01689             if (scalars)
01690             {
01691                 // do the special scalar case separately for efficiency
01692                 real x1 = *(X1()->value->value.data());
01693                 real y = *obs.data();
01694                 *X0numerator.data() += posterior * y * x1;
01695                 *denom.data() += posterior * x1 * x1;
01696                 // propagate EMBprop to X1
01697                 if (!X1()->isNonRandom())
01698                 {
01699                     real x0 = *(X0()->value->value.data());
01700                     if (x0==0.0)
01701                         PLERROR("ProductRandomVariable: can't divide by X0==0");
01702                     *tmp3.data() = y / x0;
01703                     X1()->EMBprop(vtmp3,posterior);
01704                 }
01705             }
01706             else
01707             {
01708                 Mat matObs = obs.toMat(m,l);
01709                 Mat& x1 = X1()->value->matValue;
01710                 // numerator += posterior * obs * x1'
01711                 productTranspose(tmp1, matObs,x1);
01712                 multiplyAcc(X0numerator, tmp1,posterior);
01713                 // denominator += posterior * x1 * x1'
01714                 productTranspose(tmp2, x1,x1);
01715                 multiplyAcc(denom, tmp2,posterior);
01716                 // propagate EMBprop to X1
01717                 if (!X1()->isNonRandom())
01718                 {
01719                     Mat& x0 = X0()->value->matValue;
01720                     // solve x0 * tmp3 = matObs
01721                     solveLinearSystem(x0,matObs,tmp3); 
01722                     X1()->EMBprop(vtmp3,posterior);
01723                 }
01724             }
01725         }
01726         else // learn_X1()
01727         {
01728             if (scalars)
01729             {
01730                 // do the special scalar case separately for efficiency
01731                 real x0 = *(X0()->value->value.data());
01732                 *X1numerator.data() += posterior * *obs.data() * x0;
01733                 *denom.data() += posterior * x0 * x0;
01734                 // propagate EMBprop to X0
01735                 if (!X0()->isNonRandom())
01736                 {
01737                     real x1 = *(X1()->value->value.data());
01738                     if (x1==0.0)
01739                         PLERROR("ProductRandomVariable: can't divide by X1==0");
01740                     *tmp3.data() = obs[0] / x1;
01741                     X0()->EMBprop(vtmp3,posterior);
01742                 }
01743             }
01744             else
01745             {
01746                 Mat matObs = obs.toMat(m,l);
01747                 Mat& x0 = X0()->value->matValue;
01748                 // numerator += posterior * x0' * obs
01749                 transposeProduct(tmp1, x0,matObs);
01750                 multiplyAcc(X1numerator, tmp1,posterior);
01751                 // denominator += posterior * x0' * x0
01752                 transposeProduct(tmp2, x0,x0);
01753                 multiplyAcc(denom, tmp2,posterior);
01754                 // propagate EMBprop to X0
01755                 if (!X0()->isNonRandom())
01756                 {
01757                     Mat& x1 = X1()->value->matValue;
01758                     // solve tmp3 * x1 = matObs
01759                     solveTransposeLinearSystem(x1,matObs,tmp3); 
01760                     X1()->EMBprop(vtmp3,posterior);
01761                 }
01762             }
01763         }
01764     }
01765     else
01766     {
01767         if (scalars)
01768         {
01769             if (!X1()->isNonRandom())
01770             {
01771                 real x0 = *(X0()->value->value.data());
01772                 if (x0==0.0)
01773                     PLERROR("ProductRandomVariable: can't divide by X0==0");
01774                 *tmp3.data() = obs[0] / x0;
01775                 X1()->EMBprop(vtmp3,posterior);
01776             }
01777             if (!X0()->isNonRandom())
01778             {
01779                 real x1 = *(X1()->value->value.data());
01780                 if (x1==0.0)
01781                     PLERROR("ProductRandomVariable: can't divide by X1==0");
01782                 *tmp3.data() = obs[0] / x1;
01783                 X0()->EMBprop(vtmp3,posterior);
01784             }
01785         }
01786         else
01787         {
01788             if (!X1()->isConstant())
01789             {
01790                 Mat matObs = obs.toMat(m,l);
01791                 Mat& x0 = X0()->value->matValue;
01792                 solveLinearSystem(x0,matObs,tmp3); // solve x0 * tmp3 = matObs
01793                 X1()->EMBprop(vtmp3,posterior);
01794             }
01795             if (!X0()->isConstant())
01796             {
01797                 Mat matObs = obs.toMat(m,l);
01798                 Mat& x1 = X1()->value->matValue;
01799                 // solve tmp3 * x1 = matObs
01800                 solveTransposeLinearSystem(x1,matObs,tmp3); 
01801                 X1()->EMBprop(vtmp3,posterior);
01802             }
01803         }
01804     }
01805 }
01806 
01807 void ProductRandomVariable::EMUpdate()
01808 {
01809     if (EMmark) return;
01810     EMmark=true;
01811     if (learn_something)
01812     {
01813         if (learn_X0())
01814         {
01815             if (scalars)
01816             {
01817                 if (denom(0,0)>0)
01818                     X0()->value->value[0] = X0numerator(0,0)/denom(0,0);
01819             }
01820             else
01821                 solveTransposeLinearSystemByCholesky(denom,X0numerator,
01822                                                      X0()->value->matValue,
01823                                                      &tmp2,&tmp4);
01824             if (!X1()->isConstant())
01825                 X1()->EMUpdate();
01826         }
01827         else // learn_X1()
01828         {
01829             if (scalars)
01830             {
01831                 if (denom(0,0)>0)
01832                     X1()->value->value[0] = X1numerator(0,0)/denom(0,0);
01833             }
01834             else
01835                 solveLinearSystemByCholesky(denom,X1numerator,
01836                                             X1()->value->matValue,
01837                                             &tmp2,&tmp4);
01838             if (!X0()->isConstant())
01839                 X0()->EMUpdate();
01840         }
01841     }
01842     else 
01843     {
01844         if (!X0()->isConstant())
01845             X0()->EMUpdate();
01846         if (!X1()->isConstant())
01847             X1()->EMUpdate();
01848     }
01849 }
01850 
01851 
01854 DiagonalNormalRandomVariable::DiagonalNormalRandomVariable
01855 (const RandomVar& mean, const RandomVar& log_var, 
01856  real minimum_standard_deviation)
01857     :StochasticRandomVariable(mean & log_var, mean->length()),
01858      minimum_variance(minimum_standard_deviation*minimum_standard_deviation),
01859      normfactor(mean->length()*Log2Pi), shared_variance(log_var->length()==1),
01860      mu_num(mean->length()), sigma_num(log_var->length())
01861 {
01862 }
01863 
01864 Var DiagonalNormalRandomVariable::logP(const Var& obs, 
01865                                        const RVInstanceArray& RHS, RVInstanceArray* parameters_to_learn)
01866 {
01867     if (mean()->isMarked() && log_variance()->isMarked())
01868     {
01869         if (log_variance()->value->getName()[0]=='#') 
01870             log_variance()->value->setName("log_variance");
01871         if (mean()->value->getName()[0]=='#') 
01872             mean()->value->setName("mean");
01873         Var variance = minimum_variance+exp(log_variance()->value);
01874         variance->setName("variance");
01875         if (shared_variance)
01876             return (-0.5)*(sum(square(obs-mean()->value))/variance+
01877                            (mean()->length())*log(variance) + normfactor);
01878         else
01879             return (-0.5)*(sum(square(obs-mean()->value)/variance)+
01880                            sum(log(variance))+ normfactor);
01881     }
01882     // else
01883     // probably not feasible..., but try in case we know a trick
01884     if (mean()->isMarked())
01885         return 
01886             PLearn::logP(ConditionalExpression(RVInstance(marginalize(this,
01887                                                                       log_variance()),
01888                                                           obs),RHS),true,parameters_to_learn); 
01889     else
01890         return PLearn::logP(ConditionalExpression(RVInstance(marginalize(this,mean()),
01891                                                              obs),RHS),true,parameters_to_learn); 
01892 }
01893 
01894 void DiagonalNormalRandomVariable::setValueFromParentsValue()
01895 {
01896     value = 
01897         Var(new DiagonalNormalSampleVariable(mean()->value,
01898                                              sqrt(minimum_variance+
01899                                                   exp(log_variance()->value))));
01900 }
01901 
01902 void DiagonalNormalRandomVariable::EMEpochInitialize()
01903 {
01904     if (EMmark) return;
01905     RandomVariable::EMEpochInitialize();
01906     if (learn_the_mean())
01907         mu_num.clear();
01908     if (learn_the_variance())
01909         sigma_num.clear();
01910     denom = 0.0;
01911 }
01912 
01913 void DiagonalNormalRandomVariable::EMBprop(const Vec obs, real posterior)
01914 {
01915     if (learn_the_mean())
01916         multiplyAcc(mu_num, obs,posterior);
01917     else if (!mean()->isConstant())
01918     {
01919         if (!shared_variance)
01920             PLERROR("DiagonalNormalRandomVariable: don't know how to EMBprop "
01921                     "into mean if variance is not shared");
01922         mean()->EMBprop(obs,posterior/
01923                         (minimum_variance+exp(log_variance()->value->value[0])));
01924     }
01925     if (learn_the_variance())
01926     {
01927         if (learn_the_mean())
01928         {
01929             // sigma_num[i] += obs[i]*obs[i]*posterior
01930             if (shared_variance)
01931                 sigma_num[0] += posterior*pownorm(obs)/mean()->length();
01932             else
01933                 squareMultiplyAcc(sigma_num, obs,posterior);
01934         }
01935         else
01936         {
01937             // sigma_num[i] += (obs[i]-mean[i])^2*posterior
01938             if (shared_variance)
01939                 sigma_num[0] += posterior*powdistance(obs,mean()->value->value)
01940                     /mean()->length();
01941             else
01942                 diffSquareMultiplyAcc(sigma_num, obs,
01943                                       mean()->value->value,
01944                                       posterior);
01945         }
01946     }
01947     else if (!log_variance()->isConstant())
01948     {
01949         // use sigma_num as a temporary for log_var's observation
01950         if (shared_variance)
01951             log_variance()->EMBprop(Vec(1,powdistance(obs,mean()->value->value)
01952                                         /mean()->length()),
01953                                     posterior);
01954         else
01955         {
01956             substract(obs,mean()->value->value,sigma_num);
01957             apply(sigma_num,sigma_num,square_f);
01958             log_variance()->EMBprop(sigma_num,posterior);
01959         }
01960     }
01961     if (learn_the_mean() || learn_the_variance()) denom += posterior;
01962 }
01963 
01964 void DiagonalNormalRandomVariable::EMUpdate()
01965 {
01966     if (EMmark) return;
01967     EMmark=true;
01968     // maybe we should issue a warning if
01969     // (learn_the_mean || learn_the_variance) && denom==0
01970     // (it means that all posteriors reaching EMBprop were 0)
01971     //
01972     if (denom>0 && (learn_the_mean() || learn_the_variance()))
01973     {
01974         Vec lv = log_variance()->value->value;
01975         Vec mv = mean()->value->value;
01976         if (learn_the_mean())
01977             multiply(mu_num,real(1.0/denom),mv);
01978         if (learn_the_variance())
01979         {
01980             if (learn_the_mean())
01981             {
01982                 // variance = sigma_num/denom - squared(mean)
01983                 sigma_num /= denom;
01984                 multiply(mv,mv,mu_num); // use mu_num as a temporary vec
01985                 if (shared_variance)
01986                     lv[0] = sigma_num[0] - PLearn::mean(mu_num);
01987                 else
01988                     substract(sigma_num,mu_num,lv);
01989                 // now lv really holds variance
01990 
01991                 // log_variance = log(max(0,variance-minimum_variance))
01992                 substract(lv,minimum_variance,lv);
01993                 max(lv,real(0.),lv);
01994                 apply(lv,lv,safeflog);
01995             }
01996             else
01997             {
01998                 multiply(sigma_num,1/denom,lv);
01999                 // now log_variance really holds variance
02000 
02001                 // log_variance = log(max(0,variance-minimum_variance))
02002                 substract(lv,minimum_variance,lv);
02003                 max(lv,real(0.),lv);
02004                 apply(lv,lv,safeflog);
02005             }
02006         }
02007     }
02008     if (!learn_the_mean() && !mean()->isConstant())
02009         mean()->EMUpdate();
02010     if (!learn_the_variance() && !log_variance()->isConstant())
02011         log_variance()->EMUpdate();
02012 }
02013 
02014 
02017 MixtureRandomVariable::MixtureRandomVariable
02018 (const RVArray& the_components,const RandomVar& logweights)
02019     :StochasticRandomVariable(logweights,the_components[0]->length()),
02020      components(the_components), posteriors(logweights->length()),
02021      sum_posteriors(logweights->length()), 
02022      componentsLogP(logweights->length()),
02023      lw(logweights->length())
02024 {
02025 }
02026 
02027 Var MixtureRandomVariable::logP(const Var& obs, const RVInstanceArray& RHS,
02028                                 RVInstanceArray* parameters_to_learn)
02029 {
02030     if (parameters_to_learn!=0) return ElogP(obs,*parameters_to_learn,RHS);
02031     if (log_weights()->isMarked())
02032     {
02033         int n=posteriors.length();
02034         if (log_weights()->value->getName()[0]=='#') 
02035             log_weights()->value->setName("log_weights");
02036         Var weights = softmax(log_weights()->value);
02037         weights->setName("weights");
02038         lw = log(weights);
02039         for (int i=0;i<n;i++)
02040             componentsLogP[i] = components[i]->logP(obs,RHS) + lw->subVec(i,1);
02041         logp = logadd(vconcat(componentsLogP));
02042         return logp;
02043     }
02044     // else
02045     // probably not feasible..., but try in case we know a trick
02046     return PLearn::logP(ConditionalExpression
02047                         (RVInstance(marginalize(this,log_weights()),obs),RHS),true,0);
02048 }
02049 
02050 // compute symbolically 
02051 //   E[ logP(obs,i|new_params) | old_params,obs ]
02052 // where the expectation is over the mixture component index i,
02053 // with "posterior" probabilities P(i|obs,old_params).
02054 // The parameters_to_learn[i].v hold the "new parameters" (to be optimized)
02055 // while the parameters_to_learn[i].V->value hold the "current value"
02056 // of the parameters for the EM algorithm.
02057 Var MixtureRandomVariable::ElogP(const Var& obs, 
02058                                  RVInstanceArray& parameters_to_learn,
02059                                  const RVInstanceArray& RHS)
02060 {
02061     if (log_weights()->isMarked())
02062     {
02063         int n=posteriors.length();
02064 
02065         // (1) using the "current" value of the parameters
02066         Var weights = softmax(log_weights()->value);
02067         lw = log(weights);
02068         for (int i=0;i<n;i++)
02069             // componentsLogP[i] = log(P(obs|i)*P(i)) = log(P(obs,i))
02070             componentsLogP[i] = components[i]->logP(obs,RHS) + lw->subVec(i,1);
02071         // logp = log(P(obs)) = log(sum_i P(obs,i))
02072         logp = logadd(vconcat(componentsLogP));
02073         // now compute log-posteriors by normalization
02074         for (int i=0;i<n;i++)
02075             // componentsLogP[i] = log(P(i|obs))=log(P(obs,i)/P(obs))
02076             componentsLogP[i] = componentsLogP[i] - logp;
02077 
02078         // (2) now put the "new" value of the parameters (swap with v fields)
02079         parameters_to_learn.swap_v_and_Vvalue();
02080         // unmark parents and re-compute value's in terms of ancestors' values
02081         unmarkAncestors();
02082         markRHSandSetKnownValues(RHS);
02083      
02084         // (3) and compute the  logP of each component weighted by its posterior
02085         weights = softmax(log_weights()->value);
02086         for (int i=0;i<n;i++)
02087             componentsLogP[i] = exp(components[i]->logP(obs,RHS,&parameters_to_learn) 
02088                                     + componentsLogP[i]);
02089         logp = sum(vconcat(componentsLogP));
02090 
02091         // (4) now put back the "current" value of parameters in their value field
02092         parameters_to_learn.swap_v_and_Vvalue();
02093         // unmark parents and re-compute value's in terms of ancestors' values
02094         unmarkAncestors();
02095         markRHSandSetKnownValues(RHS);
02096 
02097         return logp;
02098     }
02099     // else
02100     // probably not feasible..., but try in case we know a trick
02101     return PLearn::logP(ConditionalExpression
02102                         (RVInstance(marginalize(this,log_weights()),obs),
02103                          RHS),true,&parameters_to_learn);
02104 }
02105 
02106 void MixtureRandomVariable::setValueFromParentsValue()
02107 {
02108     Var index = new MultinomialSampleVariable(softmax(log_weights()->value));
02109     value = components.values()[index]; 
02110 }
02111 
02112 void MixtureRandomVariable::
02113 EMTrainingInitialize(const RVArray& parameters_to_learn)
02114 {
02115     if (EMmark) return;
02116     EMmark=true;
02117     learn_the_weights() = parameters_to_learn.contains(log_weights())
02118         && log_weights()->isConstant();
02119     for (int i=0;i<components.size();i++)
02120         components[i]->EMTrainingInitialize(parameters_to_learn);
02121 }
02122 
02123 void MixtureRandomVariable::EMEpochInitialize()
02124 {
02125     if (EMmark) return;
02126     RandomVariable::EMEpochInitialize();
02127     if (learn_the_weights())
02128         sum_posteriors.clear();
02129     for (int i=0;i<components.size();i++)
02130         components[i]->EMEpochInitialize();
02131 }
02132 
02133 void MixtureRandomVariable::EMBprop(const Vec obs, real posterior)
02134 {
02135     // ASSUME THAT AN FPROP HAS BEEN PERFORMED
02136     // so that weights and componentsLogP hold appropriate value
02137     //
02138     // compute posterior vector for this observation
02139     // posteriors = posterior*(components[i]->logP(obs)*weights/normalize)
02140     real log_p = logp->value[0];
02141     real *p = posteriors.data();
02142     int n = lw->value.length();
02143     for (int i=0;i<n;i++)
02144         p[i] = componentsLogP[i]->value[0] - log_p;
02145 #ifdef _MSC_VER
02146     apply(posteriors,posteriors,(tRealFunc)exp);
02147 #else
02148     apply(posteriors,posteriors,safeexp);
02149 #endif
02150     if (fabs(sum(posteriors)-1)>1e-5)
02151     {
02152         cout << "sum(posteriors) = " << sum(posteriors) << "!" << endl;
02153     }
02154     posteriors *= posterior;
02155 
02156     if (learn_the_weights())
02157         sum_posteriors+=posteriors; 
02158 
02159     // propagate to components
02160     for (int i=0;i<n;i++)
02161         components[i]->EMBprop(obs,posteriors[i]);
02162 }
02163 
02164 void MixtureRandomVariable::EMUpdate()
02165 {
02166     if (EMmark) return;
02167     EMmark=true;
02168     // update weights
02169     if (learn_the_weights())
02170     {
02171         real denom = sum(sum_posteriors);
02172         if (denom>0)
02173         {
02174             multiply(sum_posteriors,real(1.0/denom),posteriors);
02175             apply(posteriors,log_weights()->value->value,safeflog);
02176         }
02177     }
02178     // propagate to components
02179     for (int i=0;i<components.size();i++)
02180         components[i]->EMUpdate();
02181 }
02182 
02183 bool MixtureRandomVariable::canStopEM()
02184 {
02185     // propagate to components
02186     bool can=log_weights()->canStopEM();
02187     for (int i=0;i<components.size() && !can;i++)
02188         can = components[i]->canStopEM();
02189     return can;
02190 }
02191 
02192 bool MixtureRandomVariable::isDiscrete()
02193 {
02194     return components[0]->isDiscrete();
02195 }
02196 
02197 void MixtureRandomVariable::setKnownValues()
02198 {
02199     if (!pmark && !marked)
02200     {
02201         pmark = true;
02202         log_weights()->setKnownValues();
02203         for (int i=0;i<components.size();i++)
02204             components[i]->setKnownValues();
02205         setValueFromParentsValue();
02206     }
02207 }
02208  
02209 void MixtureRandomVariable::unmarkAncestors() 
02210 { 
02211     if (pmark)
02212     {
02213         marked=false; 
02214         pmark=false;
02215         log_weights()->unmarkAncestors();
02216         for (int i=0;i<components.size();i++)
02217             components[i]->unmarkAncestors();
02218     }
02219 }
02220 
02221 void MixtureRandomVariable::clearEMmarks() 
02222 { 
02223     if (EMmark)
02224     {
02225         EMmark=false; 
02226         log_weights()->clearEMmarks();
02227         for (int i=0;i<components.size();i++)
02228             components[i]->clearEMmarks();
02229     }
02230 }
02231 
02234 MultinomialRandomVariable::
02235 MultinomialRandomVariable(const RandomVar& log_probabilities)
02236     :StochasticRandomVariable(log_probabilities,1),
02237      sum_posteriors(log_probabilities->length())
02238 {
02239 }
02240 
02241 Var MultinomialRandomVariable::logP(const Var& obs, const RVInstanceArray& RHS,
02242                                     RVInstanceArray* parameters_to_learn)
02243 {
02244     if (log_probabilities()->isMarked())
02245         return log(softmax(log_probabilities()->value))[obs];
02246     // else
02247     // probably not feasible..., but try in case we know a trick
02248     return PLearn::logP(ConditionalExpression
02249                         (RVInstance(marginalize(this,log_probabilities()),obs),
02250                          RHS),true,parameters_to_learn); 
02251 }
02252 
02253 void MultinomialRandomVariable::EMEpochInitialize()
02254 {
02255     if (EMmark) return;
02256     RandomVariable::EMEpochInitialize();
02257     if (learn_the_probabilities())
02258         sum_posteriors.clear();
02259 }
02260 
02261 void MultinomialRandomVariable::EMBprop(const Vec obs, real posterior)
02262 {
02263     if (learn_the_probabilities())
02264     {
02265         real *p = sum_posteriors.data();
02266         p[(int)obs[0]] += posterior;
02267     }
02268 }
02269 
02270 void MultinomialRandomVariable::EMUpdate()
02271 {
02272     if (EMmark) return;
02273     EMmark=true;
02274     if (learn_the_probabilities())
02275     {
02276         real denom = sum(sum_posteriors);
02277         if (denom>0)
02278             // update probabilities
02279         {
02280             multiply(sum_posteriors,real(1.0/denom),sum_posteriors);
02281             apply(sum_posteriors,log_probabilities()->value->value,safeflog);
02282         }
02283         // maybe should WARN the user if denom==0 here
02284     }
02285 }
02286 
02287 bool MultinomialRandomVariable::isDiscrete()
02288 {
02289     return true;
02290 }
02291 
02292 void MultinomialRandomVariable::setValueFromParentsValue()
02293 {
02294     value = 
02295         Var(new MultinomialSampleVariable(softmax(log_probabilities()->value)));
02296 }
02297 
02300 SubVecRandomVariable::SubVecRandomVariable
02301 (const RandomVar& v, int the_start, int the_len)
02302     :FunctionalRandomVariable(v, the_len), start(the_start)
02303 {
02304     if (v->length() < the_len)
02305         PLERROR("new SubVecRandomVariable: input should have length at least %d, has %d",
02306                 the_len, v->length());
02307 }
02308 
02309 void SubVecRandomVariable::setValueFromParentsValue()
02310 {
02311     if (marked) return;
02312     value = parents[0]->value->subVec(start,value->length());
02313 }
02314 
02315 bool SubVecRandomVariable::
02316 invertible(const Var& obs, RVInstanceArray& unobserved_parents,
02317            Var** JacobianCorrection)
02318 {
02319     if (length()!=parents[0]->length())
02320         return false;
02321     // the only invertible case is when this RandomVar is a copy of its parent
02322     unobserved_parents[0].v = obs;
02323     return true;
02324 }
02325 
02326 void SubVecRandomVariable::EMBprop(const Vec obs, real posterior)
02327 {
02328     if (length()!=parents[0]->length())
02329         PLERROR("SubVecRandomVariable: can't EMBprop unless length==parent.length()");
02330     // the only "invertible" case is when this RandomVar is a copy of its parent
02331     parents[0]->EMBprop(obs,posterior);
02332 }
02333 
02336 ExtendedRandomVariable::ExtendedRandomVariable
02337 (const RandomVar& v, real fillvalue, int nextend)
02338     :FunctionalRandomVariable(v, v->length()+nextend), 
02339      n_extend(nextend), fill_value(fillvalue)
02340 {
02341 }
02342 
02343 void ExtendedRandomVariable::setValueFromParentsValue()
02344 {
02345     if (marked) return;
02346     value = extend(parents[0]->value,fill_value,n_extend);
02347 }
02348 
02349 bool ExtendedRandomVariable::
02350 invertible(const Var& obs, RVInstanceArray& unobserved_parents,
02351            Var** JacobianCorrection)
02352 {
02353     int n = n_extend*parents[0]->width();
02354     unobserved_parents[0].v = obs->subVec(parents[0]->value->length(),n);
02355     return true;
02356 }
02357 
02358 void ExtendedRandomVariable::EMBprop(const Vec obs, real posterior)
02359 {
02360     int n = n_extend*parents[0]->width();
02361     parents[0]->EMBprop(obs.subVec(parents[0]->value->length(),n),posterior);
02362 }
02363 
02366 // *** A REVOIR *** Pascal
02367 
02368 ConcatColumnsRandomVariable::ConcatColumnsRandomVariable(const RVArray& a)
02369     :FunctionalRandomVariable(a, a.length())
02370 {
02371     setValueFromParentsValue(); // just to check compatibility
02372     // for (int i=0;i<a.size();i++)
02373     // int n_rows = a[0]->value->matValue.length();
02374     // Je commente ca parce que la methode n'existe plus, mais ca avait surement son utilite... 
02375     // seeAsMatrix(n_rows,length()/n_rows);
02376 }
02377 
02378 void ConcatColumnsRandomVariable::setValueFromParentsValue()
02379 {
02380     if (marked) return;
02381     value = hconcat(parents.values());
02382 }
02383 
02384 bool ConcatColumnsRandomVariable::
02385 invertible(const Var& obs, RVInstanceArray& unobserved_parents,
02386            Var** JacobianCorrection)
02387 {
02388     PLERROR("ConcatColumnsRandomVariable::invertible not yet implemented");
02389     return true;
02390 }
02391 
02392 void ConcatColumnsRandomVariable::EMBprop(const Vec obs, real posterior)
02393 {
02394     PLERROR("ConcatColumnsRandomVariable::EMBprop not yet implemented");
02395 }
02396 
02397 /*** RandomVarVMatrix ***/
02398 
02399 // DEPRECATED: this class should be rewritten entirely or erased.
02400 // It probably won't work in its current state.
02401 
02402 RandomVarVMatrix::
02403 RandomVarVMatrix(ConditionalExpression conditional_expression)
02404     :VMatrix(-1,-1), instance(Sample(conditional_expression)) // extract the "sampling algorithm"
02405 {
02406     // make sure all non-random dependencies are computed
02407     instance->fprop_from_all_sources();
02408     // extract the path of dependencies from all stochastically sampled Vars to instance
02409     instance->random_sources().setMark(); // mark the random sources
02410     instance->markPath(); // mark successors of the random sources
02411     instance->buildPath(prop_path); // extract path from the random sources to instance
02412     // and clear marks
02413 }
02414 
02415 //  template <>
02416 //  void deepCopyField(RandomVar& field, CopiesMap& copies)
02417 //  {
02418 //    if (field)
02419 //      field = static_cast<RandomVariable*>(field->deepCopy(copies));
02420 //  }
02421 
02422 } // end of namespace PLearn
02423 
02424 
02425 /*
02426   Local Variables:
02427   mode:c++
02428   c-basic-offset:4
02429   c-file-style:"stroustrup"
02430   c-file-offsets:((innamespace . 0)(inline-open . 0))
02431   indent-tabs-mode:nil
02432   fill-column:79
02433   End:
02434 */
02435 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=79 :
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines