#include "DataAdapters/ncfilealt.h"
#include "mdatetime.h"
#include "simple2ddata.h"
#include "specfunc.h"
#include "vartype.h"
#include <regex.h>
#include <stdlib.h>
#include <sys/types.h>

#if !defined(M__NEMO)
#define M__NEMO

using michlib::Ceil;
using michlib::Floor;
using michlib::int2;
using michlib::MDateTime;
using michlib::NCFileA;
using michlib::ToGeoDomainNeg;

class NEMOData
{
 enum Type
 {
  TYPE_UNKNOWN,
  TYPE_DT,
  TYPE_NRT,
  TYPE_NRT6
 };

 NCFileA                nc, nct, ncs;
 size_t                 xb, yb, xe, ye, layer;
 std::vector<real>      depths;
 std::vector<MDateTime> times;
 size_t                 nx;
 Type                   type = TYPE_UNKNOWN;

 class EnvVar
 {
  MString name, oldvalue;
  bool    activated, saved;

  public:
  EnvVar(): activated(false) {}
  ~EnvVar() { Deactivate(); }

  void Activate(const MString& var, const MString& val)
  {
   if(activated) Deactivate();
   name         = var;
   char* curval = getenv(name.Buf());
   if(nullptr == curval)
    saved = false;
   else
   {
    oldvalue = curval;
    saved    = true;
   }
   setenv(name.Buf(), val.Buf(), 1);
  }
  void Deactivate()
  {
   if(!activated) return;
   if(saved)
    setenv(name.Buf(), oldvalue.Buf(), 1);
   else
    unsetenv(name.Buf());
   activated = false;
  }
 };

 EnvVar proxy;

 public:
 using Data = Simple2DData;

 Data ReadVar(const NCFileA& f, const MString& name, size_t i) const
 {
  using DataType             = int2;
  constexpr DataType fill    = -32767;
  real               unitmul = 1.0;

  auto offset = f.A<double>(name, "add_offset");
  auto scale  = f.A<double>(name, "scale_factor");
  if(!offset || !scale) return Data();

  auto unit = f.A<MString>(name, "units");
  if(unit && unit.Get() == "m s-1") unitmul = 100.0;

  Data data((xb < xe) ? (xe - xb + 1) : (nx + xe - xb + 1), ye - yb + 1, Lonb(), Latb());

  if(xb < xe)
  {
   auto var = f.V<DataType>(name, {"longitude", xb, xe - xb + 1}, {"latitude", yb, ye - yb + 1}, {"time", i, 1}, {"depth", layer, 1});
   if(!var) return Data();
   if(var.DimLen(0) != data.Nx() || var.DimLen(1) != data.Ny()) return Data();

   for(size_t ix = 0; ix < var.DimLen(0); ix++)
    for(size_t iy = 0; iy < var.DimLen(1); iy++)
    {
     DataType v   = var(ix, iy);
     data(ix, iy) = (v == fill) ? Data::Fillval() : ((v * scale + offset) * unitmul);
    }
  }
  else
  {
   auto var1 = f.V<DataType>(name, {"longitude", xb}, {"latitude", yb, ye - yb + 1}, {"time", i, 1}, {"depth", layer, 1});
   auto var2 = f.V<DataType>(name, {"longitude", 0, xe + 1}, {"latitude", yb, ye - yb + 1}, {"time", i, 1}, {"depth", layer, 1});
   if(!(var1 && var2)) return Data();
   if((var1.DimLen(0) + var2.DimLen(0)) != data.Nx() || var1.DimLen(1) != data.Ny() || var2.DimLen(1) != data.Ny()) return Data();
   for(size_t ix = 0; ix < var1.DimLen(0); ix++)
    for(size_t iy = 0; iy < var1.DimLen(1); iy++)
    {
     DataType v   = var1(ix, iy);
     data(ix, iy) = (v == fill) ? Data::Fillval() : ((v * scale + offset) * unitmul);
    }
   for(size_t ix = 0; ix < var2.DimLen(0); ix++)
    for(size_t iy = 0; iy < var2.DimLen(1); iy++)
    {
     DataType v                    = var2(ix, iy);
     data(ix + var1.DimLen(0), iy) = (v == fill) ? Data::Fillval() : ((v * scale + offset) * unitmul);
    }
  }
  return data;
 }

