From e2c8cba6786fa56b2dc0d1bcdcf59a76ac1e67b4 Mon Sep 17 00:00:00 2001
From: Tim Peterson <tim.peterson@monash.edu>
Date: Fri, 12 May 2023 13:13:43 +1000
Subject: [PATCH] Calibration user exp. bug fixed. Doc updated

krige.head.calib outputs changed to a list variable defining the entire mapping method and inputs. Inputting this to krige.head will map the head without any of the previously required handling of calibrated parameter values.  See ?krige.head.calib for details.
---
 .Rproj.user/19E94E50/pcs/source-pane.pper     |   2 +-
 .../19E94E50/pcs/windowlayoutstate.pper       |   4 +-
 .Rproj.user/19E94E50/persistent-state         |   2 +-
 .Rproj.user/19E94E50/sources/per/t/D19F9BF0   |  26 --
 .Rproj.user/19E94E50/sources/prop/0072A15F    |   4 +-
 .Rproj.user/19E94E50/sources/prop/0243BC54    |   6 +
 .Rproj.user/19E94E50/sources/prop/06EF12A0    |   2 +-
 .Rproj.user/19E94E50/sources/prop/57B75005    |   6 +
 .Rproj.user/19E94E50/sources/prop/AF65AB58    |   2 +-
 .Rproj.user/19E94E50/sources/prop/BD952B7F    |   2 +-
 .Rproj.user/19E94E50/sources/prop/D3423677    |   4 +-
 .Rproj.user/19E94E50/sources/prop/INDEX       |   2 +
 .../session-4dda64d8/0FF087FE-contents        |   0
 .../session-4dda64d8/221ACDA6-contents        |   0
 .../{per/t => session-4dda64d8}/4C0D952C      |  10 +-
 .../t => session-4dda64d8}/4C0D952C-contents  |  56 ++--
 .../5B688876-contents}                        | 125 +-------
 .../{per/t => session-4dda64d8}/5C989810      |   4 +-
 .../t => session-4dda64d8}/5C989810-contents  |   0
 .../session-4dda64d8/648DFF61-contents        |   0
 .../session-4dda64d8/65610FB2-contents        |   0
 .../t/F6CA6EB6 => session-4dda64d8/842B6FEA}  |  16 +-
 .../session-4dda64d8/842B6FEA-contents        | 275 ++++++++++++++++++
 .../t => session-4dda64d8}/8974836D-contents  |   0
 .../session-4dda64d8/8C88B17A-contents        | 183 ++++++++++++
 .../session-4dda64d8/9E018CFD-contents        |   0
 .../session-4dda64d8/A6CEE51C-contents        |   0
 .../session-4dda64d8/A9B84F7A-contents        |   0
 .../{per/u => session-4dda64d8}/B06A395E      |   8 +-
 .../u => session-4dda64d8}/B06A395E-contents  |   5 +-
 .../session-4dda64d8/B094F956-contents        |   0
 .../{per/t => session-4dda64d8}/B29BDB7C      |   8 +-
 .../t => session-4dda64d8}/B29BDB7C-contents  |   2 +-
 .../{per/t => session-4dda64d8}/D0D0B5EC      |  12 +-
 .../t => session-4dda64d8}/D0D0B5EC-contents  | 138 ++++++---
 .../t => session-4dda64d8}/D19F9BF0-contents  |   0
 .../{per/t => session-4dda64d8}/F1F4D846      |   2 +-
 .../t => session-4dda64d8}/F1F4D846-contents  |   0
 .../t/8974836D => session-4dda64d8/F6CA6EB6}  |  18 +-
 .../session-4dda64d8/F6CA6EB6-contents        | 129 ++++++++
 .../sources/session-4dda64d8/lock_file        |   0
 .Rproj.user/shared/notebooks/paths            |   5 +
 DESCRIPTION                                   |   2 +-
 R/get.objFunc.R                               | 125 +-------
 R/get.variogramModel.R                        | 129 ++++++++
 R/krige.head.R                                |  56 ++--
 R/krige.head.calib.R                          | 138 ++++++---
 man/krige.head.Rd                             |  28 +-
 man/krige.head.calib.Rd                       |  59 ++--
 49 files changed, 1102 insertions(+), 493 deletions(-)
 delete mode 100644 .Rproj.user/19E94E50/sources/per/t/D19F9BF0
 create mode 100644 .Rproj.user/19E94E50/sources/prop/0243BC54
 create mode 100644 .Rproj.user/19E94E50/sources/prop/57B75005
 create mode 100644 .Rproj.user/19E94E50/sources/session-4dda64d8/0FF087FE-contents
 create mode 100644 .Rproj.user/19E94E50/sources/session-4dda64d8/221ACDA6-contents
 rename .Rproj.user/19E94E50/sources/{per/t => session-4dda64d8}/4C0D952C (75%)
 rename .Rproj.user/19E94E50/sources/{per/t => session-4dda64d8}/4C0D952C-contents (96%)
 rename .Rproj.user/19E94E50/sources/{per/t/F6CA6EB6-contents => session-4dda64d8/5B688876-contents} (75%)
 rename .Rproj.user/19E94E50/sources/{per/t => session-4dda64d8}/5C989810 (92%)
 rename .Rproj.user/19E94E50/sources/{per/t => session-4dda64d8}/5C989810-contents (100%)
 create mode 100644 .Rproj.user/19E94E50/sources/session-4dda64d8/648DFF61-contents
 create mode 100644 .Rproj.user/19E94E50/sources/session-4dda64d8/65610FB2-contents
 rename .Rproj.user/19E94E50/sources/{per/t/F6CA6EB6 => session-4dda64d8/842B6FEA} (64%)
 create mode 100644 .Rproj.user/19E94E50/sources/session-4dda64d8/842B6FEA-contents
 rename .Rproj.user/19E94E50/sources/{per/t => session-4dda64d8}/8974836D-contents (100%)
 create mode 100644 .Rproj.user/19E94E50/sources/session-4dda64d8/8C88B17A-contents
 create mode 100644 .Rproj.user/19E94E50/sources/session-4dda64d8/9E018CFD-contents
 create mode 100644 .Rproj.user/19E94E50/sources/session-4dda64d8/A6CEE51C-contents
 create mode 100644 .Rproj.user/19E94E50/sources/session-4dda64d8/A9B84F7A-contents
 rename .Rproj.user/19E94E50/sources/{per/u => session-4dda64d8}/B06A395E (80%)
 rename .Rproj.user/19E94E50/sources/{per/u => session-4dda64d8}/B06A395E-contents (83%)
 create mode 100644 .Rproj.user/19E94E50/sources/session-4dda64d8/B094F956-contents
 rename .Rproj.user/19E94E50/sources/{per/t => session-4dda64d8}/B29BDB7C (78%)
 rename .Rproj.user/19E94E50/sources/{per/t => session-4dda64d8}/B29BDB7C-contents (95%)
 rename .Rproj.user/19E94E50/sources/{per/t => session-4dda64d8}/D0D0B5EC (70%)
 rename .Rproj.user/19E94E50/sources/{per/t => session-4dda64d8}/D0D0B5EC-contents (90%)
 rename .Rproj.user/19E94E50/sources/{per/t => session-4dda64d8}/D19F9BF0-contents (100%)
 rename .Rproj.user/19E94E50/sources/{per/t => session-4dda64d8}/F1F4D846 (96%)
 rename .Rproj.user/19E94E50/sources/{per/t => session-4dda64d8}/F1F4D846-contents (100%)
 rename .Rproj.user/19E94E50/sources/{per/t/8974836D => session-4dda64d8/F6CA6EB6} (51%)
 create mode 100644 .Rproj.user/19E94E50/sources/session-4dda64d8/F6CA6EB6-contents
 create mode 100644 .Rproj.user/19E94E50/sources/session-4dda64d8/lock_file
 create mode 100644 R/get.variogramModel.R

diff --git a/.Rproj.user/19E94E50/pcs/source-pane.pper b/.Rproj.user/19E94E50/pcs/source-pane.pper
index ddca97d..902cc6f 100644
--- a/.Rproj.user/19E94E50/pcs/source-pane.pper
+++ b/.Rproj.user/19E94E50/pcs/source-pane.pper
@@ -1,3 +1,3 @@
 {
-    "activeTab": 2
+    "activeTab": 0
 }
\ No newline at end of file
diff --git a/.Rproj.user/19E94E50/pcs/windowlayoutstate.pper b/.Rproj.user/19E94E50/pcs/windowlayoutstate.pper
index ecf223a..a74678f 100644
--- a/.Rproj.user/19E94E50/pcs/windowlayoutstate.pper
+++ b/.Rproj.user/19E94E50/pcs/windowlayoutstate.pper
@@ -1,12 +1,12 @@
 {
     "left": {
-        "splitterpos": 313,
+        "splitterpos": 320,
         "topwindowstate": "NORMAL",
         "panelheight": 1079,
         "windowheight": 1117
     },
     "right": {
-        "splitterpos": 490,
+        "splitterpos": 501,
         "topwindowstate": "NORMAL",
         "panelheight": 1079,
         "windowheight": 1117
diff --git a/.Rproj.user/19E94E50/persistent-state b/.Rproj.user/19E94E50/persistent-state
index 107899c..74a7c6b 100644
--- a/.Rproj.user/19E94E50/persistent-state
+++ b/.Rproj.user/19E94E50/persistent-state
@@ -1,6 +1,6 @@
 build-last-errors="[]"
 build-last-errors-base-dir="C:/Users/tpet0008/Documents/HydroMap/"
-build-last-outputs="[{\"type\":0,\"output\":\"==> Rcmd.exe INSTALL --preclean --no-multiarch --with-keep.source HydroMap\\n\\n\"},{\"type\":1,\"output\":\"* installing to library 'C:/Program Files/R/R-4.0.3/library'\\r\\n\"},{\"type\":1,\"output\":\"* installing *source* package 'HydroMap' ...\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** using staged installation\\r\\n\"},{\"type\":1,\"output\":\"** R\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** data\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"*** moving datasets to lazyload DB\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** byte-compile and prepare package for lazy loading\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** help\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"  converting help for package 'HydroMap'\\r\\n\"},{\"type\":1,\"output\":\"    import.DEM                              html  \"},{\"type\":1,\"output\":\"*** installing help indices\\r\\n\"},{\"type\":1,\"output\":\"    finding HTML links ... done\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"\\r\\n\"},{\"type\":1,\"output\":\"    import.pointData                        html  \\r\\n\"},{\"type\":1,\"output\":\"    krige.head                              html  \"},{\"type\":1,\"output\":\"\\r\\n\"},{\"type\":1,\"output\":\"    krige.head.calib                        html  \"},{\"type\":1,\"output\":\"\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** building package indices\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** testing if installed package can be loaded from temporary location\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** testing if installed package can be loaded from final location\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** testing if installed package keeps a record of temporary installation path\\r\\n\"},{\"type\":1,\"output\":\"* DONE (HydroMap)\\r\\n\"},{\"type\":1,\"output\":\"\"}]"
+build-last-outputs="[{\"type\":0,\"output\":\"==> Rcmd.exe INSTALL --preclean --no-multiarch --with-keep.source HydroMap\\n\\n\"},{\"type\":1,\"output\":\"* installing to library 'C:/Program Files/R/R-4.0.3/library'\\r\\n\"},{\"type\":1,\"output\":\"* installing *source* package 'HydroMap' ...\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** using staged installation\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** R\\r\\n\"},{\"type\":1,\"output\":\"** data\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"*** moving datasets to lazyload DB\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** byte-compile and prepare package for lazy loading\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** help\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"  converting help for package 'HydroMap'\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"*** installing help indices\\r\\n\"},{\"type\":1,\"output\":\"    finding HTML links ...\"},{\"type\":1,\"output\":\"    import.DEM                              html  \\r\\n\"},{\"type\":1,\"output\":\"    import.pointData                        html  \\r\\n\"},{\"type\":1,\"output\":\"    krige.head                              html  \"},{\"type\":1,\"output\":\" done\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"\\r\\n\"},{\"type\":1,\"output\":\"    krige.head.calib                        html  \"},{\"type\":1,\"output\":\"\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** building package indices\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** testing if installed package can be loaded from temporary location\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** testing if installed package can be loaded from final location\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** testing if installed package keeps a record of temporary installation path\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"* DONE (HydroMap)\\r\\n\"},{\"type\":1,\"output\":\"\"}]"
 compile_pdf_state="{\"tab_visible\":false,\"running\":false,\"target_file\":\"\",\"output\":\"\",\"errors\":[]}"
 files.monitored-path=""
 find-in-files-state="{\"handle\":\"\",\"input\":\"\",\"path\":\"\",\"regex\":false,\"ignoreCase\":false,\"results\":{\"file\":[],\"line\":[],\"lineValue\":[],\"matchOn\":[],\"matchOff\":[],\"replaceMatchOn\":[],\"replaceMatchOff\":[]},\"running\":false,\"replace\":false,\"preview\":false,\"gitFlag\":false,\"replacePattern\":\"\"}"
diff --git a/.Rproj.user/19E94E50/sources/per/t/D19F9BF0 b/.Rproj.user/19E94E50/sources/per/t/D19F9BF0
deleted file mode 100644
index e9cd595..0000000
--- a/.Rproj.user/19E94E50/sources/per/t/D19F9BF0
+++ /dev/null
@@ -1,26 +0,0 @@
-{
-    "id": "D19F9BF0",
-    "path": "C:/Users/tpet0008/Documents/HydroMap/R/get.allData.R",
-    "project_path": "R/get.allData.R",
-    "type": "r_source",
-    "hash": "625844390",
-    "contents": "",
-    "dirty": false,
-    "created": 1683615909780.0,
-    "source_on_save": false,
-    "relative_order": 7,
-    "properties": {
-        "source_window_id": "",
-        "Source": "Source",
-        "cursorPosition": "99,1",
-        "scrollLine": "0"
-    },
-    "folds": "",
-    "lastKnownWriteTime": 1683421624,
-    "encoding": "UTF-8",
-    "collab_server": "",
-    "source_window": "",
-    "last_content_update": 1683421624,
-    "read_only": false,
-    "read_only_alternatives": []
-}
\ No newline at end of file
diff --git a/.Rproj.user/19E94E50/sources/prop/0072A15F b/.Rproj.user/19E94E50/sources/prop/0072A15F
index 1fa2bfe..286e64b 100644
--- a/.Rproj.user/19E94E50/sources/prop/0072A15F
+++ b/.Rproj.user/19E94E50/sources/prop/0072A15F
@@ -1,6 +1,6 @@
 {
     "source_window_id": "",
     "Source": "Source",
-    "cursorPosition": "202,24",
-    "scrollLine": "184"
+    "cursorPosition": "26,60",
+    "scrollLine": "11"
 }
\ No newline at end of file
diff --git a/.Rproj.user/19E94E50/sources/prop/0243BC54 b/.Rproj.user/19E94E50/sources/prop/0243BC54
new file mode 100644
index 0000000..112488e
--- /dev/null
+++ b/.Rproj.user/19E94E50/sources/prop/0243BC54
@@ -0,0 +1,6 @@
+{
+    "source_window_id": "",
+    "Source": "Source",
+    "cursorPosition": "26,2",
+    "scrollLine": "0"
+}
\ No newline at end of file
diff --git a/.Rproj.user/19E94E50/sources/prop/06EF12A0 b/.Rproj.user/19E94E50/sources/prop/06EF12A0
index 4c571e9..f0add8c 100644
--- a/.Rproj.user/19E94E50/sources/prop/06EF12A0
+++ b/.Rproj.user/19E94E50/sources/prop/06EF12A0
@@ -1,6 +1,6 @@
 {
     "source_window_id": "",
     "Source": "Source",
-    "cursorPosition": "4,4",
+    "cursorPosition": "6,90",
     "scrollLine": "0"
 }
\ No newline at end of file
diff --git a/.Rproj.user/19E94E50/sources/prop/57B75005 b/.Rproj.user/19E94E50/sources/prop/57B75005
new file mode 100644
index 0000000..b67e905
--- /dev/null
+++ b/.Rproj.user/19E94E50/sources/prop/57B75005
@@ -0,0 +1,6 @@
+{
+    "source_window_id": "",
+    "Source": "Source",
+    "cursorPosition": "2,2",
+    "scrollLine": "0"
+}
\ No newline at end of file
diff --git a/.Rproj.user/19E94E50/sources/prop/AF65AB58 b/.Rproj.user/19E94E50/sources/prop/AF65AB58
index eda262d..c0e00d1 100644
--- a/.Rproj.user/19E94E50/sources/prop/AF65AB58
+++ b/.Rproj.user/19E94E50/sources/prop/AF65AB58
@@ -2,5 +2,5 @@
     "source_window_id": "",
     "Source": "Source",
     "cursorPosition": "109,71",
-    "scrollLine": "88"
+    "scrollLine": "0"
 }
\ No newline at end of file
diff --git a/.Rproj.user/19E94E50/sources/prop/BD952B7F b/.Rproj.user/19E94E50/sources/prop/BD952B7F
index 5f0afbb..bfed50f 100644
--- a/.Rproj.user/19E94E50/sources/prop/BD952B7F
+++ b/.Rproj.user/19E94E50/sources/prop/BD952B7F
@@ -1,5 +1,5 @@
 {
-    "cursorPosition": "2,12",
+    "cursorPosition": "2,16",
     "scrollLine": "0",
     "Source": "Source"
 }
\ No newline at end of file
diff --git a/.Rproj.user/19E94E50/sources/prop/D3423677 b/.Rproj.user/19E94E50/sources/prop/D3423677
index 0606a08..cc29640 100644
--- a/.Rproj.user/19E94E50/sources/prop/D3423677
+++ b/.Rproj.user/19E94E50/sources/prop/D3423677
@@ -1,6 +1,6 @@
 {
     "source_window_id": "",
     "Source": "Source",
-    "cursorPosition": "292,2",
-    "scrollLine": "0"
+    "cursorPosition": "93,8",
+    "scrollLine": "69"
 }
\ No newline at end of file
diff --git a/.Rproj.user/19E94E50/sources/prop/INDEX b/.Rproj.user/19E94E50/sources/prop/INDEX
index e0c6746..25b5cf8 100644
--- a/.Rproj.user/19E94E50/sources/prop/INDEX
+++ b/.Rproj.user/19E94E50/sources/prop/INDEX
@@ -4,6 +4,8 @@ C%3A%2FUsers%2Ftpet0008%2FDocuments%2FHydroMap%2FR%2Fget.allData.R="E8D41511"
 C%3A%2FUsers%2Ftpet0008%2FDocuments%2FHydroMap%2FR%2Fget.objFunc.R="D3423677"
 C%3A%2FUsers%2Ftpet0008%2FDocuments%2FHydroMap%2FR%2Fget.pointData.R="00FD9ABF"
 C%3A%2FUsers%2Ftpet0008%2FDocuments%2FHydroMap%2FR%2Fget.smoothedDEM.R="FBFA4863"
+C%3A%2FUsers%2Ftpet0008%2FDocuments%2FHydroMap%2FR%2Fget.variogram.R="0243BC54"
+C%3A%2FUsers%2Ftpet0008%2FDocuments%2FHydroMap%2FR%2Fget.variogramModel.R="57B75005"
 C%3A%2FUsers%2Ftpet0008%2FDocuments%2FHydroMap%2FR%2Fkrige.head.R="0072A15F"
 C%3A%2FUsers%2Ftpet0008%2FDocuments%2FHydroMap%2FR%2Fkrige.head.calib.R="06EF12A0"
 C%3A%2FUsers%2Ftpet0008%2FDocuments%2FHydroMap%2FR%2Fkrige.head.crossval.R="AF65AB58"
diff --git a/.Rproj.user/19E94E50/sources/session-4dda64d8/0FF087FE-contents b/.Rproj.user/19E94E50/sources/session-4dda64d8/0FF087FE-contents
new file mode 100644
index 0000000..e69de29
diff --git a/.Rproj.user/19E94E50/sources/session-4dda64d8/221ACDA6-contents b/.Rproj.user/19E94E50/sources/session-4dda64d8/221ACDA6-contents
new file mode 100644
index 0000000..e69de29
diff --git a/.Rproj.user/19E94E50/sources/per/t/4C0D952C b/.Rproj.user/19E94E50/sources/session-4dda64d8/4C0D952C
similarity index 75%
rename from .Rproj.user/19E94E50/sources/per/t/4C0D952C
rename to .Rproj.user/19E94E50/sources/session-4dda64d8/4C0D952C
index 4bb0069..d577755 100644
--- a/.Rproj.user/19E94E50/sources/per/t/4C0D952C
+++ b/.Rproj.user/19E94E50/sources/session-4dda64d8/4C0D952C
@@ -3,7 +3,7 @@
     "path": "C:/Users/tpet0008/Documents/HydroMap/R/krige.head.R",
     "project_path": "R/krige.head.R",
     "type": "r_source",
-    "hash": "2758867917",
+    "hash": "1487510892",
     "contents": "",
     "dirty": false,
     "created": 1683426418056.0,
@@ -12,15 +12,15 @@
     "properties": {
         "source_window_id": "",
         "Source": "Source",
-        "cursorPosition": "202,24",
-        "scrollLine": "184"
+        "cursorPosition": "26,60",
+        "scrollLine": "11"
     },
     "folds": "",
-    "lastKnownWriteTime": 1683775888,
+    "lastKnownWriteTime": 1683860925,
     "encoding": "UTF-8",
     "collab_server": "",
     "source_window": "",
-    "last_content_update": 1683775888647,
+    "last_content_update": 1683860925272,
     "read_only": false,
     "read_only_alternatives": []
 }
