Testing linear analisys

actions/actionuv.h

@@ -237,6 +237,7 @@ template<class D> MString ActionUV::DoAction(const CLArgs& args, D& ds)
     if(!err.Exist() && !headwrited) err = fw.AddVariable("div", "", "Velocity divergence", "(" + velunit + ")/" + distunit, "");
     if(!err.Exist() && !headwrited) err = fw.AddVariable("rot", "", "Velocity rotor", "(" + velunit + ")/" + distunit, "");
     if(!err.Exist() && !headwrited) err = fw.AddVariable("ow", "", "Okubo-Weiss parameter", "(" + velunit + ")2/" + distunit + "2", "");
+    if(!err.Exist() && !headwrited) err = fw.AddVariable("discr", "", "", "(" + velunit + ")2/" + distunit + "2", "");
     if(!err.Exist() && !headwrited) err = fw.AddVariable("ke", "", "Squared velocity module, u^2+v^2", "(" + velunit + ")2", "");
     if(!err.Exist() && !headwrited) err = fw.AddVariable("eke", "", "Squared velocity dispersion aka eddy kinetic energy, <u^2+v^2>-<u>^2-<v>^2", "(" + velunit + ")2", "");
     if(!err.Exist() && !headwrited) err = fw.WriteGrid(data);
@@ -246,6 +247,7 @@ template<class D> MString ActionUV::DoAction(const CLArgs& args, D& ds)
     if(!err.Exist()) err = fw.WriteVariable(data, "div", [&data = std::as_const(data)](size_t i, size_t j) { return data.Div(i, j); }, it);
     if(!err.Exist()) err = fw.WriteVariable(data, "rot", [&data = std::as_const(data)](size_t i, size_t j) { return data.Rot(i, j); }, it);
     if(!err.Exist()) err = fw.WriteVariable(data, "ow", [&data = std::as_const(data)](size_t i, size_t j) { return data.OW(i, j); }, it);