 Data ReadVar(const NCFileA& f, const MString& name, const std::vector<size_t>& tindex) const
 {
  Data out;
  if(tindex.size() == 0) return out;

  std::vector<size_t> count;
  for(size_t i = 0; i < tindex.size(); i++)
  {
   Data dat = ReadVar(f, name, tindex[i]);
   if(!dat) return Data();
   out.Add(dat, count);
  }
  out.Div(count);
  return out;
 }

 public:
 NEMOData() = default;

 // TODO: RetVal
 bool Open(const MString& stype, const MString& cred, const MString& proxyurl = "")
 {
  if(proxyurl.Exist()) proxy.Activate("all_proxy", proxyurl);
  MString url, urlt, urls;
  if(stype == "DT")
  {
   url  = "https://" + cred + "@my.cmems-du.eu/thredds/dodsC/cmems_mod_glo_phy_my_0.083_P1D-m";
   type = TYPE_DT;
  }
  if(stype == "NRT")
  {
   url  = "https://" + cred + "@nrt.cmems-du.eu/thredds/dodsC/global-analysis-forecast-phy-001-024";
   type = TYPE_NRT;
  }
  if(stype == "NRT6")
  {
   url  = "https://" + cred + "@nrt.cmems-du.eu/thredds/dodsC/global-analysis-forecast-phy-001-024-3dinst-uovo";
   urlt = "https://" + cred + "@nrt.cmems-du.eu/thredds/dodsC/global-analysis-forecast-phy-001-024-3dinst-thetao";
   urls = "https://" + cred + "@nrt.cmems-du.eu/thredds/dodsC/global-analysis-forecast-phy-001-024-3dinst-so";
   type = TYPE_NRT6;
  }
  nc.Reset(url);
  if(!nc) return false;

  {
   auto nx_ = nc.D("longitude");
   if(!nx_)
   {
    nc.Reset();
    return false;
   }
   nx = nx_;
  }

  if(urlt.Exist())
  {
   nct.Reset(urlt);
   if(!nct)
   {
    nc.Reset();
    return false;
   }
  }
  if(urls.Exist())
  {
   ncs.Reset(urls);
   if(!ncs)
   {
    nc.Reset();
    nct.Reset();
    return false;
   }
  }

  auto rdepths = nc.V<float>("depth");
  if(!rdepths)
  {
   nc.Reset();
   nct.Reset();
   ncs.Reset();
   return false;
  }
  depths.resize(rdepths.DimLen(0));
  for(size_t i = 0; i < depths.size(); i++) depths[i] = rdepths(i);

  auto timeD = nc.V<double>("time");
  auto timeF = nc.V<float>("time");
  if(!(timeD || timeF))
  {
   nc.Reset();
   nct.Reset();
   ncs.Reset();
   return false;
  }
  MDateTime refdate("1950-01-01");
  timeD ? times.resize(timeD.DimLen(0)) : times.resize(timeF.DimLen(0));
  for(size_t i = 0; i < times.size(); i++) times[i] = refdate + static_cast<time_t>(timeD ? timeD(i) : timeF(i)) * 3600;

  return true;
 }

 MString Info() const
 {
  MString d;
  for(size_t i = 0; i < NDepths(); i++) d += MString(" ") + "(" + i + " " + Depth(i) + ")";

  // clang-format off
  return
   "Dataset: " + Title() + "\n" +
   " Begin date: " + Time(0).ToString() + "\n" +
   " End date: " + Time(NTimes()-1).ToString() + "\n" +
   " Time step: " + Timestep() + " seconds\n" +
   " Time moments: " + NTimes() + "\n" +
   " Depths:" + d;
  // clang-format on
 }

 MString Dump() const
 {
  // clang-format off
  return
   "Current settings:\n" + MString() +
   " Longitudes: from " + Lonb() + " (" + xb + ") to "+ Lone() + " (" + xe + ")\n" +
   " Latitudes: from " + Latb() + " (" + yb + ") to "+ Late() + " (" + ye + ")\n" +
   " Depth: layer " + layer + ", depth " + Depth() + " m\n";
  // clang-format on
 }