\ No newline at end of file
diff --git a/.Rproj.user/19E94E50/sources/per/t/4C0D952C-contents b/.Rproj.user/19E94E50/sources/session-4dda64d8/4C0D952C-contents
similarity index 96%
rename from .Rproj.user/19E94E50/sources/per/t/4C0D952C-contents
rename to .Rproj.user/19E94E50/sources/session-4dda64d8/4C0D952C-contents
index fdcd247..b5c38cc 100644
--- a/.Rproj.user/19E94E50/sources/per/t/4C0D952C-contents
+++ b/.Rproj.user/19E94E50/sources/session-4dda64d8/4C0D952C-contents
@@ -10,18 +10,12 @@
 #' these features can be individually controlled by the user.
 #'
 #' Importantly, if the mapping parameters are not specified by the user, then this function estimates the parameters using
-#' a mixed data-type (i.e. real and integer parameters) split-sample maximum likelihood global optimisation. The optimisation by default
-#' includes the variogram parameters (e.g. range, sill and nugget) and the search parameters for local kriging (e.g. radius, minimum and
-#' maximum number of observations to use). Optimising these parameters is not common in kriging. It is done herein because trials for Victoria,
-#' Australia, showed that calibrating these parameters produced significantly lower cross-validation errors (i.e. the error in predicting the observations
-#' removed from the optimisation) compared to the standard approach of graphical estimation from an experimental variogram. The optimisation is
-#' numerically challenging and the following factors should be considered before use:
-#'
-#' \itemize{
-#'  \item{Optimisation of the parameters \code{mrvbf.pslope}, \code{mrvbf.ppctl} and \code{smooth.std} often required the creating of raster grids for every parameter combination. To ease the computation burden, these parameters should be treated as discrete, not continuous, numbers.}
-#'  \item{The optimisation package \code{rgeoud} is used herein. For control the optimisation process, consider directly using \code{\link{krige.head.calib}}.}
-#'  \item{Trials have established default calibration parameters and settings that were effective for Victoria, Australia. There is no guarantee they will be effective for other regions.}
-#' }
+#' a mixed data-type (i.e. real and integer parameters) split-sample maximum likelihood global optimisation. See \code{\link{krige.head.calib}} for details.
+#' 
+#' However, it is advised to separately calibrate the mapping parameters using \code{\link{krige.head.calib}}. The output from 
+#' the calibration can then be passed into \code{krige.head} via the input \code{calibration.results}. This will use the calibrated setting 
+#' and the training observed data set to map the head. To use all of the observed data, or alternate data, input \code{data}. 
+#' See \code{\link{krige.head.calib}} Examples for details.
 #'
 #' In using this function, the primary user decisions are:
 #' \itemize{
@@ -30,6 +24,9 @@
 #'  \item{The type of variogram model, defined by the input \code{model}}
 #'  }
 #'
+#' @param \code{calibration.results} List variable returned from \code{\link{krige.head.calib}}. Inputting this applies the calibration results here. All other inputs are ignores, except \code{data}, \code{newdata} and \code{grid}. 
+#' Defauli is \code{NULL}.
+#' 
 #' @param \code{formula} defines the R formula (as a character or formula data type) to be used to interpolate the heads. The left hand side of the formula must be \code{head}.
 #' The right hand side can contain any or all of the following terms: \code{elev} for the land surface elevation;
 #' \code{MrVBF} for the Multiresolution Index of Valley Bottom Flatness as a measure of valley-ness at each DEM grd cell;
@@ -165,17 +162,10 @@
 #' 
 #' # Export grids to an ARCMAP *.asc file
 #' write.asciigrid(heads,'DBNS.asc',3,na.value = -999)
