Performance improvement, calibration tool refactoring, extended regression tests, Julia 1.12.6

DDDEcco/RunDDDv2_simplified.ipynb

@@ -0,0 +1,291 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "\"\"\"\n",
+    "Runs DDD model either in single run model or for calibration\n",
+    "The model itself is called as a function which calls on several functions\n",
+    "This one is for ECCO use not MRT and not operative or HSO\n",
+    "\"\"\"\n",
+    "\n",
+    "using DataFrames\n",
+    "using Plots\n",
+    "using CSV\n",
+    "using BlackBoxOptim\n",
+    "include(joinpath(dirname(@__DIR__), \"DDDFunctions\", \"DDDAllTerrain22012024.jl\"))\n",
+    "\n",
+    "catchment = \"55.4\"  # stationnumber\n",
+    "navn = \"Royknes\"\n",
+    "\n",
+    "TR = \"3h\"         # this is just a marker for naming files, does NOT set the temporal resolution\n",
+    "\n",
+    "ptqfile = string(\"\\\\\\\\nve.no\\\\fil\\\\h\\\\HB\\\\HB-modellering\\\\DDD3h\\\\inndataV2\\\\\",catchment,\"\\\\\",\n",
+    "    catchment,\"_\",TR,\"_ptq_DDDv2_kal.csv\") \n",
+    "#ptqfile = string(\"\\\\\\\\nve.no\\\\fil\\\\h\\\\HM\\\\Interne Prosjekter\\\\UrbanDesignFlom\\\\utdata\\\\\",navn,\"\\\\SimTimeseries\\\\FraHSO\\\\\",\n",
+    "#    \"Royknes_noDatetime_7a565a2c-f20d-459c-ae55-a631d4a6dc33_3H_ptq.csv\")\n",
+    "#ptqfile = joinpath(\"/hdata/h/HB/HB-modellering/DDD3h/inndataV2\", catchment, \"$(catchment)_$(TR)_ptq_DDDv2_kal.csv\") # Linux\n",
+    "\n",
+    "r2fil = string(\"\\\\\\\\nve.no\\\\fil\\\\h\\\\HM\\\\Interne Prosjekter\\\\UrbanDesignFlom\\\\utdata\\\\\",navn,\"\\\\r2_\",\n",
+    "    catchment,\"_\",TR,\"_test.csv\")\n",
+    "#r2fil = \"/hdata/hmdata01/DDD_calibration/test_notebook/r2_$(catchment)_$(TR).csv\" # Linux\n",
+    "\n",
+    "utfile = string(\"\\\\\\\\nve.no\\\\fil\\\\h\\\\HM\\\\Interne Prosjekter\\\\UrbanDesignFlom\\\\utdata\\\\\",navn,\"\\\\simres_\",\n",
+    "    catchment,\"_\",TR,\"_DDDv2.csv\")\n",
+    "#utfile = \"/hdata/hmdata01/DDD_calibration/test_notebook/simres_$(catchment)_$(TR).csv\"\n",
+    "\n",
+    "paramfile = string(\"\\\\\\\\nve.no\\\\fil\\\\h\\\\HM\\\\Interne Prosjekter\\\\UrbanDesignFlom\\\\DDDurbanParameters\\\\\",catchment,\"\\\\ParDDDv2_\",\n",
+    "    catchment,\"_3h.csv\")\n",
+    "#paramfile = \"/hdata/h/HM/HydrologiskeModeller/DDDJulia/KalibreringDDDv2_2025/Parameters/BestPrimo2025/Best_par_$(catchment)_$(TR)_DDDv2.csv\" # Linux\n",
+    "\n",
+    "spinup = (31*4) #days used to spin up the model. \n",
+    "\n",
+    "prm = CSV.read(paramfile,DataFrame,header=[\"Name\", \"val\"], delim=';')\n",
+    "#prm = CSV.read(paramfile,header=[\"Name\", \"val\"], delim=';')\n",
+    "\n",
+    "#            u,          pro           TX,         Pkorr        skorr,     GscInt      OVP          \n",
+    "#         OVIP       Lv            rv        \n",
+    "tprm = [prm.val[20], prm.val[21], prm.val[22], prm.val[18], prm.val[19],prm.val[33], prm.val[34], \n",
+    "    prm.val[35],prm.val[36],prm.val[37]]\n",
+    "println(tprm)\n",
+    "\n",
+    "\n",
+    "Gshape, Gscale = Big2SmallLambda(prm.val[32], prm.val[33]) # Coverting integrated celerity to layers takes too long in calibration: preprocessing\n",
+    "Gpar = [Gshape, Gscale]\n",
+    "\n",
+    "println(prm.val[32],\" \", prm.val[33])\n",
+    "\n",
+    "startsim = 1 \n",
+    "kal = 0\n",
+    "modstate = 0\n",
+    "savestate = 0\n",
+    "\n",
+    "t1= time_ns()\n",
+    "\n",
+    "function calib_wrapper_model(Gpar,startsim, tprm, prm, ptqfile, utfile, r2fil, modstate, savestate, kal, spinup)\n",
+    " qobs, qberegn, KGE, NSE, bias = DDDAllTerrain(Gpar,startsim, tprm, prm, ptqfile, utfile, r2fil, modstate, savestate,\n",
+    "        kal, spinup)  \n",
+    " return qobs,qberegn, KGE,NSE,bias \n",
+    "end\n",
+    "\n",
+    "function calib_single_wsh(Gpar,startsim, tprm, prm, ptqfile, utfile, r2fil, modstate, savestate, kal, spinup)\n",
+    " qobs, qberegn, KGE, NSE, bias = DDDAllTerrain(Gpar,startsim, tprm, prm, ptqfile, utfile, r2fil, modstate, savestate,\n",
+    "        kal, spinup)    \n",
+    " return (1.0 - KGE)\n",
+    "end\n",
+    "\n",
+    "if(kal == 0)\n",
+    "    qobs,qberegn,KGE,NSE, bias = calib_wrapper_model(Gpar,startsim, tprm, prm, ptqfile, utfile, r2fil,\n",
+    "        modstate, savestate,kal, spinup) # a single run \n",
+    "    \n",
+    "    println(catchment)\n",
+    "    println(\"KGE=\",round(KGE,digits=3))\n",
+    "    println(\"NSE=\",round(NSE,digits=3))\n",
+    "    println(\"bias=\",round(bias,digits=3))\n",
+    "end\n",
+    "\n",
+    "if(kal == 1) # calibrate\n",
+    "    #                   u,        pro,         TX,        Pkorr,    skorr,          GscInt,         OVP     OVIP \n",
+    "    param_range = [(1.0,3.0), (0.05,0.05), (-0.5, 0.5), (0.5, 2.0), (0.5,2.0), (0.065,0.075), (tprm[7],tprm[7]),\n",
+    "        (tprm[8],tprm[8]), (tprm[9],tprm[9]),(tprm[10],tprm[10])] # \n",
+    "    \n",
+    "    println(param_range)\n",
+    "    calib_single_wsh_tmp(param) = calib_single_wsh(Gpar,startsim, param, prm, ptqfile, utfile, r2fil,\n",
+    "                                           modstate, savestate, kal, spinup)\n",
+    "    res = bboptimize(calib_single_wsh_tmp; SearchRange = param_range, MaxSteps = 1000, TraceMode = :verbose)\n",
+    "    param_hydro = best_candidate(res)\n",
+    "    println(param_hydro)\n",
+    "end\n",
+    "\n",
+    "t2 = time_ns()\n",
+    "println(\"Pkorr=\", round(tprm[4],digits=3))\n",
+    "println(\"Time elapsed[s]= \",(t2-t1)/1.0e9)\n",
+    "\n",
+    "if(kal==0)\n",
+    " plot(qobs[6500:8000], color=\"black\",label = \"Observed\",lw =1)\n",
+    " plot!(qberegn[6500:8000],col= \"blue\",label = \"Simulated\", lw = 2)\n",
+    "end \n",
+    " "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#function CeleritySubSurface(NoL, Gshape, Gscale, midDL, Timeresinsec) \n",
+    "\n",
+    " using Distributions  \n",
+    " using Plots\n",
+    "  NoL = 5\n",
+    "  Gshape = 4.0\n",
+    "  meanVel =0.00012\n",
+    "  midDL = 129.84 \n",
+    "  Timeresinsec = 10800\n",
+    "  a = 3 # rapid layer multiplier for DF aquifers\n",
+    "  Gscale = meanVel*Timeresinsec/(Gshape*midDL) \n",
+    "\n",
+    "k = zeros(Float64,NoL)                   #celerity of subsurface (and overland) flow\n",
+    "probvec = zeros(Float64,NoL)            #all leves and overland flow level\n",
+    "\n",
+    "dp = 1/(NoL-1)                              #Overland flow level (nol=1], fixed celerity and extremely high capacity(2000 mm)\n",
+    "\n",
+    "for i in reverse(1:(NoL-1))\n",
+    "    probvec[i+1] = i*dp - dp/2                #celerities are estimated at center of level, hence dp/2\n",
+    "end\n",
+    "\n",
+    "probvec = 1 .- probvec\n",
+    "probvec[1] = 0.99                            #Quantile in celerity distribution for overland flow fixed at 0.99\n",
+    "\n",
+    "g = Gamma(Gshape, Gscale)    \n",
+    "k[1:NoL] = quantile.(g,probvec[1:NoL])*midDL/Timeresinsec\n",
+    "##Dupuit Forschheimer \n",
+    "k[NoL] = 4*meanVel/(3+a) \n",
+    "k[4] = 4*meanVel/(3+a) \n",
+    "k[3] = 4*meanVel/(3+a) \n",
+    "k[2] = k[3]*a\n",
+    "k[1] = 0.0005\n",
+    "\n",
+    "#plot(k[1:NoL], color=\"black\",label = string(\"Shape=\",Gshape),lw =1)\n",
+    "#plot!(k[1:NoL], color=\"green\",label = string(\"Shape=\",Gshape),lw =1)\n",
+    "plot!(k[1:NoL], color=\"red\",label = string(\"DF\"),lw =3)\n",
+    "# plot!(qberegn[5000:7000],col= \"blue\",label = \"Simulated\", lw = 2)\n",
+    "\n",
+    "#return k                        #Celerities [m/s] \n",
+    "#end"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    " using Distributions\n",
+    " include(joinpath(dirname(@__DIR__), \"DDDFunctions\", \"SingleUH.jl\"))\n",
+    "#Røykenes\n",
+    "  GshInt = 0.8\n",
+    "  GscInt = 0.055 \n",
+    "  midDL = 156.07\n",
+    "  maxDl = 928\n",
+    "  MAD = 5.24\n",
+    "  Timeresinsec = 10800\n",
+    "  NoL = 2\n",
+    "  area2 = 50090000 \n",
+    "  gtcel = 0.99\n",
+    " \n",
+    "mLam = GshInt*GscInt\n",
+    "varLam = GshInt*(GscInt)^2                           #Yevjevich p.145\n",
+    "meanIntk = mLam*midDL/Timeresinsec                   #mean celerity estimated through Integrated Celerity\n",
+    "antBox = Int(trunc(maxDl/(meanIntk*Timeresinsec)))+1 #Temporal length UH_MAD\n",
+    "UH_MAD = zeros(Float64,antBox)\n",
+    "sRes = zeros(Float64,antBox) # saturation sum\n",
+    "\n",
+    "#Unit hydrograph for MAD\n",
+    "UH_MAD = SingleUH(meanIntk,Timeresinsec, midDL, maxDl, 0)\n",
+    "\n",
+    "StSt = (1000*MAD*Timeresinsec)/(area2)         # Steady state Input eq. output in mm\n",
+    "sRes[1] = 0\n",
+    "sRes[2:antBox] .= StSt.*UH_MAD[2:antBox]\n",
+    "\n",
+    "for i in 3: antBox\n",
+    "  sRes[i:antBox] .= sRes[i:antBox] + StSt.*UH_MAD[i:antBox]\n",
+    "end\n",
+    "\n",
+    "mRes = sum(sRes)\n",
+    "Fact = mLam/mRes\n",
+    "stdRes = (varLam/Fact^2)^0.5                   # see Haan p.51\n",
+    "\n",
+    "GshRes =  1.83 #mRes^2/stdRes^2\n",
+    "GscRes =  11.43 #stdRes^2/mRes\n",
+    "\n",
+    "MLev = [1/(NoL-1):1/(NoL-1):1.0;]              # (sequence)Quantiles  to calculate reservoir levels [0.1:0.1:0.9;]\n",
+    "\n",
+    "MLev[NoL-1] = gtcel                            # quantile for start overland flow\n",
+    "Res_prob = zeros(Float64,(NoL-1))\n",
+    "Magkap = zeros(Float64,NoL)\n",
+    "g = Gamma(GshRes,GscRes) \n",
+    "#calculates the reservoir levels associated with quantiles. Mean is GshRes*GscRes\n",
+    "Res_prob .= quantile.(g,MLev)\n",
+    "\n",
+    "#Capasity of Layers\n",
+    "ssRes1 = zeros(Float64,NoL)\n",
+    "ssRes1[1] = 2000                               # capacity of overland flow level\n",
+    "\n",
+    "for i in 2:(NoL-1)\n",
+    "  ssRes1[i] = Res_prob[NoL-i+1]-Res_prob[(NoL-i)]\n",
+    "end\n",
+    "\n",
+    "ssRes1[NoL] = Res_prob[1]                     # capasity for the first slowest level         \n",
+    "                       \n",
+    "Magkap = ssRes1                                 # capasity for Layers\n",
+    "M = Res_prob[(NoL-1)]                           # Total groundwater reservoir\n",
+    "println(GshRes) \n",
+    "println(GscRes)\n",
+    "println(\"Magkap fra Subrutine \",Magkap)\n",
+    "println(\"M fra Subrutine \", M)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "NoL = 2\n",
+    "GshRes =  1.83 #mRes^2/stdRes^2\n",
+    "GscRes =  11.43 #stdRes^2/mRes\n",
+    "gtcel = 0.99\n",
+    "\n",
+    "MLev = [1/(NoL-1):1/(NoL-1):1.0;]              # (sequence)Quantiles  to calculate reservoir levels [0.1:0.1:0.9;]\n",
+    "\n",
+    "MLev[NoL-1] = gtcel                            # quantile for start overland flow\n",
+    "Res_prob = zeros(Float64,(NoL-1))\n",
+    "Magkap = zeros(Float64,NoL)\n",
+    "g = Gamma(GshRes,GscRes) \n",
+    "#calculates the reservoir levels associated with quantiles. Mean is GshRes*GscRes\n",
+    "Res_prob .= quantile.(g,MLev)\n",
+    "\n",
+    "#Capasity of Layers\n",
+    "ssRes1 = zeros(Float64,NoL)\n",
+    "ssRes1[1] = 2000                               # capacity of overland flow level\n",
+    "\n",
+    "for i in 2:(NoL-1)\n",
+    "  ssRes1[i] = Res_prob[NoL-i+1]-Res_prob[(NoL-i)]\n",
+    "end\n",
+    "\n",
+    "ssRes1[NoL] = Res_prob[1]                     # capasity for the first slowest level         \n",
+    "                       \n",
+    "Magkap = ssRes1                                 # capasity for Layers\n",
+    "M = Res_prob[(NoL-1)]                           # Total groundwater reservoir\n",
+    "println(GshRes) \n",
+    "println(GscRes)\n",
+    "println(\"Magkap fra Subrutine \",Magkap)\n",
+    "println(\"M fra Subrutine \", M)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Julia DDD 1.12",
+   "language": "julia",
+   "name": "julia-ddd-1.12"
+  },
+  "language_info": {
+   "file_extension": ".jl",
+   "mimetype": "application/julia",
+   "name": "julia",
+   "version": "1.12.7"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

DDDFunctions/AlbedoUEB_debug04072022.jl

@@ -33,7 +33,7 @@ r2 = min(r1^10,1)
 r3 = 0.03
 d_tau =  ((r1+r2+r3)/tau0)*Timeresinsec
 
-taux = taux + d_tau
+taux += d_tau
 F_age = taux/(1+taux)
 
 alfa_vd = (1-C_v*F_age)*alfa_v0

DDDFunctions/Big2SmallLambda.jl

@@ -136,11 +136,9 @@ for lag in 1:10 #Loop for 10 Layers
 end #loop for levels
 
 #Estimate parameters for trlam (lambda)
-dx = trlam
-dy = arg 
-     @.modelgamma(x,p) = quantile(Gamma(p[1],p[2]),x)
-     p0=[GshInt*0.5, 2*GscInt]
-     gammafit = curve_fit(modelgamma,dy, dx , p0)
+modelgamma(x, p) = quantile.(Gamma(p[1], p[2]), x)
+gammafit = curve_fit(modelgamma, arg, trlam, [GshInt*0.5, 2*GscInt],
+                     autodiff=AutoFiniteDiff(fdjtype=Val(:central)))
 
 #using Plots
 #plot(quantile.(Gamma(coef(gammafit)[1],coef(gammafit)[2]),dy),dy)

DDDFunctions/BogLayerUpdate.jl

@@ -8,13 +8,12 @@
 #--------------------------------------------------------------------------
 
 function BogLayerUpdate!(outbog, BogLayers, UHBog, nodaysvector)
-    qlayer = outbog * UHBog    #finds response  in mm!!!!  for bog, a vector     
-    for t in 2:nodaysvector
-        BogLayers[t-1] = BogLayers[t] + qlayer[t] # shifts the level of the matrix one timestep ahead
-    end
     if nodaysvector == 1
         BogLayers .= 0.
-#        BogLayers .= qlayer
+    else
+        @inbounds @simd for t in 2:nodaysvector
+            BogLayers[t-1] = BogLayers[t] + outbog * UHBog[t] # finds response (mm) and shifts the matrix one timestep ahead
+        end
     end
 end            
 

DDDFunctions/CloudCoverGlac_debug04072022.jl

@@ -1,5 +1,3 @@
-
-
 function CloudCover(regn,Timeresinsec)
 
  grense = 0.98 
@@ -10,16 +8,12 @@ function CloudCover(regn,Timeresinsec)
  if regn > pgrense # 1 for 1 time funker bra, 3 for 3 timer?
    Cl = 1.0
    Wa = 1.0
- end      
-
- if (regn > pgrense2 && regn <=pgrense)
-   Cl = rand(Uniform(grense,0.99))
-   Wa = (Cl+1.667)/2.667
- end
-
- if (regn <= pgrense2)
-   Cl = rand(Uniform(0.01,grense)) 
-   Wa = (Cl+1.667)/2.667
+ elseif regn > pgrense2
+   Cl = grense + rand() * (0.99 - grense)
+   Wa = min((Cl+1.667)/2.667, 1.0)
+ else
+   Cl = 0.01 + rand() * (grense - 0.01)
+   Wa = min((Cl+1.667)/2.667, 1.0)
  end       
 
  #if (Cl >= grense)
@@ -34,9 +28,6 @@ function CloudCover(regn,Timeresinsec)
  #   Wa = (Cl+1.167)/2.667  # parametrization of Relative Humidity from Herrero and Polo(HESS, 2012)
  #end
 
- if (Wa > 1.0)
-    Wa = 1.0# 0.8349  #1.0
- end
  #println(Cl, " ", Wa)  
 
  return Cl, Wa