00001 // @(#)root/mathmore:$Id: GSLMCIntegrationWorkspace.h 36806 2010-11-20 11:09:14Z moneta $ 00002 // Author: Magdalena Slawinska 08/2007 00003 00004 /********************************************************************** 00005 * * 00006 * Copyright (c) 2007 ROOT Foundation, CERN/PH-SFT * 00007 * * 00008 * This library is free software; you can redistribute it and/or * 00009 * modify it under the terms of the GNU General Public License * 00010 * as published by the Free Software Foundation; either version 2 * 00011 * of the License, or (at your option) any later version. * 00012 * * 00013 * This library is distributed in the hope that it will be useful, * 00014 * but WITHOUT ANY WARRANTY; without even the implied warranty of * 00015 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * 00016 * General Public License for more details. * 00017 * * 00018 * You should have received a copy of the GNU General Public License * 00019 * along with this library (see file COPYING); if not, write * 00020 * to the Free Software Foundation, Inc., 59 Temple Place, Suite * 00021 * 330, Boston, MA 02111-1307 USA, or contact the author. * 00022 * * 00023 **********************************************************************/ 00024 00025 // Header file for class GSLIntegratorWorkspace 00026 // 00027 // Author: Magdalena Slawinska 00028 // 00029 00030 00031 00032 #ifndef ROOT_Math_GSLMCIntegrationWorkspace 00033 #define ROOT_Math_GSLMCIntegrationWorkspace 00034 00035 #include "gsl/gsl_math.h" 00036 #include "gsl/gsl_monte.h" 00037 #include "gsl/gsl_monte_vegas.h" 00038 #include "gsl/gsl_monte_miser.h" 00039 #include "gsl/gsl_monte_plain.h" 00040 00041 #include "Math/MCParameters.h" 00042 #include "Math/MCIntegrationTypes.h" 00043 00044 namespace ROOT { 00045 namespace Math { 00046 00047 00048 00049 class GSLMCIntegrationWorkspace { 00050 00051 public : 00052 00053 GSLMCIntegrationWorkspace() {} 00054 00055 virtual ~GSLMCIntegrationWorkspace() { Clear(); } 00056 00057 virtual MCIntegration::Type Type() const = 0; 00058 00059 virtual size_t NDim() const { return 0; } 00060 00061 /// initialize the workspace creating the GSL pointer if it is not there 00062 virtual bool Init(size_t dim) = 0; 00063 00064 /// re-initialize an existing the workspace 00065 virtual bool ReInit() = 0; 00066 00067 /// free the workspace deleting the GSL pointer 00068 virtual void Clear() {} 00069 00070 /// retrieve option pointer corresponding to parameters 00071 /// create a new object to be managed by the user 00072 virtual ROOT::Math::IOptions * Options() const = 0; 00073 00074 private: 00075 00076 00077 }; 00078 00079 /** 00080 workspace for VEGAS 00081 */ 00082 class GSLVegasIntegrationWorkspace : public GSLMCIntegrationWorkspace { 00083 00084 public : 00085 00086 GSLVegasIntegrationWorkspace(size_t dim = 0) : 00087 fWs(0) 00088 { 00089 if (dim > 0) Init(dim); 00090 } 00091 00092 bool Init(size_t dim) { 00093 fWs = gsl_monte_vegas_alloc( dim); 00094 if (fWs) SetVegasParameters(); 00095 return (fWs != 0); 00096 } 00097 00098 bool ReInit() { 00099 // according to the code - reinit just reset default GSL values 00100 if (!fWs) return false; 00101 int iret = gsl_monte_vegas_init( fWs ); 00102 SetVegasParameters(); 00103 return (iret == 0); 00104 } 00105 00106 void Clear() { 00107 if (fWs) gsl_monte_vegas_free( fWs); 00108 fWs = 0; 00109 } 00110 00111 gsl_monte_vegas_state * GetWS() { return fWs; } 00112 00113 void SetParameters(const struct VegasParameters &p) { 00114 fParams = p; 00115 if (fWs) SetVegasParameters(); 00116 } 00117 00118 size_t NDim() const { return (fWs) ? fWs->dim : 0; } 00119 00120 double Result() const { return (fWs) ? fWs->result : -1;} 00121 00122 double Sigma() const { return (fWs) ? fWs->sigma : 0;} 00123 00124 double Chisq() const { return (fWs) ? fWs->chisq: -1;} 00125 00126 MCIntegration::Type Type() const { return MCIntegration::kVEGAS; } 00127 00128 const VegasParameters & Parameters() const { return fParams; } 00129 VegasParameters & Parameters() { return fParams; } 00130 00131 virtual ROOT::Math::IOptions * Options() const { 00132 return fParams(); 00133 } 00134 00135 00136 private: 00137 00138 void SetVegasParameters() { 00139 fWs->alpha = fParams.alpha; 00140 fWs->iterations = fParams.iterations; 00141 fWs->stage = fParams.stage; 00142 fWs->mode = fParams.mode; 00143 fWs->verbose = fParams.verbose; 00144 } 00145 00146 00147 gsl_monte_vegas_state * fWs; 00148 VegasParameters fParams; 00149 00150 }; 00151 00152 00153 /** 00154 Workspace for MISER 00155 */ 00156 class GSLMiserIntegrationWorkspace : public GSLMCIntegrationWorkspace { 00157 00158 public : 00159 00160 GSLMiserIntegrationWorkspace(size_t dim = 0) : 00161 fHaveNewParams(false), 00162 fWs(0) 00163 { 00164 if (dim > 0) Init(dim); 00165 } 00166 00167 00168 bool Init(size_t dim) { 00169 fWs = gsl_monte_miser_alloc( dim); 00170 // need this to set parameters according to dimension 00171 if (!fHaveNewParams) fParams = MiserParameters(dim); 00172 if (fWs) SetMiserParameters(); 00173 return (fWs != 0); 00174 } 00175 00176 bool ReInit() { 00177 // according to the code - reinit just reset default GSL values 00178 if (!fWs) return false; 00179 int iret = gsl_monte_miser_init( fWs ); 00180 SetMiserParameters(); 00181 return (iret == 0); 00182 } 00183 00184 void Clear() { 00185 if (fWs) gsl_monte_miser_free( fWs); 00186 fWs = 0; 00187 } 00188 00189 gsl_monte_miser_state * GetWS() { return fWs; } 00190 00191 void SetParameters(const MiserParameters &p) { 00192 fParams = p; 00193 fHaveNewParams = true; 00194 if (fWs) SetMiserParameters(); 00195 } 00196 00197 size_t NDim() const { return (fWs) ? fWs->dim : 0; } 00198 00199 MCIntegration::Type Type() const { return MCIntegration::kMISER; } 00200 00201 00202 const MiserParameters & Parameters() const { return fParams; } 00203 MiserParameters & Parameters() { return fParams; } 00204 00205 virtual ROOT::Math::IOptions * Options() const { 00206 return fParams(); 00207 } 00208 00209 private: 00210 00211 void SetMiserParameters() 00212 { 00213 fWs->estimate_frac = fParams.estimate_frac; 00214 fWs->min_calls = fParams.min_calls; 00215 fWs->min_calls_per_bisection = fParams.min_calls_per_bisection; 00216 fWs->alpha = fParams.alpha; 00217 fWs->dither = fParams.dither; 00218 } 00219 00220 00221 bool fHaveNewParams; 00222 gsl_monte_miser_state * fWs; 00223 MiserParameters fParams; 00224 00225 }; 00226 00227 00228 00229 00230 class GSLPlainIntegrationWorkspace : public GSLMCIntegrationWorkspace{ 00231 00232 public : 00233 00234 GSLPlainIntegrationWorkspace() : 00235 fWs(0) 00236 { } 00237 00238 bool Init(size_t dim) { 00239 fWs = gsl_monte_plain_alloc( dim); 00240 // no parameter exists for plain 00241 return (fWs != 0); 00242 } 00243 00244 bool ReInit() { 00245 if (!fWs) return false; 00246 int iret = gsl_monte_plain_init( fWs ); 00247 return (iret == GSL_SUCCESS); 00248 } 00249 00250 void Clear() { 00251 if (fWs) gsl_monte_plain_free( fWs); 00252 fWs = 0; 00253 } 00254 00255 gsl_monte_plain_state * GetWS() { return fWs; } 00256 00257 //void SetParameters(const struct PlainParameters &p); 00258 00259 MCIntegration::Type Type() const { return MCIntegration::kPLAIN; } 00260 00261 size_t NDim() const { return (fWs) ? fWs->dim : 0; } 00262 00263 virtual ROOT::Math::IOptions * Options() const { 00264 return 0; 00265 } 00266 00267 00268 private: 00269 00270 gsl_monte_plain_state * fWs; 00271 00272 00273 }; 00274 00275 00276 } // namespace Math 00277 } // namespace ROOT 00278 00279 00280 00281 00282 #endif /* ROOT_Math_GSLMCIntegrationWorkspace */