-#' 
-#' # Re-calibrate the parameters and map using the default settings.
-#' # Note, here only those parameters that are a vector are calibrated. This is in addition
-#' # to the variogram parameters. 
-#' varigram.model = vgm(psill=10, model='Mat', range= 10000 , nugget=1, kappa=0.1);
-#' heads <- krige.head(formula=f, grid=DEM, data=obs.data, data.errvar.colname='total_err_var',
-#' nmax=Inf, maxdist=Inf, smooth.std = seq(0.5,1.5,length.out=11), 
-#' model = 'Mat',  fit.variogram.type=1, debug.level=1)
 #'
 #' @export
 krige.head <- function(
+  calibration.results = NULL,  
   formula = as.formula("head ~ elev + MrVBF + MrRTF + smoothing"),
   grid=NULL,
   grid.landtype.colname = NULL,
@@ -231,6 +221,32 @@ krige.head <- function(
   if (!do.grid.est && nsim!=0)
     stop('Point estimation cannot be undertaken for geostatistical simulations. To undertake gridded simulations, set newdata=NULL and nsim>0.')
 
+  # Use calibration results if provided. Else use individual input.
+  if (!is.null(calibration.results)) {
+    
+    formula = calibration.results$inputs$formula
+    if (is.null(grid))
+      grid = calibration.results$inputs$grid
+    grid.landtype.colname = calibration.results$inputs$grid.landtype.colname;
+    if (is.null(data))
+      data = calibration.results$inputs$data
+    data.fixedHead = calibration.results$inputs$data.fixedHead
+    if (is.null(newdata))
+      newdata = calibration.results$inputs$newdata
+    data.errvar.colname = calibration.results$inputsdata.errvar.colname
+    model= calibration.results$variogramModel
+    fit.variogram.type = 3; 
+    
+    mrvbf.pslope = calibration.results$parameter.Values$mrvbf.pslope
+    mrvbf.ppctl = calibration.results$parameter.Values$mrvbf.ppctl
+    smooth.std = calibration.results$parameter.Values$smooth.std
+    nmax = calibration.results$parameter.Values$nmax
+    nmax.fixedHead = calibration.results$parameter.Values$nmax.fixedHead
+    maxdist = calibration.results$parameter.Values$maxdist
+    nmin = calibration.results$parameter.Values$nmin 
+    omax = calibration.results$parameter.Values$omax
+  }
+  
   # Get variable names and check if using MrVBF or MrRTF or DEM
   if (debug.level>0)
     message('... Checking formula terms.');
diff --git a/.Rproj.user/19E94E50/sources/per/t/F6CA6EB6-contents b/.Rproj.user/19E94E50/sources/session-4dda64d8/5B688876-contents
similarity index 75%
rename from .Rproj.user/19E94E50/sources/per/t/F6CA6EB6-contents
rename to .Rproj.user/19E94E50/sources/session-4dda64d8/5B688876-contents
index 62c608c..d8f8c69 100644
--- a/.Rproj.user/19E94E50/sources/per/t/F6CA6EB6-contents
+++ b/.Rproj.user/19E94E50/sources/session-4dda64d8/5B688876-contents
@@ -121,140 +121,25 @@ get.objFunc <- function(
       return(Inf)
     } 
   }
-    
-  # Build variogram model
-  if (fit.variogram.type == 1) {
   
-    # Check the format of the model input.
-    if (class(model)[1] == "variogramModel") {
-
-      # Get model types
-      variogram.type = as.character(model$model)
-      
-      # Remove 'nug' from list of model types and add to 'model' variable
-      filt = variogram.type != 'Nug'      
-      modelType = variogram.type[filt]            
-    } else if (is.character(model)) {
-      modelType = model    
-    } else {
-      stop('When fit.variogram.type==1, the model should be a gstat variogram model or a vector of model types, eg c(\'Mat\')')
-    }
-
-    # Initialise variogram settings
-    variogram.nugg = NULL
-    
-    variogram.psill = c()
-    
-    variogram.range = c()
-    
-    variogram.kappa = c()
-    
-    variogram.ang1 = c()
-    
-    variogram.anis1 = c()
+  # Build variogram model and turn off future fitting via get.variogram()
+  if (fit.variogram.type == 1) {
+    model = get.variogramModel(params, param.names, model)
+    fit.variogram.type = 3
     
-    nModel = 0
     
+    # Get variogram parameter names and values - for reporting
     # Extract variogram settings.
     for (i in 1:length(param.names)) {
-      
       if (!is.na(pmatch('psill', param.names[i])) ||
           !is.na(pmatch('range', param.names[i])) ||
           !is.na(pmatch('kappa', param.names[i])) ||
           !is.na(pmatch('nug', param.names[i]))  ||
           !is.na(pmatch('ang1', param.names[i])) ||
           !is.na(pmatch('anis1', param.names[i]))) {
-        
-        if (is.na(pmatch('nug', param.names[i]))) {
-          iModel = as.numeric(substr(param.names[i], nchar(param.names[i]), nchar(param.names[i])))
-  
-          # Initialise some extra variogram parameters in case they're not used.        
-          if (iModel>nModel) {
-            variogram.kappa[iModel] = 0.5
-            variogram.ang1[iModel] = 0.0
-            variogram.anis1[iModel] = 1.0
-          }
-          
-        }
-          
-        if (!is.na(pmatch('nug', param.names[i]))) {
-          variogram.nug = max(0, params[i])
-                     
-        } else if (!is.na(pmatch('psill', param.names[i]))) {
-          variogram.psill[iModel] = max(0, params[i])
-          
-          if (iModel > nModel)
-            nModel = iModel
-          
-        } else if (!is.na(pmatch('range', param.names[i]))) {
-          variogram.range[iModel] = max(1e-6, params[i])
-          
-          if (iModel > nModel)
-            nModel = iModel
-          
-        } else if (!is.na(pmatch('kappa', param.names[i]))) {
-          variogram.kappa[iModel] = max(0, params[i])
-          
-          if (iModel > nModel)
-            nModel = iModel
-          
-        } else if (!is.na(pmatch('ang1', param.names[i]))) {
-          variogram.ang1[iModel] = params[i]
-          
-          if (iModel > nModel)
-            nModel = iModel
-          
-        } else if (!is.na(pmatch('anis1', param.names[i]))) {
-          variogram.anis1[iModel] = params[i]
-          
-          if (iModel > nModel)
-            nModel = iModel
-          
-        }
-        
         params2Print[i] = params[i]
       }
     }
-    
-    # Build variogram
-    for (i in 1:nModel) {
-      
-      #message(paste('... DEBGUGGING: i, model type, nug, psill, range, kappa ',i, model[i], variogram.nug, variogram.psill[i], variogram.range[i], variogram.kappa[i]))
-      
-      if (i == 1 && !is.null(variogram.nug)) {
-        vgm.model = vgm(
-          nugget = variogram.nug,
-          psill = variogram.psill[i],
-          range = variogram.range[i],
-          model = modelType[i],
-          kappa = variogram.kappa[i],
-          anis = c(variogram.ang1, variogram.anis1)
-        )
-      } else if (i == 1) {
-        vgm.model = vgm(
-          psill = variogram.psill[i],
-          range = variogram.range[i],
-          model = modelType[i],
-          kappa = variogram.kappa[i],
-          anis = c(variogram.ang1, variogram.anis1)
-        )
-      } else {
-        vgm.model = vgm(
-          psill = variogram.psill[i],
-          range = variogram.range[i],
-          model = modelType[i],
-          add.to = vgm.model,
-          kappa = variogram.kappa[i],
-          anis = c(variogram.ang1, variogram.anis1)
-        )
-      }
-    }
-    model = vgm.model
-    
-    
-    # Turn off fitting within get.variogram()
-    fit.variogram.type = 3
-    
   }
 
   # Call krige.head with the input parameter names.
diff --git a/.Rproj.user/19E94E50/sources/per/t/5C989810 b/.Rproj.user/19E94E50/sources/session-4dda64d8/5C989810
similarity index 92%
rename from .Rproj.user/19E94E50/sources/per/t/5C989810
rename to .Rproj.user/19E94E50/sources/session-4dda64d8/5C989810
index e408986..d686c31 100644
--- a/.Rproj.user/19E94E50/sources/per/t/5C989810
+++ b/.Rproj.user/19E94E50/sources/session-4dda64d8/5C989810
@@ -8,12 +8,12 @@
     "dirty": false,
     "created": 1683611401988.0,
     "source_on_save": false,
-    "relative_order": 4,
+    "relative_order": 5,
     "properties": {
         "source_window_id": "",
         "Source": "Source",
         "cursorPosition": "109,71",
-        "scrollLine": "88"
+        "scrollLine": "0"
     },
     "folds": "",
     "lastKnownWriteTime": 1683611445,
diff --git a/.Rproj.user/19E94E50/sources/per/t/5C989810-contents b/.Rproj.user/19E94E50/sources/session-4dda64d8/5C989810-contents
similarity index 100%
rename from .Rproj.user/19E94E50/sources/per/t/5C989810-contents
rename to .Rproj.user/19E94E50/sources/session-4dda64d8/5C989810-contents
diff --git a/.Rproj.user/19E94E50/sources/session-4dda64d8/648DFF61-contents b/.Rproj.user/19E94E50/sources/session-4dda64d8/648DFF61-contents
new file mode 100644
index 0000000..e69de29
diff --git a/.Rproj.user/19E94E50/sources/session-4dda64d8/65610FB2-contents b/.Rproj.user/19E94E50/sources/session-4dda64d8/65610FB2-contents
new file mode 100644
index 0000000..e69de29
diff --git a/.Rproj.user/19E94E50/sources/per/t/F6CA6EB6 b/.Rproj.user/19E94E50/sources/session-4dda64d8/842B6FEA
similarity index 64%
rename from .Rproj.user/19E94E50/sources/per/t/F6CA6EB6
rename to .Rproj.user/19E94E50/sources/session-4dda64d8/842B6FEA
index 8a77753..d8d7459 100644
--- a/.Rproj.user/19E94E50/sources/per/t/F6CA6EB6
+++ b/.Rproj.user/19E94E50/sources/session-4dda64d8/842B6FEA
@@ -1,26 +1,26 @@
 {
-    "id": "F6CA6EB6",
+    "id": "842B6FEA",
     "path": "C:/Users/tpet0008/Documents/HydroMap/R/get.objFunc.R",
     "project_path": "R/get.objFunc.R",
     "type": "r_source",
-    "hash": "4111900834",
+    "hash": "0",
     "contents": "",
     "dirty": false,
-    "created": 1683611463845.0,
+    "created": 1683855550751.0,
     "source_on_save": false,
-    "relative_order": 5,
+    "relative_order": 8,
     "properties": {
         "source_window_id": "",
         "Source": "Source",
-        "cursorPosition": "292,2",
-        "scrollLine": "0"
+        "cursorPosition": "93,8",
+        "scrollLine": "69"
     },
     "folds": "",
-    "lastKnownWriteTime": 1683611478,
+    "lastKnownWriteTime": 1683849812,
     "encoding": "UTF-8",
     "collab_server": "",
     "source_window": "",
-    "last_content_update": 1683611478454,
+    "last_content_update": 1683857723006,
     "read_only": false,
     "read_only_alternatives": []
 }
\ No newline at end of file
diff --git a/.Rproj.user/19E94E50/sources/session-4dda64d8/842B6FEA-contents b/.Rproj.user/19E94E50/sources/session-4dda64d8/842B6FEA-contents
new file mode 100644
index 0000000..d8f8c69
--- /dev/null
+++ b/.Rproj.user/19E94E50/sources/session-4dda64d8/842B6FEA-contents
@@ -0,0 +1,275 @@
+get.objFunc <- function(
+  params, 
+  param.names, 
+  formula, 
+  grid, 
+  grid.landtype.colname=NULL, 
+  data, 
+  data.fixedHead=NULL,
+  newdata=NULL, 
+  data.errvar.colname = NULL, 
+  model=c('Mat'),
+  nmin=0.2, 
+  omax=NULL,
+  fit.variogram.type=1, 
+  objFunc.type=1, 
+  lookup.table=NULL, 
+  return.transformed.params = T, 
+  return.predictions = F, 
+  use.cluster=TRUE, 
+  debug.level=0)  {
+
+  # DEBUGGING
+  op <- options("warn")
+  on.exit(options(op))
+  options(warn=1)  
+  
+  # Check if any parameters are nan or inf. If so, then later return obj function value inf
+  hasImplauibleParams = is.nan(params) | is.infinite(params)
+
+  # Loop through the input parameter names to update the values within the complete list of param names
+  # If the lookup table is not empty, then back transform the parameters from an integer to true parameter value.
+  paramsList = list(mrvbf.pslope=NULL,mrvbf.ppctl=NULL,smooth.std=NULL,nmax=NULL,nmax.fixedHead=NULL,maxdist=NULL)
+  params2Print = c()
+  param.beyond.bounds = F;
+  if (is.null(lookup.table)) {
+    for (i in 1:length(param.names)) {
+      # Skip if variogram variable and the variogram parameters are continuous and to be fitted.
+      if (fit.variogram.type == 1 &&
+          (
+            !is.na(pmatch('psill', param.names[i])) ||
+            !is.na(pmatch('range', param.names[i])) ||
+            !is.na(pmatch('kappa', param.names[i])) ||
+            !is.na(pmatch('nug', param.names[i])) ||
+            !is.na(pmatch('ang1', param.names[i])) ||
+            !is.na(pmatch('anis1', param.names[i]))
+          )) {
+        next
+      }
+      paramsList[[param.names[i]]] = params[i]
+      params2Print[i] = paramsList[[param.names[i]]]
+      
+    }
+  } else {
+    lookup.table.names = names(lookup.table)
+    # Loop through the loopup table and add default (ie non-calibrated) parameters to the list.
+    # This is required to allow calibration of a partial subset of the parameters.
+    for (i in 1:length(lookup.table.names)) {
+      if (nrow(lookup.table[[lookup.table.names[i]]]) == 1) {
+        paramsList[[lookup.table.names[i]]] = lookup.table[[lookup.table.names[i]]][1, 2];  
+      }
+    }
+    
+    # Handle the parameters to be calibrated.
+    for (i in 1:length(param.names)) {
+      # Skip if param is implausible. Note, this loop is still undertaken for the
+      # plausible parameters to return the parameter transformation for integer parameters
+      if (hasImplauibleParams[i]) {
+        params2Print[i] = params[i]
+        
+        next
+        
+      }
+      
+      # Skip if variogram variable and the variogram parameters are continuous and to be fitted.
+      if (fit.variogram.type == 1 &&
+          (
+            !is.na(pmatch('psill', param.names[i])) ||
+            !is.na(pmatch('range', param.names[i])) ||
+            !is.na(pmatch('kappa', param.names[i])) ||
+            !is.na(pmatch('nug', param.names[i])) ||
+            !is.na(pmatch('ang1', param.names[i])) ||
+            !is.na(pmatch('anis1', param.names[i]))
+          )) {
+        next
+      }
+      
+      # Assess if the parameter name is listed in the lookup table. If so, round the input
+      # parameter to an integer and add the rounded value to the vector of transformed parameter values
+      # to be returned to the optimisation routine
+      isIntegerParam =  param.names[i] %in% lookup.table.names
+      
+      # Skip if the parameter is not in the lookup table (eg is a real parameter)
+      if (!isIntegerParam) {
+        paramsList[[param.names[i]]] = params[i]
+        
+        params2Print[i] = params[i]
+        
+        next
+      }
+      
+      # Round parameter to an integer
+      params[i] = round(params[i], 0)
+      
+      if (params[i] < 1 || params[i] > dim(lookup.table[[param.names[i]]])[1])
+        param.beyond.bounds = T
+      else {
+          # Transform parameter from integer lookup value
+          paramsList[[param.names[i]]] = lookup.table[[param.names[i]]][params[i], 2]
+          
+          params2Print[i] = paramsList[[param.names[i]]]
+      }
+    }
+  }
+  
+  # Return Inf if any integer paramater is out of bounds or implausible.
+  if (param.beyond.bounds || any(hasImplauibleParams)) {
+    message('... Integer parameters are out of bounds. Returning Inf.');
+    if (return.transformed.params) { 
+      return(c(Inf, params))
+    } else {
+      return(Inf)
+    } 
+  }
+  
+  # Build variogram model and turn off future fitting via get.variogram()
+  if (fit.variogram.type == 1) {
+    model = get.variogramModel(params, param.names, model)
+    fit.variogram.type = 3
+    
+    
+    # Get variogram parameter names and values - for reporting
+    # Extract variogram settings.
+    for (i in 1:length(param.names)) {
+      if (!is.na(pmatch('psill', param.names[i])) ||
+          !is.na(pmatch('range', param.names[i])) ||
+          !is.na(pmatch('kappa', param.names[i])) ||
+          !is.na(pmatch('nug', param.names[i]))  ||
+          !is.na(pmatch('ang1', param.names[i])) ||
+          !is.na(pmatch('anis1', param.names[i]))) {
+        params2Print[i] = params[i]
+      }
+    }
+  }
+
+  # Call krige.head with the input parameter names.
+  heads = do.call(krige.head, c(paramsList, list(formula=formula, grid=grid, grid.landtype.colname=grid.landtype.colname, data=data, data.fixedHead=data.fixedHead, newdata=newdata, data.errvar.colname=data.errvar.colname, model=model,
+                nmin=nmin, omax=omax, fit.variogram.type=fit.variogram.type, use.cluster=use.cluster, debug.level=debug.level)))  
+
+  if (return.predictions) {
+    heads.to.return =  heads;
+  }
+  
+  # Create filters for erroneous estimates
+  notNA = !is.na(heads$head.pred) & !is.na(heads$head.var)
+  notNAN = !is.nan(heads$head.pred) & !is.nan(heads$head.var)
+  notNegVar = heads$head.var>0
+  
+  # Report on NAs filter
+  if (any(!notNA))
+    message(paste('... The following number of newdata estimates had a NA predictions or variances and have been removed from objective function calculation: ',sum(!notNA)),immediate.=T)
+
+  # Report on NANs filter
+  if (any(!notNAN))
+    message(paste('... The following number of newdata estimates had a NAN predictions or variances and have been removed from objective function calculation: ',sum(!notNAN)),immediate.=T)
+
+  # Report on variance filter
+  if (any(!notNegVar))
+    message(paste('... The following number of newdata estimates had a kriging variance <0 and have been removed from objective function calculation: ',sum(!notNegVar)),immediate.=T)
+
+  # Apply filters
+  heads = heads[notNA & notNAN & notNegVar,]
+  newdata = newdata[notNA & notNAN & notNegVar,]
+  
+  # Add filter to returned data
+  if (return.predictions) {
+    heads.to.return$is.valid.est =  notNA & notNAN & notNegVar;
+  }
+  
+  
+  # Filter residuals to remove 1%ile and 99%ile results
+  threshold_percentile = 1/100;
+  filt = heads$resid/sqrt(heads$head.var) > quantile(heads$resid/sqrt(heads$head.var),threshold_percentile) & heads$resid/sqrt(heads$head.var) < quantile(heads$resid/sqrt(heads$head.var),1-threshold_percentile);
+  newdata.at.tails = newdata[!filt,]
+  message('... Number of outlier newdata estimates removed from objective function: ',sum(!filt),'.')
+  message('    The mean and variance of the removed obs. newdata are ',mean(newdata.at.tails$head),' and ',var(newdata.at.tails$head) )  
+  heads = heads[filt,];
+  newdata = newdata[filt,];
+  nObs = length(heads); 
+
+  # Report summmary of est  
+  if (debug.level>0) {
+      message(paste('    The mean, min & max meters depth to water level is :',
+                    mean(heads$DEM-heads$head.pred), min(heads$DEM-heads$head.pred), max(heads$DEM-heads$head.pred) ))    
+    message(paste('    The mean, min & max meters head is :',
+                  mean(heads$head.pred), min(heads$head.pred), max(heads$head.pred) ))    
+    message(paste('    The mean, min & max meters DEM elevation is :',
+                  mean(heads$DEM), min(heads$DEM), max(heads$DEM) ))    
+    
+  }              
+  
+  # Report the number of artesian points.
+  is.artesian = heads$head.pred > heads$DEM & newdata$head < heads$DEM
+  message(paste('... The number of artesian newdata estimates = ',sum(is.artesian)))
+  
+  # Calculate the requested objective functions
+  objEst = vector(mode='numeric',length(objFunc.type))
+  for (i in 1:length(objFunc.type)) {
+      if (objFunc.type[i]==1) {
+        # Calculate liklihood est divided by the number of obs (MLKV:mean log kriging variance). 
+        # Taken from Samper and Neumen 1989 and Pardo-Iguzquiza & Dowd 2013.
+        NSSE = sum((heads$resid^2)/heads$head.var);
+        
+        MLKV = sum(log(heads$head.var))
+        
+        objEst[i]= (log(2*pi) + MLKV + NSSE)/nObs
+        
+      } else if (objFunc.type[i]==2) {
+        # As for type 1 BUT with an added penality for violoating a physical constraint.
+        # For simplicity, below the constraint is simply the head being below the land surface.
+        # This could be easily extended to include multiple spatially varying constraints.
+        resid = heads$resid;
+        var  = heads$head.var;    
+        if (any(is.artesian)) {
+          penalty.resid = heads$head.pred[is.artesian] - heads$DEM[is.artesian];
+          resid[is.artesian] = abs(resid[is.artesian]) + abs(penalty.resid);
+          message(paste('    The mean, min & max meters above the land surface from the artesian points is :',
+                        sum(is.artesian), mean(penalty.resid), min(penalty.resid) ))
+        }
+        NSSE = sum((resid^2)/var);
+        MLKV = sum(log(var))
+        objEst[i]= (log(2*pi) + MLKV + NSSE)/nObs;
+        
+        
+      } else if (objFunc.type[i]==3) {
+        # RMSE
+        objEst[i]= sqrt(mean(heads$resid^2));    
+        
+      } else if (objFunc.type[i]==4) {
+        # As for type 3 BUT with an added penality for violoating a physical constraint.
+        # For simplicity, below the constraint is simple the head being below the land surface.
+        # This could be easily extended to include multiple spatially varying constraints.
+        resid = heads$resid;
+        if (any(is.artesian)) {
+          penalty.resid = heads$head.pred[is.artesian] - heads$DEM[is.artesian];
+          resid[is.artesian] = abs(resid[is.artesian]) + abs(penalty.resid);
+          message(paste('    The mean, min & max meters above the land surface from the artesian points is :',
+                        sum(is.artesian), mean(penalty.resid), min(penalty.resid) ))
+        }    
+        objEst[i]= sqrt(sum(resid^2)/nObs);    
+        
+      } else
+        stop('objective function type unknown.')    
+  }
+
+    
+  # Print paramaters to be analysed
+  message('... Objective function results:')
+  for (i in 1:length(objFunc.type)) {
+    message(cat(format(c('Obj. Type',param.names,'Obj. func. val.'),width=16,justify='right')))
+    message(cat(format(c(objFunc.type[i],params2Print,objEst[i]),width=16,justify='right')))  
+  }
+  message('                               ')
+  
+  if (return.transformed.params && return.predictions) { 
+    return(list(obj.est =objEst, params.transformed = params, est = heads.to.return))
+  } else if (!return.transformed.params && return.predictions) { 
+    return(list(obj.est =objEst, est = heads.to.return))
+  } else if (return.transformed.params && !return.predictions) { 
+      return(c(objEst, params))
+  } else {
+    return(objEst)
+  } 
+  
+    
+}
diff --git a/.Rproj.user/19E94E50/sources/per/t/8974836D-contents b/.Rproj.user/19E94E50/sources/session-4dda64d8/8974836D-contents
similarity index 100%
rename from .Rproj.user/19E94E50/sources/per/t/8974836D-contents
rename to .Rproj.user/19E94E50/sources/session-4dda64d8/8974836D-contents
diff --git a/.Rproj.user/19E94E50/sources/session-4dda64d8/8C88B17A-contents b/.Rproj.user/19E94E50/sources/session-4dda64d8/8C88B17A-contents
new file mode 100644
index 0000000..a48959e
--- /dev/null
+++ b/.Rproj.user/19E94E50/sources/session-4dda64d8/8C88B17A-contents
@@ -0,0 +1,183 @@
+get.variogram <- function(formula, data, cutoff=120000, width=500, model=c('Mat'), 
+                          nmin=0, nmax=Inf, maxdist=Inf,  tol.RMSE = 1e-4, max.Its=10, 
+                          fit.variogram.starts = 10, debug.level=0 ) {
+
+  if (debug.level>0) 
+    message('Building variogram model:')
+
+  # Check model type is acceptable
+  if (is.character(model)) {
+    # get list of acceptable model 
+    model.types = vgm();
+    for (i in 1:length(model)) { 
+      if (!model[i] %in% model)
+	      stop(paste('The following model type is not an accepted by gstat vgm():',model[i]))
+    }
+  }
+  
+  # Check if head of depth to be est.
+  if (debug.level>0) 
+    message('... Checking formula terms.');
+  var.names <- all.vars(formula);
+  do.head.est <- any(match(var.names, 'head'),na.rm=TRUE);
+  do.depth.est <- any(match(var.names, 'depth'),na.rm=TRUE);	   
+  if (do.depth.est == FALSE && do.head.est ==FALSE)
+    stop('   Variogram formula LHS must be "depth" or "head".')
+     
+  
+  # Estimate initial resdiuals
+  if (debug.level>0) 
+    message('... Building initial regression model for residuals.');
+  fit <- lm(formula, data=data);
+  resid = residuals(fit);
+  data$resid = resid;
+  
+
+  # Filter out NAs  
+  filt_na = is.na(data$resid);
+  if (any(filt_na)) {
+    if (debug.level>0) 
+      message('... Filtering out ', sum(filt_na), ' NA values from point data residuals.');
+    dataFilt =  data[!filt_na,];    
+  } else
+    dataFilt =  data
+  
+  # Define/get initial variogram
+  if (debug.level>0) 
+    message('... Extracting input variogram settings.');
+  if (class(model)[1] =="variogramModel") {
+    filt = model$model!='Nug';
+    variogram.type = model$model[filt] 
+    variogram.range = model$range[filt] 
+    variogram.psill = model$psill[filt]
+    variogram.nugget = model$psill[!filt]
+    variogram.kappa= model$kappa[filt]
+    variogram.ang1= model$ang1[filt]
+    variogram.anis1= model$anis1[filt]
+  } else {
+    if (is.character(model)) {
+      variogram.type = model;
+    } else
+      variogram.type ='Mat';
+    
+    variogram.range = 25000;
+    variogram.psill = 20;
+    variogram.nugget = 1;
+    variogram.kappa= 0.5;
+  }
+  
+  # Fit variogam to initial residuals
+  if (debug.level>0) 
+    message('... Building initial variogram.');
+  variogram_exp = variogram(resid ~ 1, dataFilt, cressie=FALSE, cutoff=cutoff, width=width);
+  model = vgm(nugget=variogram.nugget,psill=variogram.psill, model=variogram.type, range=variogram.range, kappa=variogram.kappa);
+  fit.kappa=F
+  if (variogram.type=='Mat')
+	  fit.kappa=T	
+  if (debug.level>0) 
+    message('... Fitting initial variogram.');
+  model  = fit.variogram(variogram_exp, fit.sills=TRUE,  fit.ranges=TRUE, fit.kappa = fit.kappa, model = model);       
+
+  # Update initial variogram
+  if (debug.level>0) 
+    message('... Extracting fitted initial variogram values.');
+  variogram.type = model$model[2] 
+  variogram.nugget = if(is.finite(model$psill[1])){model$psill[1]}else{1}
+  variogram.range = if(is.finite(model$range[2])){model$range[2]}else{25000}
+  variogram.psill = if(is.finite(model$psill[2])){model$psill[2]}else{20}
+  variogram.kappa= if(is.finite(model$kappa[2])){model$kappa[2]}else{0.5}
+  variogram.ang1= model$ang1[2]
+  variogram.anis1= model$anis1[2]
+  if (debug.level>0){
+    message(paste('    type =',variogram.type))
+    message(paste('    nugget =',variogram.nugget))
+    message(paste('    psill =',variogram.psill))
+    message(paste('    range =',variogram.range))
+    message(paste('    kappa =',variogram.kappa))
+  }
+  
+  # Return if the model is not to be refined.
+  if (max.Its==0) 
+    return(model)
+                     
+  # Filter out NAs  
+  filt_na = is.na(data$resid);
+  if (any(filt_na)) {
+     dataFilt =  data[!filt_na,];    
+  } else
+     dataFilt =  data;	
+    
+    
+  # Fit the model variogram using many initial starts and select the
+  # model producing the lowest weights sum of squared error
+  # --------------------------------------------------------------           
+
+    # Initialisng data
+    SSE_varg = matrix(Inf,nrow=fit.variogram.starts, ncol=5)
+
+    # Cycle through each start. 
+    if (debug.level>0) 
+      message('... Doing multistart variogram fitting...');  
+    for (j in 1:fit.variogram.starts) {
+
+    	# Randomly sample all model parameters. The samples are 
+    	# are +- one order from the parameter input.
+    	sill = 0.1*variogram.psill + (10.0*variogram.psill - 0.1*variogram.psill) * runif(1);
+    	nug = 0.1*variogram.nugget + (10.0*variogram.nugget - 0.1*variogram.nugget) * runif(1);
+    	range = 0.1*variogram.range + (10.0*variogram.range - 0.1*variogram.range) * runif(1);
+    	if (fit.kappa)
+    	   kappa = 0.1*variogram.kappa + (10.0*variogram.kappa - 0.1*variogram.kappa) * runif(1);
+    
+    	# Build experimental and fit model variogram;	
+    	if (fit.kappa) {
+     	   model_sample = vgm(nugget=nug,psill=sill, model=variogram.type, range=range, kappa=kappa);
+    	} else {
+     	   model_sample = vgm(nugget=nug,psill=sill, model=variogram.type, range=range);
+    	}
+    	model_sample  = fit.variogram(variogram_exp, fit.sills=TRUE,  fit.ranges=TRUE, fit.kappa = fit.kappa, model = model_sample,  warn.if.neg = FALSE, debug.level=0);           
+    
+    	# Get the weighted sum of squares error and report to the user.
+    	SSEerr = attr(model_sample, 'SSErr')
+    	isSingular = attr(model_sample, 'singular')
+    
+    	if (!isSingular) {
+    	  SSE_varg[j,1] = SSEerr;
+    	  SSE_varg[j,2] = model_sample$psill[1];
+    	  SSE_varg[j,3] =  model_sample$psill[2];
+    	  SSE_varg[j,4] = model_sample$range[2];
+    	  SSE_varg[j,5] = model_sample$kappa[2];
+    	  
+    	  if (debug.level>0) 
+    	    message(paste('    Fit number ',j,' SSE =',SSEerr));
+    	  
+    	} else {
+    	  if (debug.level>0) 
+    	    message(paste('    Fit number ',j,' failed because of singularity.'));
+    	  
+    	}
+  }
+      
+  # Find variogram with lowest SSE
+  if (debug.level>0) 
+    message('... Finding lowest SSE variogram fit.');
+  filt = SSE_varg[,1] == min(SSE_varg[,1]);
+  SSE_best = SSE_varg[filt,1]
+  variogram.type = model$model[2] 
+  variogram.nugget = SSE_varg[filt,2]
+  variogram.psill = SSE_varg[filt,3]
+  variogram.range = SSE_varg[filt,4]
+  variogram.kappa = SSE_varg[filt,5]
+  if (debug.level>0) {
+    message(paste('    type =',variogram.type))
+    message(paste('    nugget =',variogram.nugget))
+    message(paste('    psill =',variogram.psill))
+    message(paste('    range =',variogram.range))
+    message(paste('    kappa =',variogram.kappa))
+    message(paste('    SSE =',SSE_best))
+    message('... Building lowest SSE variogram.');
+  }
+  
+  model = vgm(nugget=variogram.nugget,psill=variogram.psill, model=variogram.type, range=variogram.range, kappa=variogram.kappa);
+
+  return(model)
+}
diff --git a/.Rproj.user/19E94E50/sources/session-4dda64d8/9E018CFD-contents b/.Rproj.user/19E94E50/sources/session-4dda64d8/9E018CFD-contents
new file mode 100644
index 0000000..e69de29
diff --git a/.Rproj.user/19E94E50/sources/session-4dda64d8/A6CEE51C-contents b/.Rproj.user/19E94E50/sources/session-4dda64d8/A6CEE51C-contents
new file mode 100644
index 0000000..e69de29
diff --git a/.Rproj.user/19E94E50/sources/session-4dda64d8/A9B84F7A-contents b/.Rproj.user/19E94E50/sources/session-4dda64d8/A9B84F7A-contents
new file mode 100644
index 0000000..e69de29
diff --git a/.Rproj.user/19E94E50/sources/per/u/B06A395E b/.Rproj.user/19E94E50/sources/session-4dda64d8/B06A395E
similarity index 80%
rename from .Rproj.user/19E94E50/sources/per/u/B06A395E
rename to .Rproj.user/19E94E50/sources/session-4dda64d8/B06A395E
index 5f6b1a2..70948fa 100644
--- a/.Rproj.user/19E94E50/sources/per/u/B06A395E
+++ b/.Rproj.user/19E94E50/sources/session-4dda64d8/B06A395E
@@ -3,17 +3,17 @@
     "path": null,
     "project_path": null,
     "type": "r_source",
-    "hash": "3925697018",
+    "hash": "0",
     "contents": "",
     "dirty": true,
     "created": 1683713112391.0,
     "source_on_save": false,
-    "relative_order": 10,
+    "relative_order": 9,
     "properties": {
         "tempName": "Untitled1",
         "source_window_id": "",
         "Source": "Source",
-        "cursorPosition": "33,65",
+        "cursorPosition": "36,92",
         "scrollLine": "0"
     },
     "folds": "",
@@ -21,7 +21,7 @@
     "encoding": "",
     "collab_server": "",
     "source_window": "",
-    "last_content_update": 1683776039788,
+    "last_content_update": 1683858735062,
     "read_only": false,
     "read_only_alternatives": []
 }
\ No newline at end of file
diff --git a/.Rproj.user/19E94E50/sources/per/u/B06A395E-contents b/.Rproj.user/19E94E50/sources/session-4dda64d8/B06A395E-contents
similarity index 83%
rename from .Rproj.user/19E94E50/sources/per/u/B06A395E-contents
rename to .Rproj.user/19E94E50/sources/session-4dda64d8/B06A395E-contents
index f460426..01d2382 100644
--- a/.Rproj.user/19E94E50/sources/per/u/B06A395E-contents
+++ b/.Rproj.user/19E94E50/sources/session-4dda64d8/B06A395E-contents
@@ -31,4 +31,7 @@ obs.data$total_err_var = pmax(obs.data$total_err_var, 0.05)
 # Calibrate the mapping parameters with 25% of the data randomly selected at 2 cores.
 calib.results <- krige.head.calib(formula=f, grid=DEM, data=obs.data, newdata=0.25, nmin=0, nmax=Inf, maxdist=Inf, omax=0,
                                   data.errvar.colname='total_err_var', model = variogram.model,  fit.variogram.type=1,
-                                  debug.level=0, use.cluster = 2)
\ No newline at end of file
+                                  pop.size.multiplier=1, debug.level=0, use.cluster = 2)
+
+# Map calibration results using all of the observed data abd all cores. 
+head.grid <- krige.head(calibration.results = calib.results, data=obs.data, use.cluster = T)
\ No newline at end of file
diff --git a/.Rproj.user/19E94E50/sources/session-4dda64d8/B094F956-contents b/.Rproj.user/19E94E50/sources/session-4dda64d8/B094F956-contents
new file mode 100644
index 0000000..e69de29
diff --git a/.Rproj.user/19E94E50/sources/per/t/B29BDB7C b/.Rproj.user/19E94E50/sources/session-4dda64d8/B29BDB7C
similarity index 78%
rename from .Rproj.user/19E94E50/sources/per/t/B29BDB7C
rename to .Rproj.user/19E94E50/sources/session-4dda64d8/B29BDB7C
index d7f0db0..24f47f7 100644
--- a/.Rproj.user/19E94E50/sources/per/t/B29BDB7C
+++ b/.Rproj.user/19E94E50/sources/session-4dda64d8/B29BDB7C
@@ -3,23 +3,23 @@
     "path": "C:/Users/tpet0008/Documents/HydroMap/DESCRIPTION",
     "project_path": "DESCRIPTION",
     "type": "dcf",