+    if(!err.Exist()) err = fw.WriteVariable(data, "discr", [&data = std::as_const(data)](size_t i, size_t j) { return data.Discr(i, j); }, it);
     if(!err.Exist()) err = fw.WriteVariable(data, "ke", [&data = std::as_const(data)](size_t i, size_t j) { return data.U2(i, j); }, it);
     if(!err.Exist())
      err = fw.WriteVariable(

include/cache.h

@@ -1,149 +1,14 @@
 #pragma once
-#include "GPL.h"
-#include <functional>
+#include "merrors.h"
 #include <libpq-fe.h>
 #include <sqlite3.h>
 #include <time.h>
 #include <variant>
 
-using michlib::GPL;
 using michlib::int_cast;
 using michlib::MString;
 using michlib::pointer_cast;
 
-class SQLiteConnection
-{
- public:
- using DBType   = sqlite3*;
- using FuncType = std::function<void(DBType)>;
-
- private:
- static DBType                db;
- static size_t                count;
- static std::vector<FuncType> destructs;
-
- public:
- SQLiteConnection()
- {
-  count++;
-  if(db == nullptr)
-  {
-   MString oldprefix = GPL.UsePrefix("SQLITE");
-   MString name      = GPL.ParameterSValue("db", "");
-   GPL.UsePrefix(oldprefix);
-
-   auto ret = sqlite3_open(name.Buf(), &db);
-   if(ret != SQLITE_OK)
-   {
-    sqlite3_close(db);
-    db = nullptr;
-   }
-  }
- }
-
- SQLiteConnection([[maybe_unused]] const SQLiteConnection& sq): SQLiteConnection() {}
- SQLiteConnection(SQLiteConnection&&) = default;
-
- SQLiteConnection& operator=([[maybe_unused]] const SQLiteConnection& sq)
- {
-  *this = {};
-  return *this;
- }
- SQLiteConnection& operator=(SQLiteConnection&&) = default;
-
- ~SQLiteConnection()
- {
-  if(count == 0) michlib::errmessage("Destructor of SQLiteConnection called on count==0");
-  if(count > 1)
-   count--;
-  else
-  {
-   count = 0;
-   if(db != nullptr)
-   {
-    for(const auto& f: destructs) f(db);
-    sqlite3_close(db);
-   }
-   db = nullptr;
-  }
- }
-
- static void AddDestructor(FuncType&& f) { destructs.emplace_back(std::move(f)); }
-
- operator DBType() const { return db; }
-
- static DBType GetDB() { return db; }
-
- explicit operator bool() const { return db != nullptr; }
-};
-
-class PostgreSQLConnection
-{
- public:
- using DBType   = PGconn*;
- using FuncType = std::function<void(DBType)>;
-
- private:
- static DBType                conn;
- static size_t                count;
- static std::vector<FuncType> destructs;
-
- public:
- PostgreSQLConnection()
- {
-  count++;
-  if(conn == nullptr)
-  {
-   MString oldprefix = GPL.UsePrefix("POSTGRES");
-   MString name      = GPL.ParameterSValue("connection", "");
-   GPL.UsePrefix(oldprefix);
-
-   conn = PQconnectdb(name.Buf());
-   if(PQstatus(conn) != CONNECTION_OK)
-   {
-    michlib::errmessage(PQerrorMessage(conn));
-    PQfinish(conn);
-    conn = nullptr;
-   }
-  }
- }
-
- PostgreSQLConnection([[maybe_unused]] const PostgreSQLConnection& pq): PostgreSQLConnection() {}
- PostgreSQLConnection(PostgreSQLConnection&&) = default;
-
- PostgreSQLConnection& operator=([[maybe_unused]] const PostgreSQLConnection& pq)
- {
-  *this = {};
-  return *this;
- }
- PostgreSQLConnection& operator=(PostgreSQLConnection&&) = default;
-
- ~PostgreSQLConnection()
- {
-  if(count == 0) michlib::errmessage("Destructor of PostgreSQLConnection called on count==0");
-  if(count > 1)
-   count--;
-  else
-  {
-   count = 0;
-   if(conn != nullptr)
-   {
-    for(const auto& f: destructs) f(conn);
-    PQfinish(conn);
-   }
-   conn = nullptr;
-  }
- }
-
- static void AddDestructor(FuncType&& f) { destructs.emplace_back(std::move(f)); }
-
- operator DBType() const { return conn; }
-
- static DBType GetDB() { return conn; }
-
- explicit operator bool() const { return conn != nullptr; }
-};
-
 class GenericCache
 {
  public:
@@ -162,17 +27,18 @@ class FakeCache: public GenericCache
 
 class SQLiteCache: public GenericCache
 {
- static bool regdest;
-
- SQLiteConnection db;
+ sqlite3* db = nullptr;
 
  public:
- bool Init()
+ bool Init(const MString& name)
  {
-  if(!db) return false;
-
-  if(!regdest)
+  Close();
+  auto ret = sqlite3_open(name.Buf(), &db);
+  if(ret != SQLITE_OK)
   {
+   Close();
+   return false;
+  }
   // Create table
   sqlite3_stmt* sqst;
   int           i;
@@ -185,25 +51,18 @@ class SQLiteCache: public GenericCache
   if(i != SQLITE_DONE)
   {
    sqlite3_finalize(sqst);
+   Close();
    return false;
   }
   sqlite3_finalize(sqst);
   sqlite3_busy_timeout(db, 1000);
+  return true;
+ }
 
-   db.AddDestructor(
-       [](SQLiteConnection::DBType db)
+ void Close()
  {
-        sqlite3_stmt* sqst = nullptr;
-        int           i    = SQLITE_OK;
-
-        if(i == SQLITE_OK) i = sqlite3_prepare_v2(db, "DELETE from `cache` WHERE exptime<?1;", -1, &sqst, 0);
-        if(i == SQLITE_OK) i = sqlite3_bind_int64(sqst, 1, time(nullptr));
-        if(i == SQLITE_OK) i = sqlite3_step(sqst);
-        sqlite3_finalize(sqst);
-       });
-   regdest = true;
-  }
-  return true;
+  if(db != nullptr) sqlite3_close(db);
+  db = nullptr;
  }
 
  virtual bool Put(const MString& key, const MString& value, size_t ttl) const override
@@ -248,16 +107,26 @@ class SQLiteCache: public GenericCache
   return {"", false};
  }
 