 bool SetRegion(real lonbin, real lonein, real latb, real late, real depth)
 {
  if(!nc) return false;
  real lonb = ToGeoDomainNeg(lonbin), lone = ToGeoDomainNeg(lonein);

  size_t xb_, xe_, yb_, ye_, layer_ = 0;

  yb_ = static_cast<size_t>(Floor((latb + 80.0) * 12.0));
  ye_ = static_cast<size_t>(Ceil((late + 80.0) * 12.0));
  if(ye_ > 2040) ye_ = 2040;
  if(yb_ >= ye_) return false;

  xb_ = static_cast<size_t>(Floor((lonb + 180.0) * 12.0));
  xe_ = static_cast<size_t>(Ceil((lone + 180.0) * 12.0));
  if(xb_ == xe_) return false;

  if(depth < 0.0 || depth > depths.back())
   layer_ = (depth < 0.0) ? 0 : (depths.size() - 1);
  else
   for(size_t i = 0; i < depths.size() - 1; i++)
   {
    if(depth >= depths[i] && depth <= depths[i + 1])
    {
     layer_ = (depth - depths[i] <= depths[i + 1] - depth) ? i : (i + 1);
     break;
    }
   }
  xb    = xb_;
  xe    = xe_;
  yb    = yb_;
  ye    = ye_;
  layer = layer_;
  return true;
 }

 real      Lonb() const { return nc ? (-180.0 + xb / 12.0) : -1000.0; }
 real      Lone() const { return nc ? (-180.0 + xe / 12.0) : -1000.0; }
 real      Latb() const { return nc ? (-80.0 + yb / 12.0) : -1000.0; }
 real      Late() const { return nc ? (-80.0 + ye / 12.0) : -1000.0; }
 size_t    Xb() const { return nc ? xb : 0; }
 size_t    Xe() const { return nc ? xe : 0; }
 size_t    Yb() const { return nc ? yb : 0; }
 size_t    Ye() const { return nc ? ye : 0; }
 size_t    Nx() const { return nc ? ((xb < xe) ? (xe - xb + 1) : (nx + xe - xb + 1)) : 0; }
 size_t    Ny() const { return nc ? (ye - yb + 1) : 0; }
 real      Depth(size_t l) const { return nc ? depths[l] : -1000.0; }
 real      Depth() const { return Depth(layer); }
 size_t    Layer() const { return layer; }
 time_t    Timestep() const { return nc ? (times[1] - times[0]) : 0; }
 MDateTime Time(size_t i) const
 {
  if((!nc) || i >= times.size()) return MDateTime();
  return times[i];
 }

 bool     isOk() const { return nc; }
 explicit operator bool() const { return nc; }

 size_t NDepths() const { return depths.size(); }
 size_t NTimes() const { return times.size(); }

 template<vartype::Vartype vt> Data Read(size_t it) const = delete;

 template<vartype::Vartype vt> Data Read(const std::vector<size_t>& tindex) const = delete;

 MString Title() const
 {
  switch(type)
  {
  case(TYPE_DT): return "NEMO Delayed time, daily mean (DT)";
  case(TYPE_NRT): return "NEMO Near-real time, daily mean (NRT)";
  case(TYPE_NRT6): return "NEMO Near-real time, 6h resolution (NRT6)";
  default: return "No title";
  }
 }

 static real Fillval() { return Data::Fillval(); }
};

template<> inline NEMOData::Data NEMOData::Read<vartype::Vartype::U>(size_t it) const { return ReadVar(nc, "uo", it); }
template<> inline NEMOData::Data NEMOData::Read<vartype::Vartype::V>(size_t it) const { return ReadVar(nc, "vo", it); }
template<> inline NEMOData::Data NEMOData::Read<vartype::Vartype::TEMP>(size_t it) const { return ReadVar(nct ? nct : nc, "thetao", it); }
template<> inline NEMOData::Data NEMOData::Read<vartype::Vartype::SAL>(size_t it) const { return ReadVar(ncs ? ncs : nc, "so", it); }

template<> inline NEMOData::Data NEMOData::Read<vartype::Vartype::U>(const std::vector<size_t>& tindex) const { return ReadVar(nc, "uo", tindex); }
template<> inline NEMOData::Data NEMOData::Read<vartype::Vartype::V>(const std::vector<size_t>& tindex) const { return ReadVar(nc, "vo", tindex); }
template<> inline NEMOData::Data NEMOData::Read<vartype::Vartype::TEMP>(const std::vector<size_t>& tindex) const { return ReadVar(nct ? nct : nc, "thetao", tindex); }
template<> inline NEMOData::Data NEMOData::Read<vartype::Vartype::SAL>(const std::vector<size_t>& tindex) const { return ReadVar(ncs ? ncs : nc, "so", tindex); }

#endif