-    "hash": "2087192970",
+    "hash": "1037317559",
     "contents": "",
     "dirty": false,
     "created": 1683421701566.0,
     "source_on_save": false,
     "relative_order": 1,
     "properties": {
-        "cursorPosition": "2,12",
+        "cursorPosition": "2,16",
         "scrollLine": "0",
         "Source": "Source"
     },
     "folds": "",
-    "lastKnownWriteTime": 1683789589,
+    "lastKnownWriteTime": 1683860980,
     "encoding": "UTF-8",
     "collab_server": "",
     "source_window": "",
-    "last_content_update": 1683789589561,
+    "last_content_update": 1683860980016,
     "read_only": false,
     "read_only_alternatives": []
 }
\ No newline at end of file
diff --git a/.Rproj.user/19E94E50/sources/per/t/B29BDB7C-contents b/.Rproj.user/19E94E50/sources/session-4dda64d8/B29BDB7C-contents
similarity index 95%
rename from .Rproj.user/19E94E50/sources/per/t/B29BDB7C-contents
rename to .Rproj.user/19E94E50/sources/session-4dda64d8/B29BDB7C-contents
index 6eeb069..84feeb4 100644
--- a/.Rproj.user/19E94E50/sources/per/t/B29BDB7C-contents
+++ b/.Rproj.user/19E94E50/sources/session-4dda64d8/B29BDB7C-contents
@@ -1,6 +1,6 @@
 Package: HydroMap
 Title: HydroMap
