smcc.cpp
#include "capd/capdlib.h"
#include "capd/mpcapdlib.h"
using namespace std;
using namespace capd;
#include "misc.h"
#include "eps.h"
#include "symb.h"
enum StabilityResults {stablemaybe=0,stableyes=1,stableno=2};
string StabilityResultsInfo[3]={"STABLEMAYBE","STABLEYES","STABLENO"};
enum VersionType {Version3D=0,Version5D=1,Version1DMap=2,Version4D=3};
string VersionInfo[4]={"3D","5D","1DMap","4D"};
struct Options {
string Action="";
double TimeStep=0.0; // set to a nonpositive value to use automatic TimeStep
double MaxTimeStep=0.0;
double MinTimeStep=0.0;
double MaxReturnTime=0.0;
double SkipTime=0.0;
double IntegrationTime=100.0;
int IterNum=100;
int SkipNum=0;
double radius=0.0;
string InitCondStr="3 0.05 0.03 0.0";
double C1r,C2r,Lr,Rr,Gr,Gar,gr,I0r,g0r,V0r,gamma0r,betar,g1r,Eagr,Earr,qr,ar,Vg0r,Lgr,Cthr,tthr,T0r,Kr;
DInterval C1i,C2i,Li,Ri,Gi,Gai,gi,I0i,g0i,V0i,gamma0i,betai,g1i,Eagi,Eari,qi,ai,Vg0i,Lgi,Cthi,tthi,T0i,Ki;
string C1s,C2s,Ls,Rs,Gs,Gas,gs,I0s,g0s,V0s,gamma0s,betas,g1s,Eags,Ears,qs,as,Vg0s,Lgs,Cths,tths,T0s,Ks;
double gminr,gmaxr;
DInterval gmini,gmaxi;
string gmins,gmaxs;
int PoincVar=1;
double PoincVal=0.0;
// int PoincWhich=poincare::PlusMinus; // +1
// int PoincWhich=poincare::Both; // 0
int PoincWhich=poincare::MinusPlus; // -1;
string BiffParName="g";
double BiffParMin=0.85e-9,BiffParMax=1.2e-9;
double ReadParMin=-1e10,ReadParMax=1e10;
int BiffParNum=101;
bool BiffChangeStart=false;
int XDisp=1;
int YDisp=2;
int Period=1;
int MaxPeriod=1;
int BinNum=1000;
double DeltaPseudoOrbit=0.01;
double XMin=0.0,XMax=0.5;
double YMin=0.0,YMax=0.5;
string InFile="";
string InFile2="";
string QuadrFile="";
string OutFile="";
string BaseFile="";
string Info="";
string EpsFileName="";
int CircuitVersion=0;
bool WriteInitialPoint=true;
bool WriteTrajectory=true;
int Order=30;
double AbsTol=-1.0;
double RelTol=-1.0;
IMatrix CoordChange;
int CoutPrec=8;
int MpPrec=128; // 64 - slightly more than double precision (53 bits)
double DeltaMin=1e-13;
bool Verbose=false;
EPSOptions EPSOpt;
double MaxPODist=1e-14; // maximum distance between P(x_k) and x_{k+1} to assume the convergence
double MpMaxPODist=1e-20; // maximum distance between P(x_k) and x_{k+1} to assume the convergence, for MpFloat computations
double MinPOSelfDist=1e-13; // minimum distance from another point in the same orbit to consider the orbit a subperiod orbit
int NewtonMaxIterNum=500;
double NewtonMaxCorrection=0.1;
bool NewtonImproveApprox=true;
double MaxTestSize=0.0;
bool SingleColor=true;
double v10=0.0;
Options();
} Opt;
Options::Options(){
/*
SMCCinitcond3.resize(3);
double X0=0.05; double Y0=0.03; double Z0=0.0;
SMCCinitcond3(1)=X0;
SMCCinitcond3(2)=Y0;
SMCCinitcond3(3)=Z0;
SMCCinitcond5.resize(5);
SMCCinitcond5(1)=SMCCinitcond3(1);
SMCCinitcond5(2)=SMCCinitcond3(2);
SMCCinitcond5(3)=SMCCinitcond3(3);
SMCCinitcond5(4)=1.18e-9;
SMCCinitcond5(5)=298.0;*/
double eta=20;
C1r=4.7e-9/eta; C1i=C1r; C1s="0.235e-9";
C2r=47e-9/eta; C2i=C2r; C2s="2.35e-9";
Rr=1.6e3*eta; Ri=Rr; Rs="32e3";
Lr=8.5e-3*eta; Li=Lr; Ls="170e-3";
Gr=1/Rr; Gi=Gr; Gs="3.125e-5";
Gar=3.8137e-5; Gai=Gar; Gas="3.8137e-5";
gr=1e-9; gi=gr; gs="1e-9";
I0r=6.14e-5; I0i=I0r; I0s="6.14e-5";
g0r=2.75e-10; g0i=g0r; g0s="2.75e-10";
V0r=0.43; V0i=V0r; V0s="0.43";
gamma0r=16.5*(1.1); gamma0i=gamma0r; gamma0s="18.15";
betar=1.25; betai=betar; betas="1.25";
g1r = 1e-9; g1i=g1r; g1s = "1e-9";
Eagr=1.501*1.6e-19; Eagi=Eagr; Eags="2.4016e-19";
Earr=1.5*1.6e-19; Eari=Earr; Ears="2.4e-19";
qr=1.6e-19; qi=qr; qs="1.6e-19";
ar=0.25e-9; ai=ar; as="0.25e-9";
Vg0r=300.0; Vg0i=Vg0r; Vg0s="300.0";
Lgr=5e-9; Lgi=Lgr; Lgs="5e-9";
Cthr=3.18e-16; Cthi=Cthr; Cths="3.18e-16";
tthr=2.3e-10; tthi=tthr; tths="2.3e-10";
T0r=298.0; T0i=T0r; T0s="298.0";
Kr=1.38e-23; Ki=Kr; Ks="1.38e-23";
gminr=1e-10; gmini=gminr; gmins="1e-10";
gmaxr=1.7e-9; gmaxi=gmaxi; gmaxs="1.7e-9";
EPSOpt.Diameter=0.0001;
EPSOpt.DiameterSymbol=1.0;
EPSOpt.LineWidth=0.0;
EPSOpt.bbXMax=600;
EPSOpt.bbYMax=300;
EPSOpt.XLabel="y";
EPSOpt.YLabel="z";
EPSOpt.XLabelShow=true;
EPSOpt.YLabelShow=true;
EPSOpt.XYLabelFontSize=18;
EPSOpt.DrawGrid=true;
EPSOpt.DrawGridX=true;
EPSOpt.DrawGridY=true;
EPSOpt.RMarg=30;
// EPSOpt.Clip=false;
EPSOpt.DrawBackground=false;
EPSOpt.BackgroundColor="1 1 1";
EPSOpt.DefineMinMax(0.0,1.0,0.0,1.0);
}
double MaxPODist(const double & fmax){ return Opt.MaxPODist;}
MpFloat MaxPODist(const MpFloat & fmax){ return Opt.MpMaxPODist;}
/************** SAVE CONFIGURATION ******************/
string ConfigFileName="smcc.cfg";
class ApplConfigFile : public ConfigFile {
public:
ApplConfigFile():ConfigFile(ConfigFileName){;};
void Init();
void Read(string FN);
void Save(string FN);
};
void ApplConfigFile::Read(string FN){
ConfigFile::Read(FN);
}
void ApplConfigFile::Save(string FN){
ConfigFile::Save(FN);
}
void ApplConfigFile::Init(){
AddComment("# Memristor, Configuration File");
AddString("Action",&(Opt.Action));
AddDouble("TimeStep",&(Opt.TimeStep));
AddDouble("MaxTimeStep",&(Opt.MaxTimeStep));
AddDouble("MinTimeStep",&(Opt.MinTimeStep));
AddDouble("MaxReturnTime",&(Opt.MaxReturnTime));
AddDouble("SkipTime",&(Opt.SkipTime));
AddDouble("IntegrationTime",&(Opt.IntegrationTime));
AddInt("IterNum",&(Opt.IterNum));
AddInt("SkipNum",&(Opt.SkipNum));
AddDouble("radius",&(Opt.radius));
AddString("InitCondStr",&(Opt.InitCondStr));
AddDouble("C1r",&(Opt.C1r));
AddDouble("C2r",&(Opt.C2r));
AddDouble("gr",&(Opt.gr));
AddString("gs",&(Opt.gs));
AddDouble("gminr",&(Opt.gminr));
AddString("gmins",&(Opt.gmins));
AddDouble("gmaxr",&(Opt.gmaxr));
AddString("gmaxs",&(Opt.gmaxs));
AddInt("PoincVar",&(Opt.PoincVar));
AddDouble("PoincVal",&(Opt.PoincVal));
AddInt("PoincWhich",&(Opt.PoincWhich));
AddString("BiffParName",&(Opt.BiffParName));
AddDouble("BiffParMin",&(Opt.BiffParMin));
AddDouble("BiffParMax",&(Opt.BiffParMax));
AddDouble("ReadParMin",&(Opt.ReadParMin));
AddDouble("ReadParMax",&(Opt.ReadParMax));
AddInt("BiffParNum",&(Opt.BiffParNum));
AddBool("BiffChangeStart",&(Opt.BiffChangeStart));
AddInt("XDisp",&(Opt.XDisp));
AddInt("YDisp",&(Opt.YDisp));
AddInt("Period",&(Opt.Period));
AddInt("MaxPeriod",&(Opt.MaxPeriod));
AddInt("BinNum",&(Opt.BinNum));
AddDouble("DeltaPseudoOrbit",&(Opt.DeltaPseudoOrbit));
AddString("InFile",&(Opt.InFile));
AddString("InFile2",&(Opt.InFile2));
AddString("QuadrFile",&(Opt.QuadrFile));
AddString("OutFile",&(Opt.OutFile));
AddString("BaseFile",&(Opt.BaseFile));
AddString("Info",&(Opt.Info));
AddString("EpsFileName",&(Opt.EpsFileName));
AddDouble("XMin",&(Opt.XMin));
AddDouble("XMax",&(Opt.XMax));
AddDouble("YMin",&(Opt.YMin));
AddDouble("YMax",&(Opt.YMax));
AddInt("Order",&(Opt.Order));
AddBool("WriteInitialPoint",&(Opt.WriteInitialPoint));
AddBool("WriteTrajectory",&(Opt.WriteTrajectory));
AddInt("CoutPrec",&(Opt.CoutPrec));
AddDouble("AbsTol",&(Opt.AbsTol));
AddDouble("RelTol",&(Opt.RelTol));
AddInt("MpPrec",&(Opt.MpPrec));
AddDouble("DeltaMin",&(Opt.DeltaMin));
AddInt("CircuitVersion",&(Opt.CircuitVersion));
AddBool("Verbose",&(Opt.Verbose));
AddDouble("MaxPODist",&(Opt.MaxPODist));
AddDouble("MpMaxPODist",&(Opt.MpMaxPODist));
AddDouble("MinPOSelfDist",&(Opt.MinPOSelfDist));
AddInt("NewtonMaxIterNum",&(Opt.NewtonMaxIterNum));
AddDouble("NewtonMaxCorrection",&(Opt.NewtonMaxCorrection));
AddBool("NewtonImproveApprox",&(Opt.NewtonImproveApprox));
AddDouble("MaxTestSize",&(Opt.MaxTestSize));
AddBool("SingleColor",&(Opt.SingleColor));
AddDouble("v10",&(Opt.v10));
AddDouble("EPSOpt.XMin",&(Opt.EPSOpt.XMin));
AddDouble("EPSOpt.XMax",&(Opt.EPSOpt.XMax));
AddDouble("EPSOpt.YMin",&(Opt.EPSOpt.YMin));
AddDouble("EPSOpt.YMax",&(Opt.EPSOpt.YMax));
AddString("EPSOpt.XLabel",&(Opt.EPSOpt.XLabel));
AddString("EPSOpt.YLabel",&(Opt.EPSOpt.YLabel));
AddDouble("EPSOpt.Diameter",&(Opt.EPSOpt.Diameter));
AddDouble("EPSOpt.LineWidth",&(Opt.EPSOpt.LineWidth));
AddInt("EPSOpt.bbXMax",&(Opt.EPSOpt.bbXMax));
AddInt("EPSOpt.bbYMax",&(Opt.EPSOpt.bbYMax));
}
void Usage(){
cout << "smcc.x options " << endl;
exit(1);
}
void CommandLineOptions(int argc, char *argv[]){
ApplConfigFile CF;
CF.Init();
for (int i=1;i<argc;i++){
char *arg=argv[i];
if (arg[0]=='-'){
switch (arg[1]){
default:
if (++i>=argc) Usage();
if (!CF.GetValue(argv[i-1],argv[i])) Usage();
}
}
}
// copy parameters to interval values:
Opt.C1i=Opt.C1r; Opt.C2i=Opt.C2r; Opt.Li=Opt.Lr;
Opt.Ri=Opt.Rr; Opt.Gi=Opt.Gr; Opt.Gai=Opt.Gar; Opt.gi=Opt.gr; Opt.I0i=Opt.I0r; Opt.g0i=Opt.g0r; Opt.V0i=Opt.V0r;
}
double Dist(const DVector & x, const DVector & y){
double res=0.0;
if (x.dimension()!=y.dimension()) return 1e100;
for (int i=0;i<x.dimension();++i){
double d=fabs(x[i]-y[i]);
if (d>res) res=d;
}
return res;
}
template <class TVector>
typename TVector::ScalarType Dist(const TVector & x, const TVector & y){
typename TVector::ScalarType res=0.0;
if (x.dimension()!=y.dimension()) return 1e100;
for (int i=0;i<x.dimension();++i){
typename TVector::ScalarType d=abs(x[i]-y[i]);
if (d>res) res=d;
}
return res;
}
template <class TVector>
typename TVector::ScalarType VectorMaxAbsVal(const TVector & v){
typename TVector::ScalarType res=0.0;
for (int d=0;d<v.dimension();++d){
if (abs(v[d])>res) res=abs(v[d]);
}
return res;
}
template <class TMatrix>
typename TMatrix::ScalarType MatrixMaxAbsVal(const TMatrix & x){
typename TMatrix::ScalarType res=0.0;
for (int i=1;i<=x.numberOfRows();++i){
for (int j=1;j<=x.numberOfColumns();++j){
if (capd::abs(x(i,j))>res) res=capd::abs(x(i,j));
}
}
return res;
}
template <class TInterval>
typename TInterval::ScalarType IntervalMaxDiam(const TInterval & x){
return diam(x).rightBound();
}
template <class TIVector>
typename TIVector::ScalarType VectorMaxDiam(const TIVector & x){
typename TIVector::ScalarType res(0.0);
for (int i=0;i<x.dimension();++i)
if (diam(x[i]).rightBound()>res) res=diam(x[i]).rightBound();
return res;
}
template <class TIMatrix>
typename TIMatrix::ScalarType MatrixMaxDiam(const TIMatrix & x){
typename TIMatrix::ScalarType res(0.0);
for (int i=1;i<=x.numberOfRows();++i)
for (int j=1;j<=x.numberOfColumns();++j)
if (diam(x(i,j)).rightBound()>res) res=diam(x(i,j)).rightBound();
return res;
}
template <class TIVector>
typename TIVector::ScalarType VectorMid(const TIVector & v){
typename TIVector::ScalarType res=0.0;
for (int d=0;d<v.dimension();++d)
if (abs(v[d])>res) res=0.5*(v[d].leftBound()+v[d].rightBound());
return res;
}
template <class TMatrix>
TMatrix MatrixRemoveRowColumn(TMatrix & M, int pos){ // remove the row with the index "pos" and the column with the index "pos", attention: indexes are counted from 0
int dim=M.numberOfRows();
if (dim!=M.numberOfColumns()){
cerr << "MatrixRemoveRowColumn: " << "wrong dimensions." << endl;
return M;
}
TMatrix MPart((dim-1),(dim-1),0.0);
for (int i=0, posi=0;i<dim;++i){
if (i==pos) continue;
for (int j=0, posj=0;j<dim;++j){
if (j==pos) continue;
MPart[posi][posj]=M[i][j];
posj++;
}
posi++;
}
return MPart;
}
template <class TIVector>
bool EmptyIntersection(const TIVector & V1, const TIVector & V2){
TIVector res;
try{
res=intersection(V1,V2);
}
catch(exception& e){
// cout << e.what() << endl;
return true;
}
return false;
}
template <class TVector, class TScalar>
int POMinPer(const TVector & xfull, int period, TScalar MaxAllowedDist){
int dim=xfull.dimension()/period;
TVector x0(dim);
for (int i=0;i<dim;++i) x0[i]=xfull[i];
for (int p=1;p<period;++p){
TVector x(dim);
for (int i=0;i<dim;++i) x[i]=xfull[p*dim+i];
if (MaxDist(x0,x)<MaxAllowedDist) return p;
}
return period;
}
string SystemName(int Version){
if (Version==Version3D) return "smcc3";
if (Version==Version5D) return "smcc5";
cout << "Unknown version" << Version << endl;
return "";
}
//im = I0*exp(-g/g0)*0.5*(exp(V1/V0)-exp(V1/V0));
//gamma = (gamma0-beta*(g/g1)*(g/g1)*(g/g1));
string VectorFieldString(int Version){
if (Version==Version3D) return "par:C1,C2,L,G,Ga,g,I0,g0,V0;var:V1,V2,IL;fun:"
"(Ga*V1+G*(V2-V1)-I0*exp(-g/g0)*0.5*(exp(V1/V0)-exp(-V1/V0)))/C1,"
"(G*(V1-V2)+IL)/C2,"
"-V2/L;";
if (Version==Version5D) return "par:C1,C2,L,G,Ga,I0,g0,V0,gamma0,beta,g1,Eag,Ear,q,a,Vg0,Lg,Cth,tth,T0,K;var:V1,V2,IL,g,T;fun:"
"(Ga*V1+G*(V2-V1)-I0*exp(-g/g0)*0.5*(exp(V1/V0)-exp(-V1/V0)))/C1,"
"(G*(V1-V2)+IL)/C2,"
"-V2/L,"
"-Vg0*(exp(-Eag/(K*T)+q*a*(gamma0-beta*(g/g1)*(g/g1)*(g/g1))*V1/(Lg*K*T))-exp(-Ear/(K*T)-q*a*(gamma0-beta*(g/g1)*(g/g1)*(g/g1))*V1/(Lg*K*T)))/2,"
"V1*I0*exp(-g/g0)*0.5*(exp(V1/V0)-exp(-V1/V0))/Cth-(T-T0)/tth;";
if (Version==Version1DMap) return "par:G,Ga,I0,g0,V0,g;var:V1;fun:(Ga-G)*V1-I0*exp(-g/g0)*0.5*(exp(V1/V0)-exp(-V1/V0));";
if (Version==Version4D) return "par:C1,C2,L,G,Ga,I0,g0,V0,Cth,tth,T0,g;var:V1,V2,IL,T;fun:"
"(Ga*V1+G*(V2-V1)-I0*exp(-g/g0)*0.5*(exp(V1/V0)-exp(-V1/V0)))/C1,"
"(G*(V1-V2)+IL)/C2,"
"-V2/L,"
"V1*I0*exp(-g/g0)*0.5*(exp(V1/V0)-exp(-V1/V0))/Cth-(T-T0)/tth;";
cout << "Unknown version" << Version << endl;
return "";
}
template <class TVector>
string VectorToString(const TVector & x){
ostringstream oss;
oss.precision(Opt.CoutPrec);
int dim=x.dimension();
oss << dim;
for (int d=1;d<=dim;d++) oss << ' ' << x(d);
return oss.str();
}
template <class TVector>
bool StringToVector(string s, TVector & X){
istringstream is(s);
int dim=0;
is >> dim;
if (dim<=0) return false;
X.resize(dim);
for (int d=1;d<=dim;d++) is >> X(d);
return true;
}
template <class TVector>
bool InitialConditions(int Version, string InitCondString, TVector & X){
if (Version!=Version3D && Version!=Version5D){cout << "Unknown version" << Version << endl; return false;}
if (!StringToVector(InitCondString,X)) return false;
int dim=X.dimension();
if (Version==Version3D && dim!=3 || Version==Version5D && dim!=5) return false;
cout << "InitCond: " << X << endl;
return true;
}
template <class TMap>
void WriteParameters(string s, TMap m){
// cout << "WriteParameters" << endl;
int pos=s.find(";");
string sub=(pos==string::npos?s:s.substr(0, pos));
pos=sub.find(":");
sub=(pos==string::npos?sub:sub.substr(pos+1));
cout << "Parameters: " << sub << endl;
for (int i=0;;++i){
if (sub=="") break;
pos=sub.find(",");
// cout << "sub " << sub << " " << endl;
if (pos==string::npos){
cout << sub << ": " << m.getParameter(i);
break;
}
cout << sub.substr(0, pos) << ": " << m.getParameter(i) << ' ' << endl;
sub=sub.substr(pos+1);
}
cout << endl;
}
template <class TMap>
void SetParameters(int Version, double x, TMap & m){
//Version==Version3D, 3D system, "par:C1,C2,L,G,Ga,g,I0,g0,V0;var:V1,V2,IL;"
//Version==Version5D, 5D system, "par:C1,C2,L,G,Ga,I0,g0,V0,gamma0,beta,g1,Eag,Ear,q,a,Vg0,Lg,Cth,tth,T0,K;var:V1,V2,IL,g,T;fun:"
//Version==Version1DMap, 1D function, "par:G,Ga,I0,g0,V0,g;var:V1;fun:(Ga-G)*V1-I0*exp(-g/g0)*0.5*(exp(V1/V0)-exp(-V1/V0));"
//Version==Version4D, 5D system with fixed g, "par:C1,C2,L,G,Ga,I0,g0,V0,Cth,tth,T0,g;var:V1,V2,IL,T;fun:"
typedef typename TMap::ScalarType TScalar;
m.setParameter("G",TScalar(Opt.Gr));
m.setParameter("Ga",TScalar(Opt.Gar));
m.setParameter("I0",TScalar(Opt.I0r));
m.setParameter("g0",TScalar(Opt.g0r));
m.setParameter("V0",TScalar(Opt.V0r));
if (Version==Version3D || Version==Version5D || Version==Version4D){
m.setParameter("C1",TScalar(Opt.C1r));
m.setParameter("C2",TScalar(Opt.C2r));
m.setParameter("L",TScalar(Opt.Lr));
}
if (Version==Version5D || Version==Version4D){
m.setParameter("Cth",TScalar(Opt.Cthr));
m.setParameter("tth",TScalar(Opt.tthr));
m.setParameter("T0",TScalar(Opt.T0r));
}
if (Version==Version3D || Version==Version1DMap || Version==Version4D){
m.setParameter("g",TScalar(Opt.gr));
if (Opt.Verbose) WriteParameters(VectorFieldString(Version),m);
return;
}
if (Version==Version5D){
m.setParameter("gamma0",TScalar(Opt.gamma0r));
m.setParameter("beta",TScalar(Opt.betar));
m.setParameter("g1",TScalar(Opt.g1r));
m.setParameter("Eag",TScalar(Opt.Eagr));
m.setParameter("Ear",TScalar(Opt.Earr));
m.setParameter("q",TScalar(Opt.qr));
m.setParameter("a",TScalar(Opt.ar));
m.setParameter("Vg0",TScalar(Opt.Vg0r));
m.setParameter("Lg",TScalar(Opt.Lgr));
m.setParameter("K",TScalar(Opt.Kr));
if (Opt.Verbose) WriteParameters(VectorFieldString(Version),m);
return;
}
cout << "Unknown version" << Version << endl;
}
template <class TMap>
void SetParameters(int Version, MpFloat x, TMap & m){
//Version==Version3D, 3D system, "par:C1,C2,L,G,Ga,g,I0,g0,V0;var:V1,V2,IL;"
//Version==Version5D, 5D system, "par:C1,C2,L,G,Ga,I0,g0,V0,gamma0,beta,g1,Eag,Ear,q,a,Vg0,Lg,Cth,tth,T0,K;var:V1,V2,IL,g,T;fun:"
//Version==Version1DMap, 1D function, "par:G,Ga,I0,g0,V0,g;var:V1;fun:(Ga-G)*V1-I0*exp(-g/g0)*0.5*(exp(V1/V0)-exp(-V1/V0));"
//Version==Version4D, 5D system with fixed g, "par:C1,C2,L,G,Ga,I0,g0,V0,Cth,tth,T0,g;var:V1,V2,IL,T;fun:"
typedef typename TMap::ScalarType TScalar;
m.setParameter("G",TScalar(Opt.Gs));
m.setParameter("Ga",TScalar(Opt.Gas));
m.setParameter("I0",TScalar(Opt.I0s));
m.setParameter("g0",TScalar(Opt.g0s));
m.setParameter("V0",TScalar(Opt.V0s));
if (Version==Version3D || Version==Version5D || Version==Version4D){
m.setParameter("C1",TScalar(Opt.C1s));
m.setParameter("C2",TScalar(Opt.C2s));
m.setParameter("L",TScalar(Opt.Ls));
}
if (Version==Version5D || Version==Version4D){
m.setParameter("Cth",TScalar(Opt.Cths));
m.setParameter("tth",TScalar(Opt.tths));
m.setParameter("T0",TScalar(Opt.T0s));
}
if (Version==Version3D || Version==Version1DMap || Version==Version4D){
m.setParameter("g",TScalar(Opt.gs));
if (Opt.Verbose) WriteParameters(VectorFieldString(Version),m);
return;
}
if (Version==Version5D){
m.setParameter("gamma0",TScalar(Opt.gamma0s));
m.setParameter("beta",TScalar(Opt.betas));
m.setParameter("g1",TScalar(Opt.g1s));
m.setParameter("Eag",TScalar(Opt.Eags));
m.setParameter("Ear",TScalar(Opt.Ears));
m.setParameter("q",TScalar(Opt.qs));
m.setParameter("a",TScalar(Opt.as));
m.setParameter("Vg0",TScalar(Opt.Vg0s));
m.setParameter("Lg",TScalar(Opt.Lgs));
m.setParameter("K",TScalar(Opt.Ks));
if (Opt.Verbose) WriteParameters(VectorFieldString(Version),m);
return;
}
cout << "Unknown version" << Version << endl;
}
template <typename TScalar>
struct DynSyst{
typedef capd::vectalg::Matrix<TScalar,0,0> TMatrix;
typedef capd::map::Map<TMatrix> TMap;
typedef capd::dynsys::BasicOdeSolver<TMap> TOdeSolver;
typedef capd::poincare::TimeMap<TOdeSolver> TTimeMap;
typedef capd::poincare::CoordinateSection<TMatrix> TCoordinateSection;
typedef capd::poincare::BasicPoincareMap<TOdeSolver> TPoincareMap;
int PoincWhich;
int PoincVar;
double PoincVal;
TMap m;
TOdeSolver solver;
TTimeMap tm;
TCoordinateSection section;
TPoincareMap pm;
DynSyst(int Order, double TimeStep, double MaxTimeStep, double AbsTol, double RelTol, int aPoincWhich, int aPoincVar, double aPoincVal):
PoincWhich(aPoincWhich),
PoincVar(aPoincVar),
PoincVal(aPoincVal),
m(VectorFieldString(Opt.CircuitVersion)),
solver(m,Order),
tm(solver),
section(m.dimension(),PoincVar,PoincVal), // the Poincare section is y=PoincVal, i.e. index PoincVar of {0,1,2,...,m.dimension()-1}
pm(solver,section,poincare::CrossingDirection(PoincWhich)){
SetParameters(Opt.CircuitVersion,TScalar(0.0),m);
tm.stopAfterStep(false);
if (AbsTol>0.0) solver.setAbsoluteTolerance(AbsTol);
if (RelTol>0.0) solver.setRelativeTolerance(RelTol);
if (Opt.MinTimeStep>0) solver.setStepControl(capd::dynsys::DLastTermsStepControl(2,Opt.MinTimeStep));
if (MaxTimeStep>0) solver.setMaxStep(MaxTimeStep);
if (TimeStep>0) solver.setStep(TimeStep);
if (Opt.MaxReturnTime>0) pm.setMaxReturnTime(Opt.MaxReturnTime);
};
};
typedef DynSyst<double> DDynSyst;
typedef DynSyst<MpFloat> MpDynSyst;
template <typename TInterval>
struct IDynSyst{
typedef typename TInterval::BoundType TScalar;
typedef capd::vectalg::Matrix<TInterval,0,0> TIMatrix;
typedef capd::map::Map<TIMatrix> TIMap;
typedef capd::dynsys::OdeSolver<TIMap> TITaylor;
typedef capd::dynsys::OdeSolver<TIMap> TIOdeSolver;
typedef capd::poincare::TimeMap<TIOdeSolver> TITimeMap;
typedef capd::poincare::CoordinateSection<TIMatrix> TICoordinateSection;
typedef capd::poincare::PoincareMap<TIOdeSolver> TIPoincareMap;
int PoincWhich;
int PoincVar;
double PoincVal;
TIMap m;
TITaylor solver;
TITimeMap tm;
TICoordinateSection section;
TIPoincareMap pm;
IDynSyst(int Order, double TimeStep, double MaxTimeStep, double AbsTol, double RelTol, int aPoincWhich, int aPoincVar, double aPoincVal):
PoincWhich(aPoincWhich),
PoincVar(aPoincVar),
PoincVal(aPoincVal),
m(VectorFieldString(Opt.CircuitVersion)),
solver(m,Order),
tm(solver),
section(m.dimension(),PoincVar,TInterval(PoincVal)), // the Poincare section is x[PoincVar]=PoincVal
pm(solver,section,poincare::CrossingDirection(PoincWhich)){
SetParameters(Opt.CircuitVersion,TScalar(0.0),m);
tm.stopAfterStep(false);
if (AbsTol>0.0) solver.setAbsoluteTolerance(AbsTol);
if (RelTol>0.0) solver.setRelativeTolerance(RelTol);
if (MaxTimeStep>0) solver.setMaxStep(TScalar(MaxTimeStep));
if (TimeStep>0) solver.setStep(TScalar(TimeStep));
if (Opt.MaxReturnTime>0) pm.setMaxReturnTime(Opt.MaxReturnTime);
};
};
typedef IDynSyst<Interval> DIDynSyst;
typedef IDynSyst<MpInterval> MpIDynSyst;
template <class TVector, class TScalar>
bool ComputeTrajectory(double TimeStep, double SkipTime, double IntegrationTime, TVector x, vector <TVector> & xres, vector <TScalar> & tres){
DynSyst<TScalar> syst(Opt.Order,TimeStep,Opt.MaxTimeStep,Opt.AbsTol, Opt.RelTol, Opt.PoincWhich, Opt.PoincVar, Opt.PoincVal);
TVector Fx=syst.m(x);
cout << "x: " << x << " Fx " << Fx << endl;
TScalar t=0.0;
cout << "t: " << t << " x: " << x << " F(x): " << syst.m(x) << endl;
TScalar finalTime=SkipTime;
syst.tm.stopAfterStep(false);
do{
x=syst.tm(finalTime,x,t);
} while(!syst.tm.completed());
syst.tm.stopAfterStep(true);
t=0.0;
finalTime=IntegrationTime;
xres.clear();
tres.clear();
do{
xres.push_back(x);
tres.push_back(t);
x = syst.tm(finalTime,x,t);
cout << "t: " << t << "x: " << x << " F(x): " << syst.m(x) << endl;
// cout << t << " " << x << endl;
} while(!syst.tm.completed());
return true;
}
bool Trajectory(double TimeStep, double SkipTime, double IntegrationTime){
cout << "Trajectory" << " TimeStep: " << TimeStep << endl;
DVector x;
if (!InitialConditions(Opt.CircuitVersion,Opt.InitCondStr,x)) return false;
vector <DVector> xres;
vector <double> tres;
if (!ComputeTrajectory(TimeStep,SkipTime,IntegrationTime,x,xres,tres)){
cerr << "Trajectory: Error" << endl;
return false;
}
for (int p=0;p<xres.size();++p){
cout << "P: " << tres[p];
for (int i=0;i<xres[p].dimension();++i) cout << " " << xres[p][i];
cout << endl;
}
return true;
}
bool TrajectoryMp(double TimeStep, double SkipTime, double IntegrationTime){
cout << "TrajectoryMp" << " TimeStep: " << TimeStep << endl;
MpVector x;
if (!InitialConditions(Opt.CircuitVersion,Opt.InitCondStr,x)) return false;
vector <MpVector> xres;
vector <MpFloat> tres;
if (!ComputeTrajectory(TimeStep,SkipTime,IntegrationTime,x,xres,tres)){
cerr << "Trajectory: Error" << endl;
return false;
}
for (int p=0;p<xres.size();++p){
cout << "P: " << tres[p];
for (int i=0;i<xres[p].dimension();++i) cout << " " << xres[p][i];
cout << endl;
}
return true;
}
bool ComputeTrajectoryInterval(double TimeStep, double IntegrationTime){
try{
DIDynSyst systi(Opt.Order,Opt.TimeStep,Opt.MaxTimeStep,Opt.AbsTol, Opt.RelTol, Opt.PoincWhich, Opt.PoincVar, Opt.PoincVal);
DVector x0;
if (!InitialConditions(Opt.CircuitVersion,Opt.InitCondStr,x0)) return false;
IVector c(x0.dimension());
for (int d=0;d<x0.dimension();++d) c[d]=x0[d];
// take some box around c
c += interval(-1,1)*Opt.radius;
// define a doubleton representation of the interval vector c
C0HORect2Set s(c);
// we integrate the set s over the time IntegrationTime
interval T(IntegrationTime);
cout << "InitialSet: " << c << " diameter: " << diam(c) << endl;
IVector result;
if (Opt.WriteTrajectory){
systi.tm.stopAfterStep(true);
interval time_done(0.0);
do{
result = systi.tm(IntegrationTime,s);
interval stepMade = systi.solver.getStep();
time_done+=stepMade;
IVector resultlin=Opt.CoordChange*result;
Interval r=sqrt(sqr(resultlin[1])+sqr(resultlin[2]));
cout << "P: " << time_done << " " << result[0] << " " << result[1] << " " << result[2] << " r: " << r << endl;
} while(!systi.tm.completed());
} else {
systi.tm.stopAfterStep(false);
result = systi.tm(IntegrationTime,s);
}
cout << "IntegrationTime: " << T << endl;
cout << "Image: " << result << " diameter: " << diam(result) << endl;
}catch(exception& e)
{
cout << "Exception caught!" << e.what() << endl;
}
return true;
}
template <class DynSyst, class TVector, class TScalar>
bool PoincMapDS(DynSyst & ds, const TVector & x, TVector & y, TScalar & t){
t=0.0;
try{
y=ds.pm(x,t);
y[ds.PoincVar]=ds.PoincVal; // !!! assumption: poincare map is defined by the plane x[Opt.PoincVar]=Opt.PoincVal
} catch(exception& e) {
// cout << "PoincMap: Exception caught!" << e.what() << endl;
if (Opt.Verbose) cout << "Poincmap Error: x: " << x << endl;
return false;
}
return true;
}
template <class DynSyst, class TMatrix, class TVector, class TScalar>
bool PoincMapDS(DynSyst & ds, const TVector & x, TVector & y, TMatrix & DP, TScalar & t){
t=0.0;
try{
y=ds.pm(x,DP,t);
DP=ds.pm.computeDP(y,DP);
y[ds.PoincVar]=ds.PoincVal; // !!! assumption: poincare map is defined by the plane x[Opt.PoincVar]=Opt.PoincVal
} catch(exception& e) {
cout << "Exception caught!" << e.what() << endl;
return false;
}
return true;
}
bool PoincMapDS(IPoincareMap & pm, const IVector & x, IVector & y, Interval & t){
t=0.0;
try{
C0HOTripletonSet s1(x);
// C0Rect2Set s1(x); // works worse than C0HOTripletonSet
y=pm(s1,t);
// y[ds.PoincVar]=ds.PoincVal; // !!! we assume that poincare map is defined by the plane x=r
} catch(exception& e) {
cout << "Exception caught!" << e.what() << endl;
return false;
}
return true;
}
bool PoincMapDS(DIDynSyst & ds, const IVector & x, IVector & y, Interval & t){
t=0.0;
try{
C0HOTripletonSet s1(x);
// C0Rect2Set s1(x); // works worse than C0HOTripletonSet
y=ds.pm(s1,t);
y[ds.PoincVar]=ds.PoincVal; // !!! we assume that poincare map is defined by the plane x=r
} catch(exception& e) {
cout << "Exception caught!" << e.what() << endl;
return false;
}
return true;
}
bool PoincMapDS(MpIDynSyst & ds, const MpIVector & x, MpIVector & y, MpInterval & t){
t=0.0;
try{
// typedef capd::dynset::C0HOSet<MpC0TripletonSet> capd::MpC0HOTripletonSet;
MpC0Rect2Set s1(x);
// MpC0HOTripletonSet s1(x);
y=ds.pm(s1,t);
y[ds.PoincVar]=ds.PoincVal; // !!! we assume that poincare map is defined by the plane x=r
} catch(exception& e) {
cout << "Exception caught!" << e.what() << endl;
return false;
}
return true;
}
bool PoincMapDS(DIDynSyst & ds, const IVector & x, IVector & y, IMatrix & DP, Interval & t){
t=0.0;
try{
C1HORect2Set s2(x);
// typedef capd::dynset::C1HOSet<capd::C1Rect2Set> capd::C1HORect2Set
// C1Rect2Set s2(x);
y = ds.pm(s2,DP,t);
DP = ds.pm.computeDP(y,DP);
y[ds.PoincVar]=ds.PoincVal; // !!! we assume that poincare map is defined by the plane x=r
} catch(exception& e) {
cout << "Exception caught!" << e.what() << endl;
return false;
}
return true;
}
bool PoincMapDS(MpIDynSyst & ds, const MpIVector & x, MpIVector & y, MpIMatrix & DP, MpInterval & t){
t=0.0;
try{
MpC1Rect2Set s2(x);
// typedef capd::dynset::C1HOSet<MpC1Rect2Set> capd::MpC1HORect2Set
// MpC1HORect2Set s2(x); // compilation error
y=ds.pm(s2,DP,t);
DP=ds.pm.computeDP(y,DP);
y[ds.PoincVar]=ds.PoincVal; // !!! we assume that poincare map is defined by the plane x=r
} catch(exception& e) {
cout << "Exception caught!" << e.what() << endl;
return false;
}
return true;
}
bool TestPoincMapI(){
DIDynSyst isyst(Opt.Order,Opt.TimeStep,Opt.MaxTimeStep,Opt.AbsTol,Opt.RelTol,Opt.PoincWhich,Opt.PoincVar,Opt.PoincVal);
DVector x0;
if (!InitialConditions(Opt.CircuitVersion,Opt.InitCondStr,x0)) return false;
DVector xr=x0;
IVector x(x0.dimension());
for (int d=0;d<x0.dimension();++d) x[d]=x0[d];
// take some box around c
x += interval(-1,1)*Opt.radius;
IVector y;
IMatrix DP;
Interval t;
PoincMapDS(isyst,x,y,DP,t);
cout << "x: " << x << " diam: " << diam(x) << endl;
cout << "y: " << y << " diam: " << diam(y) << endl;
cout << "t: " << t << " diam: " << diam(t) << endl;
return true;
}
template <class DynSyst, class TVector, class TScalar>
bool ComputeIterates(DynSyst & ds, TVector & x, int SkipNum, int IterNum, vector <TVector> & xres, vector <TScalar> & tres){
xres.clear();
tres.clear();
TVector y;
TScalar t=0.0;
TScalar tfull=0.0;
bool StopOnPoincMapFailure=false;
for (int iter=0;iter<SkipNum;++iter){
if (!PoincMapDS(ds,x,y,t)){
if (StopOnPoincMapFailure) return false;
ds.tm.stopAfterStep(true);
do{
y=ds.tm(Opt.MaxReturnTime,x,t);
// cout << t << " " << x << endl;
x=y;
} while(!ds.tm.completed());
// cout << y << " " << t << endl;
}
x=y;
tfull+=t;
ds.tm.stopAfterStep(false);
}
for (int iter=0;iter<IterNum;++iter){
if (!PoincMapDS(ds,x,y,t)){
if (StopOnPoincMapFailure) return false;
// cout << "ab " << iter << endl;
ds.tm.stopAfterStep(true);
do{
y=ds.tm(Opt.MaxReturnTime,x,t);
// cout << t << " " << y << endl;
x=y;
} while(!ds.tm.completed());
// cout << y << " " << t << endl;
}
x=y;
tfull+=t;
ds.tm.stopAfterStep(false);
xres.push_back(x);
tres.push_back(tfull);
}
return true;
}
template <class DynSyst, class TVector>
bool ComputeIterates(DynSyst & ds, TVector & x, int SkipNum, int IterNum, TVector & popos, TVector & tvec){
vector <TVector> xres;
vector <typename TVector::ScalarType> tres;
if (!ComputeIterates(ds,x,SkipNum,IterNum,xres,tres)) return false;
popos.resize(x.dimension()*IterNum);
tvec.resize(IterNum);
for (int p=0;p<IterNum;++p){
for (int i=0;i<x.dimension();++i) popos[p*x.dimension()+i]=xres[p][i];
tvec[p]=tres[p];
}
return true;
}
bool ReturnMapTrajectory(double TimeStep, int SkipNum, int IterNum, bool WriteTrajectory, bool WritePar){
vector <DVector> xres;
vector <double> tres;
DVector x;
if (!InitialConditions(Opt.CircuitVersion,Opt.InitCondStr,x)) return false;
DDynSyst syst(Opt.Order,TimeStep,Opt.MaxTimeStep,Opt.AbsTol, Opt.RelTol, Opt.PoincWhich, Opt.PoincVar, Opt.PoincVal);
if (!ComputeIterates(syst,x,SkipNum,IterNum,xres,tres)){
cerr << "ReturnMapTrajectory: Error" << endl;
return false;
}
if (xres.size()==0) return false;
for (int p=0;p<xres.size();++p){
if (!WriteTrajectory) continue;
cout << "P: " ;
if (WritePar) {
if (Opt.BiffParName=="g") cout << Opt.gr << " ";
if (Opt.BiffParName=="G") cout << Opt.Gr << " ";
}
cout << p+SkipNum+1 << " " << tres[p] << " " << xres[p][0] << " " << xres[p][1] << " " << xres[p][2];
cout << endl;
}
if (Opt.BiffChangeStart) Opt.InitCondStr=VectorToString(xres[xres.size()-1]);
// if (Opt.BiffChangeStart && Opt.CircuitVersion==Version3D) Opt.SMCCinitcond3=xres[xres.size()-1];
// if (Opt.BiffChangeStart && Opt.CircuitVersion==Version5D) Opt.SMCCinitcond5=xres[xres.size()-1];
if (Opt.OutFile!="") VectorSetSave(xres,Opt.OutFile,Opt.CoutPrec);
return true;
}
bool MpReturnMapTrajectory(double TimeStep, int SkipNum, int IterNum, bool WriteTrajectory, bool WritePar){
DVector x0;
if (!InitialConditions(Opt.CircuitVersion,Opt.InitCondStr,x0)) return false;
MpVector x(x0.dimension());
for (int d=0;d<x0.dimension();++d) x[d]=x0[d];
vector <MpVector> xres;
vector <MpFloat> tres;
MpDynSyst syst(Opt.Order,TimeStep,Opt.MaxTimeStep,Opt.AbsTol,Opt.RelTol,Opt.PoincWhich,Opt.PoincVar,Opt.PoincVal);
if (!ComputeIterates(syst,x,SkipNum,IterNum,xres,tres)){
cerr << "ReturnMapTrajectory: Error" << endl;
return false;
}
for (int p=0;p<xres.size();++p){
if (!WriteTrajectory) continue;
cout << "P: " ;
if (WritePar) {
if (Opt.BiffParName=="g") cout << Opt.gr << " ";
if (Opt.BiffParName=="G") cout << Opt.Gr << " ";
}
cout << p+SkipNum+1 << " " << tres[p] << " " << xres[p][0] << " " << xres[p][1] << " " << xres[p][2];
cout << endl;
}
return true;
}
int POMinPeriod(const DVector & xfull, int period, double MaxAllowedDist){
int dim=xfull.dimension()/period;
DVector x0(dim);
for (int i=0;i<dim;++i) x0[i]=xfull[i];
for (int p=1;p<period;++p){
DVector x(dim);
for (int i=0;i<dim;++i) x[i]=xfull[p*dim+i];
if (MaxDist(x0,x)<MaxAllowedDist) return p;
}
return period;
}
template <class TMatrix, class TVector>
bool ComputeNewtonCorrection(const TMatrix & J, const TVector & F, TVector & h){
try{
h=capd::matrixAlgorithms::gauss(J,F);
}
catch(exception& e){
if (Opt.Verbose) cout << e.what() << endl;
return false;
}
return true;
}
template <class TIMatrix, class TIVector>
bool ComputeKrawczykCorrection(const TIVector & xfull, const TIVector & xfullcenter, const TIMatrix & J, const TIVector & F, TIVector & h){
typedef typename TIMatrix::ScalarType TInterval;
typedef typename TInterval::BoundType TScalar;
typedef capd::vectalg::Matrix<TScalar,0,0> TMatrix;
int dim=J.numberOfRows();
TIMatrix I=TIMatrix::Identity(dim);
try{
TMatrix JScalar=Mid(J);
TMatrix CScalar=matrixAlgorithms::inverseMatrix(JScalar);
TIMatrix C=Convert(CScalar);
h=-C*F+(I-C*J)*(xfull-xfullcenter);
}
catch(exception& e){
if (Opt.Verbose) cout << e.what() << endl;
return false;
}
return true;
}
template <class DynSyst, class TVector, class TScalar, class TMatrix>
bool GlobMap(DynSyst & ds, TVector & XFull, TScalar & return_time, TVector & F, TMatrix & J, TMatrix & J2Full){
int dim=ds.m.dimension();
int period=XFull.dimension()/dim;
F.resize((dim-1)*period);
J.resize((dim-1)*period,(dim-1)*period);
// f_k=P(x_{k})-x_{k+1}
return_time=0.0;
J2Full=TMatrix::Identity(dim-1);
for (int p=0;p<period;++p){
TVector x(dim);
XFull[p*dim+Opt.PoincVar]=Opt.PoincVal;
for (int i=0;i<dim;++i) x[i]=XFull[p*dim+i];
TScalar t=0.0;
TVector y(dim,0.0);
TMatrix DP(dim,dim);
if (!PoincMapDS(ds,x,y,DP,t)) return false;
return_time+=t;
// update F
int pos=(p-1+period)%period;
TVector x2=VectorRemovePos(x,Opt.PoincVar);
TVector y2=VectorRemovePos(y,Opt.PoincVar);
for (int i=0;i<dim-1;++i) F[(dim-1)*p+i]+=y2[i];
for (int i=0;i<dim-1;++i) F[(dim-1)*pos+i]-=x2[i];
// update Jacobian matrix
TMatrix J2=MatrixRemoveRowCol(DP,Opt.PoincVar,Opt.PoincVar);
J2Full=J2*J2Full;
for (int i1=0;i1<dim-1;++i1){
for (int i2=0;i2<dim-1;++i2) J[(dim-1)*p+i1][(dim-1)*p+i2]+=J2[i1][i2];
}
for (int i=0;i<dim-1;++i) J[(dim-1)*pos+i][(dim-1)*p+i]-=1.0;
}
return true;
}
template <class DynSyst, class TScalar, class TVector, class TMatrix>
bool NewtonReal(DynSyst & ds, TVector & popos, TVector & xres, TScalar & t, TVector & F, TMatrix & J, TVector & h, TMatrix & J2Full, bool verbose = false){
// apply the Newton operator to approximate zeros of f defined by f_k=P(x_{k})-x_{k+1}
int dim=ds.m.dimension();
int period=popos.dimension()/dim;
if (!GlobMap(ds,popos,t,F,J,J2Full)) return false;
if (!ComputeNewtonCorrection(J,F,h)) return false;
xres=popos;
for (int p=0;p<period;++p){
for (int i=0,ip=0;i<dim;++i){
if (i==Opt.PoincVar) continue;
xres[dim*p+i]-=h[(dim-1)*p+ip];
ip++;
}
}
return true;
}
template <class IDynSyst, class TIVector, class TIMatrix>
bool NewtonInterval(IDynSyst & ids, TIVector & popos, typename TIVector::ScalarType & return_time, TIVector & N, TIVector & x2full, TIMatrix & J2Full){
typedef typename TIVector::ScalarType TScalar;
int dim=ids.m.dimension();
int period=popos.dimension()/dim;
x2full.resize((dim-1)*period);
for (int p=0;p<period;++p){
TIVector x(dim);
for (int i=0;i<dim;++i) x[i]=popos[dim*p+i];
TIVector x2=VectorRemovePos(x,Opt.PoincVar);
for (int i=0;i<dim-1;++i) x2full[(dim-1)*p+i]=x2[i]; // skip one dimension
}
J2Full=TIMatrix::Identity(dim-1);
TIVector x2fullcenter=midVector(x2full);
TIMatrix J((dim-1)*period,(dim-1)*period,0.0);
TIVector F((dim-1)*period,0.0);
// f_k=P(x_{k})-x_{k+1}
return_time=0.0;
for (int p=0;p<period;++p){
TIVector x=VectorPart(popos,dim,p*dim);
TIVector xcenter = midVector(x);
typename TIVector::ScalarType t(0.0);
TIVector ycenter;
if (!PoincMapDS(ids,xcenter,ycenter,t)) return false;
TIMatrix DP(dim,dim);
TIVector y(dim);
if (!PoincMapDS(ids,x,y,DP,t)) return false;
return_time+=t;
// update F
int pos=(p-1+period)%period;
TIVector x2center=VectorRemovePos(xcenter,Opt.PoincVar);
TIVector y2center=VectorRemovePos(ycenter,Opt.PoincVar);
for (int i=0;i<dim-1;++i) F[(dim-1)*p+i]+=y2center[i];
for (int i=0;i<dim-1;++i) F[(dim-1)*pos+i]-=x2center[i];
// update Jacobian matrix
TIMatrix J2=MatrixRemoveRowCol(DP,Opt.PoincVar,Opt.PoincVar);
J2Full=J2*J2Full;
for (int i1=0;i1<dim-1;++i1){
for (int i2=0;i2<dim-1;++i2) J[(dim-1)*p+i1][(dim-1)*p+i2]+=J2[i1][i2];
}
// TScalar One(1.0);
for (int i=0;i<dim-1;++i) J[(dim-1)*pos+i][(dim-1)*p+i]-=TScalar(1.0);
}
TIVector h;
// if (Opt.IntervalOperator==NewtonOperator){
if (!ComputeNewtonCorrection(J,F,h)) return false;
N=x2fullcenter-h;
// }
/*
if (Opt.IntervalOperator==KrawczykOperator){
if (!ComputeKrawczykCorrection(x2full,x2fullcenter,J,F,h)) return false;
N=x2fullcenter+h;
}
*/
return true;
}
template <class DynSyst, class TVector, class TScalar, class TMatrix>
bool NewtonRealIterate(DynSyst & ds, TVector & popos, TScalar & t, TScalar & MaxF, TMatrix & J2Full, bool verbose = false){
// apply the Newton method to find an approximate position of the periodic orbit close to popos
// periodic orbits are zeros of f defined by f_k=P(x_{k})-x_{k+1}
// typedef typename TVector::ScalarType TScalar;
int dim=ds.m.dimension();
int period=popos.dimension()/dim;
TScalar lastfmax=1e100;
for (int iter=0;iter<Opt.NewtonMaxIterNum;++iter){
if (verbose) cout << "Newtoniter: " << iter << endl;
TVector xres,F,h;
TMatrix J,J2Full;
if (!NewtonReal(ds,popos,xres,t,F,J,h,J2Full,verbose)){
cout << "NewtonPO error" << endl;
return false;
}
popos=xres;
TScalar fmax=VectorMaxAbsVal(F);
TScalar maxh=VectorMaxAbsVal(h);
if (verbose) cout << "fmax: " << fmax << " maxh: " << maxh << endl;
if (fmax<MaxPODist(fmax)){
int min_period=POMinPer(popos,period,Opt.MinPOSelfDist);
if (min_period<period){
if (verbose) cout << "Newton stopped" << " period: " << period << " minperiod: " << min_period << endl;
return false;
}
MaxF=fmax;
break;
}
if (maxh>Opt.NewtonMaxCorrection){
if (verbose) cout << "Newton stopped" << " maxh: " << maxh << " > " << Opt.NewtonMaxCorrection << endl;
return false;
}
if (fmax>0.5*lastfmax && fmax>100*MaxPODist(fmax)){
if (verbose) cout << "Newton stopped" << " fmax: " << fmax << " " << lastfmax << " " << iter << endl;
return false;
}
lastfmax=fmax;
}
return true;
}
template <class IDynSyst, class TVector, class TIVector, class TInterval>
bool ProveNewton(IDynSyst & ids, TVector & xfullr, TIVector & xfull, TInterval & return_time, int & isstable){
typedef capd::vectalg::Matrix<TInterval,0,0> TIMatrix;
int dim=ids.m.dimension();
int period=xfullr.dimension()/dim;
if (Opt.Verbose) cout << "xfullr: " << xfullr << endl;
xfull.resize(xfullr.dimension());
for (int i=0;i<xfull.dimension();++i) xfull[i]=xfullr[i];
for (int iter=0;iter<Opt.NewtonMaxIterNum;++iter){
TIVector N;
TIVector x2full((dim-1)*period);
TIMatrix J2Full=TIMatrix::Identity(dim-1);
if (!NewtonInterval(ids,xfull,return_time,N,x2full,J2Full)) return false;
if(subsetInterior(N,x2full)){
if (Opt.Verbose) cout << "xfullr: " << xfullr << endl;
if (Opt.Verbose) cout << "xfull: " << xfull << endl;
if (Opt.Verbose) cout << "iter: " << iter << " proof completed" << endl;
isstable=stablemaybe;
bool CheckStability=true;
if (dim==3 && CheckStability){
bool stable=true;
bool unstable=false;
TInterval determinant=J2Full(1,1)*J2Full(2,2)-J2Full(1,2)*J2Full(2,1);
TInterval trace=J2Full(1,1)+J2Full(2,2);
TInterval a1=-trace;
TInterval a0=determinant;
TInterval t1=a0+a1+1;
TInterval t2=a0-a1+1;
if (a0.leftBound()>1) unstable=true;
if (t1.rightBound()<0) unstable=true;
if (t2.rightBound()<0) unstable=true;
if (a0.rightBound()>1) stable=false;
if (t1.leftBound()<0) stable=false;
if (t2.leftBound()<0) stable=false;
if (Opt.Verbose) cout << "a0: " << a0 << " a1: " << a1 << " t1: " << t1 << " t2: " << t2 << endl;
isstable=stablemaybe;
if (stable==true && unstable==false) isstable=stableyes;
if (stable==false && unstable==true) isstable=stableno;
TInterval delta=sqr(a1)-4*a0;
if (Opt.Verbose) cout << "delta: " << delta << endl;
TInterval l1,l2;
if (delta>=0){
if (Opt.Verbose) cout << "Positive delta: " << delta << endl;
l1=(-a1-sqrt(delta))/2;
l2=(-a1+sqrt(delta))/2;
if (Opt.Verbose) cout << "Real: delta: " << delta << " l1: " << l1 << " l2: " << l2 << endl;
} else {
if (delta<=0){
if (Opt.Verbose) cout << "Negative delta: " << delta << endl;
l2=l1=sqrt(a0); // absolute values of complex eigenvalues
TInterval omega;
omega=sqrt(-delta)/2;
if (Opt.Verbose) cout << "Complex: delta: " << delta << " l1: " << l1 << " omega: " << omega << endl;
} else {
if (Opt.Verbose) cout << "PositiveOrNegative delta: " << delta << endl;
TInterval deltaplus=TInterval(0,delta.rightBound());
TInterval deltaminus=TInterval(delta.leftBound(),0);
TInterval l1plus=(-a1-sqrt(deltaplus))/2;
TInterval l2plus=(-a1+sqrt(deltaplus))/2;
if (Opt.Verbose) cout << "Real: delta: " << delta << " l1+: " << l1plus << " l2-: " << l2plus << endl;
TInterval a0plus=TInterval(0,a0.rightBound());
TInterval lminus=sqrt(a0plus); // absolute values of complex eigenvalues
TInterval omega;
omega=sqrt(-deltaminus)/2;
if (Opt.Verbose) cout << "Complex: delta: " << delta << " l-: " << lminus << " omega: " << omega << endl;
l1=intervalHull(abs(l1plus),abs(l2plus));
l1=intervalHull(l1,lminus);
l2=l1;
if (Opt.Verbose) cout << "RealComplex: delta: " << delta << " labs: " << l1 << endl;
}
}
}
return true;
}
if (Opt.Verbose) cout << "iter: " << iter << " N is not a subset of X: Diam(x): " << VectorMaxDiam(x2full) << " Diam(N): " << VectorMaxDiam(N) << " inflating" << endl;
if (Opt.Verbose) cout << "x2full: " << x2full << endl;
if (Opt.Verbose) cout << "N: " << N << endl;
for (int p=0;p<period;++p){
typename TInterval::BoundType eps=0.1;
// copy and inflate
for (int i=0,ip=0;i<dim;++i){
if (i==Opt.PoincVar) continue;
xfull[dim*p+i]=(1+eps)*N[(dim-1)*p+ip]-eps*N[(dim-1)*p+ip];
ip++;
}
}
}
return false;
}
bool FindPO(int period){
DVector x;
if (!InitialConditions(Opt.CircuitVersion,Opt.InitCondStr,x)) return false;
DDynSyst syst(Opt.Order,Opt.TimeStep,Opt.MaxTimeStep,Opt.AbsTol,Opt.RelTol,Opt.PoincWhich,Opt.PoincVar,Opt.PoincVal);
cout << "dimension: " << syst.m.dimension() << endl;
cout << "g: " << syst.m.getParameter("g") << endl;
DVector PO;
DVector tvec;
if (!ComputeIterates(syst,x,0,period,PO,tvec)) return false;
double t(0.0);
for (int i=0;i<tvec.dimension();++i) t+=tvec[i];
cout << "PO: " << period << " " << t << " " << PO << endl;
if (Opt.NewtonImproveApprox){
DMatrix J2Full;
double MaxF(0.0);
if (NewtonRealIterate(syst,PO,t,MaxF,J2Full,Opt.Verbose)){
cout << "PO: " << period << " " << t << " " << PO << " MaxF: " << MaxF << endl;
} else {
cout << "ConvError " << period << " " << t << " " << x << endl;
return false;
}
}
return true;
}
bool FindPOMp(int period){
DVector x0;
if (!InitialConditions(Opt.CircuitVersion,Opt.InitCondStr,x0)) return false;
MpDynSyst syst(Opt.Order,Opt.TimeStep,Opt.MaxTimeStep,Opt.AbsTol,Opt.RelTol,Opt.PoincWhich,Opt.PoincVar,Opt.PoincVal);
cout << "dimension: " << syst.m.dimension() << endl;
cout << "g: " << syst.m.getParameter("g") << endl;
MpVector x(x0.dimension());
for (int d=0;d<x0.dimension();++d) x[d]=x0[d];
MpVector PO;
MpVector tvec;
if (!ComputeIterates(syst,x,0,period,PO,tvec)) return false;
MpFloat t(0.0);
for (int i=0;i<tvec.dimension();++i) t+=tvec[i];
cout << "PO: " << period << " " << t << " " << PO << endl;
if (Opt.NewtonImproveApprox){
MpMatrix J2Full;
MpFloat MaxF(0.0);
if (NewtonRealIterate(syst,PO,t,MaxF,J2Full,Opt.Verbose)){
cout << "PO: " << period << " " << t << " " << PO << " MaxF: " << MaxF << endl;
} else {
if (Opt.Verbose) cout << "ConvError " << period << " " << t << " " << x << endl;
return false;
}
}
return true;
}
bool ProvePO(int period){
DVector x;
if (!InitialConditions(Opt.CircuitVersion,Opt.InitCondStr,x)) return false;
DDynSyst syst(Opt.Order,Opt.TimeStep,Opt.MaxTimeStep,Opt.AbsTol,Opt.RelTol,Opt.PoincWhich,Opt.PoincVar,Opt.PoincVal);
cout << "dimension: " << syst.m.dimension() << endl;
cout << "g: " << syst.m.getParameter("g") << endl;
DVector PO;
DVector tvec;
if (!ComputeIterates(syst,x,0,period,PO,tvec)) return false;
double t(0.0);
for (int i=0;i<tvec.dimension();++i) t+=tvec[i];
cout << "PO: " << period << " " << t << " " << PO << endl;
if (Opt.NewtonImproveApprox){
DMatrix J2Full;
double MaxF(0.0);
if (NewtonRealIterate(syst,PO,t,MaxF,J2Full,Opt.Verbose)){
cout << "PO: " << period << " " << t << " " << PO << " MaxF: " << MaxF << endl;
} else {
if (Opt.Verbose) cout << "ConvError " << period << " " << t << " " << x << endl;
return false;
}
}
DIDynSyst isyst(Opt.Order,Opt.TimeStep,Opt.MaxTimeStep,Opt.AbsTol,Opt.RelTol,Opt.PoincWhich,Opt.PoincVar,Opt.PoincVal);
IVector xres;
DInterval tbound;
TimeCount TC;
int isstable;
if (ProveNewton(isyst,PO,xres,tbound,isstable)){
cout << "PO " << period << " " << tbound << " " << xres << " Diam: " << VectorMaxDiam(xres) << " " << StabilityResultsInfo[isstable] << " " << " t: " << TC << endl;
return true;
}
cout << "Proof failed" << endl;
return false;
}
bool ProvePOMp(int period){
DVector x0;
if (!InitialConditions(Opt.CircuitVersion,Opt.InitCondStr,x0)) return false;
MpDynSyst syst(Opt.Order,Opt.TimeStep,Opt.MaxTimeStep,Opt.AbsTol,Opt.RelTol,Opt.PoincWhich,Opt.PoincVar,Opt.PoincVal);
cout << "dimension: " << syst.m.dimension() << endl;
cout << "g: " << syst.m.getParameter("g") << endl;
MpVector x(x0.dimension());
for (int d=0;d<x0.dimension();++d) x[d]=x0[d];
MpVector PO;
MpVector tvec;
if (!ComputeIterates(syst,x,0,period,PO,tvec)) return false;
MpFloat t(0.0);
for (int i=0;i<tvec.dimension();++i) t+=tvec[i];
cout << "PO: " << period << " " << t << " " << PO << endl;
if (Opt.NewtonImproveApprox){
MpMatrix J2Full;
MpFloat MaxF(0.0);
if (NewtonRealIterate(syst,PO,t,MaxF,J2Full,Opt.Verbose)){
cout << "PO: " << period << " " << t << " " << PO << " MaxF: " << MaxF << endl;
} else {
if (Opt.Verbose) cout << "ConvError " << period << " " << t << " " << x << endl;
return false;
}
}
MpIDynSyst isyst(Opt.Order,Opt.TimeStep,Opt.MaxTimeStep,Opt.AbsTol,Opt.RelTol,Opt.PoincWhich,Opt.PoincVar,Opt.PoincVal);
MpIVector xres;
MpInterval tbound;
TimeCount TC;
int isstable;
if (ProveNewton(isyst,PO,xres,tbound,isstable)){
cout << "PO " << period << " " << tbound << " " << xres << " Diam: " << VectorMaxDiam(xres) << " " << StabilityResultsInfo[isstable] << " t: " << TC << endl;
return true;
}
cout << "Proof failed" << endl;
return false;
}
bool BiffDiag(){
cout << "PoincWhich: " << Opt.PoincWhich << endl;
for (int p=0;p<Opt.BiffParNum;++p){
double BiffParVal=Opt.BiffParMin+(Opt.BiffParMax-Opt.BiffParMin)*p/(Opt.BiffParNum-1);
if (Opt.BiffParName=="g") Opt.gr=BiffParVal;
if (Opt.BiffParName=="G") Opt.Gr=BiffParVal;
DVector x;
if (!InitialConditions(Opt.CircuitVersion,Opt.InitCondStr,x)) return false;
cout << "BiffPar: " << Opt.BiffParName << " BiffParVal: " << BiffParVal << " PoincVal: " << Opt.PoincVal << " Init: " << x << endl;
ReturnMapTrajectory(Opt.TimeStep,Opt.SkipNum,Opt.IterNum,true,true);
}
return true;
}
void BiffDiagEps(string InFileName, int BinNum, double ParMin, double ParMax, double XMin, double XMax, string EPSFileName, EPSOptions EPSOpt){
ofstream OutFile;
SaveEPS aEPS;
aEPS.open(EPSFileName);
cout << "SaveEPS" << endl;
EPSOpt.EPSFileCreator="smcc.x (Z. Galias)";
ostringstream oss;
oss << "SMCC circuit, bifurcation diagram";
EPSOpt.EPSFileTitle=oss.str();
EPSOpt.XLabel=Opt.BiffParName;
EPSOpt.YLabel="y";
EPSOpt.DefineMinMax(ParMin,ParMax,XMin,XMax);
cout << "Saving " << EPSFileName << endl;
aEPS.Init(EPSOpt);
aEPS.OutFile.precision(15);
aEPS.OutFile << "grestore\n";
// aEPS.OutFile << "0 0 1" << " rgb\n";
// aEPS.FillRectR(EPSOpt.bbXMin-1-EPSOpt.LMarg,EPSOpt.bbYMin-EPSOpt.BMarg,EPSOpt.bbXMax+EPSOpt.RMarg,EPSOpt.bbYMax+1+EPSOpt.TMarg);
aEPS.OutFile << "gsave\n";
aEPS.RescaleToNDigits(EPSOpt.XMin,EPSOpt.XMax,EPSOpt.YMin,EPSOpt.YMax,EPSOpt.accuracy);
aEPS.DefSizesAutomaticOut(EPSOpt.Diameter);
aEPS.DefLineWidthAutomaticOut(EPSOpt.LineWidth);
aEPS.OutFile << "0 0 1" << " rgb\n";
ifstream InFile(InFileName.c_str());
double valold=-1000.0; // some negative value which cannot appears in input data
vector <int> PointsInBin(BinNum,0);
int PWNumFound=0;
int PONumFound=0;
int BinsVisitedNumLast=-1;
while (InFile){
double tmp,val;
vector <double> x(3);
InFile >> val >> tmp >> tmp >> x[0] >> x[1] >> x[2];
if ((val<ParMin) || (val>ParMax)) continue;
if (val!=valold && valold>0.0){
double da=val-valold;
double a=valold;
int PointsInBinMax=0;
int PointNum=0;
for (int b=0;b<BinNum;++b){
PointNum+=PointsInBin[b];
if (PointsInBin[b]>PointsInBinMax) PointsInBinMax=PointsInBin[b];
}
int BinsVisitedNum=0;
for (int b=0;b<BinNum;++b){
if (PointsInBin[b]>0) BinsVisitedNum++;
}
if (BinsVisitedNum<BinNum/10){
PONumFound++;
cout << "PO possible: " << Opt.BiffParName << ": " << valold << " BinsVisitedNum: " << BinsVisitedNum << endl;
if (BinsVisitedNum!=BinsVisitedNumLast){
PWNumFound++;
cout << "PW possible: " << Opt.BiffParName << ": " << valold << " BinsVisitedNum: " << BinsVisitedNum << endl;
}
}
BinsVisitedNumLast=BinsVisitedNum;
for (int b=0;b<BinNum;++b){
if (PointsInBin[b]==0) continue;
double val1=((XMax-XMin)*b)/BinNum+XMin;
double val2=((XMax-XMin)*(b+1))/BinNum+XMin;
double darkness=0.5;
int Col=0;
double ratio=1.0*(PointsInBinMax-PointsInBin[b])/PointsInBinMax;
if (Col==0) aEPS.OutFile << darkness*ratio << " " << darkness*ratio << " " << 1.0 << " rgb\n"; // blue levels
if (Col==1) aEPS.OutFile << 1.0 << " " << 0.8*ratio << " " << 0.8*ratio << " rgb\n"; // red levels
if (Col==2) aEPS.OutFile << 0.8*ratio << " " << 0.8*ratio << " " << 0.8*ratio << " rgb\n"; // grey levels
if (Col==3) aEPS.OutFile << 0.8*ratio << " " << 1.0 << " " << 0.8*ratio << " rgb\n"; // green levels
if (Col==4) aEPS.OutFile << PointsInBin[b]/PointsInBinMax << " " << 0.0 << " " << PointsInBin[b]/PointsInBinMax << " rgb\n"; // magenta 1 0 1
if (Col==5) aEPS.OutFile << 0.0 << " " << 0.8*ratio << " " << 0.8*ratio << " rgb\n"; // cyan 0 1 1
// aEPS.OutFile << 0.8*(PointsInBinMax-PointsInBin[b])/PointsInBinMax << " setgray\n";
aEPS.FillRectR(a-0.49*da,val1,a+0.49*da,val2);
// aEPS.FillRectR(a-0.4*da,val1,a+0.4*da,val2);
PointsInBin[b]=0;
}
}
int b=int(((x[Opt.XDisp]-XMin)/(XMax-XMin))*BinNum);
if (b>=0 && b<BinNum) PointsInBin[b]++;
valold=val;
}
aEPS.OutFile << "grestore\n";
aEPS.Finish();
cout << "gv " << EPSFileName << endl;
cout << "PONumFound: " << PONumFound << " PWNumFound: " << PWNumFound << endl;
}
void BiffDiagPGM(string InFileName, int BinNum, double ParMin, double ParMax, double XMin, double XMax, string OutFileName){
cout << "Saving " << OutFileName << endl;
ifstream InFile(InFileName.c_str());
double valold=-1.0;
vector <int> PointsInBin(BinNum,0);
vector < vector <int> > res;
int xpos=0;
while (InFile){
double tmp,val;
vector <double> x(3);
InFile >> val >> tmp >> tmp >> x[0] >> x[1] >> x[2];
if (val<ParMin || val>ParMax) continue;
if (val!=valold && valold>0.0){
double da=val-valold;
double a=valold;
int PointsInBinMax=0;
int PointNum=0;
for (int b=0;b<BinNum;++b){
PointNum+=PointsInBin[b];
if (PointsInBin[b]>PointsInBinMax) PointsInBinMax=PointsInBin[b];
}
vector <int> col(BinNum,255); // initialize colors to be white
for (int b=0;b<BinNum;++b){
if (PointsInBin[b]==0) continue;
col[b]= int(255*0.8*(PointsInBinMax-PointsInBin[b])/PointsInBinMax);
// OutFile << xpos << " " << b << " " << 0.8*(PointsInBinMax-PointsInBin[b])/PointsInBinMax << endl;
PointsInBin[b]=0;
}
res.push_back(col);
cout << xpos << endl;
xpos++;
}
int b=int(((x[Opt.XDisp]-XMin)/(XMax-XMin))*BinNum);
if (b>=0 && b<BinNum) PointsInBin[b]++;
valold=val;
}
InFile.close();
ofstream OutFile(OutFileName.c_str());
OutFile << "P2" << endl;
OutFile << res.size() << ' ' << BinNum << endl;
OutFile << 255 << endl;
for (int y=BinNum-1;y>=0;--y){
OutFile << res[0][y];
for (int x=1;x<res.size();++x) OutFile << ' ' << int(res[x][y]);
OutFile << endl;
}
OutFile.close();
cout << "File saved: " << OutFileName << endl;
}
void BiffDiagPGM2(string InFileName, int BinNum, double ParMin, double ParMax, int ParNum, double XMin, double XMax, string OutFileName){
cout << "Saving " << OutFileName << endl;
ifstream InFile(InFileName.c_str());
double valold=-1.0;
vector <int> PointsInBin(BinNum,0);
vector < vector <int> > res(ParNum,vector<int>(BinNum,0)); // res[p][b] stores the number of entries with the parameter p in the bin b
int xpos=0;
double ParDelta=(ParMax-ParMin)/(ParNum-1);
int pnum=0;
while (InFile){
double tmp;
double val=ParMin-1.0;
vector <double> x(3);
InFile >> val >> tmp >> tmp >> x[0] >> x[1] >> x[2];
if (val<Opt.ReadParMin || val>Opt.ReadParMax) continue;
if (valold!=val){
cout << pnum++ << endl;
valold=val;
}
int b=int(((x[Opt.XDisp]-XMin)/(XMax-XMin))*BinNum);
int p=int((val+0.1*ParDelta-ParMin)/(ParMax-ParMin)*ParNum);
// cout << val << " " << x[Opt.XDisp]<< " " << p << " " << b << endl;
if (p>=0 && p<ParNum && b>=0 && b<BinNum) res[p][b]++;
/*
if (val!=valold && valold>0.0){
double da=val-valold;
double a=valold;
int PointsInBinMax=0;
int PointNum=0;
for (int b=0;b<BinNum;++b){
PointNum+=PointsInBin[b];
if (PointsInBin[b]>PointsInBinMax) PointsInBinMax=PointsInBin[b];
}
vector <int> col(BinNum,255); // initialize colors to be white
for (int b=0;b<BinNum;++b){
if (PointsInBin[b]==0) continue;
col[b]= int(255*0.8*(PointsInBinMax-PointsInBin[b])/PointsInBinMax);
// OutFile << xpos << " " << b << " " << 0.8*(PointsInBinMax-PointsInBin[b])/PointsInBinMax << endl;
PointsInBin[b]=0;
}
res.push_back(col);
cout << xpos << endl;
xpos++;
}*/
}
InFile.close();
ofstream OutFile(OutFileName.c_str());
OutFile << "P2" << endl;
OutFile << res.size() << ' ' << BinNum << endl;
OutFile << 255 << endl;
vector < vector <int> > colors(ParNum,vector<int>(BinNum,255)); // initialize colors to be white
for (int p=1;p<ParNum;++p){
int PointsInBinMax=0;
int PointNum=0;
for (int b=0;b<BinNum;++b){
PointNum+=res[p][b];
if (res[p][b]>PointsInBinMax) PointsInBinMax=res[p][b];
}
for (int b=0;b<BinNum;++b){
if (res[p][b]>0) colors[p][b]=int(255*0.8*(PointsInBinMax-res[p][b])/PointsInBinMax);
}
}
for (int y=BinNum-1;y>=0;--y){
OutFile << colors[0][y];
for (int x=1;x<res.size();++x) OutFile << ' ' << colors[x][y];
OutFile << endl;
}
OutFile.close();
cout << "File saved: " << OutFileName << endl;
}
bool PlotTrajectoryEps(double TimeStep, double SkipTime, double IntegrationTime, string FileName, EPSOptions EPSOpt){
cout << "PlotTrajectoryEps" << endl;
int XDisp=Opt.XDisp;
int YDisp=Opt.YDisp;
EPSOpt.LineWidth=0.0002;
EPSOpt.DiameterSymbol=10;
EPSOpt.Diameter=0.002;
EPSOpt.BoxLineWidth=0.01;
EPSOpt.TicsLineWidth=0.01;
EPSOpt.GridLineWidth=0.01;
EPSOpt.EPSFileCreator="smcc.x (Z. Galias)";
ostringstream oss;
oss << SystemName(Opt.CircuitVersion);
EPSOpt.EPSFileTitle=oss.str();
cout << "Saving " << FileName << endl;
SaveEPS aEPS(FileName.c_str());
aEPS.Init(EPSOpt);
aEPS.OutFile << "grestore\n";
aEPS.OutFile.precision(10);
aEPS.OutFile << "gsave\n";
aEPS.RescaleToBoundingBox(EPSOpt.XMin,EPSOpt.XMax,EPSOpt.YMin,EPSOpt.YMax);
// aEPS.SetClip(EPSOpt.bbXMin,EPSOpt.bbXMax,EPSOpt.bbYMin,EPSOpt.bbYMax);
aEPS.DefSizesAutomaticOut(EPSOpt.Diameter);
aEPS.DefLineWidthAutomaticOut(EPSOpt.LineWidth);
// Plot Trajectory
aEPS.OutFile << "0 0 1 rgb" << endl;
DVector x;
if (!InitialConditions(Opt.CircuitVersion,Opt.InitCondStr,x)) return false;
vector <DVector> xres;
vector <double> tres;
if (ComputeTrajectory(TimeStep,SkipTime,IntegrationTime,x,xres,tres)){
aEPS.OutFile.precision(6);
for (int p=0;p<xres.size();++p){
DVector x=xres[p];
aEPS.OutFile << aEPS.FunX(x(XDisp)) << ' ' << aEPS.FunY(x(YDisp)) << (p==0?" m ":" l ") << endl;
}
aEPS.OutFile << "s" << endl;
} else {
cerr << "ComputeTrajectory: Error" << endl;
}
aEPS.OutFile << "grestore\n";
aEPS.Finish();
cout << "gv " << FileName << endl;
return true;
}
/*
void PlotPoincTrajectoryEps(int SkipNum, int IterNum, string FileName, EPSOptions EPSOpt){
cout << "PlotPoincTrajectoryEps" << endl;
int XDisp=Opt.XDisp;
int YDisp=Opt.YDisp;
EPSOpt.LineWidth=0.002;
EPSOpt.DiameterSymbol=10;
EPSOpt.BoxLineWidth=0.01;
EPSOpt.TicsLineWidth=0.01;
EPSOpt.GridLineWidth=0.01;
EPSOpt.EPSFileCreator="ccpw (Z. Galias)";
ostringstream oss;
oss << Opt.SystemName << " return map trajectory";
EPSOpt.EPSFileTitle=oss.str();
cout << "Saving " << FileName << endl;
SaveEPS aEPS(FileName.c_str());
aEPS.Init(EPSOpt);
aEPS.OutFile << "grestore\n";
aEPS.OutFile.precision(10);
aEPS.OutFile << "gsave\n";
aEPS.RescaleToBoundingBox(EPSOpt.XMin,EPSOpt.XMax,EPSOpt.YMin,EPSOpt.YMax);
// aEPS.SetClip(EPSOpt.bbXMin,EPSOpt.bbXMax,EPSOpt.bbYMin,EPSOpt.bbYMax);
aEPS.DefSizesAutomaticOut(EPSOpt.Diameter);
aEPS.DefLineWidthAutomaticOut(EPSOpt.LineWidth);
// Plot Trajectory
DVector x(3);
x[0]=Opt.X0;
x[1]=Opt.Y0;
x[2]=Opt.Z0;
vector <DVector> xres;
vector <double> tres;
if (ComputeReturnMapTrajectory(SkipNum,IterNum,x,xres,tres,"",0.0)){
aEPS.OutFile.precision(5);
aEPS.OutFile << "0 0 1 rgb" << endl;
bool UseColors=Opt.ReturnMapColors;
if (UseColors){
if (XDisp!=YDisp){
for (int p=0;p<xres.size();++p) if (xres[p](YDisp)<0.0) aEPS.DrawPointR(xres[p](XDisp),xres[p](YDisp));
} else {
for (int p=0;p<xres.size()-1;++p) if (xres[p](YDisp)<0.0) aEPS.DrawPointR(xres[p](XDisp),xres[p+1](YDisp));
}
aEPS.OutFile << "1 0 0 rgb" << endl;
if (XDisp!=YDisp){
for (int p=0;p<xres.size();++p){
cout << "p: " << p << " x: " << xres[p] << endl;
if (xres[p](YDisp)>0.0) aEPS.DrawPointR(xres[p](XDisp),xres[p](YDisp));
}
} else {
for (int p=0;p<xres.size()-1;++p) if (xres[p](YDisp)>0.0) aEPS.DrawPointR(xres[p](XDisp),xres[p+1](YDisp));
}
} else {
if (XDisp!=YDisp){
for (int p=0;p<xres.size();++p) aEPS.DrawPointR(xres[p](XDisp),xres[p](YDisp));
} else {
for (int p=0;p<xres.size()-1;++p) aEPS.DrawPointR(xres[p](XDisp),xres[p+1](YDisp));
}
}
} else {
cerr << "ComputeReturnMapTrajectory: Error" << endl;
}
aEPS.OutFile << "grestore\n";
aEPS.Finish();
cout << "gv " << FileName << endl;
}
*/
template <class TMatrix, class TVector>
void CharPoly(const TMatrix & J, TVector & p){
// compute characteristic polynomial p of the matrix J
// works for real and interval matrices
int n=J.numberOfRows();
bool UseOldVersion=false;
if (UseOldVersion && (n==3)){
p=TVector(4);
p[3]=1.0; // coefficient at \lambda^3
p[2]=-J(1,1)-J(2,2)-J(3,3); // -tr(A)
p[1]=J(1,1)*J(2,2)+J(2,2)*J(3,3)+J(3,3)*J(1,1)-J(1,2)*J(2,1)-J(2,3)*J(3,2)-J(3,1)*J(1,3);
p[0]=-(J(1,1)*J(2,2)*J(3,3)+J(2,1)*J(3,2)*J(1,3)+J(3,1)*J(1,2)*J(2,3)-J(3,1)*J(2,2)*J(1,3)-J(3,2)*J(2,3)*J(1,1)-J(2,1)*J(1,2)*J(3,3)); // -det(J)
return;
}
// compute characteristic polynomial using Faddeev-LeVerrier algorithm
p=TVector(n+1,0.0);
p[n]=1.0;
TMatrix M(n,n,0.0);
for (int k=1;k<=n;++k){
M=J*M+p[n-k+1]*TMatrix::Identity(n);
p[n-k]=-Trace(J*M)/k;
}
}
template <class TVector>
bool IsHurwitz(TVector & p){
// verify whether polynomial p is Hurwitz (has all zeros with negative real part)
// works for real polynowmial coefficients (nonrigorous) and interval polynowmial coefficients (rigorous)
typedef typename TVector::ScalarType TScalar;
typedef capd::vectalg::Matrix<TScalar,0,0> TMatrix;
int n=p.dimension()-1;
bool UseOldVersion=false;
if (UseOldVersion && (n==3)){
bool stable=true;
bool unstable=false;
for (int k=0;k<3;++k){
if (!(p[k]>0)) stable=false;
if (p[k]<0) unstable=true;
}
if (!(p[2]*p[1]-p[0]*p[3]>0)) stable=false;
if (p[2]*p[1]-p[0]*p[3]<0) unstable=true;
if (stable==unstable) cout << "STABILITY UNKNOWN: Stable: " << (stable?"YES":"NO") << " Unstable: " << (unstable?"YES":"NO") << endl; // should rarely happen
return stable;
}
if (p[n]<0) p=-p; // make sure that the coefficient at p[n] is positive
if (!(p[n]<0) && !(p[n]>0)){
cout << "IsHurwitz, coefficient p[n] contains zero" << endl;
return false;
}
bool stable=true;
bool unstable=false;
for (int k=0;k<=n;++k){ // verify whether all coefficients are positive
if (!(p[k]>0)) stable=false;
if (p[k]<0) unstable=true;
}
if (stable==unstable){
cout << "STABILITY UNKNOWN: some coefficients contain zero" << endl;
return false;
}
if (!stable){
cout << "UNSTABLE, some coefficients are negative" << endl;
return stable;
}
TMatrix H(n,n,0.0); // Hurwitz matrix
MatrixSetValue(H,0.0);
for (int k=1;k<=n;++k){
int col=n-2*k+1;
for (int j=1;j<=n;++j){
if (col>=0 && col<=n) H(k,j)=p[col];
col++;
}
}
if (Opt.Verbose) cout << "Hurwitz: " << H << endl;
for (int k=1;k<=n;++k){
// skip n, n-2, n-4, ..., Hurwitz Lienard Criterion,
// it has been verified before that all coefficients are positive
if (k%2==n%2) continue;
TMatrix T(k,k,0.0);
for (int i=1;i<=k;++i){
for (int j=1;j<=k;++j) T(i,j)=H(i,j);
}
TScalar DetT;
try{
DetT=capd::matrixAlgorithms::det(T);
}
catch(exception& e){
cout << e.what() << endl;
return false;
}
if (Opt.Verbose) cout << "Minor k: " << k << " Det: " << DetT << endl;
if (DetT<=0.0) stable=false;
if (!(DetT<0.0) && !(DetT>0.0)) stable=false; //DetT contains zero
if (DetT<0.0) unstable=true;
}
if (stable==unstable) cout << "STABILITY UNKNOWN: Stable: " << (stable?"YES":"NO") << " Unstable: " << (unstable?"YES":"NO") << endl; // should rarely happen
return stable;
}
template <class TScalar, class TInterval>
bool NontrivialV1(const TInterval & g, TScalar & v1r, TInterval & v1){
// solve (Ga-G)v-I0*exp(-g/g0)sinh(v/V0)=0 in v using Newton method starting from v=v1
typedef capd::vectalg::Vector<TScalar,0> TVector;
typedef capd::vectalg::Matrix<TScalar,0,0> TMatrix;
typedef capd::vectalg::Vector<TInterval,0> TIVector;
typedef capd::vectalg::Matrix<TInterval,0,0> TIMatrix;
TInterval mu=(Opt.I0r/Opt.V0r)*exp(-g/Opt.g0r);
if (Opt.Gar-Opt.Gr<=mu) return false; // no nontrivial solutions
capd::map::Map<TMatrix> m(VectorFieldString(Version1DMap)); // Create real valued 1D map
SetParameters(Version1DMap,TScalar(0.0),m);
m.setParameter("g",g.leftBound());
// Find approximate position of zero, use v1r as initial point for the Newton method
// v1r=v1.leftBound();
for (int k=0;k<20;++k){
TMatrix J(1,1,0.0);
TVector F=m(TVector({v1r}),J);
if (Opt.Verbose) cout << "v1: " << v1r << " F: " << F << endl;
v1r=v1r-F(1)/J(1,1); // Newton method
}
// Prove the existence of zero
capd::map::Map<TIMatrix> mi(VectorFieldString(Version1DMap)); // Create interval valued 1D map
SetParameters(Version1DMap,TScalar(0.0),mi);
mi.setParameter("g",g);
v1=TInterval(v1r);
bool proved=false;
for (int k=0;k<20;++k){
TInterval v1mid=mid(v1);
TIMatrix J(1,1,0.0);
TIVector Fmid=mi(TIVector({v1mid}));
TIVector F=mi(TIVector({v1}),J);
TInterval Nv1=v1mid-Fmid(1)/J(1,1); // Newton method
if (Opt.Verbose) cout << "IntervalNewton iter: " << k << " v1: " << v1 << " diam(v1): " << DiamInterval(v1) << " Nv1: " << Nv1 << " diam(Nv1): " << DiamInterval(Nv1) << endl;
if (subsetInterior(Nv1,v1)){
v1=Nv1;
proved=true;
break;
}
TScalar eps=0.1;
v1=(1+eps)*Nv1-eps*Nv1; // inflate
}
if (Opt.Verbose) cout << "v1: " << v1 << " proved: " << (proved?"TRUE":"FALSE") << endl;
return proved;
}
template <class TInterval>
void smcc3FPs(){
typedef typename TInterval::BoundType TScalar;
Opt.CircuitVersion=0;
IDynSyst<TInterval> systi(Opt.Order,Opt.TimeStep,Opt.MaxTimeStep,Opt.AbsTol, Opt.RelTol, Opt.PoincWhich, Opt.PoincVar, Opt.PoincVal);
cout << "g: " << systi.m.getParameter("g") << endl;
TInterval g=systi.m.getParameter("g");
TInterval g0=systi.m.getParameter("g0");
TInterval I0=systi.m.getParameter("I0");
TInterval V0=systi.m.getParameter("V0");
TInterval G=systi.m.getParameter("G");
TInterval Ga=systi.m.getParameter("Ga");
TInterval C1=systi.m.getParameter("C1");
TInterval C2=systi.m.getParameter("C2");
TInterval L=systi.m.getParameter("L");
TInterval mu=(I0/V0)*exp(-g/g0);
if (Opt.Verbose){
cout << "mu: " << mu << endl;
if (Ga-G<=mu) cout << "Single fixed point" << endl;
if (Ga-G>mu) cout << "Three fixed points" << endl;
}
typedef capd::vectalg::Matrix<TInterval,0,0> TIMatrix;
typedef capd::vectalg::Vector<TInterval,0> TIVector;
TInterval v1(0.0);
TIVector FP({v1,TInterval(0.0),-G*v1});
TIMatrix J(3,3);
// Jacobian computed by hand
// J(1,1)=(Ga-G)/C1-(mu/C1)*(exp(v1/V0)+exp(-v1/V0))/2; J(1,2)=G/C1; J(1,3)=0;
// J(2,1)=G/C2; J(2,2)=-G/C2; J(2,3)=1/C2;
// J(3,1)=0; J(3,2)=-1/L; J(3,3)=0;
// cout << "J: " << J << endl;
TIVector F=systi.m(FP,J); // compute Jacobian using CAPD
cout << "FP: " << FP << endl;
cout << "F(FP): " << F << endl;
if (Opt.Verbose) cout << "J: " << J << endl;
// characteristic equation: c3\lambda^3+c2\lambda^2+c1\lambda+c0=0
TIVector C;
CharPoly(J,C); // compute characteristic polynomial p of the matrix J
cout << "CharPoly: " << C << endl;
bool stable=IsHurwitz(C);
cout << "Stable: "<< (stable?"TRUE":"FALSE") << endl;
if (Ga-G<=mu) return; // only trivial fixed point
// solve (Ga-G)v-I0*exp(-g/g0)sinh(v/V0)=0
for (int sign=0;sign<=1;++sign){ // find negative and positive nontrivial solutions v1
TScalar v1r=sqrt(6*(Opt.Gar-Opt.Gr)*power(Opt.V0r,3)*exp(g.leftBound()/Opt.g0r)/Opt.I0r-6*sqr(Opt.V0r));
if (sign==0) v1r=-v1r;
bool proved=NontrivialV1(g,v1r,v1);
if (!proved) continue;
TIMatrix J(3,3);
TIVector FP({v1,TInterval(0.0),-G*v1});
TIVector F=systi.m(FP,J); // compute Jacobian using CAPD
cout << "FP: " << FP << endl;
cout << "F(FP): " << F << endl;
if (Opt.Verbose) cout << "J: " << J << endl;
TIVector C;
CharPoly(J,C); // compute characteristic polynomial p of the matrix J
cout << "CharPoly: " << C << endl;
bool stable=IsHurwitz(C);
cout << "Stable: " << (stable?"TRUE":"FALSE") << endl;
}
}
template <class TIMatrix>
void CheckStability(const TIMatrix & J){
typedef typename TIMatrix::ScalarType TInterval;
typedef capd::vectalg::Vector<TInterval,0> TIVector;
int n=J.numberOfRows();
TIVector p;
CharPoly(J,p); // compute characteristic polynomial p of the matrix J
cout << "CharPoly: " << p << endl;
bool stable=IsHurwitz(p); // verify whether p is Hurwitz
cout << "Stable: " << (stable?"YES":"NO") << endl;
}
void ComputeEVs(DMatrix J){
int n=J.numberOfRows();
DVector rV(n),iV(n);
try{
capd::alglib::computeEigenvalues(J,rV,iV);
}
catch(exception& e){
cout << "ComputeEVs: " << e.what() << endl;
return;
}
cout << "EigenValues = " << capd::alglib::eigenvaluesToString(rV, iV, ", ") << endl;
CheckStability(J);
}
void ComputeEVs(const MpMatrix & J){
CheckStability(J); // for multiprecision computeEigenvalues does not work, stability is verified using Hurwitz criterion
}
template <class TInterval>
void smcc5FPs(){
// find a fixed point of smcc5 starting from Opt.v10
typedef typename TInterval::BoundType TScalar;
typedef capd::vectalg::Vector<TScalar,0> TVector;
typedef capd::vectalg::Matrix<TScalar,0,0> TMatrix;
typedef capd::vectalg::Vector<TInterval,0> TIVector;
typedef capd::vectalg::Matrix<TInterval,0,0> TIMatrix;
Opt.CircuitVersion=1;
DynSyst<TScalar> syst(Opt.Order,Opt.TimeStep,Opt.MaxTimeStep,Opt.AbsTol, Opt.RelTol, Opt.PoincWhich, Opt.PoincVar, Opt.PoincVal);
TScalar gamma0=syst.m.getParameter("gamma0");
TScalar g0=syst.m.getParameter("g0");
TScalar g1=syst.m.getParameter("g1");
TScalar Ga=syst.m.getParameter("Ga");
TScalar G=syst.m.getParameter("G");
TScalar I0=syst.m.getParameter("I0");
TScalar V0=syst.m.getParameter("V0");
TScalar T0=syst.m.getParameter("T0");
TScalar tth=syst.m.getParameter("tth");
TScalar Cth=syst.m.getParameter("Cth");
TScalar beta=syst.m.getParameter("beta");
TScalar a=syst.m.getParameter("a");
TScalar q=syst.m.getParameter("q");
TScalar Ear=syst.m.getParameter("Ear");
TScalar Eag=syst.m.getParameter("Eag");
TScalar Lg=syst.m.getParameter("Lg");
TScalar c1(beta*power(g0/g1,3));
TScalar c2(Lg*(Ear-Eag)/(2*a*q));
TScalar v10=-c2/gamma0;
if (Opt.v10!=0) v10=Opt.v10;
cout << "v10: " << v10 << endl;
TScalar v1=v10;
// Solve v1*(gamma0+c1*power(log((Ga-G)*v1/(I0*sinh(v1/V0))),3))+c2=0 using the Newton method
for (int i=0;i<20;++i){
TScalar F=v1*(gamma0+c1*power(log((Ga-G)*v1/(I0*sinh(v1/V0))),3))+c2;
TScalar FPrim=gamma0+c1*power(log((Ga-G)*v1/(I0*sinh(v1/V0))),3)+v1*c1*3*power(log((Ga-G)*v1/(I0*sinh(v1/V0))),2)*(1/v1-(1/sinh(v1/V0))*(cosh(v1/V0))*(1/V0));
TScalar h=-F/FPrim;
v1=v1+h;
F=v1*(gamma0+c1*power(log((Ga-G)*v1/(I0*sinh(v1/V0))),3))+c2;
if (Opt.Verbose) cout << "iter: " << i << " v1: " << v1 << " F: " << F << endl;
}
// Fixed Point
TVector xr({v1,TScalar(0.0),-G*v1,-g0*log(((Ga-G)*v1)/(I0*sinh(v1/V0))),T0+(tth/Cth)*sqr(v1)*(Ga-G)});
TMatrix J(5,5); // Jacobian
TVector F=syst.m(xr,J);
cout << "FP: " << xr << endl;
cout << "F(FP): " << F << endl;
cout << "J: " << J << endl;
ComputeEVs(J);
v10=v1;
bool prove=true;
if (prove) {
IDynSyst<TInterval> syst(Opt.Order,Opt.TimeStep,Opt.MaxTimeStep,Opt.AbsTol, Opt.RelTol, Opt.PoincWhich, Opt.PoincVar, Opt.PoincVal);
TInterval gamma0=syst.m.getParameter("gamma0");
TInterval g0=syst.m.getParameter("g0");
TInterval Ga=syst.m.getParameter("Ga");
TInterval G=syst.m.getParameter("G");
TInterval I0=syst.m.getParameter("I0");
TInterval V0=syst.m.getParameter("V0");
TInterval T0=syst.m.getParameter("T0");
TInterval c1(beta*power(g0/g1,3));
TInterval c2(Lg*(Ear-Eag)/(2*a*q));
// Prove the existence of a zero of v1*(gamma0+c1*power(log((Ga-G)*v1/(I0*sinh(v1/V0))),3))+c2=0 using the interval Newton method
bool proved=false;
TInterval v1=v10;
for (int k=0;k<20;++k){
TInterval v1mid=mid(v1);
TInterval F=v1mid*(gamma0+c1*power(log((Ga-G)*v1mid/(I0*sinh(v1mid/V0))),3))+c2;
TInterval FPrim=gamma0+c1*power(log((Ga-G)*v1/(I0*sinh(v1/V0))),3)+v1*c1*3*power(log((Ga-G)*v1/(I0*sinh(v1/V0))),2)*(1/v1-(1/sinh(v1/V0))*(cosh(v1/V0))*(1/V0));
TInterval Nv1=v1mid-F/FPrim;
if (Opt.Verbose) cout << "IntervalNewton iter: " << k << " v1: " << v1 << " diam(v1): " << diam(v1) << " Nv1: " << Nv1 << " diam(Nv1): " << diam(Nv1) << endl;
if (subsetInterior(Nv1,v1)){
v1=Nv1;
proved=true;
break;
}
TScalar eps=0.1;
v1=(1+eps)*Nv1-eps*Nv1; // inflate
}
cout << "proved: " << (proved?"true":"false") << endl;
if (!proved) return;
TIVector FP({v1,TInterval(0.0),-G*v1,-g0*log(((Ga-G)*v1)/(I0*sinh(v1/V0))),T0+(TInterval(tth)/TInterval(Cth))*sqr(v1)*(Ga-G)});
cout << "FP: " << FP << endl;;
for (int i=0;i<FP.dimension();++i) cout << diam(FP[i]).rightBound() << " ";
cout << endl;
TIMatrix J(FP.dimension(),FP.dimension());
TIVector F=syst.m(FP,J);
cout << "F(FP): " << F << endl;
CheckStability(J);
}
}
void setgmin(double & gmin){ gmin=Opt.gminr;}
void setgmax(double & gmax){ gmax=Opt.gmaxr;}
void setgmin(MpFloat & gmin){ gmin=MpFloat(Opt.gmins);}
void setgmax(MpFloat & gmax){ gmax=MpFloat(Opt.gmaxs);}
template <class TInterval, class TScalar>
void smcc5FRBorderv1(TScalar v1r, TInterval v1, int d, TScalar gr){
if ((d<0) || (d>1)){cout << "d out of range" << endl; return;}
typedef capd::vectalg::Vector<TScalar,0> TVector;
typedef capd::vectalg::Matrix<TScalar,0,0> TMatrix;
typedef capd::vectalg::Vector<TInterval,0> TIVector;
typedef capd::vectalg::Matrix<TInterval,0,0> TIMatrix;
Opt.CircuitVersion=1;
bool RealValuedComputations=true;
if (RealValuedComputations){
cout << "-------------REAL COMPUTATIONS-------------------------------" << endl;
DynSyst<TScalar> systr(Opt.Order,Opt.TimeStep,Opt.MaxTimeStep,Opt.AbsTol, Opt.RelTol, Opt.PoincWhich, Opt.PoincVar, Opt.PoincVal);
TVector FP({v1r,0*v1r,-systr.m.getParameter("G")*v1r,gr,systr.m.getParameter("T0")+(systr.m.getParameter("tth")/systr.m.getParameter("Cth"))*sqr(v1r)*(systr.m.getParameter("Ga")-systr.m.getParameter("G"))});
cout << "FPr: " << FP << endl;
TVector F=systr.m(FP);
cout << "F(FPr): " << F << endl;
if (F(4)>0.0) cout << "g: " << gr << " g grows" << endl;
if (F(4)<0.0) cout << "g: " << gr << " g decreases" << endl;
TScalar gamma=systr.m.getParameter("gamma0")-systr.m.getParameter("beta")*power(gr/systr.m.getParameter("g1"),3);
TScalar tmp=2*systr.m.getParameter("a")*systr.m.getParameter("q")*gamma*v1r-systr.m.getParameter("Lg")*(systr.m.getParameter("Eag")-systr.m.getParameter("Ear"));
bool FixedPointExists=false;
if ((d==0) && (tmp>=0.0)) FixedPointExists=true;
if ((d==1) && (tmp<=0.0)) FixedPointExists=true;
if (FixedPointExists){
cout << "Fixed point exists: " << FP << endl;
} else {
cout << "No fixed point at: " << FP << endl;
}
// define 4D dynamical system (g is fixed)
capd::map::Map<TMatrix> m(VectorFieldString(Version4D));
SetParameters(Version4D,TScalar(0.0),m);
m.setParameter("g",gr);
TVector FP4({FP(1),FP(2),FP(3),FP(5)});
TMatrix J(4,4);
F=m(FP4,J);
cout << "FP4r: " << FP4 << endl;
cout << "F(FP4r): " << F << endl;
if (Opt.Verbose) cout << "Jr: " << J << endl;
ComputeEVs(J);
}
cout << "-------------INTERVAL COMPUTATIONS-------------------------------" << endl;
IDynSyst<TInterval> systi(Opt.Order,Opt.TimeStep,Opt.MaxTimeStep,Opt.AbsTol, Opt.RelTol, Opt.PoincWhich, Opt.PoincVar, Opt.PoincVal);
TInterval g(gr);
TIVector FP({v1,0*v1,-systi.m.getParameter("G")*v1,g,systi.m.getParameter("T0")+(systi.m.getParameter("tth")/systi.m.getParameter("Cth"))*sqr(v1)*(systi.m.getParameter("Ga")-systi.m.getParameter("G"))});
cout << "FP: " << FP << endl;
TIVector F=systi.m(FP);
cout << "F(FP): " << F << endl;
if (F(4)>0.0) cout << "g: " << g << " g grows" << endl;
if (F(4)<0.0) cout << "g: " << g << " g decreases" << endl;
TInterval gamma=systi.m.getParameter("gamma0")-systi.m.getParameter("beta")*power(g/systi.m.getParameter("g1"),3);
TInterval tmp=2*systi.m.getParameter("a")*systi.m.getParameter("q")*gamma*v1-systi.m.getParameter("Lg")*(systi.m.getParameter("Eag")-systi.m.getParameter("Ear"));
bool FixedPointExists=false;
if ((d==0) && (tmp>=0.0)) FixedPointExists=true;
if ((d==1) && (tmp<=0.0)) FixedPointExists=true;
if (FixedPointExists){
cout << "Fixed point exists: " << FP << endl;
} else {
cout << "No fixed point at: " << FP << endl;
}
// define 4D dynamical system (g is fixed)
capd::map::Map<TIMatrix> m(VectorFieldString(Version4D));
SetParameters(Version4D,TScalar(0.0),m);
m.setParameter("g",g);
TIVector FP4({FP(1),FP(2),FP(3),FP(5)});
TIMatrix J(4,4);
F=m(FP4,J);
cout << "FP4: " << FP4 << endl;
cout << "F(FP4): " << F << endl;
if (Opt.Verbose) cout << "J: " << J << endl;
CheckStability(J);
}
template <class TInterval>
void smcc5FPBorders(){
// find fixed points at borders g=gmin and g=gmax
typedef typename TInterval::BoundType TScalar;
for (int d=0;d<=1;++d){ // test g=gmin (d==0) and g=gmax (d==1)
TScalar gr;
if (d==0) setgmin(gr);
if (d==1) setgmax(gr);
cout << "--------------------------------------------" << endl;
cout << "g: " << gr << " trivial fixed point." << endl;
TInterval v1(0.0);
smcc5FRBorderv1(TScalar(0.0),v1,d,gr); // test whether v1 generates a fixed point
cout << "--------------------------------------------" << endl;
cout << "g: " << gr << " nontrivial fixed point." << endl;
TScalar mu=(Opt.I0r/Opt.V0r)*exp(-gr/Opt.g0r);
if (Opt.Gar-Opt.Gr<=mu){
cout << "No nontrivial solutions for g: " << gr << endl;
continue;
}
for (int sign=0;sign<=1;++sign){ // find negative and positive nontrivial solutions v1
TScalar v1r=sqrt(6*(Opt.Gar-Opt.Gr)*power(Opt.V0r,3)*exp(gr/Opt.g0r)/Opt.I0r-6*sqr(Opt.V0r));
if (sign==0) v1r=-v1r;
if (Opt.Verbose) cout << "v1r: " << v1r << endl;
TInterval g(gr);
if (!NontrivialV1(g,v1r,v1)){
cout << "Solution not found for g: " << gr << " v1r: " << v1r << endl;
continue;
}
if (Opt.Verbose) cout << "Solution found for g: " << gr << " v1r: " << v1r << " v1: " << v1 << endl;
smcc5FRBorderv1(v1r,v1,d,gr); // test whether v1 generates a fixed point
}
}
}
void QuadrSet2QuadrSet(string InFile, string OutFile){
cout << "QuadrSet2QuadrSet" << endl;
QuadrSet QS;
QS.Read(InFile);
QS.Save(OutFile,Opt.CoutPrec);
}
void QuadrSet2PolySet(string InFile, string OutFile){
QuadrSet QS;
QS.Read(InFile);
PolygonSet Polygons=QS.GetPolygonSet(false);
PolygonSetSave(Polygons,OutFile,Opt.CoutPrec);
}
void QuadrSetSaveEPS(string InFile, string EpsFile, bool singlecolor){
QuadrSet QS;
QS.Read(InFile);
QS.PrintEPS(EpsFile,Opt.EPSOpt,Opt.XDisp,Opt.YDisp,false,false,singlecolor,false);
}
bool BrokenLineImageInOutPolygonSet(const BrokenLine & S, const PolygonSet & PS, IPoincareMap & pm, bool CheckInside, int XDisp, int YDisp, int & TestNum, vector <IVector> & Cover, vector <IVector> & CoverIm, Interval & tbound){
// CheckInside=true -> check if S is inside PS
// CheckInside=false -> check is S has empty intersection with PS
static double MinOKSize=1e10;
static double MaxOKSize=0.0;
TimeCount TC;
bool TestPoints=true;
for (int e=0;e<S.size()-1;++e){ // process each edge
list <IVector> xsqu,xequ;
int TestPointsNum=2;
if (TestPoints){
TestPointsNum=100;
for (int k=0;k<TestPointsNum;++k){
xsqu.push_back(S[e]*(1.0-1.0*k/TestPointsNum)+S[e+1]*(1.0*k/TestPointsNum));
xequ.push_back(S[e]*(1.0-1.0*(k+1)/TestPointsNum)+S[e+1]*(1.0*(k+1)/TestPointsNum));
}
} else {
xsqu.push_back(S[e]);
xequ.push_back(S[e+1]);
}
double ToDo=EuclDist(S[e],S[e+1]).rightBound();
double Done=0.0;
for (;;){
if (xsqu.empty()) break;
IVector xs=xsqu.front(); xsqu.pop_front();
IVector xe=xequ.front(); xequ.pop_front();
// if (TestPoints) xe=xs;
double TestSize=EuclDist(xs,xe).rightBound();
IVector x=Hull(xs,xe);
IVector y;
Interval t;
bool ok=true;
if (Opt.MaxTestSize>0.0 && TestSize>Opt.MaxTestSize){ ok=false;}
if (ok){
TestNum++;
if (!PoincMapDS(pm,x,y,t)){
cout << "PoincOKerror " << x << " " << y << " " << t << endl;
ok=false;
} else {
cout << "PoincOK" << x << " " << y << " " << t << endl;
}
}
if (ok){
if (CheckInside){
bool inside=false;
for (int pp=0;pp<PS.size();++pp){
if (PointInPolyInterior(y,PS[pp],XDisp,YDisp)==1) inside=true;
}
if (!inside) ok=false;
} else {
bool outside=true;
for (int pp=0;pp<PS.size();++pp){
if (PointInPolyInterior(y,PS[pp],XDisp,YDisp)!=0) outside=false;
}
if (!outside) ok=false;
}
}
if (Opt.Verbose) cout << "BLIInOutPS " << (CheckInside?1:0) << " e: " << e << " Done: " << Done/ToDo << " TestNum: " << TestNum << " x: " << x << " Diam: " << VectorMaxDiam(x).rightBound() << " y: " << y << " Diam: " << VectorMaxDiam(y).rightBound() << " " << (ok?"OK":"NO") << " time: " << TC << endl;
if (ok){
if (leftBound(tbound)<0.0) tbound=t;
else tbound=intervalHull(tbound,t);
Cover.push_back(x);
CoverIm.push_back(y);
Done+=TestSize;
if (TestSize<MinOKSize) MinOKSize=TestSize;
if (TestSize>MaxOKSize) MaxOKSize=TestSize;
continue;
}
IVector xm=0.5*(xs+xe); // split interval (xs,xe) into two parts and add to queues
xsqu.push_back(xs); xequ.push_back(xm);
xsqu.push_back(xm); xequ.push_back(xe);
}
cout << "Edge: " << e << " TestNum: " << TestNum << " CoverSize: " << Cover.size() << " MinOKSize: " << MinOKSize << " MaxOKSize: " << MaxOKSize << " time: " << TC << endl;
}
return true;
}
void PrintQS(SaveEPS & aEPS, const QuadrSet & QS, int XDisp, int YDisp, EPSOptions & EPSOpt, bool fill=false, bool DoSymmetric=false){
cout << "Printing QuadrSet: " << QS.QuadrNum() << endl;
for (int i=0;i<QS.QuadrNum();i++){
Quadr Q=QS.Q[i];
IVector Range=Q.GetRange();
if (Inf(Range(XDisp))>EPSOpt.XMax || Sup(Range(XDisp))<EPSOpt.XMin || Inf(Range(YDisp))>EPSOpt.YMax || Sup(Range(YDisp))<EPSOpt.YMin) continue;
aEPS.OutFile << ColorString(Opt.SingleColor?0:i) << " rgb\n";
for (int j=0;j<Q.SegmentNum();j++){
if (j==0) aEPS.OutFile << "n " << aEPS.FunX(Mid(Q.S[j].P[0](XDisp))) << " " << aEPS.FunY(Mid(Q.S[j].P[0](YDisp))) << " m ";
for (int k=1;k<Q.S[j].PointNum();k++){
aEPS.OutFile << aEPS.FunX(Mid(Q.S[j].P[k](XDisp))) << " " << aEPS.FunY(Mid(Q.S[j].P[k](YDisp))) << " l ";
}
}
if (fill) aEPS.OutFile << "c f\n";
else aEPS.OutFile << "c s\n";
if (DoSymmetric){
for (int j=0;j<Q.SegmentNum();j++){
if (j==0) aEPS.OutFile << "n " << aEPS.FunX(-Mid(Q.S[j].P[0](XDisp))) << " " << aEPS.FunY(-Mid(Q.S[j].P[0](YDisp))) << " m ";
for (int k=1;k<Q.S[j].PointNum();k++){
aEPS.OutFile << aEPS.FunX(-Mid(Q.S[j].P[k](XDisp))) << " " << aEPS.FunY(-Mid(Q.S[j].P[k](YDisp))) << " l ";
}
}
if (fill) aEPS.OutFile << "c f\n";
else aEPS.OutFile << "c s\n";
}
for (int j=0;j<Q.SegmentNum();j++){
aEPS.OutFile << ColorString(Opt.SingleColor?0:i) << " rgb\n";
aEPS.OutFile << "n " << aEPS.FunX(Mid(Q.S[j].P[0](XDisp))) << " " << aEPS.FunY(Mid(Q.S[j].P[0](YDisp))) << " m ";
for (int k=1;k<Q.S[j].PointNum();k++){
aEPS.OutFile << aEPS.FunX(Mid(Q.S[j].P[k](XDisp))) << " " << aEPS.FunY(Mid(Q.S[j].P[k](YDisp))) << " l ";
}
aEPS.OutFile << "s" << endl;
if (DoSymmetric){
aEPS.OutFile << "n " << aEPS.FunX(-Mid(Q.S[j].P[0](XDisp))) << " " << aEPS.FunY(-Mid(Q.S[j].P[0](YDisp))) << " m ";
for (int k=1;k<Q.S[j].PointNum();k++){
aEPS.OutFile << aEPS.FunX(-Mid(Q.S[j].P[k](XDisp))) << " " << aEPS.FunY(-Mid(Q.S[j].P[k](YDisp))) << " l ";
}
aEPS.OutFile << "s" << endl;
}
}
}
}
void VectorSetAndQSSaveEPS(const vector <DVector> & VSet, const QuadrSet & QS, string FileName, EPSOptions EPSOpt, bool fill=false, bool DoSymmetric=false, string info=""){
cout << "IVectorSetSaveEps" << endl;
int XDisp=Opt.XDisp;
int YDisp=Opt.YDisp;
EPSOpt.LineWidth=0.002;
EPSOpt.DiameterSymbol=10;
EPSOpt.Diameter=0.002;
EPSOpt.BoxLineWidth=0.01;
EPSOpt.TicsLineWidth=0.01;
EPSOpt.GridLineWidth=0.01;
EPSOpt.EPSFileCreator="smcc (Z. Galias)";
ostringstream oss;
oss << SystemName(Opt.CircuitVersion) << info;
EPSOpt.EPSFileTitle=oss.str();
cout << "Saving " << FileName << endl;
SaveEPS aEPS(FileName.c_str());
aEPS.Init(EPSOpt);
aEPS.OutFile << "grestore\n";
aEPS.OutFile.precision(10);
aEPS.OutFile << "gsave\n";
aEPS.RescaleToBoundingBox(EPSOpt.XMin,EPSOpt.XMax,EPSOpt.YMin,EPSOpt.YMax);
// aEPS.SetClip(EPSOpt.bbXMin,EPSOpt.bbXMax,EPSOpt.bbYMin,EPSOpt.bbYMax);
aEPS.DefSizesAutomaticOut(EPSOpt.Diameter);
aEPS.DefLineWidthAutomaticOut(EPSOpt.LineWidth);
aEPS.OutFile.precision(5);
aEPS.OutFile << "0 setlinewidth" << endl;
aEPS.OutFile << "0 0 1 rgb" << endl;
// Plot Trajectory
bool UseColors=!Opt.SingleColor;
if (VSet.size()>0){
cout << "Printing VSet: " << VSet.size() << endl;
if (UseColors){
if (XDisp!=YDisp){
for (int p=0;p<VSet.size();++p) if (VSet[p](YDisp)<0.0) aEPS.DrawPointR(VSet[p](XDisp),VSet[p](YDisp));
} else {
for (int p=0;p<VSet.size()-1;++p) if (VSet[p](YDisp)<0.0) aEPS.DrawPointR(VSet[p](XDisp),VSet[p+1](YDisp));
}
aEPS.OutFile << "1 0 0 rgb" << endl;
if (XDisp!=YDisp){
for (int p=0;p<VSet.size();++p){
cout << "p: " << p << " x: " << VSet[p] << endl;
if (VSet[p](YDisp)>0.0) aEPS.DrawPointR(VSet[p](XDisp),VSet[p](YDisp));
}
} else {
for (int p=0;p<VSet.size()-1;++p) if (VSet[p](YDisp)>0.0) aEPS.DrawPointR(VSet[p](XDisp),VSet[p+1](YDisp));
}
} else {
if (XDisp!=YDisp){
for (int p=0;p<VSet.size();++p) aEPS.DrawPointR(VSet[p](XDisp),VSet[p](YDisp));
} else {
for (int p=0;p<VSet.size()-1;++p) aEPS.DrawPointR(VSet[p](XDisp),VSet[p+1](YDisp));
}
}
}
if (QS.size()>0) PrintQS(aEPS,QS,XDisp,YDisp,EPSOpt,fill,DoSymmetric);
aEPS.OutFile << "grestore\n";
aEPS.Finish();
cout << "gv " << FileName << endl;
}
void IVectorSetsAndQSSaveEPS(const vector <IVector> & IVSet1, const vector <IVector> & IVSet2, const QuadrSet & QS, string FileName, EPSOptions EPSOpt, bool fill=false, bool DoSymmetric=false, string info=""){
cout << "IVectorSetSaveEps" << endl;
int XDisp=Opt.XDisp;
int YDisp=Opt.YDisp;
EPSOpt.LineWidth=0.002;
EPSOpt.DiameterSymbol=10;
EPSOpt.Diameter=0.002;
EPSOpt.BoxLineWidth=0.01;
EPSOpt.TicsLineWidth=0.01;
EPSOpt.GridLineWidth=0.01;
EPSOpt.EPSFileCreator="smcc (Z. Galias)";
ostringstream oss;
oss << SystemName(Opt.CircuitVersion) << info;
EPSOpt.EPSFileTitle=oss.str();
cout << "Saving " << FileName << endl;
SaveEPS aEPS(FileName.c_str());
aEPS.Init(EPSOpt);
aEPS.OutFile << "grestore\n";
aEPS.OutFile.precision(10);
aEPS.OutFile << "gsave\n";
aEPS.RescaleToBoundingBox(EPSOpt.XMin,EPSOpt.XMax,EPSOpt.YMin,EPSOpt.YMax);
// aEPS.SetClip(EPSOpt.bbXMin,EPSOpt.bbXMax,EPSOpt.bbYMin,EPSOpt.bbYMax);
aEPS.DefSizesAutomaticOut(EPSOpt.Diameter);
aEPS.DefLineWidthAutomaticOut(EPSOpt.LineWidth);
aEPS.OutFile.precision(5);
aEPS.OutFile << "0 setlinewidth" << endl;
if (IVSet1.size()>0){
cout << "Printing IVector: " << IVSet1.size() << endl;
aEPS.OutFile << "0 1 1 rgb" << endl;
for (int p=0;p<IVSet1.size();++p) aEPS.DrawRectR(Inf(IVSet1[p](XDisp)),Inf(IVSet1[p](YDisp)),Sup(IVSet1[p](XDisp)),Sup(IVSet1[p](YDisp)));
}
if (IVSet2.size()>0){
cout << "Printing IVector: " << IVSet2.size() << endl;
aEPS.OutFile << "1 0 1 rgb" << endl;
for (int p=0;p<IVSet2.size();++p) aEPS.DrawRectR(Inf(IVSet2[p](XDisp)),Inf(IVSet2[p](YDisp)),Sup(IVSet2[p](XDisp)),Sup(IVSet2[p](YDisp)));
}
if (QS.size()>0) PrintQS(aEPS,QS,XDisp,YDisp,EPSOpt,fill,DoSymmetric);
aEPS.OutFile << "grestore\n";
aEPS.Finish();
cout << "gv " << FileName << endl;
}
vector <IVector> Split(const IVector & X, int x, int y){
vector <IVector> Res;
for (int ix=0;ix<2;ix++){
for (int iy=0;iy<2;iy++){
IVector Y=X;
if (ix==0) Y[x]=Interval(Inf(X[x]),Mid(X[x]));
else Y[x]=Interval(Mid(X[x]),Sup(X[x]));
if (iy==0) Y[y]=Interval(Inf(X[y]),Mid(X[y]));
else Y[y]=Interval(Mid(X[y]),Sup(X[y]));
Res.push_back(Y);
}
}
return Res;
}
bool BrokenLineImageInsidePolygonSet(const BrokenLine & BL, const PolygonSet & PS, IPoincareMap & pm, int XDisp, int YDisp, int & TestNum, vector <IVector> & Cover, vector <IVector> & CoverIm, Interval & tbound){
return BrokenLineImageInOutPolygonSet(BL,PS,pm,true,XDisp,YDisp,TestNum,Cover,CoverIm,tbound);
}
bool BrokenLineImageOutsidePolygonSet(const BrokenLine & BL, const PolygonSet & PS, IPoincareMap & pm, int XDisp, int YDisp, int & TestNum, vector <IVector> & Cover, vector <IVector> & CoverIm, Interval & tbound){
return BrokenLineImageInOutPolygonSet(BL,PS,pm,false,XDisp,YDisp,TestNum,Cover,CoverIm,tbound);
}
bool BrokenLineImageNotInBrokenLineSet(const BrokenLine & S, const BrokenLineSet & BLS, IPoincareMap & pm, int XDisp, int YDisp, int & TestNum, vector <IVector> & Cover, vector <IVector> & CoverIm, Interval & tbound){
static double MinOKSize=1e10;
static double MaxOKSize=0.0;
TimeCount TC;
for (int e=0;e<S.size()-1;++e){ // process each edge
list <IVector> xsqu,xequ;
xsqu.push_back(S[e]);
xequ.push_back(S[e+1]);
double ToDo=EuclDist(S[e],S[e+1]).rightBound();
double Done=0.0;
for (;;){
if (xsqu.empty()) break;
IVector xs=xsqu.front(); xsqu.pop_front();
IVector xe=xequ.front(); xequ.pop_front();
double TestSize=EuclDist(xs,xe).rightBound();
IVector x=Hull(xs,xe);
IVector y;
Interval t;
bool ok=true;
if (Opt.MaxTestSize>0.0 && TestSize>Opt.MaxTestSize){ ok=false;}
if (ok){
TestNum++;
if (!PoincMapDS(pm,x,y,t)) ok=false;
}
if (ok){
for (int bl=0;bl<BLS.size();++bl){
if (!PointNotInBrokenLine(y,BLS[bl],XDisp,YDisp)) ok=false;
}
}
if (Opt.Verbose) cout << "BLINotInBLS e: " << e << " Done: " << Done/ToDo << " TestNum: " << TestNum << " x: " << x << " Diam: " << VectorMaxDiam(x).rightBound() << " y: " << x << " Diam: " << VectorMaxDiam(y).rightBound() << " " << (ok?"OK":"NO") << " time: " << TC << endl;
if (ok){
if (leftBound(tbound)<0.0) tbound=t;
else tbound=intervalHull(tbound,t);
Cover.push_back(x);
CoverIm.push_back(y);
Done+=TestSize;
if (TestSize<MinOKSize) MinOKSize=TestSize;
if (TestSize>MaxOKSize) MaxOKSize=TestSize;
continue;
}
IVector xm=0.5*(xs+xe); // split interval (xs,xe) into two parts and add to queues
xsqu.push_back(xs); xequ.push_back(xm);
xsqu.push_back(xm); xequ.push_back(xe);
}
cout << "Edge: " << e << " TestNum: " << TestNum << " CoverSize: " << Cover.size() << " MinOKSize: " << MinOKSize << " MaxOKSize: " << MaxOKSize << " time: " << TC << endl;
}
return true;
}
string Int2String(int x){
ostringstream oss;
oss << x;
return oss.str();
}
bool TrappingRegionBorder(string InFile){
TimeCount TC;
int XDisp=0;
if (XDisp==Opt.PoincVar) XDisp++;
int YDisp=XDisp+1;
if (YDisp==Opt.PoincVar) YDisp++;
PolygonSet Polygons;
PolygonSetLoad(InFile,Polygons);
IDynSyst<Interval> syst(Opt.Order,Opt.TimeStep,Opt.MaxTimeStep,Opt.AbsTol, Opt.RelTol, Opt.PoincWhich, Opt.PoincVar, Opt.PoincVal);
int TestNum=0;
vector <IVector> Cover;
vector <IVector> CoverIm;
bool res=true;
Interval tbound=-1.0;
for (int p=0;p<Polygons.size();++p){ // process each polygon
if (BrokenLineImageInsidePolygonSet(Polygons[p],Polygons,syst.pm,XDisp,YDisp,TestNum,Cover,CoverIm,tbound)){
} else {
res=false;
cout << "Error" << endl;
}
}
cout << "CoverSize: " << Cover.size() << " ReturnTimeBound: " << tbound << endl;
if (Opt.BaseFile!=""){
string CoverFileNameTxt=Opt.BaseFile+"BC"+".txt"; // Border Cover
string CoverImFileNameTxt=Opt.BaseFile+"BCIm"+".txt"; // Border Cover Image
IVectorSetSave(Cover,CoverFileNameTxt,Opt.CoutPrec);
IVectorSetSave(CoverIm,CoverImFileNameTxt,Opt.CoutPrec);
}
return res;
}
void TrappingRegionInterior(string InFile, bool CheckInside, bool BorderOnly=false){
TimeCount TC;
int XDisp=0;
if (XDisp==Opt.PoincVar) XDisp++;
int YDisp=XDisp+1;
if (YDisp==Opt.PoincVar) YDisp++;
PolygonSet Polygons;
PolygonSetLoad(InFile,Polygons);
IDynSyst<Interval> syst(Opt.Order,Opt.TimeStep,Opt.MaxTimeStep,Opt.AbsTol, Opt.RelTol, Opt.PoincWhich, Opt.PoincVar, Opt.PoincVal);
list <IVector> qu;
int TestNum=0;
int SkipNum=0;
vector <IVector> Cover;
vector <IVector> CoverIm;
double MinOKSize=1e10;
double MaxOKSize=0.0;
Interval tbound=-1.0;
for (int p=0;p<Polygons.size();++p){
IVector Range=GetRange(Polygons[p]);
cout << "Poly: " << p << " Range: " << Range << endl;
qu.push_back(Range);
vector <IVector> CoverSingle;
vector <IVector> CoverImSingle;
for (;;){
if (qu.empty()) break;
IVector x=qu.front();
qu.pop_front();
Interval t;
bool ok=true;
IVector y;
double TestSize=VectorMaxDiam(x).rightBound();
if (Opt.MaxTestSize>0.0 && TestSize>Opt.MaxTestSize){ SkipNum++; ok=false;}
if (ok){
TestNum++;
if (!PoincMapDS(syst.pm,x,y,t)) ok=false;
}
if (ok && CheckInside){
bool inside=false;
for (int pp=0;pp<Polygons.size();++pp){
if (PointInPolyInterior(y,Polygons[pp],XDisp,YDisp)==1) inside=true;
}
if (!inside) ok=false;
}
cout << "Tested: " << TestNum << " Skipped: " << SkipNum << " ToDo: " << qu.size() << " x: " << x << " Diamx: " << TestSize << " Diamy: " << VectorMaxDiam(y).rightBound() << " " << (ok?"OK":"NO") << " time: " << TC << endl;
if (ok){
if (leftBound(tbound)<0.0) tbound=t;
else tbound=intervalHull(tbound,t);
if (TestSize<MinOKSize) MinOKSize=TestSize;
if (TestSize>MaxOKSize) MaxOKSize=TestSize;
CoverSingle.push_back(x);
CoverImSingle.push_back(y);
continue;
}
vector <IVector> xsplit=Split(x,XDisp,YDisp); // split x into 4 parts and add to queues
for (int i=0;i<xsplit.size();++i){
bool OutsideAll=true;
bool InsideOne=false;
for (int p=0;p<Polygons.size();++p){
int res=PointInPolyInterior(xsplit[i],Polygons[p],XDisp,YDisp);
if (res!=0) OutsideAll=false; // xsplit[i] is not outside Polygons[p]
if (res==1) InsideOne=true; // xsplit[i] is inside Polygons[p]
}
if (OutsideAll) continue; // xsplit[i] is outside Polygons[p], do not add it to the list
if (InsideOne && BorderOnly) continue; // xsplit[i] is inside Polygons[p], do not add it to the list is BorderOnly is set
qu.push_back(xsplit[i]);
}
}
for (int p=0;p<CoverSingle.size();++p) Cover.push_back(CoverSingle[p]);
for (int p=0;p<CoverImSingle.size();++p) CoverIm.push_back(CoverImSingle[p]);
string CoverFileNameTxt=Opt.BaseFile+"IC"+Int2String(p)+".txt"; // Interior Cover
string CoverImFileNameTxt=Opt.BaseFile+"ICIm"+Int2String(p)+".txt"; // Interior Cover Image
IVectorSetSave(CoverSingle,CoverFileNameTxt,Opt.CoutPrec);
IVectorSetSave(CoverImSingle,CoverImFileNameTxt,Opt.CoutPrec);
}
cout << "Tested: " << TestNum << " Skipped: " << SkipNum << " CoverSize: " << Cover.size() << " CheckInside: " << (CheckInside?"OK":"NO") << " MinOKDiam: " << MinOKSize << " MaxOKDiam: " << MaxOKSize << " ReturnTimeBound " << tbound << " time: " << TC << endl;
if (Opt.BaseFile!=""){
string CoverFileNameTxt=Opt.BaseFile+"IC"+".txt"; // Interior Cover
string CoverImFileNameTxt=Opt.BaseFile+"ICIm"+".txt"; // Interior Cover Image
IVectorSetSave(Cover,CoverFileNameTxt,Opt.CoutPrec);
IVectorSetSave(CoverIm,CoverImFileNameTxt,Opt.CoutPrec);
}
}
bool CoveringRelation(string InFile){
TimeCount TC;
int XDisp=0;
if (XDisp==Opt.PoincVar) XDisp++;
int YDisp=XDisp+1;
if (YDisp==Opt.PoincVar) YDisp++;
QuadrSet QS;
QS.Read(InFile);
PolygonSet Polygons=QS.GetPolygonSet(false);
IDynSyst<Interval> syst(Opt.Order,Opt.TimeStep,Opt.MaxTimeStep,Opt.AbsTol, Opt.RelTol, Opt.PoincWhich, Opt.PoincVar, Opt.PoincVal);
bool res;
int TestNum=0;
int CoverSize=0;
vector <IVector> VCover;
vector <IVector> HCover;
vector <IVector> VCoverIm;
vector <IVector> HCoverIm;
vector < vector <bool> > CoverPossible(QS.size(), vector<bool>(QS.size(),false));
for (int q=0;q<QS.size();++q){
if (QS.Q[q].size()!=4){ cout << "This is not a quadrangle" << endl; continue;}
// identify which covering relations are to be checked
// at least one point along horizontal edge of q should be inside i
for (int s=0;s<QS.Q[q].size();++s){
if (s%2==0) continue; // process only horizontal edges
for (int p=0;p<QS.Q[q].S[s].size();++p){
if (p==0 || p==QS.Q[q].S[s].size()-1) continue;
IVector x=QS.Q[q].S[s].P[p];
IVector y;
Interval t;
if (!PoincMapDS(syst.pm,x,y,t)) continue;
for (int i=0;i<QS.size();++i){
if (PointInPolyInterior(y,Polygons[i],XDisp,YDisp)==1) CoverPossible[q][i]=true;
}
}
}
}
for (int q=0;q<QS.size();++q){
if (QS.Q[q].size()!=4){ cout << "This is not a quadrangle" << endl; continue;}
// identify which covering relations are to be checked
// corners of q cannot be inside i
for (int s=0;s<QS.Q[q].size();++s){
if (s%2==0) continue; // process only horizontal edges
for (int p=0;p<QS.Q[q].S[s].size();++p){
if (p>0 && p<QS.Q[q].S[s].size()-1) continue;
IVector x=QS.Q[q].S[s].P[p];
IVector y;
Interval t;
if (!PoincMapDS(syst.pm,x,y,t)) continue;
for (int i=0;i<QS.size();++i){
if (PointInPolyInterior(y,Polygons[i],XDisp,YDisp)==1) CoverPossible[q][i]=false;
}
}
}
}
cout << "CoveringPossible: " << endl;
for (int q=0;q<CoverPossible.size();++q){
// cout << q << " ";
for (int p=0;p<CoverPossible[q].size();++p) cout << (CoverPossible[q][p]?"1":"0");
cout << endl;
}
Interval tbound=-1.0;
vector < vector <bool> > CoverDone(QS.size(), vector<bool>(QS.size(),false));
for (int q=0;q<QS.size();++q){
vector <IVector> VCoverSingle;
vector <IVector> HCoverSingle;
vector <IVector> VCoverImSingle;
vector <IVector> HCoverImSingle;
// if (Opt.QuadrOne>=0 && Opt.QuadrOne!=q) continue;
cout << "Covering relations for quadrangle q: " << q << endl;
if (QS.Q[q].size()!=4){ cout << "This is not a quadrangle" << endl; continue;}
PolygonSet PSChecked;
BrokenLineSet PSCheckedHor;
for (int p=0;p<CoverPossible[q].size();++p){
cout << (CoverPossible[q][p]?"1":"0");
if (CoverPossible[q][p]){
PSChecked.push_back(Polygons[p]);
PSCheckedHor.push_back(QS.Q[p].S[1].P); // add horizontal edges
PSCheckedHor.push_back(QS.Q[p].S[3].P);
}
}
cout << endl;
// check if images of vertical edges of q are outside polygons for which covering relations are being checked
for (int s=0;s<QS.Q[q].size();++s){
if (s%2==1) continue; // skip horizontal edges
if (BrokenLineImageOutsidePolygonSet(QS.Q[q].S[s].P,PSChecked,syst.pm,XDisp,YDisp,TestNum,VCoverSingle,VCoverImSingle,tbound)){
for (int p=0;p<CoverPossible[q].size();++p) CoverDone[q][p]=CoverPossible[q][p];
} else {
for (int p=0;p<CoverPossible[q].size();++p) CoverDone[q][p]=false;
res=false;
cout << "Error" << endl;
}
}
// check if images of horizontal edges have empty intersection with horizontal edges of polygons in PSChecked
for (int s=0;s<QS.Q[q].size();++s){
if (s%2==0) continue; // skip vertical edges
if (BrokenLineImageNotInBrokenLineSet(QS.Q[q].S[s].P,PSCheckedHor,syst.pm,XDisp,YDisp,TestNum,HCoverSingle,HCoverImSingle,tbound)){
for (int p=0;p<CoverPossible[q].size();++p) CoverDone[q][p]=CoverPossible[q][p];
} else {
for (int p=0;p<CoverPossible[q].size();++p) CoverDone[q][p]=false;
res=false;
cout << "Error" << endl;
}
}
for (int p=0;p<HCoverSingle.size();++p) HCover.push_back(HCoverSingle[p]);
for (int p=0;p<VCoverSingle.size();++p) VCover.push_back(VCoverSingle[p]);
for (int p=0;p<HCoverImSingle.size();++p) HCoverIm.push_back(HCoverImSingle[p]);
for (int p=0;p<VCoverImSingle.size();++p) VCoverIm.push_back(VCoverImSingle[p]);
string VCoverFileNameTxt=Opt.BaseFile+"VC"+Int2String(q)+".txt";
string VCoverImFileNameTxt=Opt.BaseFile+"VCIm"+Int2String(q)+".txt";
string HCoverFileNameTxt=Opt.BaseFile+"HC"+Int2String(q)+".txt";
string HCoverImFileNameTxt=Opt.BaseFile+"HCIm"+Int2String(q)+".txt";
IVectorSetSave(VCoverSingle,VCoverFileNameTxt,Opt.CoutPrec);
IVectorSetSave(VCoverImSingle,VCoverImFileNameTxt,Opt.CoutPrec);
IVectorSetSave(HCoverSingle,HCoverFileNameTxt,Opt.CoutPrec);
IVectorSetSave(HCoverImSingle,HCoverImFileNameTxt,Opt.CoutPrec);
}
cout << "VCoverSize: " << VCover.size() << " HCoverSize: " << HCover.size() << " Time: " << TC << endl;
if (Opt.BaseFile!=""){
string VCoverFileNameTxt=Opt.BaseFile+"VC"+".txt";
string VCoverImFileNameTxt=Opt.BaseFile+"VCIm"+".txt";
string HCoverFileNameTxt=Opt.BaseFile+"HC"+".txt";
string HCoverImFileNameTxt=Opt.BaseFile+"HCIm"+".txt";
IVectorSetSave(VCover,VCoverFileNameTxt,Opt.CoutPrec);
IVectorSetSave(VCoverIm,VCoverImFileNameTxt,Opt.CoutPrec);
IVectorSetSave(HCover,HCoverFileNameTxt,Opt.CoutPrec);
IVectorSetSave(HCoverIm,HCoverImFileNameTxt,Opt.CoutPrec);
}
cout << "ReturnTimeBound: " << tbound << endl;
cout << "CoveringTable: " << endl;
for (int q=0;q<CoverDone.size();++q){
// cout << q << " ";
for (int p=0;p<CoverDone[q].size();++p) cout << (CoverDone[q][p]?"1":"0");
cout << endl;
}
return res;
}
void VectorSetAndQSSaveEPS(string InFile, string QuadrFile, string EpsFile, string info=""){
DVectorSet VSet;
QuadrSet QS;
if (InFile!="") VectorSetLoad(InFile,VSet);
if (QuadrFile!="") QS.Read(QuadrFile);
VectorSetAndQSSaveEPS(VSet,QS,EpsFile,Opt.EPSOpt,false,false,info);
}
void IVectorSetsAndQSSaveEPS(string InFile1, string InFile2, string QuadrFile, string EpsFile, string info=""){
IVectorSet IVSet1;
IVectorSet IVSet2;
QuadrSet QS;
if (InFile1!="") IVectorSetLoad(InFile1,IVSet1);
if (InFile2!="") IVectorSetLoad(InFile2,IVSet2);
if (QuadrFile!="") QS.Read(QuadrFile);
IVectorSetsAndQSSaveEPS(IVSet1,IVSet2,QS,EpsFile,Opt.EPSOpt,false,false,info);
}
int main(int argc, char *argv[]){
TimeCount TC;
cout << "StartDay: " << CurrentDay() << " StartHour: " << CurrentHour() << " HostName: " << HostName() << endl;
CommandLineOptions(argc,argv);
cout.precision(Opt.CoutPrec);
MpFloat::setDefaultPrecision(Opt.MpPrec);
if (Opt.Action=="Trajectory") Trajectory(Opt.TimeStep,Opt.SkipTime,Opt.IntegrationTime); // Action==1
if (Opt.Action=="TrajectoryMp") TrajectoryMp(Opt.TimeStep,Opt.SkipTime,Opt.IntegrationTime); // Action==1
if (Opt.Action=="PlotTrajectoryEps") PlotTrajectoryEps(Opt.TimeStep,Opt.SkipTime,Opt.IntegrationTime,Opt.EpsFileName,Opt.EPSOpt);
if (Opt.Action=="TrajectoryInterval") ComputeTrajectoryInterval(Opt.TimeStep,Opt.IntegrationTime); // Action==2
if (Opt.Action=="ReturnMapTrajectory") ReturnMapTrajectory(Opt.TimeStep,Opt.SkipNum,Opt.IterNum,true,false);
if (Opt.Action=="TestPoincMapI") TestPoincMapI();
// if (Opt.Action=="ReturnMapTrajectoryI") ReturnMapTrajectoryI(Opt.TimeStep,Opt.SkipNum,Opt.IterNum,true,false);
if (Opt.Action=="MpReturnMapTrajectory") MpReturnMapTrajectory(Opt.TimeStep,Opt.SkipNum,Opt.IterNum,true,false);
if (Opt.Action=="BiffDiag") BiffDiag(); // Action==5
if (Opt.Action=="BiffDiagEps") BiffDiagEps(Opt.InFile,Opt.BinNum,Opt.BiffParMin,Opt.BiffParMax,Opt.XMin,Opt.XMax,Opt.EpsFileName,Opt.EPSOpt); // Action==7
if (Opt.Action=="BiffDiagPGM") BiffDiagPGM(Opt.InFile,Opt.BinNum,Opt.BiffParMin,Opt.BiffParMax,Opt.XMin,Opt.XMax,Opt.EpsFileName); // Action==8
if (Opt.Action=="BiffDiagPGM2") BiffDiagPGM2(Opt.InFile,Opt.BinNum,Opt.BiffParMin,Opt.BiffParMax,Opt.BiffParNum,Opt.XMin,Opt.XMax,Opt.EpsFileName); // Action==8
if (Opt.Action=="FindPO") FindPO(Opt.Period);
if (Opt.Action=="FindPOMp") FindPOMp(Opt.Period);
if (Opt.Action=="ProvePO") ProvePO(Opt.Period);
if (Opt.Action=="ProvePOMp") ProvePOMp(Opt.Period);
if (Opt.Action=="smcc3FPs") smcc3FPs<Interval>();
if (Opt.Action=="smcc3FPsMp") smcc3FPs<MpInterval>();
if (Opt.Action=="smcc5FPs") smcc5FPs<Interval>();
if (Opt.Action=="smcc5FPsMp") smcc5FPs<MpInterval>();
if (Opt.Action=="smcc5FPBorders") smcc5FPBorders<Interval>();
if (Opt.Action=="smcc5FPBordersMp") smcc5FPBorders<MpInterval>();
if (Opt.Action=="QuadrSet2PolySet") QuadrSet2PolySet(Opt.QuadrFile,Opt.OutFile);
if (Opt.Action=="QuadrSet2QuadrSet") QuadrSet2QuadrSet(Opt.QuadrFile,Opt.OutFile);
if (Opt.Action=="TrappingRegionBorder") TrappingRegionBorder(Opt.QuadrFile);
if (Opt.Action=="TrappingRegionImageExists") TrappingRegionInterior(Opt.QuadrFile,false);
if (Opt.Action=="TrappingRegionImageInside") TrappingRegionInterior(Opt.QuadrFile,true);
if (Opt.Action=="CoveringRelation") CoveringRelation(Opt.QuadrFile);
if (Opt.Action=="QuadrSetSaveEPS") QuadrSetSaveEPS(Opt.QuadrFile,Opt.EpsFileName,Opt.SingleColor);
if (Opt.Action=="IVectorSetsAndQSSaveEPS") IVectorSetsAndQSSaveEPS(Opt.InFile,Opt.InFile2,Opt.QuadrFile,Opt.EpsFileName,Opt.Info);
if (Opt.Action=="VectorSetAndQSSaveEPS") VectorSetAndQSSaveEPS(Opt.InFile,Opt.QuadrFile,Opt.EpsFileName,Opt.Info);
cout << "Time: " << TC << " Elapsed: " << TC.ElapsedTimeString() << " HostName: " << HostName() << endl;
return 0;
}