diff --git a/src/ANN.h b/src/ANN.h index df0885e..78c5d1a 100644 --- a/src/ANN.h +++ b/src/ANN.h @@ -32,7 +32,6 @@ Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 US // Activation functions inline float sigmoid(float x) { - return 1.0f / (1.0f + std::exp(-x)); } @@ -44,10 +43,16 @@ inline float linear(float x) { return x; } +// Define the tanh activation function +inline float tanh_activation(float x) { + return std::tanh(x); +} + // Activation function lookup inline float (*getActivationFunction(const std::string& name))(float) { if (name == "RELU") return relu; if (name == "SIGM") return sigmoid; + if (name == "TANH") return tanh_activation; // Added tanh if (name == "LINE") return linear; return nullptr; } @@ -138,10 +143,6 @@ struct Network { if (std::string(header) != "FFANN001") { throw std::runtime_error("Invalid file format."); } - - // int inputSize = 0; // This will be set based on the first layer's weight matrix dimensions - //std::vector> layers; // Vector to store layers polymorphically - for (int i = 0; i < numLayers; ++i) { char activation[5] = {0}; int width, height; @@ -171,7 +172,6 @@ struct Network { layer->loadWeightsAndBiases(weightData, biasData); layers.push_back(std::unique_ptr(layer)); std::cout << "adding dense "<activation == relu) act = "ReLU"; if (denseLayer->activation == sigmoid) act = "Sigmoid"; if (denseLayer->activation == linear) act = "Linear"; + if (denseLayer->activation == tanh_activation) act = "TanH"; ss << "Dense Layer: Input Size = " << denseLayer->weights[0].size() << ", Output Size = " << denseLayer->neurons.size() << ", Activation Function = " << act << "\n"; diff --git a/src/FromObsModels.h b/src/FromObsModels.h new file mode 100644 index 0000000..d0b7806 --- /dev/null +++ b/src/FromObsModels.h @@ -0,0 +1,42 @@ +/* + +Copyright (C) 2012 ForeFire Team, SPE, Universit� de Corse. + +This program is free software; you can redistribute it and/or +modify it under the terms of the GNU Lesser General Public +License as published by the Free Software Foundation; either +version 2.1 of the License, or (at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +Lesser General Public License for more details. + +You should have received a copy of the GNU Lesser General Public +License along with this program; if not, write to the Free Software +Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 US + +*/ + +#ifndef FROMOBSMODEL_H_ +#define FROMOBSMODEL_H_ + +#include "FireDomain.h" + +using namespace std; + +namespace libforefire { + +struct SensibleheatFlux{ + double flaming; + double smoldering; + + SensibleheatFlux(double f, double s) : flaming(f), smoldering(s) {} +}; + +/*! \compute heat flux from local input */ +SensibleheatFlux computeHeatFLuxFromBmap(const double&, const double&, const double&, const double&, const double&, const double&, const double&); + +} /* namespace libforefire */ + +#endif /* FROMOBSMODEL_H_ */ diff --git a/src/HeatFluxFromObsModel.cpp b/src/HeatFluxFromObsModel.cpp index 55fe369..db30eaa 100644 --- a/src/HeatFluxFromObsModel.cpp +++ b/src/HeatFluxFromObsModel.cpp @@ -21,7 +21,7 @@ Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 US #include "FluxModel.h" #include "FireDomain.h" - +#include "FromObsModels.h" using namespace std; namespace libforefire { @@ -61,70 +61,6 @@ class HeatFluxFromObsModel: public FluxModel { }; -struct SensibleheatFlux{ - double flaming; - double smoldering; - - SensibleheatFlux(double f, double s) : flaming(f), smoldering(s) {} -}; - - - -FluxModel* getHeatFluxFromObsModel(const int& = 0, DataBroker* = 0); - -/*! \compute heat flux from local input */ -SensibleheatFlux computeHeatFLuxFromBmap(const double&, const double&, const double&, const double&, const double&, const double&, const double&); - - -/* name of the model */ -const string HeatFluxFromObsModel::name = "heatFluxFromObs"; - -/* registration */ -int HeatFluxFromObsModel::isInitialized = - FireDomain::registerFluxModelInstantiator(name, getHeatFluxFromObsModel ); - -/* instantiation */ -FluxModel* getHeatFluxFromObsModel(const int& index, DataBroker* db) { - return new HeatFluxFromObsModel(index, db); -} - -/* constructor */ -HeatFluxFromObsModel::HeatFluxFromObsModel( - const int & mindex, DataBroker* db) : FluxModel(mindex, db) { - - /* defining the properties needed for the model */ - evaporationTime_data = registerProperty("FromObs_evaporationTime"); - residenceTime_data = registerProperty("FromObs_residenceTime"); - burningTime_data = registerProperty("FromObs_burningTime"); - nominalHeatFlux_f_data = registerProperty("FromObs_NominalHeatFlux_flaming"); - nominalHeatFlux_s_data = registerProperty("FromObs_NominalHeatFlux_smoldering"); - - /* allocating the vector for the values of these properties */ - if ( numProperties > 0 ) properties = new double[numProperties]; - - /* registering the model in the data broker */ - dataBroker->registerFluxModel(this); - - /* Definition of the coefficients */ - /*nominalHeatFlux = 150000.; - if ( params->isValued("nominalHeatFlux") ) - nominalHeatFlux = params->getDouble("nominalHeatFlux");*/ - -} - -/* destructor (shoudn't be modified) */ -HeatFluxFromObsModel::~HeatFluxFromObsModel() { - if ( properties != 0 ) delete properties; -} - -/* accessor to the name of the model */ -string HeatFluxFromObsModel::getName(){ - return name; -} - -/* ****************** */ -/* Model for the flux */ -/* ****************** */ SensibleheatFlux computeHeatFLuxFromBmap(const double& burningTime, const double& residenceTime, const double& nominalHeatFlux_f, const double& nominalHeatFlux_s, const double& bt, const double& et, const double& at){ @@ -186,6 +122,56 @@ SensibleheatFlux computeHeatFLuxFromBmap(const double& burningTime, const doubl } +FluxModel* getHeatFluxFromObsModel(const int& = 0, DataBroker* = 0); + +/* name of the model */ +const string HeatFluxFromObsModel::name = "heatFluxFromObs"; + +/* registration */ +int HeatFluxFromObsModel::isInitialized = + FireDomain::registerFluxModelInstantiator(name, getHeatFluxFromObsModel ); + +/* instantiation */ +FluxModel* getHeatFluxFromObsModel(const int& index, DataBroker* db) { + return new HeatFluxFromObsModel(index, db); +} + +/* constructor */ +HeatFluxFromObsModel::HeatFluxFromObsModel( + const int & mindex, DataBroker* db) : FluxModel(mindex, db) { + + /* defining the properties needed for the model */ + evaporationTime_data = registerProperty("FromObs_evaporationTime"); + residenceTime_data = registerProperty("FromObs_residenceTime"); + burningTime_data = registerProperty("FromObs_burningTime"); + nominalHeatFlux_f_data = registerProperty("FromObs_NominalHeatFlux_flaming"); + nominalHeatFlux_s_data = registerProperty("FromObs_NominalHeatFlux_smoldering"); + + /* allocating the vector for the values of these properties */ + if ( numProperties > 0 ) properties = new double[numProperties]; + + /* registering the model in the data broker */ + dataBroker->registerFluxModel(this); + + /* Definition of the coefficients */ + /*nominalHeatFlux = 150000.; + if ( params->isValued("nominalHeatFlux") ) + nominalHeatFlux = params->getDouble("nominalHeatFlux");*/ + +} + +/* destructor (shoudn't be modified) */ +HeatFluxFromObsModel::~HeatFluxFromObsModel() { + if ( properties != 0 ) delete properties; +} + +/* accessor to the name of the model */ +string HeatFluxFromObsModel::getName(){ + return name; +} + + + double HeatFluxFromObsModel::getValue(double* valueOf , const double& bt, const double& et, const double& at){ diff --git a/src/RothermelAndrews2018.cpp b/src/RothermelAndrews2018.cpp index d6b26f1..17fbea1 100644 --- a/src/RothermelAndrews2018.cpp +++ b/src/RothermelAndrews2018.cpp @@ -45,17 +45,17 @@ class RothermelAndrews2018: public PropagationModel { static int isInitialized; double windReductionFactor; /*! properties needed by the model */ - size_t wv; - size_t slope; - size_t wo; - size_t fd; - size_t fpsa; - size_t mf; - size_t pp; - size_t h; - size_t st; - size_t se; - size_t me; + size_t wv_; + size_t slope_; + size_t wo_; + size_t fd_; + size_t fpsa_; + size_t mf_; + size_t pp_; + size_t h_; + size_t st_; + size_t se_; + size_t me_; /*! result of the model */ double getSpeed(double*); @@ -89,18 +89,18 @@ RothermelAndrews2018::RothermelAndrews2018(const int & mindex, DataBroker* db) /* defining the properties needed for the model */ windReductionFactor = params->getDouble("windReductionFactor"); - wv = registerProperty("normalWind"); - slope = registerProperty("slope"); - - wo = registerProperty("fuel.fl1h_tac"); - fd = registerProperty("fuel.fd_ft"); - fpsa = registerProperty("fuel.SAVcar_ftinv"); - mf = registerProperty("fuel.mdOnDry1h_r"); - pp = registerProperty("fuel.fuelDens_lbft3"); - h = registerProperty("fuel.H_BTUlb"); - me = registerProperty("fuel.Dme_pc"); - st = registerProperty("fuel.totMineral_r"); - se = registerProperty("fuel.effectMineral_r"); + wv_ = registerProperty("normalWind"); + slope_ = registerProperty("slope"); + + wo_ = registerProperty("fuel.fl1h_tac"); + fd_ = registerProperty("fuel.fd_ft"); + fpsa_ = registerProperty("fuel.SAVcar_ftinv"); + mf_ = registerProperty("fuel.mdOnDry1h_r"); + pp_ = registerProperty("fuel.fuelDens_lbft3"); + h_ = registerProperty("fuel.H_BTUlb"); + me_ = registerProperty("fuel.Dme_pc"); + st_ = registerProperty("fuel.totMineral_r"); + se_ = registerProperty("fuel.effectMineral_r"); /* allocating the vector for the values of these properties */ if ( numProperties > 0 ) properties = new double[numProperties]; @@ -124,7 +124,7 @@ string RothermelAndrews2018::getName(){ /* Model for the propagation velovity of the front */ /* *********************************************** */ -double RothermelAndrews2018::getSpeed(double* Z){ +double RothermelAndrews2018::getSpeed(double* valueOf){ // Constants double msToftmin = 196.85039; @@ -133,35 +133,43 @@ double RothermelAndrews2018::getSpeed(double* Z){ double tacTokgm2 = 0.224; double pcTor = 0.01; - double wvC = msToftmin * Z[wv]; // convert m/s to feat/min - double meC = Z[me] * pcTor; // convert percentage to ratio - double woC = Z[wo] * tacTokgm2 / lbft2Tokgm2; // convert US short tons per acre to pounds per square foot - - if (wvC < 0) wvC = 0; - - if(woC > 0){ - double Beta_op = 3.348 * pow(Z[fpsa], -0.8189); // Optimum packing ratio - double ODBD = woC / Z[fd]; // Ovendry bulk density + // Fuel properties + double tan_slope = 0; + if (valueOf[slope_]>0){ // Only defined for null or positive slope + tan_slope = valueOf[slope_]; + } + double wv = msToftmin * valueOf[wv_]; // convert m/s to feat/min + double me = valueOf[me_] * pcTor; // convert percentage to ratio + double pp = valueOf[pp_]; + double mf = valueOf[mf_]; + double fpsa = valueOf[fpsa_]; + double fd = valueOf[fd_]; + double wo = valueOf[wo_] * tacTokgm2 / lbft2Tokgm2; // convert US short tons per acre to pounds per square foot + double h = valueOf[h_]; + double st = valueOf[st_]; + double se = valueOf[se_]; + + if (wv < 0) wv = 0; + + if(wo > 0){ + double Beta_op = 3.348 * pow(fpsa, -0.8189); // Optimum packing ratio + double ODBD = wo / fd; // Ovendry bulk density double Beta = ODBD / pp; // Packing ratio - double WN = woC / (1 + Z[st]); // Net fuel loading + double WN = wo / (1 + st); // Net fuel loading double A = 133.0 / pow(fpsa, 0.7913); // updated A double T_max = pow(fpsa,1.5) * pow(495.0 + 0.0594 * pow(fpsa, 1.5),-1.0); // Maximum reaction velocity double T = T_max * pow((Beta / Beta_op), A) * exp(A * (1 - Beta / Beta_op)); // Optimum reaction velocity - double NM = 1. - 2.59 * (Z[mf] / meC) + 5.11 * pow(mf / meC, 2.) - 3.52 * pow(mf / meC,3.); // Moisture damping coeff. - double NS = 0.174 * pow(Z[se], -0.19); // Mineral damping coefficient - double RI = T * WN * Z[h] * NM * NS; + double NM = 1. - 2.59 * (mf / me) + 5.11 * pow(mf / me, 2.) - 3.52 * pow(mf / me,3.); // Moisture damping coeff. + double NS = 0.174 * pow(se, -0.19); // Mineral damping coefficient + double RI = T * WN * h * NM * NS; double PFR = pow(192.0 + 0.2595 * fpsa, -1) * exp((0.792 + 0.681 * pow(fpsa, 0.5)) * (Beta + 0.1)); // Propogating flux ratio // Wind Coefficien t double B = 0.02526 * pow(fpsa, 0.54); double C = 7.47 * exp(-0.1333 * pow(fpsa, 0.55)); double E = 0.715 * exp(-3.59 * pow(10, -4) * fpsa); - if (wvC > (0.9 * RI)) wvC = 0.9 * RI; - double WC = (C * pow(wvC, B)) * pow((Beta / Beta_op), (-E)); - - double SC = 0; - if (Z[slope] >0){ // Rothermel only for upslope, assuming no slope if negative slope is encountered - SC = 5.275*(pow(Beta, -0.3))*pow(Z[slope], 2); - } + if (wv > (0.9 * RI)) wv = 0.9 * RI; + double WC = (C * pow(wv, B)) * pow((Beta / Beta_op), (-E)); + double SC = 5.275*(pow(Beta, -0.3))*pow(tan_slope, 2); // Heat sink double EHN = exp(-138. / fpsa); // Effective Heating Number = f(surface are volume ratio) double QIG = 250. + 1116. * mf; // Heat of preignition= f(moisture content) diff --git a/src/ScalarFromObsModel.cpp b/src/ScalarFromObsModel.cpp new file mode 100644 index 0000000..9aa998b --- /dev/null +++ b/src/ScalarFromObsModel.cpp @@ -0,0 +1,159 @@ +/* + +Copyright (C) 2012 ForeFire Team, SPE, UniversitŽ de Corse. + +This program is free software; you can redistribute it and/or +modify it under the terms of the GNU Lesser General Public +License as published by the Free Software Foundation; either +version 2.1 of the License, or (at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +Lesser General Public License for more details. + +You should have received a copy of the GNU Lesser General Public +License along with this program; if not, write to the Free Software +Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 US + +*/ + + + +#include "FluxModel.h" +#include "FireDomain.h" +#include "FromObsModels.h" + +using namespace std; + +namespace libforefire { + +class ScalarFromObsModel: public FluxModel { + + /*! name the model */ + static const string name; + + /*! boolean for initialization */ + static int isInitialized; + + /*! properties needed by the model */ + size_t nominalHeatFlux_f_data; + size_t nominalHeatFlux_s_data; + + size_t evaporationTime_data; + size_t residenceTime_data; + size_t burningTime_data; + + size_t EFScalar_f_data; + size_t EFScalar_s_data; + + size_t radiation_fraction_f_data; + size_t radiation_fraction_s_data; + + size_t conversion_factor_data; + + /*! coefficients needed by the model */ + + /*! local variables */ + + /*! result of the model */ + double getValue(double*, const double& + , const double&, const double&); + +public: + ScalarFromObsModel(const int& = 0, DataBroker* = 0); + virtual ~ScalarFromObsModel(); + + string getName(); +}; + + +FluxModel* getScalarFromObsModel(const int& = 0, DataBroker* = 0); + +/* name of the model */ +const string ScalarFromObsModel::name = "SFObs"; + +/* registration */ +int ScalarFromObsModel::isInitialized = + FireDomain::registerFluxModelInstantiator(name, getScalarFromObsModel ); + +/* instantiation */ +FluxModel* getScalarFromObsModel(const int& index, DataBroker* db) { + return new ScalarFromObsModel(index, db); +} + +/* constructor */ +ScalarFromObsModel::ScalarFromObsModel( + const int & mindex, DataBroker* db) : FluxModel(mindex, db) { + + /* defining the properties needed for the model */ + nominalHeatFlux_f_data = registerProperty("FromObs_NominalHeatFlux_flaming"); + nominalHeatFlux_s_data = registerProperty("FromObs_NominalHeatFlux_smoldering"); + + evaporationTime_data = registerProperty("FromObs_evaporationTime"); + residenceTime_data = registerProperty("FromObs_residenceTime"); + burningTime_data = registerProperty("FromObs_burningTime"); + + EFScalar_f_data = registerProperty("FromObs_EFScalar_flaming"); + EFScalar_s_data = registerProperty("FromObs_EFScalar_smoldering"); + + radiation_fraction_f_data = registerProperty("FromObs_radiationFraction_flaming"); + radiation_fraction_s_data = registerProperty("FromObs_radiationFraction_smoldering"); + + conversion_factor_data = registerProperty("FromObs_conversionFactor"); + + /* allocating the vector for the values of these properties */ + if ( numProperties > 0 ) properties = new double[numProperties]; + + /* registering the model in the data broker */ + dataBroker->registerFluxModel(this); + + /* Definition of the coefficients */ +} + +/* destructor (shoudn't be modified) */ +ScalarFromObsModel::~ScalarFromObsModel() { + if ( properties != 0 ) delete properties; +} + +/* accessor to the name of the model */ +string ScalarFromObsModel::getName(){ + return name; +} + +/* ****************** */ +/* Model for the flux */ +/* ****************** */ +double ScalarFromObsModel::getValue(double* valueOf + , const double& bt, const double& et, const double& at){ + /* Mean heat flux released between the time interval [bt, et] */ + /* The heat flux is supposed to be constant from the arrival time (at) + * and for a period of time of 'burningDuration', constant of the model */ + + double nominalHeatFlux_f = valueOf[nominalHeatFlux_f_data]; + double nominalHeatFlux_s = valueOf[nominalHeatFlux_s_data]; + + double evaporationTime = valueOf[evaporationTime_data]; + double residenceTime = valueOf[residenceTime_data]; + double burningTime = valueOf[burningTime_data]; + + double EFScalar_f = valueOf[EFScalar_f_data]; + double EFScalar_s = valueOf[EFScalar_s_data]; + + double radiation_fraction_f = valueOf[radiation_fraction_f_data]; + double radiation_fraction_s = valueOf[radiation_fraction_s_data]; + + double conversion_factor = valueOf[conversion_factor_data]; + + SensibleheatFlux sensibleheatFlux = computeHeatFLuxFromBmap(burningTime,residenceTime, nominalHeatFlux_f,nominalHeatFlux_s, bt, et, evaporationTime+at); + + /*double ScalarFlux = EFScalar * (HeatFlux/(1.-radiation_fraction)); */ + double ScalarFlux = EFScalar_f * conversion_factor * radiation_fraction_f * 1.e-6 * sensibleheatFlux.flaming / (1.-radiation_fraction_f) + + EFScalar_s * conversion_factor * radiation_fraction_s * 1.e-6 * sensibleheatFlux.smoldering / (1.-radiation_fraction_s); + + return ScalarFlux; + + +} + +} /* namespace libforefire */