openModeller  Version 1.4.0
Rule.cpp
Go to the documentation of this file.
00001 /* **************************************
00002  *  GARP Modeling Package
00003  *
00004  * **************************************
00005  *
00006  * Copyright (c), The Center for Research, University of Kansas, 2385 Irving Hill Road, Lawrence, KS 66044-4755, USA.
00007  * Copyright (C), David R.B. Stockwell of Symbiotik Pty. Ltd.
00008  *
00009  * This program is free software; you can redistribute it and/or modify
00010  * it under the terms of the license that is distributed with the software.
00011  *
00012  * This program is distributed in the hope that it will be useful,
00013  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00014  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015  * license.txt file included with this software for more details.
00016  */
00017 
00018 // Rule.cpp: implementation of the Rule class.
00019 //
00021 
00022 #include "Rule.h"
00023 #include "EnvCellSet.h"
00024 #include "Utilities.h"
00025 
00026 #include <openmodeller/om.hh>
00027 
00028 #include <math.h> 
00029 
00030 // ==========================================================================
00031 //  Rule implelentation
00032 // ==========================================================================
00033 Rule::Rule() 
00034 {
00035   Gene = NULL;
00036   intGenes = 0;
00037   
00038   for (int i = 0; i < 10; i++)
00039     dblPerformance[i] = 0.0;
00040 
00041   bGeneIsActive = NULL;
00042   iGeneIndex = NULL;
00043   iActiveGenes = 0;
00044 
00045   blnNeedsEvaluation = true;
00046   intGens = 0;
00047   intTrials = 0;
00048   intScreener = 0;
00049   intScreen = 0;
00050   intLength = 0;
00051   intNumber = 0;
00052   intConclusion = 0;
00053   chrOrigin = (char) 0;
00054   chrPad = (char) 0;
00055   
00056   lId = 0;
00057   iOrigGen = 0;
00058 
00059   _pXYs = 0.0;
00060   _no = 0;
00061   _dA = 0.0;
00062   _dSig = 0.0;
00063 }
00064 
00065 // ==========================================================================
00066 Rule::~Rule() 
00067 {
00068   if (Gene)
00069     delete[] Gene;
00070 }
00071 
00072 // ==========================================================================
00073 Rule * Rule::clone()
00074 {
00075   int i;
00076 
00077   Rule * newRule = objFactory();
00078 
00079   for (i = 0; i < 10; i++) 
00080     newRule->dblPerformance[i] = dblPerformance[i];
00081   
00082   newRule->intGenes = intGenes;
00083   newRule->intGens = intGens;
00084   newRule->blnNeedsEvaluation = blnNeedsEvaluation;
00085   newRule->intTrials = intTrials;
00086   newRule->intScreener = intScreener;
00087   newRule->intScreen = intScreen;
00088   newRule->intLength = intLength;
00089   newRule->chrOrigin = chrOrigin;
00090   newRule->chrPad = chrPad;
00091 
00092   if (newRule->Gene) delete[] newRule->Gene;
00093   newRule->Gene = new BYTE[newRule->intLength];
00094 
00095   for (i = 0;i < intLength; i++)
00096     { newRule->Gene[i] = Gene[i]; }
00097 
00098   // parameters related to the algorithm
00099   newRule->bGeneIsActive = bGeneIsActive;
00100   newRule->iGeneIndex = iGeneIndex;
00101   newRule->iActiveGenes = iActiveGenes;
00102 
00103   newRule->lId = lId;
00104   newRule->iOrigGen = iOrigGen;
00105 
00106   newRule->_pXYs = _pXYs;
00107   newRule->_no = _no;
00108   newRule->_dA = _dA;
00109   newRule->_dSig = _dSig;
00110 
00111   return newRule;
00112 }
00113 
00114 // ==========================================================================
00115 void Rule::copy(Rule * fromRule)
00116 {
00117   int i;
00118 
00119   if (type() != fromRule->type())
00120     throw GarpException(1, "Cannot copy rules with different types");
00121 
00122   for (i = 0; i < 10; i++) 
00123     dblPerformance[i] = fromRule->dblPerformance[i];
00124   
00125   intGenes      = fromRule->intGenes;
00126   intGens       = fromRule->intGens;
00127   blnNeedsEvaluation  = fromRule->blnNeedsEvaluation;
00128   intTrials     = fromRule->intTrials;
00129   intScreener     = fromRule->intScreener;
00130   intScreen     = fromRule->intScreen;
00131   intLength     = fromRule->intLength;
00132   chrOrigin     = fromRule->chrOrigin;
00133   chrPad        = fromRule->chrPad;
00134 
00135   if (Gene) delete[] Gene;
00136   Gene = new BYTE[intLength];
00137 
00138   for (i = 0; i < intLength; i++)
00139     { Gene[i] = fromRule->Gene[i]; }
00140 
00141   // parameters related to the algorithm
00142   bGeneIsActive = fromRule->bGeneIsActive;
00143   iGeneIndex    = fromRule->iGeneIndex;
00144   iActiveGenes  = fromRule->iActiveGenes;
00145 
00146   lId = fromRule->lId;
00147   iOrigGen = fromRule->iOrigGen;
00148 
00149   _pXYs = fromRule->_pXYs;
00150   _no = fromRule->_no;
00151   _dA = fromRule->_dA;
00152   _dSig = fromRule->_dSig;
00153 }
00154 // ==========================================================================
00155 void
00156 Rule::RestoreRule( double *perf, unsigned char *genes, int arry_len, int *gene_index )
00157 {
00158   for( int i=0; i<10; i++ ) {
00159 
00160     dblPerformance[i] = perf[i];
00161   }
00162 
00163   iGeneIndex = gene_index;
00164 
00165   if (Gene)
00166     delete [] Gene;
00167 
00168   Gene = new BYTE[arry_len];
00169 
00170   intLength = arry_len;
00171   iActiveGenes = arry_len/2;
00172   for( int i=0; i<intLength; i++ ) {
00173     Gene[i] = (BYTE) genes[i];
00174   }
00175 
00176 
00177 }
00178 
00179 // ==========================================================================
00180 void Rule::initialize(EnvCellSet * objEnvCellSet, const RuleSet * objRuleSet, 
00181             bool * bGeneIsActivePtr, int * iGeneIndexPtr, int iActGenes)
00182 {
00183   int i, p;
00184 
00185   // parameters related to the algorithm
00186   bGeneIsActive = bGeneIsActivePtr;
00187   iGeneIndex = iGeneIndexPtr;
00188   iActiveGenes = iActGenes;
00189 
00190   if (!iGeneIndex)
00191     i = 0;
00192 
00193   intTrials = 0;
00194   blnNeedsEvaluation = true;
00195   intGens = 0;
00196   chrOrigin = type();
00197   intScreen = 1;
00198   chrPad = ' ';
00199 
00200   // number of genes is equal to the number of values in a cell
00201   // dont take the mask into account
00202   intGenes = objEnvCellSet->get(0)->size() - 1;
00203   intLength = intGenes * 2;
00204 
00205   //printf("intGenes=%3d intLength=%3d ActGenes=%3d\n", intGenes, intLength, iActGenes);
00206 
00207   if (Gene) delete[] Gene;
00208   Gene = new BYTE[intLength];
00209 
00210   for (i = 1; i < intGenes; i++)
00211   {
00212     Gene[i * 2]     = 0;
00213     Gene[i * 2 + 1] = 255;
00214   }
00215 
00216   p = GarpUtil::randint(0, objEnvCellSet->count() - 1);
00217   Gene[0] = objEnvCellSet->get(p)->values[0];
00218     Gene[1] = 0;
00219 
00220   intConclusion = Gene[0];
00221 }
00222 
00223 // ==========================================================================
00224 char * Rule::toString()
00225 {
00226   int i;
00227   char typ = (char) 0;
00228   char * strText = new char[1024];
00229 
00230   // rule type
00231   typ = this->type();
00232   strcpy(strText, "");
00233   sprintf(strText, "%c ", typ);
00234 
00235   // genes
00236   for (i = 0; i < intLength; i++)
00237     sprintf(strText, "%s %3d ", strText, Gene[i]);
00238 
00239   // performance
00240   for (i = 0; i < 10; i++)
00241     sprintf(strText, "%s %-5.3f ", strText, dblPerformance[i]);
00242 
00243   // other values
00244 
00245   // check string sizes
00246   if (strlen(strText) > 1024)
00247     throw GarpException(82, "String size exceeded in Rule::toString()");
00248 
00249   // return result
00250   return strText;
00251 }
00252 
00253 // ==========================================================================
00254 char * Rule::toXML()
00255 {
00256   /* XML Sample
00257 
00258     <Rule Type="r">
00259       <Genes> 1 0 64 128 12 15 0 255 </Genes>
00260       <Performance> 0.0000 1.2500 ... </Performance>
00261     </Rule>
00262   */
00263 
00264   int i;
00265   char strAux[64];
00266   char * strXML;
00267 
00268   strXML = new char[1024];
00269 
00270   // rule type
00271   sprintf(strXML, "<Rule Type=\"%c\" Id=\"%d\" OrigGen=\"%d\">\n", this->type(), lId, iOrigGen);
00272 
00273   // genes
00274   strcat(strXML, "  <Genes>");
00275   for (i = 0; i < intLength; i++)
00276   {
00277     sprintf(strAux, " %3d", Gene[i]);
00278     strcat(strXML, strAux); 
00279   }
00280   strcat(strXML, "</Genes>\n");
00281 
00282   // performance
00283   strcat(strXML, "  <Performance>");
00284   for (i = 0; i < 10; i++)
00285   {
00286     sprintf(strAux, " %-5.3f", dblPerformance[i]);
00287     strcat(strXML, strAux); 
00288   }
00289 
00290   strcat(strXML, "</Performance>");
00291 
00292   strcat(strXML, "</Rule>");
00293 
00294   // other values
00295 
00296   // check string sizes
00297   if (strlen(strXML) > 1024)
00298     throw GarpException(82, "String size exceeded in Rule::toXML()");
00299 
00300   // return result
00301   return strXML;
00302 }
00303 
00304 // ==========================================================================
00305 double Rule::getCertainty(EnvCell * cell)
00306 { return (Gene[0] == cell->values[0]); }
00307 
00308 // ==========================================================================
00309 double Rule::getError(BYTE pred, EnvCell * cell)
00310 { return (double) abs( pred - cell->values[0] ); }
00311  
00312 // ==========================================================================
00313 void Rule::mutate(int intTemperature)
00314 {
00315   BYTE aux1, aux2, aux_g1, aux_g2;
00316   int rnd1, rnd2;
00317   int j, k;
00318 
00319     j = GarpUtil::randint(1, iActiveGenes);
00320   k = 2 * iGeneIndex[j];
00321 
00322   rnd1 = GarpUtil::randint(-intTemperature, intTemperature);
00323   rnd2 = GarpUtil::randint(-intTemperature, intTemperature);
00324 
00325   if ((k + 1 > intLength) || (k < 2))
00326     throw GarpException(100, "Array read out of bounds (Rule::mutate)");
00327 
00328   aux1 = Gene[k];
00329   aux2 = Gene[k + 1];
00330   
00331   Gene[k]     = (BYTE) (aux1 - rnd1);
00332   Gene[k + 1] = (BYTE) (aux2 + rnd2);
00333 
00334   aux_g1 = MIN(Gene[k], Gene[k + 1]);
00335   aux_g2 = MAX(Gene[k], Gene[k + 1]);
00336 
00337   Gene[k]     = aux_g1;
00338   Gene[k + 1] = aux_g2;
00339 
00340   if (type() == 'a')
00341     Gene[k + 1] = Gene[k];
00342 
00343   blnNeedsEvaluation = true;
00344   intGens = 0;
00345   chrOrigin = 'm';
00346 }
00347 
00348 // ==========================================================================
00349 bool Rule::similar(Rule * objOtherRule)
00350 {
00351   bool found;
00352   int j, k;
00353   int ng_2k, ng_2k_1, bg_2k, bg_2k_1;
00354 
00355   if (type() == objOtherRule->type())
00356   {
00357     // check rule value (presence/absence)
00358     if (Gene[0] != objOtherRule->Gene[0]) 
00359       return 0;
00360 
00361     for (j = 1, found = true; (j < iActiveGenes) && (found); j += 1)
00362     {
00363       k = iGeneIndex[j];
00364 
00365       if ((k * 2 + 1 > intLength) || (k * 2 < 2))
00366         throw GarpException(100, "Array read out of bounds (Rule::similar)");
00367 
00368       ng_2k_1 = Gene[k * 2 + 1];
00369       ng_2k   = Gene[k * 2];
00370 
00371       bg_2k_1 = objOtherRule->Gene[k * 2 + 1];
00372       bg_2k   = objOtherRule->Gene[k * 2];
00373 
00374       found = !( ((ng_2k_1 - ng_2k) == 255 && 
00375             (bg_2k_1 - bg_2k) != 255 ) || 
00376              ((ng_2k_1 - ng_2k) != 255 && 
00377               (bg_2k_1 - bg_2k) == 255 ) );
00378     }
00379 
00380     return found;
00381   }
00382 
00383   return false;
00384 }
00385 
00386 // ==========================================================================
00387 double Rule::testWithData(EnvCellSet * objTrainSet)
00388 {
00389   const int MAX_UTILS = 10;
00390   double Utility[MAX_UTILS];
00391 
00392   EnvCell * cell;
00393   int n, i, rnd, no=0;
00394   double prediction, certainty, strength;
00395   double pXs=0,pYs=0,pXYs=0,pYcXs=0,pXSs=0,pYcs=0;
00396   
00397   for (i = 1; i < MAX_UTILS; i++) 
00398     Utility[i] = 0.0;
00399 
00400   Utility[0] = 1.0;
00401 
00402   n = objTrainSet->size();
00403   for(i=0; i < n; i++)
00404     { 
00405     // Get an in range data point
00406     rnd = GarpUtil::randint(0, n - 1);
00407     cell = objTrainSet->get(rnd);
00408   
00409     strength = getStrength(cell);
00410     certainty = getCertainty(cell);
00411     prediction = getError(128, cell);
00412 
00413     pXs  += strength;
00414     pYs  += certainty;
00415     pYcs += prediction;
00416 
00417     if (strength > 0)
00418     {
00419       pXSs  += strength;
00420       pYcXs += getError((BYTE)prediction, cell);
00421       pXYs  += (MIN(certainty, strength)) / strength;
00422       no++;
00423     }
00424     }
00425 
00426   // Priors
00427   Utility[1] = pXs/n;   // proportion 
00428   Utility[2] = pYs/n;   // Prior probability
00429   Utility[3] = pYcs/n;
00430   
00431   // Posteriors
00432   if (no > 0) 
00433     Utility[4] = pXSs/no;
00434   else 
00435     Utility[4] = 0;
00436   
00437   if (no > 0)         // Posterior probability
00438     Utility[5] = pXYs/no; 
00439   else 
00440     Utility[5] = 0;
00441   
00442   if (no>0) 
00443     Utility[6] = (pYcXs/no);
00444   else 
00445     Utility[6] = 0;
00446   
00447   Utility[7] = ((double)no)/n;
00448 
00449   // Crisp Significance
00450   if (no >= MIN_SIG_NO && Utility[2] > 0 && Utility[2] < 1.0) 
00451     Utility[8]= (pXYs-Utility[2] * no) / sqrt(no * Utility[2] * (1 - Utility[2]));
00452   else Utility[8] = 0;
00453 
00454   Utility[9] = 0.0;
00455 
00456   //flags!!! not implemented yet
00457   //if (Postflag)   
00458     Utility[0] *= Utility[5];
00459 
00460   //if (Compflag) 
00461   //  Utility[0] *= Utility[7];
00462   
00463   //if (Sigflag)  
00464     Utility[0] *= Utility[8];
00465   
00466     //printf("%c] u0=%+7.3f u1=%+7.3f u2=%+7.3f u8=%+7.3f pXYs=%+7.3f no=%4d n=%4d\n",
00467     //type(), Utility[0], Utility[1], Utility[2], Utility[8], pXYs, no, n);
00468 
00469   //if (Ecoflag)
00470   //  Utility[0] *= ecoSpace();
00471   
00472   //if (Lengthflag) 
00473   //  Utility[0] /= length();
00474   //*/
00475   
00476   // Record performance in rule
00477 
00478   for (i = 0; i < 10; i++) 
00479     dblPerformance[i] = Utility[i];
00480   
00481   _pXYs = pXYs;
00482   _no = no;
00483   _dA = pYs / n;
00484   _dSig = Utility[8];
00485 
00486   return (Utility[0]);
00487 }
00488 
00489 // ==========================================================================
00490 //  GarpRule implelentation
00491 // ==========================================================================
00492 GarpRule::GarpRule() { }
00493 
00494 // ==========================================================================
00495 GarpRule::~GarpRule() { }
00496 
00497 // ==========================================================================
00498 //  RangeRule implelentation
00499 // ==========================================================================
00500 RangeRule::RangeRule() { }
00501 
00502 // ==========================================================================
00503 RangeRule::~RangeRule() { }
00504 
00505 // ==========================================================================
00506 void RangeRule::fromString(char * strRule)
00507 {
00508 }
00509 
00510 // ==========================================================================
00511 void RangeRule::initialize(EnvCellSet * objEnvCellSet, const RuleSet * objRuleSet, 
00512                bool * geneIsActivePtr, int * geneIndexPtr, int iActGenes)
00513 {
00514   int i, j, k;
00515 
00516   // call parent initialize
00517   GarpRule::initialize(objEnvCellSet, objRuleSet, geneIsActivePtr, geneIndexPtr, iActGenes);
00518 
00519   // loop iterates through variables
00520   for(k = 1; k < iActiveGenes; k++)
00521   {
00522     i = GarpUtil::randint(1, iActiveGenes);
00523   
00524     j = iGeneIndex[i];
00525     
00526     if ((j * 2 + 1 > intLength) || (j * 2 < 2))
00527       throw GarpException(100, "Array out of bounds (RangeRule::initialize)");
00528 
00529     bioclimRange(objEnvCellSet, Gene[0], GarpUtil::random() * 0.1, j);
00530   }
00531 }
00532 
00533 // ==========================================================================
00534 void RangeRule::bioclimRange(EnvCellSet * objEnvCellSet, BYTE pred, double level, int var)
00535 {
00536   BioclimHistogram * histogram;
00537   int sum, n;
00538   int UL = 0;
00539   int LL = 0;
00540     
00541   histogram = objEnvCellSet->getBioclimHistogram();
00542 
00543   sum = 0;
00544 
00545     for (n = 0; n < 256; n++)
00546     {
00547     sum += histogram->matrix[pred][var][n];
00548     if (sum>(level * histogram->matrix[pred][0][pred])) 
00549     {
00550       // printf("%d %d %f",var,sum,level);
00551       LL = n;
00552       break;
00553     }
00554   }
00555 
00556     sum = 0;    
00557     
00558   for (n = 255; n >= 0; n--)
00559     {
00560     sum += histogram->matrix[pred][var][n];
00561     if (sum > (level * histogram->matrix[pred][0][pred])) 
00562     {
00563       // printf("%d %d %f",var,sum,level);
00564       UL = n;
00565       break;
00566     } 
00567   }
00568 
00569     Gene[var * 2 + 1] = (BYTE) UL;
00570     Gene[var * 2]     = (BYTE) LL;
00571 }
00572 
00573 // ==========================================================================
00574 bool RangeRule::applyToCell(EnvCell * cell)
00575 {
00576   // visit each of the genes
00577   for (int i = 1; i < iActiveGenes; i++)
00578   {
00579     if (!((Gene[iGeneIndex[i] * 2] == 0) && (Gene[iGeneIndex[i] * 2 + 1] == 255)))
00580       if (GarpUtil::notBetween(cell->values[iGeneIndex[i]], Gene[iGeneIndex[i] * 2], Gene[iGeneIndex[i] * 2 + 1]))
00581         return false;
00582   }
00583 
00584   return true;
00585 }
00586 
00587 // ==========================================================================
00588 double RangeRule::getStrength(EnvCell * cell)
00589 {
00590   int a, b, c;
00591   int i, k; 
00592 
00593   //printf("GetStrength(%2d)\n", Gene[0]);
00594   for (k = 1; k < iActiveGenes; k++)
00595     {
00596     // get the index of the kth active gene
00597     i = iGeneIndex[k];
00598 
00599     if ((i * 2 + 1 > intLength) || (i * 2 < 2))
00600       throw GarpException(100, "Array out of bounds (RangeRule::getStrength)");
00601 
00602     a = Gene[i * 2];
00603     b = Gene[i * 2 + 1];
00604     c = cell->values[i];
00605 
00606     if (!GarpUtil::membership(a, b, c)) 
00607       { 
00608         //printf("Strength = 0\n"); 
00609         return 0; 
00610       }
00611     } 
00612 
00613   //printf("Strength = 1\n");
00614   return 1.0;
00615 }
00616 
00617 // ==========================================================================
00618 //  NegatedRangeRule implelentation
00619 // ==========================================================================
00620 NegatedRangeRule::NegatedRangeRule() {} 
00621 // ==========================================================================
00622 NegatedRangeRule::~NegatedRangeRule() {} 
00623 // ==========================================================================
00624 void NegatedRangeRule::fromString(char * strRule)
00625 {
00626 }
00627 
00628 // ==========================================================================
00629 bool NegatedRangeRule::applyToCell(EnvCell * cell)
00630 {
00631   int i, j;
00632 
00633   // visit each of the genes
00634   for (i = 1; i < iActiveGenes; i++)
00635   {
00636     j = iGeneIndex[i];
00637 
00638     if ((j * 2 + 1 > intLength) || (j * 2 < 2))
00639       throw GarpException(100, "Array out of bounds (NegatedRangeRule::applyToCell)");
00640 
00641     if (!((Gene[j * 2] == 0) && (Gene[j * 2 + 1] == 255)))
00642       if (GarpUtil::notBetween(cell->values[j], Gene[j * 2], Gene[j * 2 + 1]))
00643         return true;
00644   }
00645 
00646   return false;
00647 }
00648 
00649 // ==========================================================================
00650 double NegatedRangeRule::getStrength(EnvCell * cell)
00651 {
00652   double strength, neg_strength;
00653   
00654   strength     = RangeRule::getStrength(cell);
00655   neg_strength = 1 - strength;
00656   
00657   return neg_strength;
00658 }
00659 
00660 // ==========================================================================
00661 //  AtomicRule implelentation
00662 // ==========================================================================
00663 AtomicRule::AtomicRule() { }
00664 
00665 // ==========================================================================
00666 AtomicRule::~AtomicRule() { }
00667 
00668 // ==========================================================================
00669 void AtomicRule::initialize(EnvCellSet * objEnvCellSet, const RuleSet * objRuleSet, 
00670               bool * geneIsActivePtr, int * geneIndexPtr, int iActGenes)
00671 {
00672   int i, j, k, p;
00673 
00674   // call inherited initialize
00675   GarpRule::initialize(objEnvCellSet, objRuleSet, geneIsActivePtr, geneIndexPtr, iActGenes);
00676 
00677   p = GarpUtil::randint(0, objEnvCellSet->size());
00678   EnvCell * objEnvCell = objEnvCellSet->get(p);
00679 
00680     for (i = 1; i < iActiveGenes; i++)
00681   {
00682     j = GarpUtil::randint(1, iActiveGenes);
00683     k = iGeneIndex[j];
00684 
00685     if ((k * 2 + 1 > intLength) || (k * 2 < 2))
00686       throw GarpException(100, "Array out of bounds (AtomicRule::initialize)");
00687 
00688     Gene[2 * k] = Gene[2 * k + 1] = objEnvCell->values[k]; 
00689     }  
00690 }
00691 
00692 // ==========================================================================
00693 void AtomicRule::fromString(char * strRule)
00694 {
00695 }
00696 
00697 // ==========================================================================
00698 bool AtomicRule::applyToCell(EnvCell * cell)
00699 {
00700   int i;
00701 
00702   // visit each of the genes
00703   for (i = 1; i < iActiveGenes; i++)
00704   {
00705     if (!((Gene[iGeneIndex[i] * 2] == 0) && (Gene[iGeneIndex[i] * 2 + 1] == 255)))
00706       if (!(cell->values[iGeneIndex[i]] == Gene[iGeneIndex[i] * 2]))
00707         return false;
00708   }
00709 
00710   return true;
00711 }
00712 
00713 // ==========================================================================
00714 double AtomicRule::getStrength(EnvCell * cell)
00715 {
00716   for (int i = 1; i < iActiveGenes; i++)
00717   {
00718     if (!GarpUtil::membership(Gene[iGeneIndex[i] * 2], Gene[iGeneIndex[i] * 2 + 1], cell->values[iGeneIndex[i]])) 
00719       return 0;
00720   }
00721 
00722   return 1;
00723 }
00724 
00725 // ==========================================================================
00726 //  LogitRule implelentation
00727 // ==========================================================================
00728 LogitRule::LogitRule()
00729 {
00730 }
00731 
00732 // ==========================================================================
00733 LogitRule::~LogitRule()
00734 {
00735 }
00736 
00737 // ==========================================================================
00738 void LogitRule::fromString(char * strRule)
00739 {
00740 }
00741 
00742 // ==========================================================================
00743 void LogitRule::initialize(EnvCellSet * objEnvCellSet, const RuleSet * objRuleSet, 
00744                bool * geneIsActivePtr, int * geneIndexPtr, int iActGenes)
00745 {
00746   int i, j, k;
00747   double constant, coef[2];
00748 
00749   // call inherited initialize
00750   GarpRule::initialize(objEnvCellSet, objRuleSet, geneIsActivePtr, geneIndexPtr, iActGenes);
00751 
00752   // TO DO: must take into account the variables in use (may use some, or all)
00753     for (i = 0; i < iActiveGenes; i++)
00754   {
00755     k = GarpUtil::randint(1, iActiveGenes);
00756     j = iGeneIndex[k];
00757 
00758     regression(objEnvCellSet, j, constant, coef[0], coef[1]);
00759 
00760     Gene[1]     = (int) constant;
00761     Gene[j * 2]   = (int) coef[0]; 
00762     Gene[j * 2 + 1] = (int) coef[1];
00763     }  
00764 }
00765 
00766 // ==========================================================================
00767 int LogitRule::regression(EnvCellSet * objEnvCellSet, int dep, double& constant, double& coef1, double& coef2)
00768 {
00769   double a, b, x, y, xi, yi, xiyi, xi2, xb, xx, xxi, xxiyi, xxi2;
00770   int i, n, pred;
00771   BYTE * values;
00772 
00773   n = objEnvCellSet->count();
00774   pred = 0;
00775 
00776   a = b = x = y = xi = yi = xiyi = xi2 = xb = xx = xxi = xxiyi = xxi2 = 0.0;
00777 
00778   for (i = 0; i < n; i++) 
00779   {
00780     values = objEnvCellSet->get(i)->values;
00781 
00782     y = (double) values[pred];
00783     x = (double) values[dep];
00784     
00785     xx =  x * x;
00786     xi += x;
00787     yi += y;
00788     xiyi += x * y;
00789     xi2 += x * x;
00790 
00791     xxi += xx;
00792     xxiyi += xx * y;
00793     xxi2 += xx * xx;
00794   }
00795 
00796   b = (n * xiyi - xi * yi) / (n * xi2 - (xi * xi));
00797   coef1 = ((b * 255) + 128);
00798 
00799   a = yi / n - b * xi / n;
00800   constant = a;
00801   xb = (n * xxiyi - xxi * yi) / (n * xxi2 - (xxi * xxi));
00802   coef2 = ((xb * 255) + 128);
00803 
00804   //printf("\nb= %f %f %f %f %f %f %d",b,a,xb,coef[0],coef[1],coef[2],n);
00805   
00806   return n;
00807 }
00808 
00809 // ==========================================================================
00810 bool LogitRule::applyToCell(EnvCell * cell)
00811 { return (getStrength(cell) == 1.0); }
00812 
00813 // ==========================================================================
00814 double LogitRule::getStrength(EnvCell * cell)
00815 {
00816   BYTE * Data;
00817   int i, k;
00818   double  Sum, prob, r;
00819 
00820   Sum = 0.0;
00821   prob = 0.0;
00822 
00823   Data = cell->values;
00824 
00825   for (k = 1; k < iActiveGenes; k++)
00826   {
00827     i = iGeneIndex[k];
00828 
00829     if ((i * 2 + 1 > intLength) || (i * 2 < 2))
00830       throw GarpException(100, "Array read out of bounds (LogitRule::getStrength)");
00831 
00832     if (GarpUtil::membership(Gene[i * 2], Gene[i * 2 + 1], 1) != 255) 
00833     {
00834       r = (double) (Data[i] / 254.0);
00835 
00836       Sum += ((double) (Gene[i * 2]     - 128)) * r;
00837       Sum += ((double) (Gene[i * 2 + 1] - 128)) * r * r;
00838       
00839       //printf("\ni:%d %d C:%f %d %f", i,Rule->Gene[i*2],v,Data[i],Sum);
00840     }
00841 
00842     prob = 1.0 / (1.0 + (double) exp(-Sum));
00843 
00844     //printf("\nprob %f",prob);
00845   }
00846 
00847   if (prob > 0.5)
00848     return 1.0;
00849   else
00850     return 0.0;
00851 }
00852 
00853 // ==========================================================================
00854 void LogitRule::mutate(int intTemperature)
00855 { 
00856   int j, k;
00857 
00858     k = GarpUtil::randint(1, iActiveGenes);
00859   j = 2 * iGeneIndex[k];
00860 
00861   if ((j > intLength) || (j < 0))
00862     throw GarpException(100, "Array read out of bounds (LogitRule::mutate)");
00863 
00864   Gene[j]     = Gene[j]     + GarpUtil::randint(-intTemperature, intTemperature); 
00865   Gene[j + 1] = Gene[j + 1] + GarpUtil::randint(-intTemperature, intTemperature); 
00866 
00867   blnNeedsEvaluation = true;
00868   intGens = 0;
00869   chrOrigin = 'm';
00870 }
00871 
00872 // ==========================================================================
00873 bool LogitRule::similar(Rule * objRule)
00874 {
00875   bool found;
00876   int k;
00877 
00878   LogitRule * objOtherRule = (LogitRule *) objRule;
00879 
00880   if (type() == objOtherRule->type())
00881   {
00882     // check rule value (presence/absence)
00883     if (Gene[0] != objOtherRule->Gene[0]) 
00884       return 0;
00885 
00886     for (k = 2, found = true; (k < intLength) && (found); k ++)
00887     {
00888       if ((k > intLength) || (k < 2))
00889         throw GarpException(100, "Array read out of bounds (LogitRule::similar)");
00890 
00891       found = !( ((abs(Gene[k] - 128) < 10) && 
00892             (abs(objOtherRule->Gene[k] - 128) > 10)) || 
00893              ((abs(Gene[k] - 128)>10) && 
00894               (abs(objOtherRule->Gene[k] - 128)<10)) );
00895     }
00896 
00897     return found;
00898   }
00899 
00900   return false;
00901 }
00902 
00903 // ==========================================================================
00904 void Rule::log()
00905 {
00906   printf("<%c> ", type());
00907   for (int i = 0; i < intGenes * 2; i += 2)
00908     {
00909       //if (fabs(Genes[i] - Gene[i + 1]) >= 2.0)
00910       //Log::instance()->info( "******** ******** ");
00911       //else
00912       printf( "%3d %3d ", Gene[i], Gene[i + 1] );
00913     }
00914 
00915   printf( "- (%2d) : %f\n", Gene[0], dblPerformance[0]);
00916 }
00917 
00918 // ==========================================================================
00919