- virtual ~SQLiteCache() override = default;
+ virtual ~SQLiteCache() override
+ {
+  if(!*this) return;
+
+  sqlite3_stmt* sqst = nullptr;
+  int           i    = SQLITE_OK;
+
+  if(i == SQLITE_OK) i = sqlite3_prepare_v2(db, "DELETE from `cache` WHERE exptime<?1;", -1, &sqst, 0);
+  if(i == SQLITE_OK) i = sqlite3_bind_int64(sqst, 1, time(nullptr));
+  if(i == SQLITE_OK) i = sqlite3_step(sqst);
+  sqlite3_finalize(sqst);
+  sqlite3_close_v2(db);
+ }
 
  explicit operator bool() const { return db != nullptr; }
 };
 
 class PostgreSQLCache: public GenericCache
 {
- static bool regdest;
-
- PostgreSQLConnection conn;
+ PGconn* conn = nullptr;
 
  bool CheckCon() const
  {
@@ -279,41 +148,39 @@ class PostgreSQLCache: public GenericCache
  }
 
  public:
- bool Init()
+ bool Init(const MString& name)
  {
-  if(!conn) return false;
+  Close();
+  conn = PQconnectdb(name.Buf());
 
-  if(!regdest)
-  {
-   // Create table
+  if(PQstatus(conn) != CONNECTION_OK)
   {
-    auto* res = PQexec(conn, "SET client_min_messages=WARNING;");
-    PQclear(res);
+   michlib::errmessage(PQerrorMessage(conn));
+   Close();
+   return false;
   }
+
+  // Create table
+  if(false)
   {
    auto* res = PQexec(conn, "CREATE TABLE IF NOT EXISTS cache(key TEXT PRIMARY KEY NOT NULL, value BYTEA, exptime TIMESTAMP(0) NOT NULL);");
    if(PQresultStatus(res) != PGRES_COMMAND_OK)
    {
     michlib::errmessage(PQresStatus(PQresultStatus(res)));
     michlib::errmessage(PQerrorMessage(conn));
-    }
     PQclear(res);
+    Close();
    }
-   {
-    auto* res = PQexec(conn, "SET client_min_messages=NOTICE;");
+   else
     PQclear(res);
   }
+  return true;
+ }
 
-   conn.AddDestructor(
-       [](PostgreSQLConnection::DBType conn)
+ void Close()
  {
-        auto* res = PQexec(conn, "DELETE FROM cache WHERE exptime<localtimestamp;");
-        PQclear(res);
-       });
-
-   regdest = true;
-  }
-  return true;
+  if(conn != nullptr) PQfinish(conn);
+  conn = nullptr;
  }
 
  virtual bool Put(const MString& key, const MString& value, size_t ttl) const override
@@ -368,7 +235,14 @@ class PostgreSQLCache: public GenericCache
   return {std::move(val), true};
  }
 
- virtual ~PostgreSQLCache() override = default;
+ virtual ~PostgreSQLCache() override
+ {
+  if(!CheckCon()) return;
+
+  auto* res = PQexec(conn, "DELETE FROM cache WHERE exptime<localtimestamp;");
+  PQclear(res);
+  Close();
+ }
 
  explicit operator bool() const { return conn != nullptr; }
 };
