Skip to content

Commit

Permalink
calibration coeff in ini
Browse files Browse the repository at this point in the history
  • Loading branch information
gantolini committed Jan 24, 2024
1 parent aa9e956 commit 903fd16
Show file tree
Hide file tree
Showing 4 changed files with 149 additions and 113 deletions.
22 changes: 14 additions & 8 deletions DATA/PROJECT/testFrostForecast/testFrostForecastSettings.ini
Original file line number Diff line number Diff line change
Expand Up @@ -7,19 +7,25 @@ utm_zone=32

[project]
csvPath=./outputData/
id="4,10,2303,6257,6265"
id=4,10,2303,6257,6265
meteoGrid=../../METEOGRID/DBGridXML_ERG5_forecast.xml
meteo_points=../../METEOPOINT/frost.db
name=testFrostForecast
var="TAVG,RHAVG,W_SCAL_INT,W_VEC_DIR,PREC,RAD"
var=TAVG,RHAVG,W_SCAL_INT,W_VEC_DIR,PREC,RAD

[reuter_param]
SE_intercept="0.182930, 0.249670, 0.169418, 0.332520, 0.274008"
SE_parRHss="0.001691, 0.002694, 0.001781, 0.003540, 0.002943"
SE_parTss="0.012427, 0.014025, 0.009732, 0.018650, 0.014905"
intercept="1.969822, 1.775462, 3.016488, 2.271830, 2.367134"
parRHss="-0.014889, -0.013332, -0.014434, -0.020210, -0.011935"
parTss="0.093746, 0.129308, 0.079078, 0.133310, 0.106645"
SE_intercept=0.182930, 0.249670, 0.169418, 0.332520, 0.274008
SE_parRHss=0.001691, 0.002694, 0.001781, 0.003540, 0.002943
SE_parTss=0.012427, 0.014025, 0.009732, 0.018650, 0.014905
intercept=1.969822, 1.775462, 3.016488, 2.271830, 2.367134
parRHss=-0.014889, -0.013332, -0.014434, -0.020210, -0.011935
parTss=0.093746, 0.129308, 0.079078, 0.133310, 0.106645

[calibration]
monthIni=3
monthFin=4
thresholdTmin=5
thresholdTrange=6

[software]
software=frostForecast
190 changes: 113 additions & 77 deletions bin/frostForecast/frost.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,11 @@ void Frost::initialize()
SE_intercept.clear();
SE_parTss.clear();
SE_parRHss.clear();

monthIni = 1;
monthFin = 12;
thresholdTmin = 0;
thresholdTrange = 10;
}

void Frost::setSettingsFileName(const QString &value)
Expand Down Expand Up @@ -143,92 +148,101 @@ int Frost::readSettings()
csvFilePath = path + QDir::cleanPath(csvFilePath);
}

QString idTemp = projectSettings->value("id","").toString();
if (idTemp.isEmpty())
idList = projectSettings->value("id").toStringList();
if (idList.isEmpty())
{
logger.writeError ("missing id station");
return ERROR_MISSINGPARAMETERS;
}
else
{
idList = idTemp.split(",");
}

QString varTemp = projectSettings->value("var","").toString();
if (varTemp.isEmpty())
varList = projectSettings->value("var").toStringList();
if (varList.isEmpty())
{
logger.writeError ("missing var");
return ERROR_MISSINGPARAMETERS;
}
else
{
varList = varTemp.split(",");
}

projectSettings->endGroup();

QStringList myList;

// reuter_param
projectSettings->beginGroup("reuter_param");
QString interceptTemp = projectSettings->value("intercept","").toString();
if (interceptTemp.isEmpty())

myList = projectSettings->value("intercept").toStringList();
if (myList.isEmpty())
{
logger.writeError ("missing reuter param intercept");
return ERROR_MISSINGPARAMETERS;
}
else
{
intercept = interceptTemp.split(",");
intercept = StringListToFloat(myList);
}
QString parTssTemp = projectSettings->value("parTss","").toString();
if (parTssTemp.isEmpty())

myList = projectSettings->value("parTss").toStringList();
if (myList.isEmpty())
{
logger.writeError ("missing reuter param parTss");
return ERROR_MISSINGPARAMETERS;
}
else
{
parTss = parTssTemp.split(",");
parTss = StringListToFloat(myList);
}
QString parRHssTemp = projectSettings->value("parRHss","").toString();
if (parRHssTemp.isEmpty())

myList = projectSettings->value("parRHss").toStringList();
if (myList.isEmpty())
{
logger.writeError ("missing reuter param parRHss");
return ERROR_MISSINGPARAMETERS;
}
else
{
parRHss = parRHssTemp.split(",");
parRHss = StringListToFloat(myList);
}
QString SEinterceptTemp = projectSettings->value("SE_intercept","").toString();
if (SEinterceptTemp.isEmpty())