-Version: 0.2
+Version: 0.2.1.0
 Authors@R: person("Tim", "Peterson", email = "timjp@unimelb.edu.au", role = c("aut", "cre"))
 Description: Interpolates groundwater level observations giving realistic flow lines.
 Depends: R (>= 3.2.2)
diff --git a/.Rproj.user/19E94E50/sources/per/t/D0D0B5EC b/.Rproj.user/19E94E50/sources/session-4dda64d8/D0D0B5EC
similarity index 70%
rename from .Rproj.user/19E94E50/sources/per/t/D0D0B5EC
rename to .Rproj.user/19E94E50/sources/session-4dda64d8/D0D0B5EC
index fc5e486..7000260 100644
--- a/.Rproj.user/19E94E50/sources/per/t/D0D0B5EC
+++ b/.Rproj.user/19E94E50/sources/session-4dda64d8/D0D0B5EC
@@ -3,24 +3,24 @@
     "path": "C:/Users/tpet0008/Documents/HydroMap/R/krige.head.calib.R",
     "project_path": "R/krige.head.calib.R",
     "type": "r_source",
-    "hash": "3896465155",
+    "hash": "0",
     "contents": "",
     "dirty": false,
     "created": 1683611155091.0,
     "source_on_save": false,
-    "relative_order": 3,
+    "relative_order": 4,
     "properties": {
         "source_window_id": "",
         "Source": "Source",
-        "cursorPosition": "4,4",
+        "cursorPosition": "6,90",
         "scrollLine": "0"
     },
-    "folds": "465|28|554|4|\n562|33|577|4|\n",
-    "lastKnownWriteTime": 1683789567,
+    "folds": "458|28|547|4|\n555|33|570|4|\n",
+    "lastKnownWriteTime": 1683859970,
     "encoding": "UTF-8",
     "collab_server": "",
     "source_window": "",
-    "last_content_update": 1683789567135,
+    "last_content_update": 1683859970018,
     "read_only": false,
     "read_only_alternatives": []
 }
\ No newline at end of file
diff --git a/.Rproj.user/19E94E50/sources/per/t/D0D0B5EC-contents b/.Rproj.user/19E94E50/sources/session-4dda64d8/D0D0B5EC-contents
similarity index 90%
rename from .Rproj.user/19E94E50/sources/per/t/D0D0B5EC-contents
rename to .Rproj.user/19E94E50/sources/session-4dda64d8/D0D0B5EC-contents
index 23c7ebe..995198e 100644
--- a/.Rproj.user/19E94E50/sources/per/t/D0D0B5EC-contents
+++ b/.Rproj.user/19E94E50/sources/session-4dda64d8/D0D0B5EC-contents
@@ -3,10 +3,8 @@
 #' \code{krige.head.calib} calibrates parameters to minimise the interpolation error.
 #'
 #' @details This function optimises the parameters for the spatial interpolation so as to minimise the
-#' interpolation error - which by default is defined using a formal likelihood function. Once this function has return
-#' the optimal parameter set, then the mapping can be undertaken using \code{\link{krige.head}}.
-#' Importantly, the mapping approach used for the calibration (e.g. the formula etc) must be identical to that
-#' input to \code{\link{krige.head}}.
+#' interpolation error, which by default is defined using a formal likelihood function. Once this function has return
+#' the optimal parameter set, the mapping can be undertaken using \code{\link{krige.head}} (see examples below).
 #'
 #' The mapping parameters are estimated using a mixed data-type (i.e. real and integer parameters) split-sample maximum likelihood global optimisation. The optimisation by default
 #' includes the variogram parameters (e.g. range, sill and nuggest) and the search parameters for local kriging (e.g. radius, minimum and
@@ -22,8 +20,8 @@
 #'  \item{The available input point data must be split up into \code{data} and \code{newdata}. The point observations within \code{data} are used to estimate the water level at the locations defined
 #'  within \code{newdata}. The difference between the observed and estimated values defines the interpolation prediction error. This process also gives the kriging variance at
 #'  each \code{newdata} location. The points used for \code{data} and \code{newdata} should both cover all types of terrain, geology, landuse and the full mapping extend so as to avoid bias. Also, changing either may
-#'  change the calibration solution. }
-#'  \item{The \code{objFunc.type} allows for maximum likelihood estimation - that is optimisation that accounts for the expected kriging error. See \code{objFunc.type} below for details.
+#'  change the calibration solution.}
+#'  \item{The \code{objFunc.type} allows for maximum likelihood estimation - that is optimisation that accounts for the expected kriging error. See \code{objFunc.type} below for details.}
 #' }
 #'
 #' In using this function, the primary user decisions are:
@@ -131,7 +129,7 @@
 #'
 #' # Load a model variogram and mapping parameters found to be effective.
 #' data('mapping.parameters')
-#
+#'
 #' # Define a simple kriging formula without MrVBF terms that does not require the package RSAGA.
 #' f <- as.formula('head ~ elev + smoothing')
 #'
@@ -140,33 +138,17 @@
 #'
 #' # Calibrate the mapping parameters with 25% of the data randomly selected at 2 cores.
 #' calib.results <- krige.head.calib(formula=f, grid=DEM, data=obs.data, newdata=0.25, nmin=0, nmax=Inf, maxdist=Inf, omax=0,
-#'    data.errvar.colname='total_err_var', model = variogram.model,  fit.variogram.type=1,
-#'    debug.level=0, use.cluster = 2)
-#'
-#' # Reformat the parameter to a named list. This is just for ease of use and the need for
-#' this may be removed in future versions.
-#' params.all = c(mrvbf.pslope=NA,mrvbf.ppctl=NA,smooth.std=NA,nmax=NA,maxdist=NA,
-#'                kappa_model_1 = NA, ang1_model_1 = NA, anis1_model_1 = NA,
-#'                nmax.fixedHead = NA);
-#' for (i in 1:length(head.calib$param.Names)) {
-#'  params.all[[calib.results$param.Names[i]]] = calib.results$par[i]
-#' }
-#' params.all.names=names(params.all);
-#'
-#' # Rebuild the variogram with the calibrated parameters.
-#' varigram.model = vgm(psill=params.all[['psill_model_1']], model='Mat',
-#'                    range= params.all[['range_model_1']] , nugget=params.all[['nug']],
-#'                    kappa=params.all[['kappa_model_1']] );
+#'                       data.errvar.colname='total_err_var', model = variogram.model,  fit.variogram.type=1,
+#'                       pop.size.multiplier=1, debug.level=0, use.cluster = 2)
 #'
-#' # Grid the observed head data.
-#' head.map <- krige.head(formula=f, grid=DEM, data=obs.data, data.errvar.colname='total_err_var',
-#' model=variogram.model, smooth.std=params.all[['smooth.std']], maxdist=params.all[['maxdist']],
-#' nmax=params.all[['nmax']], fit.variogram.type=3, debug.level=1)
+#' # Grid the observed head using ONLY the training data from the calibration (and all cores) and then map..
+#' head.grid.training <- krige.head(calibration.results = calib.results, use.cluster = T)
+#' spplot(head.grid.training)
+#' 
+#' # Grid the observed head using ALL of the obsrved data from the calibration (and all cores).
+#' head.grid.all <- krige.head(calibration.results = calib.results, data=obs.data, use.cluster = T)
+#' spplot(head.grid.all)
 #'
-#' # Write the grids to ARCMAP ASCII grids
-#' write.asciigrid(head.map, 'head.asc',1);
-#' write.asciigrid(head.map, 'KrigingVariance.asc',2);
-
 #' @export
 krige.head.calib <-
   function(formula = as.formula("head ~ elev + MrVBF + MrRTF + smoothing"),
@@ -197,6 +179,17 @@ krige.head.calib <-
     if (debug.level>0)
       message('Calibrating mapping parameters.')
 
+    # Store input parameter values (exluding variogram params). These are updated if they're calibrated.
+    solution = list()
+    solution$parameter.Values$mrvbf.pslope = mrvbf.pslope
+    solution$parameter.Values$mrvbf.ppctl = mrvbf.ppctl
+    solution$parameter.Values$smooth.std = smooth.std
+    solution$parameter.Values$nmax = nmax
+    solution$parameter.Values$nmax.fixedHead = nmax.fixedHead
+    solution$parameter.Values$maxdist = maxdist
+    solution$parameter.Values$nmin = nmin
+    solution$parameter.Values$omax = omax
+    
     # Assess formula componants
     var.names = all.vars(as.formula(formula));
     use.MrVBF = any(match(var.names, 'MrVBF'));
@@ -757,7 +750,7 @@ krige.head.calib <-
     } else {
       message('... Starting continuous value calibration of the following parameters:',paste( param.Names,sep="",collapse=', '))
     }
-    solution = genoud(
+    genoud.solution = genoud(
       HydroMap:::get.objFunc,
       nvars = nParams,
       max = doObjMaximisation,
@@ -806,7 +799,7 @@ krige.head.calib <-
 
     # Store the returned value for use when calling get.objFunc directly.
     # This is required if some of the parameters have discrete parameters (which are used by the lookup table to transform to meaningful values)
-    solution$par.pretransformed = solution$par;
+    genoud.solution$par.pretransformed = genoud.solution$par;
 
     # Back transform parameters
     lookup.table.names = names(param.IntegerLookupTable)
@@ -817,26 +810,79 @@ krige.head.calib <-
 
       # Get non-transformed value from table.
       if (isIntegerParam) {
-        solution$par[i] = param.IntegerLookupTable[[param.Names[i]]][solution$par[i], 2]
+        genoud.solution$param.Values[i] = param.IntegerLookupTable[[param.Names[i]]][genoud.solution$par[i], 2]
       }
     }
 
-
     # Print summary
     message(paste('    Best objective function value:', solution$value))
     for (i in 1:nParams)
-      message(paste('    Parameter ', param.Names[i], ' = ', solution$par[i]))
-
-    # Add parameter names to the solution list variable.
-    solution$param.Names = param.Names
-
+      message(paste('    Parameter ', param.Names[i], ' = ', genoud.solution$param.Values[i]))
+    
+    # Add formula
+    solution$inputs$formula = formula;     
+    
+    # Add data and new data
+    solution$inputs$data = data
+    solution$inputs$newdata = newdata
+    solution$inputs$grid = grid
+    
+    # Add remaining inputs 
+    solution$inputs$grid.landtype.colname = grid.landtype.colname
+    solution$inputs$data.fixedHead = data.fixedHead
+    solution$inputs$data.errvar.colname = data.errvar.colname  
+
+    ind = match(param.Names,'mrvbf.pslope')
+    ind = na.omit(ind)
+    if (length(ind)==1 && is.finite(ind))
+      solution$parameter.Values$mrvbf.pslope = genoud.solution$param.Values[ind]
+    
+    ind = match(param.Names,'mrvbf.ppctl')
+    ind = na.omit(ind)
+    if (length(ind)==1 && is.finite(ind))
+      solution$parameter.Values$mrvbf.ppctl = genoud.solution$param.Values[ind]
+    
+    ind = match(param.Names,'smooth.std')
+    ind = na.omit(ind)
+    if (length(ind)==1 && is.finite(ind))
+      solution$parameter.Values$smooth.std = genoud.solution$param.Values[ind]
+    
+    ind = match(param.Names,'nmax')
+    ind = na.omit(ind)
+    if (length(ind)==1 && is.finite(ind))
+      solution$parameter.Values$nmax = genoud.solution$param.Values[ind]
+    
+    ind = match(param.Names,'nmax.fixedHead')
+    ind = na.omit(ind)
+    if (length(ind)==1 && is.finite(ind))
+      solution$parameter.Values$nmax.fixedHead = genoud.solution$param.Values[ind]
+    
+    ind = match(param.Names,'maxdist')
+    ind = na.omit(ind)
+    if (length(ind)==1 && is.finite(ind))
+      solution$parameter.Values$maxdist = genoud.solution$param.Values[ind]
+    
+    ind = match(param.Names,'nmin')
+    ind = na.omit(ind)
+    if (length(ind)==1 && is.finite(ind))
+      solution$parameter.Values$nmin = genoud.solution$param.Values[ind]
+    
+    ind = match(param.Names,'omax')
+    ind = na.omit(ind)
+    if (length(ind)==1 && is.finite(ind))
+      solution$parameter.Values$omax = genoud.solution$param.Values[ind]
+    
+    # Add variogram model
+    if (fit.variogram.type == 1)
+      solution$variogramModel = get.variogramModel(genoud.solution$par, param.Names, model)
+    
+    
     # Add the lookup table to the solution. This is required if the the objectiveget.objFunc() is to be directly called.
     solution$lookup.table = param.IntegerLookupTable;
-
-    # Add data and new data
-    solution$data = data
-    solution$newdata = newdata
-
+    
+    # Add calibration settings
+    solution$genoud.output = genoud.solution
+    
     return(solution)
 
   }