@@ -376,23 +250,26 @@ class PostgreSQLCache: public GenericCache
 inline GenericCache* CreateCache(const MString& cachedesc)
 {
  auto i = cachedesc.GetPos(':');
- auto name = i == 0 ? cachedesc : cachedesc.SubStr(1, i - 1);
- auto par  = i == 0 ? "" : cachedesc.SubStr(i + 1, cachedesc.Len() - i);
+ if(i == 0)
+ {
+  if(cachedesc == "no") return new FakeCache;
+  return nullptr;
+ }
 
- if(name == "no") return new FakeCache;
+ auto name = cachedesc.SubStr(1, i - 1);
+ auto par  = cachedesc.SubStr(i + 1, cachedesc.Len() - i);
 
  if(name == "sqlite")
  {
   auto ret = new SQLiteCache;
-  ret->Init();
+  ret->Init(par);
   if(*ret) return ret;
   delete ret;
  }
-
  if(name == "postgre" || name == "postgres" || name == "postgresql")
  {
   auto ret = new PostgreSQLCache;
-  ret->Init();
+  ret->Init(par);
   if(*ret) return ret;
   delete ret;
  }

include/uvdata.h

@@ -293,6 +293,58 @@ template<class OneVarData> class UVData: public UVDataDims<OneVarData, internal:
   return sn * sn + ss * ss - w * w;
  }
 
+ // Discriminant of characteristic equation
+ real Discr(size_t i) const { return Discr(i % nx, i / nx); }
+ real Discr(size_t ix, size_t iy) const
+ {
+  if(!*this) return 0.0;
+  real ux = dUdX(ix, iy);
+  real uy = dUdY(ix, iy);
+  real vx = dVdX(ix, iy);
+  real vy = dVdY(ix, iy);
+
+  if(ux == Fillval() || uy == Fillval() || vx == Fillval() || vy == Fillval()) return Fillval();
+
+  real sn = ux - vy;
+  real ss = vx + uy;
+  real w  = michlib::Abs(uy - vx);
+
+  if(sn * sn + 4.0 * uy * vx >= 0)
+   return michlib::Acos(w / Sqrt(ss * ss + sn * sn)) * 2.0 / M_PI;
+  else
+  {
+   real L = (ux + vy) * 0.5;
+   real W = michlib::Sqrt(-(sn * sn + 4.0 * uy * vx)) * 0.5;
+   michlib::message(W);
+
+   real C2 = uy / W;
+   real C4 = (vy - L) / W;
+
+   real A = 2.0 * C4;
+   real B = C4 * C4 + C2 * C2 - 1.0;
+
+   real t1 = michlib::Atan2(-B, A) * 0.5;
+   real t2 = michlib::Atan2(B, -A) * 0.5;
+
+   real d1, d2;
+   {
+    real tsn  = michlib::Sin(t1);
+    real tcn  = michlib::Cos(t1);
+    real temp = tcn + C4 * tsn;
+    d1        = C2 * C2 * tsn * tsn + temp * temp;
+   }
+   {
+    real tsn  = michlib::Sin(t2);
+    real tcn  = michlib::Cos(t2);
+    real temp = tcn + C4 * tsn;
+    d2        = C2 * C2 * tsn * tsn + temp * temp;
+   }
+
+   real mind = std::min(d1, d2), maxd = std::max(d1, d2);
+   return -michlib::Sqrt(mind / maxd);
+  }
+ }
+
  auto StablePoints(size_t ix, size_t iy) const
  {
   std::vector<struct B::StPoint> points;

sources/COPERNICUS.cpp

@@ -129,9 +129,6 @@ Error COPERNICUSData::Mirror(const CLArgs& args) const
   dsets = dlist.Value();
  }
 
