openModeller
Version 1.4.0
|
00001 // 00002 // CsmKG 00003 // 00004 // Description: Csm implementation using Keiser-Gutman cutoff to discard components 00005 // 00006 // 00007 // Author: CRIA <t.sutton@reading.ac.uk>, (C) 2004 00008 // 00009 // Copyright: See COPYING file that comes with this distribution 00010 // 00011 // 00012 00013 #include <string.h> 00014 #include "csmkg.hh" 00015 00016 #include <gsl/gsl_statistics_double.h> 00017 #include <gsl/gsl_multifit_nlin.h> //remove this when we have proper covar matrix 00018 #include <gsl/gsl_math.h> 00019 #include <gsl/gsl_eigen.h> 00020 #include <gsl/gsl_sf_gamma.h> 00021 00022 #include <math.h> 00023 00024 /****************************************************************/ 00025 /********************** Algorithm's Metadata ********************/ 00026 00027 #define NUM_PARAM 0 00028 00029 00030 /************************************/ 00031 /*** Algorithm parameter metadata ***/ 00032 00033 static AlgParamMetadata *parameters = 0; 00034 /* 00035 AlgParamMetadata parameters[NUM_PARAM] = { 00036 00037 // Metadata of the first parameter. 00038 { 00039 "First parameter id", // Id. 00040 "First parameter name", // Name. 00041 "First parameter type", // Type. 00042 "First parameter overview", // Overview. 00043 "First parameter description", // Description. 00044 00045 1, // Not zero if the parameter has lower limit. 00046 0.0, // Parameter's lower limit. 00047 1, // Not zero if the parameter has upper limit. 00048 1.0, // Parameter's upper limit. 00049 "0.1" // Parameter's typical (default) value. 00050 }, 00051 }; 00052 */ 00053 00054 00055 /************************************/ 00056 /*** Algorithm's general metadata ***/ 00057 00058 static AlgMetadata metadata = { 00059 00060 "CSMKG", // Id. 00061 "Climate Space Model - Kaiser-Gutman", // Name. 00062 "0.1 alpha", // Version. 00063 00064 // Overview 00065 "Climate Space Model [CSM] is a principle components based \ 00066 algorithm developed by Dr. Neil Caithness", 00067 00068 // Description. 00069 "Climate Space Model [CSM] is a principle components based \ 00070 algorithm developed by Dr. Neil Caithness. The component \ 00071 selection process int this algorithm implementation is \ 00072 based on the Keiser-Gutman cutoff where any component with \ 00073 an eigenvalue < 1 is discarded.\ 00074 \n \ 00075 The original CSM was written as series of Matlab functions. \ 00076 ", 00077 00078 "Neil Caithness", // Author 00079 "", // Bibliography. 00080 00081 "Tim Sutton, Renato De Giovanni", // Code author. 00082 "t.sutton [at] reading.ac.uk", // Code author's contact. 00083 00084 0, // Does not accept categorical data. 00085 0, // Does not need (pseudo)absence points. 00086 00087 NUM_PARAM, // Algorithm's parameters. 00088 parameters 00089 }; 00090 00091 00092 00093 /****************************************************************/ 00094 /****************** Algorithm's factory function ****************/ 00095 00096 dllexp 00097 AlgorithmImpl * 00098 algorithmFactory() 00099 { 00100 return new CsmKG(); 00101 } 00102 00103 dllexp 00104 AlgMetadata const * 00105 algorithmMetadata() 00106 { 00107 return &metadata; 00108 } 00109 00110 /****************************************************************/ 00111 /****************************** Csm *****************************/ 00112 00117 CsmKG::CsmKG(): Csm(&metadata) 00118 { 00119 _initialized = 0; 00120 } 00121 00122 00124 CsmKG::~CsmKG() 00125 { 00126 //This is handled by the parent class I think... 00127 /* 00128 if ( _initialized ) 00129 { 00130 gsl_matrix_free (_gsl_environment_matrix); 00131 gsl_matrix_free (_gsl_covariance_matrix); 00132 gsl_vector_free (_gsl_avg_vector); 00133 gsl_vector_free (_gsl_stddev_vector); 00134 gsl_vector_free (_gsl_eigenvalue_vector); 00135 gsl_matrix_free (_gsl_eigenvector_matrix); 00136 } 00137 */ 00138 } 00139 int CsmKG::discardComponents() 00140 { 00141 printf ("vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv\n"); 00142 printf (" Starting CSM - Kaiser-Gutman \n"); 00143 printf (" Component discarding routine \n"); 00144 printf ("vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv\n"); 00145 // Discard any columns from the eigenvector where the eigenvalue is < 1 00146 // (Keiser-Gutman method) 00147 // because the eigenvalues are sorted, we can discard any columns after the 00148 // first one is encountered. 00149 00150 int myColumnNo = -1; 00151 // for debuggin print out the averages 00152 for (int j=0;j<_layer_count;j++) 00153 { 00154 float myValue = gsl_vector_get (_gsl_eigenvalue_vector,j); 00155 if (myValue < 1) 00156 { 00157 myColumnNo = j; 00158 } 00159 } 00160 if (myColumnNo > -1) 00161 { 00162 int j; 00163 _retained_components_count = myColumnNo-1; 00164 printf ("\n\nNumber of components retained: %i of %i\n", 00165 _retained_components_count, 00166 _layer_count); 00167 //so now create a local copy of the eigvect and eigval... 00168 //first clone the vector and matrix... 00169 gsl_vector * tmp_gsl_eigenvalue_vector = gsl_vector_alloc (_layer_count); 00170 gsl_matrix * tmp_gsl_eigenvector_matrix = gsl_matrix_alloc (_layer_count, _layer_count); 00171 00172 gsl_matrix_memcpy (tmp_gsl_eigenvector_matrix, _gsl_eigenvector_matrix); 00173 gsl_vector_memcpy (tmp_gsl_eigenvalue_vector,_gsl_eigenvalue_vector); 00174 00175 //now clear our current eig vec and matrix 00176 00177 gsl_vector_free (_gsl_eigenvalue_vector); 00178 gsl_matrix_free (_gsl_eigenvector_matrix); 00179 00180 //next we reassign them to the reduced size (ie without unwanted components) 00181 _gsl_eigenvalue_vector = gsl_vector_alloc (_retained_components_count); 00182 _gsl_eigenvector_matrix = gsl_matrix_alloc (_layer_count, _retained_components_count); 00183 00184 //now copy over just the components we intend to keep 00185 //...first the vector 00186 for (j=0;j<_retained_components_count;j++) 00187 { 00188 float myFloat = gsl_vector_get (tmp_gsl_eigenvalue_vector,j); 00189 gsl_vector_set (_gsl_eigenvalue_vector,j,myFloat); 00190 } 00191 //...now the matrix 00192 gsl_vector * tmp_gsl_vector = gsl_vector_alloc (_layer_count); 00193 for (j=0;j<_retained_components_count;j++) 00194 { 00195 gsl_matrix_get_col (tmp_gsl_vector, tmp_gsl_eigenvector_matrix, j); 00196 gsl_matrix_set_col (_gsl_eigenvector_matrix,j,tmp_gsl_vector); 00197 } 00198 //now clear away the temporary vars 00199 gsl_vector_free (tmp_gsl_vector); 00200 gsl_vector_free (tmp_gsl_eigenvalue_vector); 00201 gsl_matrix_free (tmp_gsl_eigenvector_matrix); 00202 00203 00204 } 00205 else 00206 { 00207 //we keep all the components! 00208 _retained_components_count = _layer_count; 00209 printf ("\n\nNumber of components retained: %i of %i\n", 00210 _retained_components_count, 00211 _layer_count); 00212 } 00213 return 1; 00214 }