diff --git a/.Rproj.user/19E94E50/sources/per/t/D19F9BF0-contents b/.Rproj.user/19E94E50/sources/session-4dda64d8/D19F9BF0-contents
similarity index 100%
rename from .Rproj.user/19E94E50/sources/per/t/D19F9BF0-contents
rename to .Rproj.user/19E94E50/sources/session-4dda64d8/D19F9BF0-contents
diff --git a/.Rproj.user/19E94E50/sources/per/t/F1F4D846 b/.Rproj.user/19E94E50/sources/session-4dda64d8/F1F4D846
similarity index 96%
rename from .Rproj.user/19E94E50/sources/per/t/F1F4D846
rename to .Rproj.user/19E94E50/sources/session-4dda64d8/F1F4D846
index a805e78..68bcf19 100644
--- a/.Rproj.user/19E94E50/sources/per/t/F1F4D846
+++ b/.Rproj.user/19E94E50/sources/session-4dda64d8/F1F4D846
@@ -8,7 +8,7 @@
     "dirty": false,
     "created": 1683621910078.0,
     "source_on_save": false,
-    "relative_order": 8,
+    "relative_order": 7,
     "properties": {
         "source_window_id": "",
         "Source": "Source",
diff --git a/.Rproj.user/19E94E50/sources/per/t/F1F4D846-contents b/.Rproj.user/19E94E50/sources/session-4dda64d8/F1F4D846-contents
similarity index 100%
rename from .Rproj.user/19E94E50/sources/per/t/F1F4D846-contents
rename to .Rproj.user/19E94E50/sources/session-4dda64d8/F1F4D846-contents
diff --git a/.Rproj.user/19E94E50/sources/per/t/8974836D b/.Rproj.user/19E94E50/sources/session-4dda64d8/F6CA6EB6
similarity index 51%
rename from .Rproj.user/19E94E50/sources/per/t/8974836D
rename to .Rproj.user/19E94E50/sources/session-4dda64d8/F6CA6EB6
index 9a8c47b..e5a21eb 100644
--- a/.Rproj.user/19E94E50/sources/per/t/8974836D
+++ b/.Rproj.user/19E94E50/sources/session-4dda64d8/F6CA6EB6
@@ -1,26 +1,26 @@
 {
-    "id": "8974836D",
-    "path": "C:/Users/tpet0008/Documents/HydroMap/R/set.env.R",
-    "project_path": "R/set.env.R",
+    "id": "F6CA6EB6",
+    "path": "C:/Users/tpet0008/Documents/HydroMap/R/get.variogramModel.R",
+    "project_path": "R/get.variogramModel.R",
     "type": "r_source",
-    "hash": "3237158787",
+    "hash": "3069732398",
     "contents": "",
     "dirty": false,
-    "created": 1683675524014.0,
+    "created": 1683611463845.0,
     "source_on_save": false,
-    "relative_order": 9,
+    "relative_order": 6,
     "properties": {
         "source_window_id": "",
         "Source": "Source",
-        "cursorPosition": "5,11",
+        "cursorPosition": "2,2",
         "scrollLine": "0"
     },
     "folds": "",
-    "lastKnownWriteTime": 1683421624,
+    "lastKnownWriteTime": 1683849805,
     "encoding": "UTF-8",
     "collab_server": "",
     "source_window": "",
-    "last_content_update": 1683421624,
+    "last_content_update": 1683849805552,
     "read_only": false,
     "read_only_alternatives": []
 }
\ No newline at end of file
diff --git a/.Rproj.user/19E94E50/sources/session-4dda64d8/F6CA6EB6-contents b/.Rproj.user/19E94E50/sources/session-4dda64d8/F6CA6EB6-contents
new file mode 100644
index 0000000..b11af11
--- /dev/null
+++ b/.Rproj.user/19E94E50/sources/session-4dda64d8/F6CA6EB6-contents
@@ -0,0 +1,129 @@
+get.variogramModel <- function(
+  params, 
+  param.names, 
+  model=c('Mat'))  {
+
+  # DEBUGGING
+  op <- options("warn")
+  on.exit(options(op))
+  options(warn=1)  
+  
+  # Check if any parameters are nan or inf. If so, then later return obj function value inf
+  hasImplauibleParams = is.nan(params) | is.infinite(params)
+
+  # Check the format of the model input.
+  if (class(model)[1] == "variogramModel") {
+    
+    # Get model types
+    variogram.type = as.character(model$model)
+    
+    # Remove 'nug' from list of model types and add to 'model' variable
+    filt = variogram.type != 'Nug'      
+    modelType = variogram.type[filt]            
+  } else if (is.character(model)) {
+    modelType = model    
+  } else {
+    stop('When calling get.variogramModel(), fit.variogram.type must equal and the model should be a gstat variogram model or a vector of model types, eg c(\'Mat\')')
+  }
+  
+  # Initialise variogram settings
+  variogram.nugg = NULL
+  variogram.psill = c()
+  variogram.range = c()
+  variogram.kappa = c()
+  variogram.ang1 = c()
+  variogram.anis1 = c()
+  nModel = 0
+  
+  # Extract variogram settings.
+  for (i in 1:length(param.names)) {
+    if (!is.na(pmatch('psill', param.names[i])) ||
+        !is.na(pmatch('range', param.names[i])) ||
+        !is.na(pmatch('kappa', param.names[i])) ||
+        !is.na(pmatch('nug', param.names[i]))  ||
+        !is.na(pmatch('ang1', param.names[i])) ||
+        !is.na(pmatch('anis1', param.names[i]))) {
+      
+      if (is.na(pmatch('nug', param.names[i]))) {
+        iModel = as.numeric(substr(param.names[i], nchar(param.names[i]), nchar(param.names[i])))
+        
+        # Initialise some extra variogram parameters in case they're not used.        
+        if (iModel>nModel) {
+          variogram.kappa[iModel] = 0.5
+          variogram.ang1[iModel] = 0.0
+          variogram.anis1[iModel] = 1.0
+        }
+      }
+      
+      if (!is.na(pmatch('nug', param.names[i]))) {
+        variogram.nug = max(0, params[i])
+        
+      } else if (!is.na(pmatch('psill', param.names[i]))) {
+        variogram.psill[iModel] = max(0, params[i])
+        
+        if (iModel > nModel)
+          nModel = iModel
+        
+      } else if (!is.na(pmatch('range', param.names[i]))) {
+        variogram.range[iModel] = max(1e-6, params[i])
+        
+        if (iModel > nModel)
+          nModel = iModel
+        
+      } else if (!is.na(pmatch('kappa', param.names[i]))) {
+        variogram.kappa[iModel] = max(0, params[i])
+        
+        if (iModel > nModel)
+          nModel = iModel
+        
+      } else if (!is.na(pmatch('ang1', param.names[i]))) {
+        variogram.ang1[iModel] = params[i]
+        
+        if (iModel > nModel)
+          nModel = iModel
+        
+      } else if (!is.na(pmatch('anis1', param.names[i]))) {
+        variogram.anis1[iModel] = params[i]
+        
+        if (iModel > nModel)
+          nModel = iModel
+        
+      }
+    }
+  }
+  
+  # Build variogram
+  for (i in 1:nModel) {
+    
+    #message(paste('... DEBGUGGING: i, model type, nug, psill, range, kappa ',i, model[i], variogram.nug, variogram.psill[i], variogram.range[i], variogram.kappa[i]))
+    
+    if (i == 1 && !is.null(variogram.nug)) {
+      vgm.model = vgm(
+        nugget = variogram.nug,
+        psill = variogram.psill[i],
+        range = variogram.range[i],
+        model = modelType[i],
+        kappa = variogram.kappa[i],
+        anis = c(variogram.ang1, variogram.anis1)
+      )
+    } else if (i == 1) {
+      vgm.model = vgm(
+        psill = variogram.psill[i],
+        range = variogram.range[i],
+        model = modelType[i],
+        kappa = variogram.kappa[i],
+        anis = c(variogram.ang1, variogram.anis1)
+      )
+    } else {
+      vgm.model = vgm(
+        psill = variogram.psill[i],
+        range = variogram.range[i],
+        model = modelType[i],
+        add.to = vgm.model,
+        kappa = variogram.kappa[i],
+        anis = c(variogram.ang1, variogram.anis1)
+      )
+    }
+  }
+  return(vgm.model)
+}
diff --git a/.Rproj.user/19E94E50/sources/session-4dda64d8/lock_file b/.Rproj.user/19E94E50/sources/session-4dda64d8/lock_file
new file mode 100644
index 0000000..e69de29
diff --git a/.Rproj.user/shared/notebooks/paths b/.Rproj.user/shared/notebooks/paths
index 3ef8679..0a6bb7c 100644
--- a/.Rproj.user/shared/notebooks/paths
+++ b/.Rproj.user/shared/notebooks/paths
@@ -1 +1,6 @@
+C:/Users/tpet0008/Documents/HydroMap/R/get.allData.R="F8B21746"
+C:/Users/tpet0008/Documents/HydroMap/R/get.objFunc.R="06297799"
 C:/Users/tpet0008/Documents/HydroMap/R/get.pointData.R="90643B0A"
+C:/Users/tpet0008/Documents/HydroMap/R/get.variogram.R="85C5643A"
+C:/Users/tpet0008/Documents/HydroMap/R/get.variogramModel.R="C55B83BF"
+C:/Users/tpet0008/Documents/HydroMap/R/set.env.R="02592D8E"
diff --git a/DESCRIPTION b/DESCRIPTION
index 6eeb069..84feeb4 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: HydroMap
 Title: HydroMap
-Version: 0.2
+Version: 0.2.1.0
 Authors@R: person("Tim", "Peterson", email = "timjp@unimelb.edu.au", role = c("aut", "cre"))
 Description: Interpolates groundwater level observations giving realistic flow lines.
 Depends: R (>= 3.2.2)
diff --git a/R/get.objFunc.R b/R/get.objFunc.R
index 62c608c..d8f8c69 100644
--- a/R/get.objFunc.R
+++ b/R/get.objFunc.R
@@ -121,140 +121,25 @@ get.objFunc <- function(
       return(Inf)
     } 
   }