- michlib::RegExpSimple filter((args.contains("filter") ? args.at("filter") : ".*").Buf());
- if(filter.Compile() != 0) return Error(pref, MString("Can't compile regular expression ") + filter.RegStr());
-
  CURLRAII chandle;
  for(const auto& dset: dsets)
  {
@@ -157,34 +154,20 @@ Error COPERNICUSData::Mirror(const CLArgs& args) const
    while(rpos != rfiles.size() || lpos != lfiles.size())
    {
     if(rpos == rfiles.size())
-     while(lpos != lfiles.size())
-     {
-      if(filter.Match(lfiles[lpos].name.Buf())) rem.push_back(lpos);
-      lpos++;
-     }
+     while(lpos != lfiles.size()) rem.push_back(lpos++);
     if(lpos == lfiles.size())
-     while(rpos != rfiles.size())
-     {
-      if(filter.Match(rfiles[rpos].name.Buf())) down.push_back(rpos);
-      rpos++;
-     }
+     while(rpos != rfiles.size()) down.push_back(rpos++);
     if(rpos == rfiles.size() || lpos == lfiles.size()) continue;
 
     if(rfiles[rpos].name < lfiles[lpos].name)
-    {
-     if(filter.Match(rfiles[rpos].name.Buf())) down.push_back(rpos);
-     rpos++;
-    }
+     down.push_back(rpos++);
     else if(lfiles[lpos].name < rfiles[rpos].name)
-    {
-     if(filter.Match(lfiles[lpos].name.Buf())) rem.push_back(lpos);
-     lpos++;
-    }
+     rem.push_back(lpos++);
     else
     {
      auto delta = rfiles[rpos].mtime.Epoch() - lfiles[lpos].mtime.Epoch();
      if(delta < 0) delta = -delta;
-     if((delta > 0 || rfiles[rpos].size != lfiles[lpos].size) && filter.Match(lfiles[lpos].name.Buf())) upd.emplace_back(rpos, lpos);
+     if(delta > 0 || rfiles[rpos].size != lfiles[lpos].size) upd.emplace_back(rpos, lpos);
      lpos++;
      rpos++;
     }

sources/NEMO.h

@@ -7,7 +7,6 @@ class NEMOData: public LayeredDataZ
  {
   TYPE_UNKNOWN,
   TYPE_DT,
-  TYPE_DT1,
   TYPE_NRT,
   TYPE_NRT6,
   TYPE_BALTICDT,
@@ -35,7 +34,6 @@ class NEMOData: public LayeredDataZ
   switch(type)
   {
   case(TYPE_DT): return "NEMO Delayed time, daily mean (DT)";
-  case(TYPE_DT1): return "NEMO Delayed time, daily mean, part 2 (DT1)";
   case(TYPE_NRT): return "NEMO Near-real time, daily mean (NRT)";
   case(TYPE_NRT6): return "NEMO Near-real time, 6h resolution (NRT6)";
   case(TYPE_BALTICDT): return "NEMO Delayed time, Baltic region, daily mean (BALTICDT)";
@@ -61,8 +59,6 @@ class NEMOData: public LayeredDataZ
   GPL.UsePrefix("NEMO");
   if(dataset == "DT")
    type = TYPE_DT;
-  else if(dataset == "DT1")
-   type = TYPE_DT1;
   else if(dataset == "NRT")
    type = TYPE_NRT;
   else if(dataset == "NRT6")

src/cache.cpp

@@ -1,14 +0,0 @@
-#define MICHLIB_NOSOURCE
-#include "cache.h"
-
-SQLiteConnection::DBType                SQLiteConnection::db        = nullptr;
-size_t                                  SQLiteConnection::count     = 0;
-std::vector<SQLiteConnection::FuncType> SQLiteConnection::destructs = {};
-
-PostgreSQLConnection::DBType                PostgreSQLConnection::conn      = nullptr;
-size_t                                      PostgreSQLConnection::count     = 0;
-std::vector<PostgreSQLConnection::FuncType> PostgreSQLConnection::destructs = {};
-
-bool SQLiteCache::regdest = false;
-
-bool PostgreSQLCache::regdest = false;

src/copcat.cpp

@@ -7,20 +7,16 @@ const MString CopernicusCatalog::caturl = "https://stac.marine.copernicus.eu/met
 
 CopernicusCatalog::CopernicusCatalog()
 {
- auto oldprefix = michlib::GPL.UsePrefix("COPERNICUS");
  // Cache
+ auto oldprefix = michlib::GPL.UsePrefix("COPERNICUS");
  cache.reset(CreateCache(michlib::GPL.ParameterSValue("Cache", "")));
+ michlib::GPL.UsePrefix(oldprefix);
  if(!cache)
  {
   michlib::errmessage("Can't init cache");
   cache.reset(new FakeCache);
  }
 
- // Proxy
- auto proxyurl = michlib::GPL.ParameterSValue("Proxy", "");
- if(proxyurl.Exist()) curl_easy_setopt(chandle, CURLOPT_PROXY, proxyurl.Buf());
-
- michlib::GPL.UsePrefix(oldprefix);
  GetCatalog();
 }
 
@@ -48,14 +44,9 @@ RetVal<std::vector<MString>> CopernicusCatalog::ProductList() const
  for(Json::ArrayIndex i = 0; i < links.size(); i++)
  {
   const auto& rel  = links[i]["rel"];
-  const auto& href = links[i]["href"];
+  const auto& titl = links[i]["title"];
 
-  if(rel.type() == Json::stringValue && href.type() == Json::stringValue && rel.asString() == "child")
-  {
-   auto str = href.asString();
-   str.erase(str.find('/'));
-   out.emplace_back(str.c_str());
-  }
+  if(rel.type() == Json::stringValue && titl.type() == Json::stringValue && rel.asString() == "child") out.emplace_back(titl.asString().c_str());
  }
  return out;
 }
@@ -71,8 +62,9 @@ RetVal<MString> CopernicusCatalog::ProductURL(const MString& prod) const
 
  for(Json::ArrayIndex i = 0; i < links.size(); i++)
  {
+  const auto& titl = links[i]["title"];
   const auto& href = links[i]["href"];
-  if(href.type() == Json::stringValue && href.asString() == (prod + "/product.stac.json").Buf()) return DirName(caturl) + "/" + MString(href.asString().c_str());
+  if(titl.type() == Json::stringValue && href.type() == Json::stringValue && titl.asString().c_str() == prod) return DirName(caturl) + "/" + MString(href.asString().c_str());
  }
  return {pref, "unknown product: " + prod};
 }
@@ -98,14 +90,9 @@ RetVal<std::vector<MString>> CopernicusCatalog::DatasetList(const MString& prod)
  for(Json::ArrayIndex i = 0; i < links.size(); i++)
  {
   const auto& rel  = links[i]["rel"];
-  const auto& href = links[i]["href"];
+  const auto& titl = links[i]["title"];
 
-  if(rel.type() == Json::stringValue && href.type() == Json::stringValue && rel.asString() == "item")
-  {
-   auto str = href.asString();
-   str.erase(str.find('/'));
-   out.emplace_back(str.c_str());
-  }
+  if(rel.type() == Json::stringValue && titl.type() == Json::stringValue && rel.asString() == "item") out.emplace_back(titl.asString().c_str());
  }
  return out;
 }
@@ -129,8 +116,9 @@ RetVal<MString> CopernicusCatalog::DatasetURL(const MString& prod, const MString
 
  for(Json::ArrayIndex i = 0; i < links.size(); i++)
  {
+  const auto& titl = links[i]["title"];
   const auto& href = links[i]["href"];
-  if(href.type() == Json::stringValue && href.asString() == (dataset + "/dataset.stac.json").Buf()) return DirName(url) + "/" + MString(href.asString().c_str());
+  if(titl.type() == Json::stringValue && href.type() == Json::stringValue && titl.asString().c_str() == dataset) return DirName(url) + "/" + MString(href.asString().c_str());
  }
  return {pref, "unknown dataset: " + dataset};
 }