Skip to content

Commit

Permalink
code for lds test for lyapunov exponents
Browse files Browse the repository at this point in the history
  • Loading branch information
@VictorHuynh
VictorHuynh committed Aug 14, 2024
1 parent b7e961a commit e79a3b6
Show file tree
Hide file tree
Showing 2 changed files with 106 additions and 0 deletions.
82 changes: 82 additions & 0 deletions demos/complexity_checker.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -92,6 +92,88 @@
"## Expected Result:\n",
" # Rossler ≈ 2.01"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"ename": "TypeError",
"evalue": "Value after * must be an iterable, not numpy.float64",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)",
"Cell \u001b[0;32mIn[4], line 19\u001b[0m\n\u001b[1;32m 16\u001b[0m x\u001b[38;5;241m=\u001b[39mnp\u001b[38;5;241m.\u001b[39marray(x)\n\u001b[1;32m 17\u001b[0m tpts \u001b[38;5;241m=\u001b[39m [tpts]\n\u001b[0;32m---> 19\u001b[0m lyapunov_spectrum \u001b[38;5;241m=\u001b[39m \u001b[43mfind_lyapunov_exponents\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtpts\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtimesteps\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmodel\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mprecomp\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mTrue\u001b[39;49;00m\u001b[43m)\u001b[49m\n\u001b[1;32m 21\u001b[0m \u001b[38;5;28mprint\u001b[39m(lyapunov_spectrum)\n",
"File \u001b[0;32m~/vincent-dynadojo/src/dynadojo/utils/complexity_measures.py:310\u001b[0m, in \u001b[0;36mfind_lyapunov_exponents\u001b[0;34m(trajectory, tpts, traj_length, model, pts_per_period, precomp, tol, min_tpts, **kwargs)\u001b[0m\n\u001b[1;32m 308\u001b[0m x \u001b[38;5;241m=\u001b[39m trajectory[i]\n\u001b[1;32m 309\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m precomp:\n\u001b[0;32m--> 310\u001b[0m jacval \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39marray(\u001b[43mmodel\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mjac\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mt\u001b[49m\u001b[43m)\u001b[49m)\n\u001b[1;32m 311\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m 312\u001b[0m rhsy \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mlambda\u001b[39;00m a: np\u001b[38;5;241m.\u001b[39marray(model\u001b[38;5;241m.\u001b[39mrhs(a, t)) \u001b[38;5;66;03m# define a function 'rhsy', using the model's right hand side diffeq\u001b[39;00m\n",
"File \u001b[0;32m~/vincent-dynadojo/src/dynadojo/systems/gilpin_lds.py:10\u001b[0m, in \u001b[0;36mDynSys.jac\u001b[0;34m(self, X, t)\u001b[0m\n\u001b[1;32m 8\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mjac\u001b[39m(\u001b[38;5;28mself\u001b[39m, X, t):\n\u001b[1;32m 9\u001b[0m \u001b[38;5;250m \u001b[39m\u001b[38;5;124;03m\"\"\"The Jacobian of the dynamical system\"\"\"\u001b[39;00m\n\u001b[0;32m---> 10\u001b[0m out \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_jac(\u001b[38;5;241m*\u001b[39mX\u001b[38;5;241m.\u001b[39mT, t, alpha\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m-\u001b[39m\u001b[38;5;241m1\u001b[39m)\n\u001b[1;32m 11\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m out\n",
"\u001b[0;31mTypeError\u001b[0m: Value after * must be an iterable, not numpy.float64"
]
}
],
"source": [
"import numpy as np\n",
"from dynadojo.systems.gilpin_lds import LDS\n",
"from dynadojo.utils.complexity_measures import find_lyapunov_exponents\n",
"\n",
"\n",
"model = LDS()\n",
"\n",
"alpha = -1 # Rate of change\n",
"x0 = 5 # Initial condition\n",
"t_start = 0 # Start time\n",
"t_end = 100 # End time\n",
"timesteps=100\n",
"\n",
"tpts = np.linspace(t_start, t_end, timesteps)\n",
"x = x0 * np.exp(alpha * tpts)\n",
"x=np.array(x)\n",
"tpts = [tpts]\n",
"\n",
"lyapunov_spectrum = find_lyapunov_exponents(x, tpts, timesteps, model, precomp=True)\n",
"\n",
"print(lyapunov_spectrum)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[ 0.06170239 -0.50872132 -12.96323588]\n"
]
}
],
"source": [
"\n",
"from dynadojo.systems.gilpin_flows import GilpinFlowsSystem\n",
"from dynadojo.wrappers import SystemChecker\n",
"from dynadojo.utils.complexity_measures import find_lyapunov_exponents\n",
"\n",
"dimension=3\n",
"system_name=\"Lorenz\"\n",
"seed=1\n",
"timesteps=1000\n",
"max_timesteps=1000\n",
"\n",
"system = SystemChecker(GilpinFlowsSystem(latent_dim=dimension, embed_dim=dimension, system_name=system_name, seed=seed))\n",
"unwrapped_system = system._system # Reach under SystemChecker to enable \"return_times\" kwarg to GilpinFlowsSystem's \"make_data\" method\n",
"model = unwrapped_system.system # Reach under GilpinFlowsSystem to access the model for Lyapunov\n",
"\n",
"# Generate in distribution trajectory data\n",
"x0 = system.make_init_conds(1)\n",
"xtpts, x = unwrapped_system.make_data(x0, timesteps=max_timesteps, return_times=True)\n",
"x0 = x0[0]\n",
"x = x[0]\n",
"xlyapunov_spectrum = find_lyapunov_exponents(x, xtpts, timesteps, model)\n",
"\n",
"print(xlyapunov_spectrum)"
]
}
],
"metadata": {
Expand Down
24 changes: 24 additions & 0 deletions src/dynadojo/systems/gilpin_lds.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,24 @@
from dysts.base import BaseDyn, staticjit

class DynSys:
def __init__(self, alpha):
self.alpha = alpha

def rhs(self, X, t):
"""The right hand side of a dynamical equation"""
out = self._rhs(*X.T, t, alpha=self.alpha)
return out
def jac(self, X, t):
"""The Jacobian of the dynamical system"""
out = self._jac(*X.T, t, alpha=self.alpha)
return out

class LDS(DynSys):
@staticjit
def _rhs(x, t, alpha):
xdot = alpha * x
return xdot
@staticjit
def _jac(x, t, alpha):
row1 = [alpha]
return row1

0 comments on commit e79a3b6

Please sign in to comment.