-    
-  # Build variogram model
-  if (fit.variogram.type == 1) {
   
-    # Check the format of the model input.
-    if (class(model)[1] == "variogramModel") {
-
-      # Get model types
-      variogram.type = as.character(model$model)
-      
-      # Remove 'nug' from list of model types and add to 'model' variable
-      filt = variogram.type != 'Nug'      
-      modelType = variogram.type[filt]            
-    } else if (is.character(model)) {
-      modelType = model    
-    } else {
-      stop('When fit.variogram.type==1, the model should be a gstat variogram model or a vector of model types, eg c(\'Mat\')')
-    }
-
-    # Initialise variogram settings
-    variogram.nugg = NULL
-    
-    variogram.psill = c()
-    
-    variogram.range = c()
-    
-    variogram.kappa = c()
-    
-    variogram.ang1 = c()
-    
-    variogram.anis1 = c()
+  # Build variogram model and turn off future fitting via get.variogram()
+  if (fit.variogram.type == 1) {
+    model = get.variogramModel(params, param.names, model)
+    fit.variogram.type = 3
     
-    nModel = 0
     
+    # Get variogram parameter names and values - for reporting
     # Extract variogram settings.
     for (i in 1:length(param.names)) {
-      
       if (!is.na(pmatch('psill', param.names[i])) ||
           !is.na(pmatch('range', param.names[i])) ||
           !is.na(pmatch('kappa', param.names[i])) ||
           !is.na(pmatch('nug', param.names[i]))  ||
           !is.na(pmatch('ang1', param.names[i])) ||
           !is.na(pmatch('anis1', param.names[i]))) {
-        
-        if (is.na(pmatch('nug', param.names[i]))) {
-          iModel = as.numeric(substr(param.names[i], nchar(param.names[i]), nchar(param.names[i])))
-  
-          # Initialise some extra variogram parameters in case they're not used.        
-          if (iModel>nModel) {
-            variogram.kappa[iModel] = 0.5
-            variogram.ang1[iModel] = 0.0
-            variogram.anis1[iModel] = 1.0
-          }
-          
-        }
-          
-        if (!is.na(pmatch('nug', param.names[i]))) {
-          variogram.nug = max(0, params[i])
-                     
-        } else if (!is.na(pmatch('psill', param.names[i]))) {
-          variogram.psill[iModel] = max(0, params[i])
-          
-          if (iModel > nModel)
-            nModel = iModel
-          
-        } else if (!is.na(pmatch('range', param.names[i]))) {
-          variogram.range[iModel] = max(1e-6, params[i])
-          
-          if (iModel > nModel)
-            nModel = iModel
-          
-        } else if (!is.na(pmatch('kappa', param.names[i]))) {
-          variogram.kappa[iModel] = max(0, params[i])
-          
-          if (iModel > nModel)
-            nModel = iModel
-          
-        } else if (!is.na(pmatch('ang1', param.names[i]))) {
-          variogram.ang1[iModel] = params[i]
-          
-          if (iModel > nModel)
-            nModel = iModel
-          
-        } else if (!is.na(pmatch('anis1', param.names[i]))) {
-          variogram.anis1[iModel] = params[i]
-          
-          if (iModel > nModel)
-            nModel = iModel
-          
-        }
-        
         params2Print[i] = params[i]
       }
     }
-    
-    # Build variogram
-    for (i in 1:nModel) {
-      
-      #message(paste('... DEBGUGGING: i, model type, nug, psill, range, kappa ',i, model[i], variogram.nug, variogram.psill[i], variogram.range[i], variogram.kappa[i]))
-      
-      if (i == 1 && !is.null(variogram.nug)) {
-        vgm.model = vgm(
-          nugget = variogram.nug,
-          psill = variogram.psill[i],
-          range = variogram.range[i],
-          model = modelType[i],
-          kappa = variogram.kappa[i],
-          anis = c(variogram.ang1, variogram.anis1)
-        )
-      } else if (i == 1) {
-        vgm.model = vgm(
-          psill = variogram.psill[i],
-          range = variogram.range[i],
-          model = modelType[i],
-          kappa = variogram.kappa[i],
-          anis = c(variogram.ang1, variogram.anis1)
-        )
-      } else {
-        vgm.model = vgm(
-          psill = variogram.psill[i],
-          range = variogram.range[i],
-          model = modelType[i],
-          add.to = vgm.model,
-          kappa = variogram.kappa[i],
-          anis = c(variogram.ang1, variogram.anis1)
-        )
-      }
-    }
-    model = vgm.model
-    
-    
-    # Turn off fitting within get.variogram()
-    fit.variogram.type = 3
-    
   }
 
   # Call krige.head with the input parameter names.
diff --git a/R/get.variogramModel.R b/R/get.variogramModel.R
new file mode 100644
index 0000000..b11af11
--- /dev/null
+++ b/R/get.variogramModel.R
@@ -0,0 +1,129 @@
+get.variogramModel <- function(
+  params, 
+  param.names, 
+  model=c('Mat'))  {
+
+  # DEBUGGING
+  op <- options("warn")
+  on.exit(options(op))
+  options(warn=1)  
+  
+  # Check if any parameters are nan or inf. If so, then later return obj function value inf
+  hasImplauibleParams = is.nan(params) | is.infinite(params)
+
+  # Check the format of the model input.
+  if (class(model)[1] == "variogramModel") {
+    
+    # Get model types
+    variogram.type = as.character(model$model)
+    
+    # Remove 'nug' from list of model types and add to 'model' variable
+    filt = variogram.type != 'Nug'      
+    modelType = variogram.type[filt]            
+  } else if (is.character(model)) {
+    modelType = model    
+  } else {
+    stop('When calling get.variogramModel(), fit.variogram.type must equal and the model should be a gstat variogram model or a vector of model types, eg c(\'Mat\')')
+  }
+  
+  # Initialise variogram settings
+  variogram.nugg = NULL
+  variogram.psill = c()
+  variogram.range = c()
+  variogram.kappa = c()
+  variogram.ang1 = c()
+  variogram.anis1 = c()
+  nModel = 0
+  
+  # Extract variogram settings.
+  for (i in 1:length(param.names)) {
+    if (!is.na(pmatch('psill', param.names[i])) ||
+        !is.na(pmatch('range', param.names[i])) ||
+        !is.na(pmatch('kappa', param.names[i])) ||
+        !is.na(pmatch('nug', param.names[i]))  ||
+        !is.na(pmatch('ang1', param.names[i])) ||
+        !is.na(pmatch('anis1', param.names[i]))) {
+      
+      if (is.na(pmatch('nug', param.names[i]))) {
+        iModel = as.numeric(substr(param.names[i], nchar(param.names[i]), nchar(param.names[i])))
+        
+        # Initialise some extra variogram parameters in case they're not used.        
+        if (iModel>nModel) {
+          variogram.kappa[iModel] = 0.5
+          variogram.ang1[iModel] = 0.0
+          variogram.anis1[iModel] = 1.0
+        }
+      }
+      
+      if (!is.na(pmatch('nug', param.names[i]))) {
+        variogram.nug = max(0, params[i])
+        
+      } else if (!is.na(pmatch('psill', param.names[i]))) {
+        variogram.psill[iModel] = max(0, params[i])
+        
+        if (iModel > nModel)
+          nModel = iModel
+        
+      } else if (!is.na(pmatch('range', param.names[i]))) {
+        variogram.range[iModel] = max(1e-6, params[i])
+        
+        if (iModel > nModel)
+          nModel = iModel
+        
+      } else if (!is.na(pmatch('kappa', param.names[i]))) {
+        variogram.kappa[iModel] = max(0, params[i])
+        
+        if (iModel > nModel)
+          nModel = iModel
+        
+      } else if (!is.na(pmatch('ang1', param.names[i]))) {
+        variogram.ang1[iModel] = params[i]
+        
+        if (iModel > nModel)
+          nModel = iModel
+        
+      } else if (!is.na(pmatch('anis1', param.names[i]))) {
+        variogram.anis1[iModel] = params[i]
+        
+        if (iModel > nModel)
+          nModel = iModel
+        
+      }
+    }
+  }
+  
+  # Build variogram
+  for (i in 1:nModel) {
+    
+    #message(paste('... DEBGUGGING: i, model type, nug, psill, range, kappa ',i, model[i], variogram.nug, variogram.psill[i], variogram.range[i], variogram.kappa[i]))
+    
+    if (i == 1 && !is.null(variogram.nug)) {
+      vgm.model = vgm(
+        nugget = variogram.nug,
+        psill = variogram.psill[i],
+        range = variogram.range[i],
+        model = modelType[i],
+        kappa = variogram.kappa[i],
+        anis = c(variogram.ang1, variogram.anis1)
+      )
+    } else if (i == 1) {
+      vgm.model = vgm(
+        psill = variogram.psill[i],
+        range = variogram.range[i],
+        model = modelType[i],
+        kappa = variogram.kappa[i],
+        anis = c(variogram.ang1, variogram.anis1)
+      )
+    } else {
+      vgm.model = vgm(
+        psill = variogram.psill[i],
+        range = variogram.range[i],
+        model = modelType[i],
+        add.to = vgm.model,
+        kappa = variogram.kappa[i],
+        anis = c(variogram.ang1, variogram.anis1)
+      )
+    }
+  }
+  return(vgm.model)
+}
diff --git a/R/krige.head.R b/R/krige.head.R
index fdcd247..b5c38cc 100644
--- a/R/krige.head.R
+++ b/R/krige.head.R
@@ -10,18 +10,12 @@
 #' these features can be individually controlled by the user.
 #'
 #' Importantly, if the mapping parameters are not specified by the user, then this function estimates the parameters using
-#' a mixed data-type (i.e. real and integer parameters) split-sample maximum likelihood global optimisation. The optimisation by default
-#' includes the variogram parameters (e.g. range, sill and nugget) and the search parameters for local kriging (e.g. radius, minimum and
-#' maximum number of observations to use). Optimising these parameters is not common in kriging. It is done herein because trials for Victoria,
-#' Australia, showed that calibrating these parameters produced significantly lower cross-validation errors (i.e. the error in predicting the observations
-#' removed from the optimisation) compared to the standard approach of graphical estimation from an experimental variogram. The optimisation is
-#' numerically challenging and the following factors should be considered before use:
-#'
-#' \itemize{
-#'  \item{Optimisation of the parameters \code{mrvbf.pslope}, \code{mrvbf.ppctl} and \code{smooth.std} often required the creating of raster grids for every parameter combination. To ease the computation burden, these parameters should be treated as discrete, not continuous, numbers.}
-#'  \item{The optimisation package \code{rgeoud} is used herein. For control the optimisation process, consider directly using \code{\link{krige.head.calib}}.}
-#'  \item{Trials have established default calibration parameters and settings that were effective for Victoria, Australia. There is no guarantee they will be effective for other regions.}
-#' }
+#' a mixed data-type (i.e. real and integer parameters) split-sample maximum likelihood global optimisation. See \code{\link{krige.head.calib}} for details.
+#' 
+#' However, it is advised to separately calibrate the mapping parameters using \code{\link{krige.head.calib}}. The output from 
+#' the calibration can then be passed into \code{krige.head} via the input \code{calibration.results}. This will use the calibrated setting 
+#' and the training observed data set to map the head. To use all of the observed data, or alternate data, input \code{data}. 
+#' See \code{\link{krige.head.calib}} Examples for details.
 #'
 #' In using this function, the primary user decisions are:
 #' \itemize{
@@ -30,6 +24,9 @@
 #'  \item{The type of variogram model, defined by the input \code{model}}
 #'  }
 #'
+#' @param \code{calibration.results} List variable returned from \code{\link{krige.head.calib}}. Inputting this applies the calibration results here. All other inputs are ignores, except \code{data}, \code{newdata} and \code{grid}. 
+#' Defauli is \code{NULL}.
+#' 
 #' @param \code{formula} defines the R formula (as a character or formula data type) to be used to interpolate the heads. The left hand side of the formula must be \code{head}.
 #' The right hand side can contain any or all of the following terms: \code{elev} for the land surface elevation;
 #' \code{MrVBF} for the Multiresolution Index of Valley Bottom Flatness as a measure of valley-ness at each DEM grd cell;
@@ -165,17 +162,10 @@
 #' 
 #' # Export grids to an ARCMAP *.asc file
 #' write.asciigrid(heads,'DBNS.asc',3,na.value = -999)
-#' 
-#' # Re-calibrate the parameters and map using the default settings.
-#' # Note, here only those parameters that are a vector are calibrated. This is in addition
-#' # to the variogram parameters. 
-#' varigram.model = vgm(psill=10, model='Mat', range= 10000 , nugget=1, kappa=0.1);
-#' heads <- krige.head(formula=f, grid=DEM, data=obs.data, data.errvar.colname='total_err_var',
-#' nmax=Inf, maxdist=Inf, smooth.std = seq(0.5,1.5,length.out=11), 
-#' model = 'Mat',  fit.variogram.type=1, debug.level=1)
 #'
 #' @export
 krige.head <- function(
+  calibration.results = NULL,  
   formula = as.formula("head ~ elev + MrVBF + MrRTF + smoothing"),
   grid=NULL,
   grid.landtype.colname = NULL,
@@ -231,6 +221,32 @@ krige.head <- function(
   if (!do.grid.est && nsim!=0)
     stop('Point estimation cannot be undertaken for geostatistical simulations. To undertake gridded simulations, set newdata=NULL and nsim>0.')
 
+  # Use calibration results if provided. Else use individual input.
+  if (!is.null(calibration.results)) {
+    
+    formula = calibration.results$inputs$formula
+    if (is.null(grid))
+      grid = calibration.results$inputs$grid
+    grid.landtype.colname = calibration.results$inputs$grid.landtype.colname;
+    if (is.null(data))
+      data = calibration.results$inputs$data
+    data.fixedHead = calibration.results$inputs$data.fixedHead
+    if (is.null(newdata))
+      newdata = calibration.results$inputs$newdata
+    data.errvar.colname = calibration.results$inputsdata.errvar.colname
+    model= calibration.results$variogramModel
+    fit.variogram.type = 3; 
+    
+    mrvbf.pslope = calibration.results$parameter.Values$mrvbf.pslope
+    mrvbf.ppctl = calibration.results$parameter.Values$mrvbf.ppctl
+    smooth.std = calibration.results$parameter.Values$smooth.std
+    nmax = calibration.results$parameter.Values$nmax
+    nmax.fixedHead = calibration.results$parameter.Values$nmax.fixedHead
+    maxdist = calibration.results$parameter.Values$maxdist
+    nmin = calibration.results$parameter.Values$nmin 
+    omax = calibration.results$parameter.Values$omax
+  }
+  
   # Get variable names and check if using MrVBF or MrRTF or DEM
   if (debug.level>0)
     message('... Checking formula terms.');
diff --git a/R/krige.head.calib.R b/R/krige.head.calib.R
index 23c7ebe..995198e 100644
--- a/R/krige.head.calib.R
+++ b/R/krige.head.calib.R
@@ -3,10 +3,8 @@
 #' \code{krige.head.calib} calibrates parameters to minimise the interpolation error.
 #'
 #' @details This function optimises the parameters for the spatial interpolation so as to minimise the
-#' interpolation error - which by default is defined using a formal likelihood function. Once this function has return
-#' the optimal parameter set, then the mapping can be undertaken using \code{\link{krige.head}}.
-#' Importantly, the mapping approach used for the calibration (e.g. the formula etc) must be identical to that
-#' input to \code{\link{krige.head}}.
+#' interpolation error, which by default is defined using a formal likelihood function. Once this function has return
+#' the optimal parameter set, the mapping can be undertaken using \code{\link{krige.head}} (see examples below).
 #'
 #' The mapping parameters are estimated using a mixed data-type (i.e. real and integer parameters) split-sample maximum likelihood global optimisation. The optimisation by default
 #' includes the variogram parameters (e.g. range, sill and nuggest) and the search parameters for local kriging (e.g. radius, minimum and
@@ -22,8 +20,8 @@
 #'  \item{The available input point data must be split up into \code{data} and \code{newdata}. The point observations within \code{data} are used to estimate the water level at the locations defined
 #'  within \code{newdata}. The difference between the observed and estimated values defines the interpolation prediction error. This process also gives the kriging variance at
 #'  each \code{newdata} location. The points used for \code{data} and \code{newdata} should both cover all types of terrain, geology, landuse and the full mapping extend so as to avoid bias. Also, changing either may
-#'  change the calibration solution. }
-#'  \item{The \code{objFunc.type} allows for maximum likelihood estimation - that is optimisation that accounts for the expected kriging error. See \code{objFunc.type} below for details.
+#'  change the calibration solution.}
+#'  \item{The \code{objFunc.type} allows for maximum likelihood estimation - that is optimisation that accounts for the expected kriging error. See \code{objFunc.type} below for details.}
 #' }
 #'
 #' In using this function, the primary user decisions are:
@@ -131,7 +129,7 @@
 #'
 #' # Load a model variogram and mapping parameters found to be effective.
 #' data('mapping.parameters')
-#
+#'
 #' # Define a simple kriging formula without MrVBF terms that does not require the package RSAGA.
 #' f <- as.formula('head ~ elev + smoothing')
 #'
@@ -140,33 +138,17 @@
 #'
 #' # Calibrate the mapping parameters with 25% of the data randomly selected at 2 cores.
 #' calib.results <- krige.head.calib(formula=f, grid=DEM, data=obs.data, newdata=0.25, nmin=0, nmax=Inf, maxdist=Inf, omax=0,
-#'    data.errvar.colname='total_err_var', model = variogram.model,  fit.variogram.type=1,
-#'    debug.level=0, use.cluster = 2)
-#'
-#' # Reformat the parameter to a named list. This is just for ease of use and the need for
-#' this may be removed in future versions.
-#' params.all = c(mrvbf.pslope=NA,mrvbf.ppctl=NA,smooth.std=NA,nmax=NA,maxdist=NA,
-#'                kappa_model_1 = NA, ang1_model_1 = NA, anis1_model_1 = NA,
-#'                nmax.fixedHead = NA);
-#' for (i in 1:length(head.calib$param.Names)) {
-#'  params.all[[calib.results$param.Names[i]]] = calib.results$par[i]
-#' }
-#' params.all.names=names(params.all);
-#'
-#' # Rebuild the variogram with the calibrated parameters.
-#' varigram.model = vgm(psill=params.all[['psill_model_1']], model='Mat',
-#'                    range= params.all[['range_model_1']] , nugget=params.all[['nug']],
-#'                    kappa=params.all[['kappa_model_1']] );
+#'                       data.errvar.colname='total_err_var', model = variogram.model,  fit.variogram.type=1,
+#'                       pop.size.multiplier=1, debug.level=0, use.cluster = 2)
 #'
-#' # Grid the observed head data.
-#' head.map <- krige.head(formula=f, grid=DEM, data=obs.data, data.errvar.colname='total_err_var',
-#' model=variogram.model, smooth.std=params.all[['smooth.std']], maxdist=params.all[['maxdist']],
-#' nmax=params.all[['nmax']], fit.variogram.type=3, debug.level=1)
+#' # Grid the observed head using ONLY the training data from the calibration (and all cores) and then map..
+#' head.grid.training <- krige.head(calibration.results = calib.results, use.cluster = T)
+#' spplot(head.grid.training)
+#' 
+#' # Grid the observed head using ALL of the obsrved data from the calibration (and all cores).
+#' head.grid.all <- krige.head(calibration.results = calib.results, data=obs.data, use.cluster = T)
+#' spplot(head.grid.all)
 #'
-#' # Write the grids to ARCMAP ASCII grids
-#' write.asciigrid(head.map, 'head.asc',1);
-#' write.asciigrid(head.map, 'KrigingVariance.asc',2);
-
 #' @export
 krige.head.calib <-
   function(formula = as.formula("head ~ elev + MrVBF + MrRTF + smoothing"),
@@ -197,6 +179,17 @@ krige.head.calib <-
     if (debug.level>0)
       message('Calibrating mapping parameters.')
 
+    # Store input parameter values (exluding variogram params). These are updated if they're calibrated.
+    solution = list()
+    solution$parameter.Values$mrvbf.pslope = mrvbf.pslope
+    solution$parameter.Values$mrvbf.ppctl = mrvbf.ppctl
+    solution$parameter.Values$smooth.std = smooth.std
+    solution$parameter.Values$nmax = nmax
+    solution$parameter.Values$nmax.fixedHead = nmax.fixedHead
+    solution$parameter.Values$maxdist = maxdist
+    solution$parameter.Values$nmin = nmin
+    solution$parameter.Values$omax = omax
+    
     # Assess formula componants
     var.names = all.vars(as.formula(formula));
     use.MrVBF = any(match(var.names, 'MrVBF'));
@@ -757,7 +750,7 @@ krige.head.calib <-
     } else {
       message('... Starting continuous value calibration of the following parameters:',paste( param.Names,sep="",collapse=', '))
     }
-    solution = genoud(
+    genoud.solution = genoud(
       HydroMap:::get.objFunc,
       nvars = nParams,
       max = doObjMaximisation,
@@ -806,7 +799,7 @@ krige.head.calib <-
 
     # Store the returned value for use when calling get.objFunc directly.
     # This is required if some of the parameters have discrete parameters (which are used by the lookup table to transform to meaningful values)
-    solution$par.pretransformed = solution$par;
+    genoud.solution$par.pretransformed = genoud.solution$par;
 
     # Back transform parameters
     lookup.table.names = names(param.IntegerLookupTable)
@@ -817,26 +810,79 @@ krige.head.calib <-
 
       # Get non-transformed value from table.
       if (isIntegerParam) {
-        solution$par[i] = param.IntegerLookupTable[[param.Names[i]]][solution$par[i], 2]
+        genoud.solution$param.Values[i] = param.IntegerLookupTable[[param.Names[i]]][genoud.solution$par[i], 2]
       }
     }
 
-
     # Print summary
     message(paste('    Best objective function value:', solution$value))
     for (i in 1:nParams)
-      message(paste('    Parameter ', param.Names[i], ' = ', solution$par[i]))
-
-    # Add parameter names to the solution list variable.
-    solution$param.Names = param.Names
-
+      message(paste('    Parameter ', param.Names[i], ' = ', genoud.solution$param.Values[i]))
+    
+    # Add formula
+    solution$inputs$formula = formula;     
+    
+    # Add data and new data
+    solution$inputs$data = data
+    solution$inputs$newdata = newdata
+    solution$inputs$grid = grid
+    
+    # Add remaining inputs 
+    solution$inputs$grid.landtype.colname = grid.landtype.colname
+    solution$inputs$data.fixedHead = data.fixedHead
+    solution$inputs$data.errvar.colname = data.errvar.colname  
+
+    ind = match(param.Names,'mrvbf.pslope')
+    ind = na.omit(ind)
+    if (length(ind)==1 && is.finite(ind))
+      solution$parameter.Values$mrvbf.pslope = genoud.solution$param.Values[ind]
+    
+    ind = match(param.Names,'mrvbf.ppctl')
+    ind = na.omit(ind)
+    if (length(ind)==1 && is.finite(ind))
+      solution$parameter.Values$mrvbf.ppctl = genoud.solution$param.Values[ind]
+    
+    ind = match(param.Names,'smooth.std')
+    ind = na.omit(ind)
+    if (length(ind)==1 && is.finite(ind))
+      solution$parameter.Values$smooth.std = genoud.solution$param.Values[ind]
+    
+    ind = match(param.Names,'nmax')
+    ind = na.omit(ind)
+    if (length(ind)==1 && is.finite(ind))
+      solution$parameter.Values$nmax = genoud.solution$param.Values[ind]
+    
+    ind = match(param.Names,'nmax.fixedHead')
+    ind = na.omit(ind)
+    if (length(ind)==1 && is.finite(ind))
+      solution$parameter.Values$nmax.fixedHead = genoud.solution$param.Values[ind]
+    
+    ind = match(param.Names,'maxdist')
+    ind = na.omit(ind)
+    if (length(ind)==1 && is.finite(ind))
+      solution$parameter.Values$maxdist = genoud.solution$param.Values[ind]
+    
+    ind = match(param.Names,'nmin')
+    ind = na.omit(ind)
+    if (length(ind)==1 && is.finite(ind))
+      solution$parameter.Values$nmin = genoud.solution$param.Values[ind]
+    
+    ind = match(param.Names,'omax')
+    ind = na.omit(ind)
+    if (length(ind)==1 && is.finite(ind))
+      solution$parameter.Values$omax = genoud.solution$param.Values[ind]
+    
+    # Add variogram model
+    if (fit.variogram.type == 1)
+      solution$variogramModel = get.variogramModel(genoud.solution$par, param.Names, model)
+    
+    
     # Add the lookup table to the solution. This is required if the the objectiveget.objFunc() is to be directly called.
     solution$lookup.table = param.IntegerLookupTable;
-
-    # Add data and new data
-    solution$data = data
-    solution$newdata = newdata
-
+    
+    # Add calibration settings
+    solution$genoud.output = genoud.solution
+    
     return(solution)
 
   }
diff --git a/man/krige.head.Rd b/man/krige.head.Rd
index 00e812f..0400e79 100644
--- a/man/krige.head.Rd
+++ b/man/krige.head.Rd
@@ -5,6 +5,7 @@
 \title{Spatially interpolates sparse groundwater level observations, and if required, estimate mapping parameters.}
 \usage{
 krige.head(
+  calibration.results = NULL,
   formula = as.formula("head ~ elev + MrVBF + MrRTF + smoothing"),
   grid = NULL,
   grid.landtype.colname = NULL,
@@ -59,6 +60,9 @@ krige.head(
 )
 }
 \arguments{
+\item{\code{calibration.results}}{List variable returned from \code{\link{krige.head.calib}}. Inputting this applies the calibration results here. All other inputs are ignores, except \code{data}, \code{newdata} and \code{grid}. 
+Defauli is \code{NULL}.}
+
 \item{\code{formula}}{defines the R formula (as a character or formula data type) to be used to interpolate the heads. The left hand side of the formula must be \code{head}.
 The right hand side can contain any or all of the following terms: \code{elev} for the land surface elevation;
 \code{MrVBF} for the Multiresolution Index of Valley Bottom Flatness as a measure of valley-ness at each DEM grd cell;
@@ -152,18 +156,12 @@ the inclusion of categorical land types and fixed head boundary conditions, such
 these features can be individually controlled by the user.
 
 Importantly, if the mapping parameters are not specified by the user, then this function estimates the parameters using
-a mixed data-type (i.e. real and integer parameters) split-sample maximum likelihood global optimisation. The optimisation by default
-includes the variogram parameters (e.g. range, sill and nugget) and the search parameters for local kriging (e.g. radius, minimum and
-maximum number of observations to use). Optimising these parameters is not common in kriging. It is done herein because trials for Victoria,
-Australia, showed that calibrating these parameters produced significantly lower cross-validation errors (i.e. the error in predicting the observations
-removed from the optimisation) compared to the standard approach of graphical estimation from an experimental variogram. The optimisation is
-numerically challenging and the following factors should be considered before use:
+a mixed data-type (i.e. real and integer parameters) split-sample maximum likelihood global optimisation. See \code{\link{krige.head.calib}} for details.
 
-\itemize{
- \item{Optimisation of the parameters \code{mrvbf.pslope}, \code{mrvbf.ppctl} and \code{smooth.std} often required the creating of raster grids for every parameter combination. To ease the computation burden, these parameters should be treated as discrete, not continuous, numbers.}
- \item{The optimisation package \code{rgeoud} is used herein. For control the optimisation process, consider directly using \code{\link{krige.head.calib}}.}
- \item{Trials have established default calibration parameters and settings that were effective for Victoria, Australia. There is no guarantee they will be effective for other regions.}
-}
+However, it is advised to separately calibrate the mapping parameters using \code{\link{krige.head.calib}}. The output from 
+the calibration can then be passed into \code{krige.head} via the input \code{calibration.results}. This will use the calibrated setting 
+and the training observed data set to map the head. To use all of the observed data, or alternate data, input \code{data}. 
+See \code{\link{krige.head.calib}} Examples for details.
 
 In using this function, the primary user decisions are:
 \itemize{
@@ -219,14 +217,6 @@ spplot(heads,3, scales = list(draw = TRUE))
 # Export grids to an ARCMAP *.asc file
 write.asciigrid(heads,'DBNS.asc',3,na.value = -999)
 
-# Re-calibrate the parameters and map using the default settings.
-# Note, here only those parameters that are a vector are calibrated. This is in addition
-# to the variogram parameters. 
-varigram.model = vgm(psill=10, model='Mat', range= 10000 , nugget=1, kappa=0.1);
-heads <- krige.head(formula=f, grid=DEM, data=obs.data, data.errvar.colname='total_err_var',
-nmax=Inf, maxdist=Inf, smooth.std = seq(0.5,1.5,length.out=11), 
-model = 'Mat',  fit.variogram.type=1, debug.level=1)
-
 }
 \references{
 Gallant, J.C., Dowling, T.I. (2003): 'A multiresolution index of valley bottom flatness for mapping depositional areas', Water Resources Research, 39/12:1347-1359
diff --git a/man/krige.head.calib.Rd b/man/krige.head.calib.Rd
index 3943fe5..f789e1c 100644
--- a/man/krige.head.calib.Rd
+++ b/man/krige.head.calib.Rd
@@ -146,7 +146,34 @@ As per \code{genoud} plus the parameter values in the transformed calibration sp
 \code{krige.head.calib} calibrates parameters to minimise the interpolation error.
 }
 \details{
+This function optimises the parameters for the spatial interpolation so as to minimise the
+interpolation error, which by default is defined using a formal likelihood function. Once this function has return
+the optimal parameter set, the mapping can be undertaken using \code{\link{krige.head}} (see examples below).
+
+The mapping parameters are estimated using a mixed data-type (i.e. real and integer parameters) split-sample maximum likelihood global optimisation. The optimisation by default
+includes the variogram parameters (e.g. range, sill and nuggest) and the search parameters for local kriging (e.g. radius, minimum and
+maximum number of observations to use). Optimising these parameters is not common in kriging. It is done herein because trials for Victoria,
+Australia, showed that calibrating these parameters produced significantly lower cross-validation errors (i.e. the error in predicting the observations
+removed from the optimisation) compared to the standard approach using fitting to an experimental variogram. The optimisation is
+numerically challenging and the following factors should be considered before use:
+
+\itemize{
+ \item{Optimisation of the parameters \code{mrvbf.pslope}, \code{mrvbf.ppctl} and \code{smooth.std} often required the creating of raster grids for every parameter combination. To ease the computation burden, these parameters should be treated as disrcete, not continuous, numbers.}
+ \item{The optimisation package \code{rgeoud} is used herein. See the \code{rgeoud} documentation for details of the optimisation scheme.}
+ \item{Trials have established default calibration parameters and settings that were effective for Victoria, Australia. There is no guarantee they will be effective for other regions.}
+ \item{The available input point data must be split up into \code{data} and \code{newdata}. The point observations within \code{data} are used to estimate the water level at the locations defined
+ within \code{newdata}. The difference between the observed and estimated values defines the interpolation prediction error. This process also gives the kriging variance at
+ each \code{newdata} location. The points used for \code{data} and \code{newdata} should both cover all types of terrain, geology, landuse and the full mapping extend so as to avoid bias. Also, changing either may
+ change the calibration solution.}
+ \item{The \code{objFunc.type} allows for maximum likelihood estimation - that is optimisation that accounts for the expected kriging error. See \code{objFunc.type} below for details.}
+}
 
+In using this function, the primary user decisions are:
+\itemize{
+ \item{The kriging with external drift formula defining the independent gridded variables deemed to predict the groundwater elevation. See the input \code{formula}.}
+ \item{The type of the objective function, defined by the input \code{objFunc.type}}
+ \item{The type of variogram model, defined by the input \code{model}}
+ }
 }
 \examples{
 # Load packages in case they have not loaded.
@@ -167,6 +194,7 @@ data('victoria.groundwater')
 
 # Load a model variogram and mapping parameters found to be effective.
 data('mapping.parameters')
+
 # Define a simple kriging formula without MrVBF terms that does not require the package RSAGA.
 f <- as.formula('head ~ elev + smoothing')
 
@@ -175,32 +203,17 @@ varigram.model <- vgm(psill=25, model="Mat", range=5000, nugget=5, kappa = 0.5)
 
 # Calibrate the mapping parameters with 25\% of the data randomly selected at 2 cores.
 calib.results <- krige.head.calib(formula=f, grid=DEM, data=obs.data, newdata=0.25, nmin=0, nmax=Inf, maxdist=Inf, omax=0,
-   data.errvar.colname='total_err_var', model = variogram.model,  fit.variogram.type=1,
-   debug.level=0, use.cluster = 2)
-
-# Reformat the parameter to a named list. This is just for ease of use and the need for
-this may be removed in future versions.
-params.all = c(mrvbf.pslope=NA,mrvbf.ppctl=NA,smooth.std=NA,nmax=NA,maxdist=NA,
-               kappa_model_1 = NA, ang1_model_1 = NA, anis1_model_1 = NA,
-               nmax.fixedHead = NA);
-for (i in 1:length(head.calib$param.Names)) {
- params.all[[calib.results$param.Names[i]]] = calib.results$par[i]
-}
-params.all.names=names(params.all);
+                      data.errvar.colname='total_err_var', model = variogram.model,  fit.variogram.type=1,
+                      pop.size.multiplier=1, debug.level=0, use.cluster = 2)
 
-# Rebuild the variogram with the calibrated parameters.
-varigram.model = vgm(psill=params.all[['psill_model_1']], model='Mat',
-                   range= params.all[['range_model_1']] , nugget=params.all[['nug']],
-                   kappa=params.all[['kappa_model_1']] );
+# Grid the observed head using ONLY the training data from the calibration (and all cores) and then map..
+head.grid.training <- krige.head(calibration.results = calib.results, use.cluster = T)
+spplot(head.grid.training)
 
-# Grid the observed head data.
-head.map <- krige.head(formula=f, grid=DEM, data=obs.data, data.errvar.colname='total_err_var',
-model=variogram.model, smooth.std=params.all[['smooth.std']], maxdist=params.all[['maxdist']],
-nmax=params.all[['nmax']], fit.variogram.type=3, debug.level=1)
+# Grid the observed head using ALL of the obsrved data from the calibration (and all cores).
+head.grid.all <- krige.head(calibration.results = calib.results, data=obs.data, use.cluster = T)
+spplot(head.grid.all)
 
-# Write the grids to ARCMAP ASCII grids
-write.asciigrid(head.map, 'head.asc',1);
-write.asciigrid(head.map, 'KrigingVariance.asc',2);
 }
 \references{
 Gallant, J.C., Dowling, T.I. (2003) A multiresolution index of valley bottom flatness for mapping depositional areas, Water Resources Research, 39/12:1347-1359