myList = projectSettings->value("SE_intercept").toStringList();
if (myList.isEmpty())
{
logger.writeError ("missing reuter param SE_intercept");
return ERROR_MISSINGPARAMETERS;
}
else
{
SE_intercept = SEinterceptTemp.split(",");
SE_intercept = StringListToFloat(myList);
}
QString SEparTssTemp = projectSettings->value("SE_parTss","").toString();
if (SEparTssTemp.isEmpty())

myList = projectSettings->value("SE_parTss").toStringList();
if (myList.isEmpty())
{
logger.writeError ("missing reuter param SE_parTss");
return ERROR_MISSINGPARAMETERS;
}
else
{
SE_parTss = SEparTssTemp.split(",");
SE_parTss = StringListToFloat(myList);
}
QString SEparRHssTemp = projectSettings->value("SE_parRHss","").toString();
if (SEparRHssTemp.isEmpty())

myList = projectSettings->value("SE_parRHss").toStringList();
if (myList.isEmpty())
{
logger.writeError ("missing reuter param SE_parRHss");
return ERROR_MISSINGPARAMETERS;
}
else
{
SE_parRHss = SEparRHssTemp.split(",");
SE_parRHss = StringListToFloat(myList);
}
projectSettings->endGroup();

// calibration
projectSettings->beginGroup("calibration");
monthIni = projectSettings->value("monthIni").toInt();
monthFin = projectSettings->value("monthFin").toInt();
thresholdTmin = projectSettings->value("thresholdTmin").toFloat();
thresholdTrange = projectSettings->value("thresholdTrange").toFloat();
projectSettings->endGroup();

return FROSTFORECAST_OK;
}
Expand Down Expand Up @@ -385,12 +399,12 @@ int Frost::getForecastData(QString id, int posIdList)
int myTmpHour;
QDateTime dateTimeTmp;

float myIntercept = intercept[posIdList].toFloat();
float myParTss = parTss[posIdList].toFloat();
float myParRHss = parRHss[posIdList].toFloat();
float mySEintercept = SE_intercept[posIdList].toFloat();
float mySEparTss = SE_parTss[posIdList].toFloat();
float mySEparRHss = SE_parRHss[posIdList].toFloat();
float myIntercept = intercept[posIdList];
float myParTss = parTss[posIdList];
float myParRHss = parRHss[posIdList];
float mySEintercept = SE_intercept[posIdList];
float mySEparTss = SE_parTss[posIdList];
float mySEparRHss = SE_parRHss[posIdList];

for (int i = 0; i <= 1; i++)
{
Expand Down Expand Up @@ -550,8 +564,51 @@ QList<QString> Frost::getIdList() const
return idList;
}

bool Frost::getRadiativeCoolingHistory(QString id, float thresholdTmin, float thresholdTRange, int monthIni, int monthFin, int timeZone,
std::vector<std::vector<float>>& outData, std::vector <std::vector <float>>& sunsetData)
void fitCoolingCoefficient(std::vector <std::vector <float>>& hourly_series, int nrIterations, float tolerance, float minValue, float maxValue, std::vector <float>& coeff)
{
float myMin, myMax;
float mySum1, mySum2;
float myFirstError = 0;
float mySecondError = 0;;
int iteration;
float Tsunset;
int time;

for(std::vector <float> hourlyT : hourly_series)
{

myMin = minValue;
myMax = maxValue;
iteration = 0;
time = 0;
do {

mySum1 = 0;
mySum2 = 0;
Tsunset = hourlyT[0];
time = 0;
for (float T : hourlyT)
{
mySum1 += pow((Tsunset - myMin * sqrt(time) - T), 2);
mySum2 += pow((Tsunset - myMax * sqrt(time) - T), 2);
time++;
}

myFirstError = sqrt(mySum1);
mySecondError = sqrt(mySum2);

if (myFirstError < mySecondError)
myMax = (myMin + myMax) / 2;
else
myMin = (myMin + myMax) / 2;

} while (fabs(myFirstError - mySecondError) > tolerance && iteration <= nrIterations);

coeff.push_back((myMin + myMax) / 2);
}
}

