01_OOPS_1yr_param_calculation_2000_2016.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Script to calculate optimal interpolation parameters\n",
    "\n",
    "* The script loops over the variable names in _varname_\n",
    "* The script calculates the optimal signal-to-noise (e) and correlation length (l) parameters for the (full) years 2000 to 2016 by cross validation.\n",
    "* These optimal parameters are stored in\n",
    "    * /_varname_/*varname*_newe.txt\n",
    "    * /_varname_/*varname*_newl.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "BLAS.set_num_threads(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Julia Version 0.6.4\n",
      "Commit 9d11f62bcb (2018-07-09 19:09 UTC)\n",
      "Platform Info:\n",
      "  OS: Linux (x86_64-pc-linux-gnu)\n",
      "  CPU: Intel(R) Xeon(R) CPU E5-2630L 0 @ 2.00GHz\n",
      "  WORD_SIZE: "
     ]
    }
   ],
   "source": [
    "versioninfo()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "64\n",
      "  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Sandybridge MAX_THREADS=16)\n",
      "  LAPACK: libopenblas64_\n",
      "  LIBM: libopenlibm\n",
      "  LLVM: libLLVM-3.9.1 (ORCJIT, sandybridge)\n"
     ]
    }
   ],
   "source": [
    "print(join([\"$p: $v\\n\" for (p,v) in Pkg.installed()]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"/home/lennerts/Diva_product_scripts_2018\""
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pwd()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "cd(\"//data/20180306_OOPS/\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"/data/20180306_OOPS\""
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pwd()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Libz: 0.4.0\n",
      "AxisAlgorithms: 0.3.0\n",
      "OffsetArrays: 0.6.0\n",
      "HTTP: 0.6.12\n",
      "CondaBinDeps: 0.1.0\n",
      "OceanPlot: 0.0.0-\n",
      "Nullables: 0.0.6\n",
      "ZMQ: 0.6.3\n",
      "DataStructures: 0.8.4\n",
      "Compat: 1.0.0\n",
      "EzXML: 0.7.1\n",
      "ShowItLikeYouBuildIt: 0.2.0\n",
      "SpecialFunctions: 0.6.0\n",
      "Blosc: 0.5.0\n",
      "ZipFile: 0.6.0\n",
      "FixedPointNumbers: 0.4.6\n",
      "SHA: 0.5.7\n",
      "Missings: 0.2.10\n",
      "MAT: 0.4.0\n",
      "HDF5: 0.9.3\n",
      "ColorTypes: 0.6.7\n",
      "BufferedStreams: 0.4.1\n",
      "MbedTLS: 0.5.12\n",
      "SortingAlgorithms: 0.2.1\n",
      "Conda: 1.0.0\n",
      "PyCall: 1.17.1\n",
      "WoodburyMatrices: 0.3.0\n",
      "JSON: 0.17.2\n",
      "StatsBase: 0.23.1\n",
      "IJulia: 1.9.1\n",
      "PyPlot: 2.6.0\n",
      "BinDeps: 0.8.8\n",
      "Parameters: 0.9.1\n",
      "DIVAnd: 2.0.0+\n",
      "Mustache: 0.3.3\n",
      "Reexport: 0.1.0\n",
      "CMakeWrapper: 0.1.0\n",
      "URIParser: 0.3.1\n",
      "Interpolations: 0.8.0\n",
      "BinaryProvider: 0.3.3\n",
      "IniFile: 0.4.0\n",
      "LaTeXStrings: 1.0.0\n",
      "DataArrays: 0.7.0\n",
      "VersionParsing: 1.1.1\n",
      "Colors: 0.8.2\n",
      "MacroTools: 0.4.2\n",
      "NCDatasets: 0.3.0\n",
      "Ratios: 0.2.0\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mRecompiling stale cache file /home/lennerts/.julia/lib/v0.6/HTTP.ji for module HTTP.\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mRecompiling stale cache file /home/lennerts/.julia/lib/v0.6/PyPlot.ji for module PyPlot.\n",
      "\u001b[39m"
     ]
    }
   ],
   "source": [
    "using DIVAnd\n",
    "using PyPlot\n",
    "using NCDatasets\n",
    "using DataStructures\n",
    "\n",
    "# Load a more efficient version of sparse matrix multiplication\n",
    "include(joinpath(Pkg.dir(\"DIVAnd\"),\"src\",\"override_ssmult.jl\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### yearlists_ function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "yearlists_ (generic function with 1 method)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function yearlists_(dataset_range, total_window_yrs)\n",
    "    # dataset_range = 2000:2012\n",
    "    # total_window_yrs = 10\n",
    "    # will return: [2000:2009, 2001:2010, 2002:2011, 2003:2012]\n",
    "    \n",
    "    n_windows = length(dataset_range) - total_window_yrs + 1\n",
    "\n",
    "    a = Array{UnitRange{Int64}}(n_windows)\n",
    "\n",
    "    for i = 1:n_windows\n",
    "        a[i] =  dataset_range[i]:(dataset_range[i] + total_window_yrs -1)\n",
    "    end\n",
    "    return(a)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "plotres (generic function with 1 method)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "function plotres(timeindex,sel,fit,erri) # fit = interpolated field\n",
    "    tmp = copy(fit)\n",
    "    tmp[erri .> .3] = NaN; # only plotting where error < 0.3\n",
    "    figure(figsize = (10,8))\n",
    "    subplot(2,1,1)\n",
    "    title(\"$(timeindex) - surface\")\n",
    "    \n",
    "    # select the data near the surface\n",
    "    selsurface = sel .& (depth .< 5)\n",
    "    vmin = minimum(value[selsurface])\n",
    "    vmax = maximum(value[selsurface])\n",
    "    \n",
    "    # plot the data\n",
    "    scatter(lon[selsurface],lat[selsurface],10,value[selsurface];\n",
    "            vmin = vmin, vmax = vmax)\n",
    "    xlim(minimum(lonr),maximum(lonr))\n",
    "    ylim(minimum(latr),maximum(latr))\n",
    "    colorbar()\n",
    "    contourf(bx,by,b', levels = [-1e5,0],colors = [[.5,.5,.5]])\n",
    "    \n",
    "    # plot the analysis\n",
    "    subplot(2,1,2)\n",
    "    pcolor(lonr,latr,tmp[:,:,1]';\n",
    "           vmin = vmin, vmax = vmax)\n",
    "    colorbar()\n",
    "    contourf(bx,by,b', levels = [-1e5,0],colors = [[.5,.5,.5]])\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "download bathymetry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mBathymetry file already downloaded\n",
      "\u001b[39m"
     ]
    }
   ],
   "source": [
    "bathname = \"gebco_30sec_4.nc\"\n",
    "\n",
    "if !isfile(bathname)\n",
    "    download(\"https://b2drop.eudat.eu/s/ACcxUEZZi6a4ziR/download\",bathname)\n",
    "else\n",
    "    info(\"Bathymetry file already downloaded\")\n",
    "end\n",
    "\n",
    "bathisglobal = true"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "resolution of grid and mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-75.0:0.1:20.0, 35.0:0.1:75.0, [2475.89 1803.55 … 290.672 483.172; 2698.02 2466.61 … 291.781 466.875; … ; 3013.61 2988.84 … 173.094 126.891; 2998.63 2965.13 … 146.625 140.734])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dx = dy = 0.1\n",
    "lonr = -75.:dx:20.; # the range of longitudes (start:step:end)\n",
    "latr = 35.:dy:75.; # the range of latitudes (start:step:end)\n",
    "depthr = [-1., 1.] # put . always! (otherwise integer and error in DIVA)\n",
    "\n",
    "# create mask\n",
    "# mask,(pm,pn),(xi,yi) = divand_rectdom(lonr,latr)\n",
    "# I changed the line below after error with diva3d\n",
    "mask,(pm,pn),(xi,yi) = DIVAnd.domain(bathname, true, lonr, latr)\n",
    "bx,by,b = load_bath(bathname,true,lonr,latr)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optimisation loop for:\n",
    "* different parameters (varname)\n",
    "* different year ranges (yearlist -> yr)\n",
    "* different seasons (monthlists -> mo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "epsilon2 = 2 # 0.5^-1 (old value: 0.5);\n",
    "len = 150_000 # IN  METERS\n",
    "\n",
    "# Months and years\n",
    "monthlists = [[1,2,3,4,5,6,7,8,9,10,11,12]]\n",
    "yearlists = yearlists_(2000:2016, 1);\n",
    "\n",
    "for varname = [\"log_tem_lon\",\"log_cal_hel\", \"log_tot_lar\",\"log_tot_sma\",\n",
    "                \"log_oit_spp\", \"log_cal_fin\",\"ratio_large_to_small\",\n",
    "                \"log_aca_spp\", \"log_met_luc\", \"log_chli\"]\n",
    "    \n",
    "    i = 1\n",
    "    newl = ones(length(yearlists)*length(monthlists))\n",
    "    newe = ones(length(yearlists)*length(monthlists))\n",
    "    \n",
    "    varDir = joinpath(pwd(),varname)\n",
    "    mkpath(varDir)\n",
    "    \n",
    "    @show varname\n",
    "    t1 = now()\n",
    "    @show t1\n",
    "    \n",
    "    # read data\n",
    "    fname = \"bigfile_$(varname).txt\"\n",
    "\n",
    "    if !isfile(fname)\n",
    "        error(\"ERROR: File not found\")\n",
    "    else\n",
    "        info(\"File found -> OK!\")\n",
    "    end\n",
    "\n",
    "    value,lon,lat,depth,time,ids = loadbigfile(fname)\n",
    "    checkobs((lon,lat,depth,time), value, ids)\n",
    "    \n",
    "    #### make selection for each timestep\n",
    "    for yr = yearlists\n",
    "        for mo = monthlists\n",
    "                print(\"start $(i)  \")\n",
    "                TS = TimeSelectorYearListMonthList([yr],[mo])\n",
    "                @show TS\n",
    "                t1 = now() # keep track of time\n",
    "                @show t1\n",
    "\n",
    "                # subset data based on timestep\n",
    "                # https://stackoverflow.com/a/29661623\n",
    "                sel = ones(length(value))\n",
    "                sel = (depth .< 1) .& (indexin(Dates.month.(time), mo) .> 0) .& (indexin(Dates.year.(time), yr) .> 0)\n",
    "\n",
    "                value_sel = value[sel]\n",
    "                lon_sel = lon[sel]\n",
    "                lat_sel = lat[sel]\n",
    "                depth_sel = depth[sel]\n",
    "                time_sel = time[sel]\n",
    "                ids_sel = ids[sel];\n",
    "                @show (size(value_sel))\n",
    "\n",
    "\n",
    "                #### run optimization for noise-to-signal ratio\n",
    "                # optimal parameter selection with mode '0' = default\n",
    "                @time bestfactorl,bestfactore, cvval,cvvalues, x2Ddata,y2Ddata,cvinter,xi2D,yi2D =\n",
    "                DIVAnd_cv(mask,(pm,pn),(xi,yi),(lon_sel,lat_sel),value_sel-mean(value_sel),len,epsilon2,2,3,0);\n",
    "\n",
    "                @show newl[i]=len*bestfactorl\n",
    "                @show newe[i]=epsilon2*bestfactore\n",
    "            \n",
    "                open(joinpath(varDir,\"$(varname)_newe.txt\"), \"a\") do io\n",
    "                    writedlm(io, newe[i])\n",
    "                end\n",
    "            \n",
    "                open(joinpath(varDir,\"$(varname)_newl.txt\"), \"a\") do io\n",
    "                    writedlm(io, newl[i])\n",
    "                end\n",
    "\n",
    "                i = i+1\n",
    "\n",
    "        end # mo\n",
    "    end # yr\n",
    "\n",
    "    #writedlm(\"$(varname)_newe.txt\", newe)\n",
    "    #writedlm(\"$(varname)_newl.txt\", newl)\n",
    "\n",
    "end # varname"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.4",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}