openModeller
Version 1.4.0
|
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