bool Frost::getRadiativeCoolingHistory(QString id, std::vector<std::vector<float>>& outData, std::vector <std::vector <float>>& sunsetData)
{
unsigned pos = 0;
bool found = false;
Expand Down Expand Up @@ -616,7 +673,7 @@ bool Frost::getRadiativeCoolingHistory(QString id, float thresholdTmin, float th
{
hourSunSetLocal = round(sunPosition.set/3600);
hourSunRiseLocal = round(sunPosition.rise/3600);
hourSunRiseUtc = hourSunRiseLocal - timeZone;
hourSunRiseUtc = hourSunRiseLocal - gisSettings.timeZone;
dateIndexSunRiseUtc = i;

if (hourSunRiseUtc < 1)
Expand All @@ -630,7 +687,7 @@ bool Frost::getRadiativeCoolingHistory(QString id, float thresholdTmin, float th
dateIndexSunRiseUtc--;
}

h = hourSunSetLocal - timeZone;
h = hourSunSetLocal - gisSettings.timeZone;
d = i-1;

isSunRise = false;
Expand Down Expand Up @@ -691,7 +748,7 @@ bool Frost::getRadiativeCoolingHistory(QString id, float thresholdTmin, float th
if (h == hourSunRiseUtc && d == dateIndexSunRiseUtc)
{
isSunRise = true;
if (dataPresent && minT < thresholdTmin && (maxT - minT >= thresholdTRange))
if (dataPresent && minT < thresholdTmin && (maxT - minT >= thresholdTrange))
{
outData.push_back(hourlyT);
sunsetData.push_back(ssData);
Expand All @@ -709,47 +766,26 @@ bool Frost::getRadiativeCoolingHistory(QString id, float thresholdTmin, float th
return (outData.size() > 0);
}


void Frost::fitCoolingCoefficient(std::vector <std::vector <float>>& hourly_series, int nrIterations, float tolerance, float minValue, float maxValue, std::vector <float>& coeff)
bool Frost::calibrateModel(QString id)
{
float myMin, myMax;
float mySum1, mySum2;
float myFirstError = 0;
float mySecondError = 0;;
int iteration;
float Tsunset;
int time;
std::vector <std::vector <float>> outData;
std::vector <std::vector <float>> sunsetData;
std::vector <float> outCoeff;
std::vector <float> weights;
float regrConst, R2, stdError;
std::vector <float> regrCoeff;

for(std::vector <float> hourlyT : hourly_series)
if (getRadiativeCoolingHistory(id, outData, sunsetData))
{
fitCoolingCoefficient(outData, 10000, float(0.001), float(0), float(10), outCoeff);

myMin = minValue;
myMax = maxValue;
iteration = 0;
time = 0;
do {
for (int j=0; j < sunsetData.size(); j++)
weights.push_back(1);

mySum1 = 0;
mySum2 = 0;
Tsunset = hourlyT[0];
time = 0;
for (float T : hourlyT)
{
mySum1 += pow((Tsunset - myMin * sqrt(time) - T), 2);
mySum2 += pow((Tsunset - myMax * sqrt(time) - T), 2);
time++;
}
statistics::weightedMultiRegressionLinearWithStats(sunsetData, outCoeff, weights, &regrConst, regrCoeff, false, true, &R2, &stdError);

myFirstError = sqrt(mySum1);
mySecondError = sqrt(mySum2);

if (myFirstError < mySecondError)
myMax = (myMin + myMax) / 2;
else
myMin = (myMin + myMax) / 2;

} while (fabs(myFirstError - mySecondError) > tolerance && iteration <= nrIterations);

coeff.push_back((myMin + myMax) / 2);
return true;
}

return false;
}
25 changes: 16 additions & 9 deletions bin/frostForecast/frost.h
Original file line number Diff line number Diff line change
Expand Up @@ -39,9 +39,8 @@ class Frost
int createCsvFile(QString id);
QList<QString> getIdList() const;

// calibration functions
bool getRadiativeCoolingHistory(QString id, float thresholdTmin, float thresholdTRange, int monthIni, int monthFin, int timeZone, std::vector<std::vector<float> > &outData, std::vector<std::vector<float> > &sunsetData);
void fitCoolingCoefficient(std::vector<std::vector<float> > &hourly_series, int nrIterations, float tolerance, float minValue, float maxValue, std::vector<float> &coeff);
bool getRadiativeCoolingHistory(QString id, std::vector<std::vector<float>>& outData, std::vector <std::vector <float>>& sunsetData);
bool calibrateModel(QString id);

private:
gis::Crit3DGisSettings gisSettings;
Expand All @@ -58,16 +57,24 @@ class Frost
QList<Crit3DMeteoPoint> meteoPointsList;
QList<QString> idList;
QList<QString> varList;
QList<QString> intercept;
QList<QString> parTss;
QList<QString> parRHss;
QList<QString> SE_intercept;
QList<QString> SE_parTss;
QList<QString> SE_parRHss;

std::vector <float> intercept;
std::vector <float> parTss;
std::vector <float> parRHss;
std::vector <float> SE_intercept;
std::vector <float> SE_parTss;
std::vector <float> SE_parRHss;

int indexSunSet;
int indexSunRise;
QString errorString;

//calibration parameters
int monthIni;
int monthFin;
float thresholdTmin;
float thresholdTrange;

QList<float> myForecast;
QList<float> myForecastMin;
QList<float> myForecastMax;
Expand Down
Loading

0 comments on commit 903fd16

Please sign in to comment.