{ "cells": [ { "cell_type": "code", "execution_count": 81, "metadata": {}, "outputs": [], "source": [ "import scipy.integrate\n", "import numpy\n", "import matplotlib.pyplot\n", "import os\n", "import json\n", "import scipy.interpolate\n", "\n", "def dfdy(t,y,system):\n", " dfdy=system.M.dot(y)+system.u(t)\n", " return dfdy\n", "\n", "def jacobi(t,y,system):\n", " return system.M\n", "\n", "#SE post calculation\n", "def dfdyS(t,S,system):\n", " #unwrap S to NxM where M is number of parameters\n", " mS=numpy.reshape(S,(system.n,system.m))\n", " mOut=system.M.dot(mS)+system.fS(t)\n", " return numpy.ravel(mOut)\n", "\n", "def jacobiSE(t,S,system):\n", " N=system.n*(system.m)\n", " fJ=numpy.zeros((N,N))\n", " #print('fJ shape {}'.format(fJ.shape))\n", " for i in range(system.m):\n", " fJ[i*system.n:(i+1)*system.n,i*system.n:(i+1)*system.n]=system.M\n", " return fJ\n", "\n", "#SE simultaeneous calculation\n", "def dfdySFull(t,S,system):\n", " #unwrap S to NxM where M is number of parameters\n", " mS=numpy.reshape(S,(system.n,system.m+1))\n", " #system.fS(y,t) is NxM matrix where M are parameters\n", " y=mS[:,0]\n", " mOut=system.M.dot(mS)+system.fSY(y,t)\n", " return numpy.ravel(mOut)\n", "\n", "def jacobiSEFull(t,S,system):\n", " N=system.n*(system.m+1)\n", " fJ=numpy.zeros((N,N))\n", " #print('fJ shape {}'.format(fJ.shape))\n", " for i in range(system.m+1):\n", " fJ[i*system.n:(i+1)*system.n,i*system.n:(i+1)*system.n]=system.M\n", " return fJ\n", "\n", "class model:\n", " def __init__(self):\n", " self.compartments={}\n", " self.seJ={}\n", " \n", " def add_source(self,compartmentName,formula):\n", " self.compartments[compartmentName]['source']=formula\n", "\n", " def add_compartment(self,compartmentName):\n", " self.compartments[compartmentName]={}\n", " self.compartments[compartmentName]['targets']={}\n", " self.compartments[compartmentName]['sensTargets']={}\n", "\n", " def bind(self,sourceCompartment,targetCompartment,k,pc,scaleToVolume=0):\n", " #establish a flow from one compartment to the other\n", " cSrc=self.compartments[sourceCompartment]\n", " cTarg=self.compartments[targetCompartment]\n", " \n", " vSPar=self.getVolumePar(sourceCompartment,scaleToVolume)\n", " vTPar=self.getVolumePar(targetCompartment,scaleToVolume)\n", " \n", " #the source equation (where we subtract the current)\n", " addValue(cSrc['targets'],sourceCompartment,-get(k)/get(pc)/get(vSPar))\n", " #the target equation (where we add the current)\n", " addValue(cTarg['targets'],sourceCompartment,get(k)/get(pc)/get(vTPar))\n", " \n", " \n", " def getDerivative(self,variable, sign, qPar, pcPar, vPar):\n", " \n", " #for flow based transfer, k=Q/vS/pc, and jacobian derivative is -Q/V/pc/pc\n", " #for diffusion dominated transfer, k=k1/pc\n", " q=get(qPar)\n", " pc=get(pcPar)\n", " v=get(vPar)\n", " if variable==\"partitionCoefficient\":\n", " return sign*(-1)*q/v/pc/pc\n", " if variable==\"diffusionCoefficient\":\n", " return sign/pc/v\n", " if variable==\"flow\":\n", " return sign/pc/v\n", " if variable==\"volume\":\n", " return sign*(-1)*q/pc/v/v\n", " \n", " def getVolumePar(self,compartment,useVolume=1):\n", " if not useVolume:\n", " return {\"value\":1}\n", " try:\n", " parVName=self.mod[\"volumes\"][compartment]\n", " parV=self.mod[\"parameters\"][parVName]\n", " except KeyError:\n", " return {\"value\":1}\n", " \n", " parV[\"name\"]=parVName\n", " return parV\n", " \n", " def bindDerivative(self,parameterName,variable,src,target,q,pc,useVolume=1):\n", " parName=parameterName\n", " srcV=self.getVolumePar(src,useVolume)\n", " doAdd=(variable!=\"volume\" or calculateDerivative(srcV))\n", " \n", " if doAdd:\n", " #the source equation (where we subtract the current)\n", " if variable==\"volume\":\n", " parName=srcV[\"name\"]\n", " \n", " csrc=self.getSEJ_comp(parName,src)\n", " addValue(csrc,src,self.getDerivative(variable,-1,q,pc,srcV))\n", " \n", " targetV=self.getVolumePar(target,useVolume)\n", " doAdd=(variable!=\"volume\" or calculateDerivative(targetV))\n", " \n", " #the target equation (where we add the current)\n", " if doAdd:\n", " if variable==\"volume\":\n", " parName=targetV[\"name\"]\n", "\n", " ctarget=self.getSEJ_comp(parName,target)\n", " addValue(ctarget,src,self.getDerivative(variable,+1,q,pc,targetV))\n", " \n", " \n", " \n", " def build(self):\n", " comps=self.compartments\n", " self.n=len(comps)\n", " self.M=numpy.zeros((self.n,self.n))\n", " self.fu=[lambda t:0]*self.n\n", " self.lut={c:i for (i,c) in zip(range(self.n),comps.keys())}\n", " for c in comps:\n", " comp=comps[c]\n", " if 'source' in comp:\n", " self.fu[self.lut[c]]=parseFunction(comp['source'])\n", " \n", " for t in comp['targets']:\n", " self.M[self.lut[c],self.lut[t]]=comp['targets'][t]\n", " #build SE part\n", " self.m=len(self.seJ)\n", " #MxNxN matrix\n", " self.lutSE={c:i for (i,c) in zip(range(self.m),self.seJ.keys())}\n", " self.fSS=numpy.zeros((self.m,self.n,self.n))\n", " for par in self.seJ:\n", " sources=self.seJ[par]\n", " for compartment in sources:\n", " targets=sources[compartment]\n", " for t in targets:\n", " #print('FSS: Adding {}/{},{}/{},{}/{}:{}'.\\\n", " # format(par,self.lutSE[par],compartment,self.lut[compartment],t,self.lut[t],targets[t]))\n", " self.fSS[self.lutSE[par],self.lut[compartment],self.lut[t]]=targets[t]\n", " \n", " \n", " \n", " #print('Done')\n", " \n", " def inspect(self):\n", " comps=self.compartments\n", " pars=self.mod['parameters']\n", " print('Compartments')\n", " for c in comps:\n", " print('{}/{}:'.format(c,self.lut[c]))\n", " comp=comps[c]\n", " if 'source' in comp:\n", " print('\\tsource\\n\\t\\t{}'.format(comp['source']))\n", " print('\\ttargets')\n", " for t in comp['targets']:\n", " print('\\t\\t{}[{},{}]: {}'.format(t,self.lut[c],self.lut[t],comp['targets'][t]))\n", " print('Flows')\n", " for f in self.flows:\n", " fName=self.flows[f]\n", " fParName=self.mod['flows'][fName]\n", " fPar=pars[fParName]\n", " print('\\t{}[{}]:{} [{}]'.format(f,fName,fParName,get(fPar)))\n", " \n", " print('Volumes')\n", " for v in self.mod['volumes']:\n", " vParName=self.mod['volumes'][v]\n", " vPar=pars[vParName]\n", " print('\\t{}:{} [{}]'.format(v,vParName,get(vPar)))\n", " \n", " print('Partition coefficients')\n", " for pc in self.mod['partitionCoefficients']:\n", " pcParName=self.mod['partitionCoefficients'][pc]\n", " pcPar=pars[pcParName]\n", " print('\\t{}:{} [{}]'.format(pc,pcParName,get(pcPar)))\n", " \n", " \n", " print('SE parameters')\n", " for p in self.seJ:\n", " print(p)\n", " sources=self.seJ[p]\n", " for compartment in sources:\n", " targets=sources[compartment]\n", " for t in targets:\n", " print('\\t SE bind {}/{}:{}'.format(compartment,t,targets[t]))\n", " #print('Done')\n", " \n", " def parse(self,file):\n", " \n", " with open(file,'r') as f:\n", " self.mod=json.load(f)\n", " for m in self.mod['compartments']:\n", " self.add_compartment(m)\n", " for s in self.mod['sources']:\n", " #src=mod['sources'][s]\n", " self.add_source(s,self.mod['sources'][s])\n", " self.flows={}\n", " pars=self.mod['parameters']\n", " for f in self.mod['flows']:\n", " #skip comments\n", " if f.find(':')<0:\n", " continue\n", " \n", " comps=f.split(':')\n", " c0=splitVector(comps[0])\n", " c1=splitVector(comps[1])\n", " for x in c0:\n", " for y in c1:\n", " pairName='{}:{}'.format(x,y)\n", " self.flows[pairName]=f\n", " \n", " for b in self.mod['bindings']['diffusion']:\n", " comps=b.split('->')\n", " try:\n", " pcParName=self.mod['partitionCoefficients'][b]\n", " pcPar=pars[pcParName]\n", " except KeyError:\n", " pcPar={\"value\":1}\n", " \n", " kParName=self.mod['bindings']['diffusion'][b]\n", " kPar=pars[kParName]\n", " self.bind(comps[0],comps[1],kPar,pcPar)\n", " \n", " #add derivatives calculateDerivative returns true\n", " if calculateDerivative(kPar):\n", " self.bindDerivative(kParName,\"diffusionParameter\",comps[0],comps[1],kPar,pcPar,0)\n", " if calculateDerivative(pcPar):\n", " self.bindDerivative(pcParName,\"partitionCoefficient\",comps[0],comps[1],kPar,pcPar,0)\n", " \n", " for q in self.mod['bindings']['flow']:\n", " comps=q.split('->')\n", " srcs=splitVector(comps[0])\n", " tgts=splitVector(comps[1])\n", " for cs in srcs:\n", " for ct in tgts:\n", " #get partition coefficient\n", " try:\n", " pcParName=self.mod['partitionCoefficients'][cs]\n", " pcPar=pars[pcParName]\n", " except KeyError:\n", " pcPar={\"value\":1}\n", " \n", " #get flow (direction could be reversed)\n", " try:\n", " qName=self.flows['{}:{}'.format(cs,ct)]\n", " except KeyError:\n", " qName=self.flows['{}:{}'.format(ct,cs)]\n", " \n", " flowParName=self.mod['flows'][qName]\n", " flowPar=pars[flowParName]\n", " \n", " self.bind(cs,ct,flowPar,pcPar,1)\n", " \n", " if calculateDerivative(pcPar):\n", " self.bindDerivative(pcParName,\"partitionCoefficient\",cs,ct,flowPar,pcPar)\n", " if calculateDerivative(flowPar):\n", " self.bindDerivative(flowParName,\"flow\",cs,ct,flowPar,pcPar)\n", " self.bindDerivative(\"x\",\"volume\",cs,ct,flowPar,pcPar)\n", " \n", " #print('bind: {}:(q={},vt={},pc={}):{}'.format(bind,q,vt,pc,q/vt/pc))\n", " \n", " self.build()\n", " \n", " def u(self,t):\n", " ub=[f(t) for f in self.fu]\n", " return numpy.array(ub)\n", " \n", " def fSY(self,y,t):\n", " #M number of sensitivity parameters\n", " #N number of equations\n", " #fSS is MxNxN\n", " \n", " #assume a tabulated solution y(t) at t spaced intervals\n", " \n", " qS=self.fSS.dot(y)\n", " #qS is MxN\n", " #but NxM is expected, so do a transpose\n", " \n", " #for simultaneous calculation, a Nx(M+1) matrix is expected\n", " tS=numpy.zeros((self.n,self.m+1))\n", " #columns from 2..M+1 are the partial derivatives \n", " tS[:,1:]=numpy.transpose(qS)\n", " #first column is the original function\n", " tS[:,0]=self.u(t)\n", " return tS\n", " \n", " def fS(self,t):\n", " #M number of sensitivity parameters\n", " #N number of equations\n", " #fSS is MxNxN\n", " \n", " #assume a tabulated solution y(t) at t spaced intervals\n", " \n", " qS=self.fSS.dot(self.getY(t))\n", " return numpy.transpose(qS)\n", " \n", " def getSEJ(self,parName):\n", " #find the sensitivity (SE) derivative of Jacobi with respect to parameter \n", " try:\n", " return self.seJ[parName]\n", " except KeyError:\n", " self.seJ[parName]={}\n", " return self.seJ[parName]\n", " \n", " def getSEJ_comp(self,parName,compartmentName):\n", " #find equation dictating concentration in compartmentName for jacobi-parameter derivative\n", " seJ=self.getSEJ(parName)\n", " \n", " try:\n", " return seJ[compartmentName]\n", " except KeyError:\n", " seJ[compartmentName]={}\n", " return seJ[compartmentName]\n", " def setY(self,t,y):\n", " self.tck=[None]*self.n\n", " for i in range(self.n):\n", " self.tck[i] = scipy.interpolate.splrep(t, y[:,i], s=0)\n", " \n", " def getY(self,t):\n", " fY=numpy.zeros(self.n)\n", " for i in range(self.n):\n", " fY[i]=scipy.interpolate.splev(t, self.tck[i], der=0)\n", " return fY\n", " \n", " def calculateUncertainty(self,sol,se):\n", " s2out=numpy.zeros(sol.shape)\n", " pars=self.mod['parameters']\n", " for parName in pars:\n", " par=pars[parName]\n", " if not calculateDerivative(par):\n", " continue\n", " v=par[\"value\"]\n", " if par['dist']=='lognormal':\n", " #this is sigma^2_lnx\n", " sln2=numpy.log(par[\"cv\"]*par[\"cv\"]+1)\n", " #have to multiplied by value to get the derivative with respect to lnx\n", " s2=sln2*v*v\n", " else:\n", " #for Gaussian, cv is sigma/value; get sigma by value multiplication\n", " s2=par[\"cv\"]*par[\"cv\"]*v*v\n", " j=self.lutSE[parName]\n", " fse=se[:,:,j]\n", " print('Calculating for {}/{}:{}'.format(parName,j,fse.shape))\n", " #se2 is nt x n\n", " se2=numpy.multiply(fse,fse)\n", " s2out+=se2*s2\n", " return numpy.sqrt(s2out)\n", " \n", "def splitVector(v):\n", " if v.find('(')<0:\n", " return [v]\n", " return v[1:-1].split(',')\n", "\n", "def parseFunction(formula):\n", " if formula['name']=='exponential':\n", " c0=formula['constant']\n", " k=formula['k']\n", " return lambda t:c0*numpy.exp(k*t)\n", " if formula['name']=='constant':\n", " c0=formula['value']\n", " return lambda t:c0\n", " if formula['name']=='Heavyside':\n", " t0=formula['limit']\n", " v=formula['value']\n", " return lambda t:v if t" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "#copy output of Thompson et al\n", "print(s1.shape)\n", "fig, axs = matplotlib.pyplot.subplots(5, 3,figsize=(15,25))\n", "name=['arterial','adipose','brain','heart','kidney','liver','lung','muscle','skin',\n", " 'splanchnic','stomach','excrement']\n", "#diazepam\n", "max=[1.5,2.6,3,4,5,2.5,6.8,1.5,1.5,4,4.2,3]\n", "#cotinine\n", "max=[9]*12\n", "max[11]=90\n", "\n", "max=[1000*x for x in max]\n", "for i in range(len(name)):\n", " row=i//3\n", " col=i%3\n", " fy=sol[:,sys.lut[name[i]]]\n", " fe=se[:,sys.lut[name[i]]]\n", " ax=axs[row,col]\n", " ax.plot(t/60,fy)\n", " ax.fill_between(t/60, fy-fe, fy + fe, color='red',alpha=0.1)\n", " ax.plot(t/60,fy-fe,color='red',linewidth=1,alpha=0.2)\n", " ax.plot(t/60,fy+fe,color='red',linewidth=1,alpha=0.2)\n", " axs[row,col].set_title(name[i])\n", " axs[row,col].set_ylim([0,max[i]])\n", " axs[row,col].set_xlim([0,250])\n" ] }, { "cell_type": "code", "execution_count": 89, "metadata": {}, "outputs": [ { "data": { "text/plain": [ " message: 'The solver successfully reached the end of the integration interval.'\n", " nfev: 1135\n", " njev: 25\n", " nlu: 25\n", " sol: None\n", " status: 0\n", " success: True\n", " t: array([0.00000000e+00, 3.20591951e-04, 6.41183902e-04, 1.28236780e-03,\n", " 1.92355170e-03, 2.56473561e-03, 8.97657462e-03, 1.53884136e-02,\n", " 2.18002527e-02, 2.82120917e-02, 5.40150345e-02, 7.98179774e-02,\n", " 1.05620920e-01, 1.31423863e-01, 1.57226806e-01, 1.99625129e-01,\n", " 2.42023452e-01, 2.84421775e-01, 3.26820098e-01, 3.69218420e-01,\n", " 4.11616743e-01, 4.79107559e-01, 5.46598375e-01, 6.14089191e-01,\n", " 6.81580007e-01, 7.49070823e-01, 8.16561639e-01, 8.84052455e-01,\n", " 9.51550514e-01, 1.01904857e+00, 1.08654663e+00, 1.15404469e+00,\n", " 1.22154275e+00, 1.28904081e+00, 1.35653887e+00, 1.44133130e+00,\n", " 1.52612374e+00, 1.61091617e+00, 1.69570861e+00, 1.78050104e+00,\n", " 1.86529348e+00, 1.96079513e+00, 2.05629678e+00, 2.15179843e+00,\n", " 2.24730008e+00, 2.34280173e+00, 2.46429320e+00, 2.58578468e+00,\n", " 2.70727615e+00, 2.82876763e+00, 2.95025910e+00, 3.07228543e+00,\n", " 3.19431176e+00, 3.31633809e+00, 3.43836441e+00, 3.56039074e+00,\n", " 3.70535909e+00, 3.85032743e+00, 3.99529578e+00, 4.14026413e+00,\n", " 4.28523247e+00, 4.61325179e+00, 4.94127111e+00, 5.26929043e+00,\n", " 5.59730975e+00, 5.92532907e+00, 6.25334839e+00, 6.49130662e+00,\n", " 6.72926484e+00, 6.96722306e+00, 7.20518128e+00, 7.44313950e+00,\n", " 7.68109773e+00, 7.91905595e+00, 8.15701417e+00, 8.38073440e+00,\n", " 8.60445462e+00, 8.82817485e+00, 9.05189507e+00, 9.27561530e+00,\n", " 9.49933552e+00, 9.72305574e+00, 9.94677597e+00, 1.00874168e+01,\n", " 1.02280576e+01, 1.03686985e+01, 1.05093393e+01, 1.07277551e+01,\n", " 1.09461710e+01, 1.11645868e+01, 1.13830026e+01, 1.16014185e+01,\n", " 1.20719265e+01, 1.25424345e+01, 1.30129425e+01, 1.34834505e+01,\n", " 1.39539585e+01, 1.40376119e+01, 1.41212653e+01, 1.42049187e+01,\n", " 1.42885721e+01, 1.43722255e+01, 1.44558790e+01, 1.45634334e+01,\n", " 1.46709877e+01, 1.47785421e+01, 1.48860965e+01, 1.49936509e+01,\n", " 1.51396824e+01, 1.52857140e+01, 1.54317455e+01, 1.55777770e+01,\n", " 1.59004276e+01, 1.62230782e+01, 1.65457289e+01, 1.68683795e+01,\n", " 1.71910301e+01, 1.77869542e+01, 1.83828783e+01, 1.85318593e+01,\n", " 1.86808404e+01, 1.88298214e+01, 1.89788024e+01, 1.91277834e+01,\n", " 1.92767645e+01, 1.93982262e+01, 1.95196879e+01, 1.96411496e+01,\n", " 1.97626114e+01, 1.98840731e+01, 2.00155725e+01, 2.01470719e+01,\n", " 2.02785713e+01, 2.04100707e+01, 2.05938577e+01, 2.07776447e+01,\n", " 2.09614317e+01, 2.11452187e+01, 2.16199944e+01, 2.20947702e+01,\n", " 2.25695459e+01, 2.30443216e+01, 2.35190974e+01, 2.43992497e+01,\n", " 2.46192878e+01, 2.48393258e+01, 2.50593639e+01, 2.52794020e+01,\n", " 2.54994401e+01, 2.57194782e+01, 2.58285538e+01, 2.59376295e+01,\n", " 2.60467052e+01, 2.61557809e+01, 2.62648566e+01, 2.64129177e+01,\n", " 2.65609789e+01, 2.67090400e+01, 2.68571012e+01, 2.70099722e+01,\n", " 2.71628433e+01, 2.73157143e+01, 2.74685853e+01, 2.80819963e+01,\n", " 2.86954072e+01, 2.93088181e+01, 2.99222290e+01, 2.99605672e+01,\n", " 2.99989054e+01, 2.99990971e+01, 2.99992887e+01, 2.99994804e+01,\n", " 2.99996721e+01, 2.99998638e+01, 2.99998744e+01, 2.99998849e+01,\n", " 2.99999061e+01, 2.99999272e+01, 2.99999483e+01, 2.99999561e+01,\n", " 2.99999640e+01, 2.99999797e+01, 2.99999953e+01, 2.99999975e+01,\n", " 2.99999996e+01, 2.99999997e+01, 2.99999998e+01, 2.99999999e+01,\n", " 3.00000000e+01, 3.00000000e+01, 3.00000002e+01, 3.00000004e+01,\n", " 3.00000006e+01, 3.00000023e+01, 3.00000040e+01, 3.00000057e+01,\n", " 3.00000228e+01, 3.00000399e+01, 3.00000570e+01, 3.00002278e+01,\n", " 3.00003987e+01, 3.00005696e+01, 3.00007405e+01, 3.00024492e+01,\n", " 3.00041580e+01, 3.00058667e+01, 3.00075755e+01, 3.00092843e+01,\n", " 3.00263719e+01, 3.00434595e+01, 3.00605471e+01, 3.00776348e+01,\n", " 3.00947224e+01, 3.01118100e+01, 3.01673762e+01, 3.02229424e+01,\n", " 3.02785086e+01, 3.03340747e+01, 3.03896409e+01, 3.04452071e+01,\n", " 3.05007733e+01, 3.05682713e+01, 3.06357694e+01, 3.07032675e+01,\n", " 3.07707655e+01, 3.08382636e+01, 3.09057616e+01, 3.09732597e+01,\n", " 3.10407577e+01, 3.11082558e+01, 3.11757538e+01, 3.12552883e+01,\n", " 3.13348228e+01, 3.14143573e+01, 3.14938918e+01, 3.15734263e+01,\n", " 3.16529608e+01, 3.17474570e+01, 3.18419532e+01, 3.19364493e+01,\n", " 3.20309455e+01, 3.21254416e+01, 3.22199378e+01, 3.23414413e+01,\n", " 3.24629448e+01, 3.25844483e+01, 3.27059518e+01, 3.28274553e+01,\n", " 3.29489589e+01, 3.30704624e+01, 3.31919659e+01, 3.33153927e+01,\n", " 3.34388196e+01, 3.35622464e+01, 3.36856733e+01, 3.38091001e+01,\n", " 3.39982846e+01, 3.41874690e+01, 3.43766535e+01, 3.45658379e+01,\n", " 3.47550224e+01, 3.51264554e+01, 3.54978884e+01, 3.58693215e+01,\n", " 3.62407545e+01, 3.66121875e+01, 3.69836205e+01, 3.72054366e+01,\n", " 3.74272527e+01, 3.76490687e+01, 3.78708848e+01, 3.80927009e+01,\n", " 3.82360839e+01, 3.83794670e+01, 3.85228501e+01, 3.86662331e+01,\n", " 3.88403885e+01, 3.90145439e+01, 3.91886992e+01, 3.93628546e+01,\n", " 3.96674704e+01, 3.99720861e+01, 4.02767019e+01, 4.05813177e+01,\n", " 4.08859335e+01, 4.14733614e+01, 4.20607893e+01, 4.24024682e+01,\n", " 4.27441471e+01, 4.30858260e+01, 4.34275049e+01, 4.35129246e+01,\n", " 4.35983443e+01, 4.36837641e+01, 4.37691838e+01, 4.38546035e+01,\n", " 4.39400232e+01, 4.40608926e+01, 4.41817619e+01, 4.43026312e+01,\n", " 4.44235006e+01, 4.45443699e+01, 4.46912355e+01, 4.48381011e+01,\n", " 4.49849667e+01, 4.51318323e+01, 4.53290529e+01, 4.55262734e+01,\n", " 4.57234940e+01, 4.59207145e+01, 4.63338678e+01, 4.67470211e+01,\n", " 4.71601743e+01, 4.75733276e+01, 4.79864809e+01, 4.87601600e+01,\n", " 4.89535798e+01, 4.91469996e+01, 4.93404193e+01, 4.95338391e+01,\n", " 4.97272589e+01, 4.99206787e+01, 5.00421128e+01, 5.01635469e+01,\n", " 5.02849810e+01, 5.04064152e+01, 5.05278493e+01, 5.06757793e+01,\n", " 5.08237093e+01, 5.09716394e+01, 5.11195694e+01, 5.12536191e+01,\n", " 5.13876688e+01, 5.15217185e+01, 5.16557682e+01, 5.21486586e+01,\n", " 5.26415490e+01, 5.31344395e+01, 5.36273299e+01, 5.41202203e+01,\n", " 5.47944536e+01, 5.54686869e+01, 5.60045516e+01, 5.65404162e+01,\n", " 5.70762808e+01, 5.76121454e+01, 5.81480100e+01, 5.86838746e+01,\n", " 5.92197392e+01, 5.97556038e+01, 6.02914684e+01, 6.08273331e+01,\n", " 6.13631977e+01, 6.18990623e+01, 6.24349269e+01, 6.29707915e+01,\n", " 6.35066561e+01, 6.40425207e+01, 6.45783853e+01, 6.51142500e+01,\n", " 6.57794325e+01, 6.64446150e+01, 6.71097975e+01, 6.77749800e+01,\n", " 6.84401625e+01, 6.91053450e+01, 6.98529979e+01, 7.06006508e+01,\n", " 7.13483037e+01, 7.20959566e+01, 7.28436094e+01, 7.35912623e+01,\n", " 7.44334931e+01, 7.52757238e+01, 7.61179545e+01, 7.69601853e+01,\n", " 7.78024160e+01, 7.86446467e+01, 7.94868775e+01, 8.03291082e+01,\n", " 8.11713389e+01, 8.21778788e+01, 8.31844187e+01, 8.41909586e+01,\n", " 8.51974985e+01, 8.62040384e+01, 8.72105782e+01, 8.85692586e+01,\n", " 8.99279389e+01, 9.12866192e+01, 9.26452995e+01, 9.40039798e+01,\n", " 9.53626601e+01, 9.68856694e+01, 9.84086787e+01, 9.99316880e+01,\n", " 1.01454697e+02, 1.02977707e+02, 1.04500716e+02, 1.06290777e+02,\n", " 1.08080838e+02, 1.09870899e+02, 1.11660959e+02, 1.13451020e+02,\n", " 1.15241081e+02, 1.17473382e+02, 1.19705683e+02, 1.21937984e+02,\n", " 1.24170285e+02, 1.26402586e+02, 1.28634887e+02, 1.31735011e+02,\n", " 1.34835135e+02, 1.37935259e+02, 1.41035384e+02, 1.44135508e+02,\n", " 1.47235632e+02, 1.50998305e+02, 1.54760979e+02, 1.58523652e+02,\n", " 1.62286326e+02, 1.66048999e+02, 1.69811673e+02, 1.74086293e+02,\n", " 1.78360913e+02, 1.82635533e+02, 1.86910153e+02, 1.91184773e+02,\n", " 1.95459393e+02, 2.00340233e+02, 2.05221073e+02, 2.10101912e+02,\n", " 2.14982752e+02, 2.19863592e+02, 2.24744432e+02, 2.29625272e+02,\n", " 2.34506112e+02, 2.39386952e+02, 2.44267792e+02, 2.49148632e+02,\n", " 2.54029472e+02, 2.58910312e+02, 2.63791151e+02, 2.68671991e+02,\n", " 2.75508564e+02, 2.82345136e+02, 2.89181708e+02, 2.96018280e+02,\n", " 3.02854852e+02, 3.09691425e+02, 3.17582683e+02, 3.25473942e+02,\n", " 3.33365201e+02, 3.41256460e+02, 3.49147719e+02, 3.57038978e+02,\n", " 3.64930237e+02, 3.72821496e+02, 3.80712755e+02, 3.88604014e+02,\n", " 3.96495273e+02, 4.04386532e+02, 4.12277790e+02, 4.20169049e+02,\n", " 4.28060308e+02, 4.35951567e+02, 4.43842826e+02, 4.51734085e+02,\n", " 4.61149736e+02, 4.70565388e+02, 4.79981039e+02, 4.89396691e+02,\n", " 4.98812342e+02, 5.08227994e+02, 5.17643645e+02, 5.27059297e+02,\n", " 5.36474948e+02, 5.47097513e+02, 5.57720077e+02, 5.68342642e+02,\n", " 5.78965206e+02, 5.89587770e+02, 6.12387555e+02, 6.35187339e+02,\n", " 6.52316298e+02, 6.69445257e+02, 6.86574217e+02, 7.03703176e+02,\n", " 7.20832135e+02, 7.37961095e+02, 7.55090054e+02, 7.72219013e+02,\n", " 7.75644805e+02, 7.76266240e+02, 7.76887675e+02, 7.77509109e+02,\n", " 7.78751979e+02, 7.79994848e+02, 7.81237718e+02, 7.82480587e+02,\n", " 7.88083218e+02, 7.93685849e+02, 7.99288481e+02, 8.04891112e+02,\n", " 8.19221953e+02, 8.33552795e+02, 8.47883636e+02, 8.62214478e+02,\n", " 8.86883807e+02, 9.11553137e+02, 9.36222466e+02, 9.60891796e+02,\n", " 9.85561125e+02, 1.01434242e+03, 1.04312371e+03, 1.07190500e+03,\n", " 1.10068629e+03, 1.12946758e+03, 1.16998037e+03, 1.21049315e+03,\n", " 1.25100594e+03, 1.29151873e+03, 1.33203151e+03, 1.39046177e+03,\n", " 1.44889202e+03, 1.50732228e+03, 1.56575253e+03, 1.62418279e+03,\n", " 1.69696214e+03, 1.76974149e+03, 1.84252084e+03, 1.91530019e+03,\n", " 1.98807954e+03, 2.06085889e+03, 2.16122300e+03, 2.26158711e+03,\n", " 2.36195123e+03, 2.46231534e+03, 2.56267945e+03, 2.66304357e+03,\n", " 2.79861731e+03, 2.93419105e+03, 3.06976480e+03, 3.20533854e+03,\n", " 3.34091228e+03, 3.47648603e+03, 3.61205977e+03, 3.74763351e+03,\n", " 3.88320726e+03, 4.01878100e+03, 4.15435474e+03, 4.28992849e+03,\n", " 4.42550223e+03, 4.56107597e+03, 4.69664972e+03, 4.84924465e+03,\n", " 5.00183958e+03, 5.15443450e+03, 5.30702943e+03, 5.45962436e+03,\n", " 5.61221929e+03, 5.76481422e+03, 5.91740915e+03, 6.07000407e+03,\n", " 6.24122942e+03, 6.41245476e+03, 6.58368011e+03, 6.75490545e+03,\n", " 6.92613079e+03, 7.09735614e+03, 7.26858148e+03, 7.43980683e+03,\n", " 7.61103217e+03, 7.80531658e+03, 7.99960098e+03, 8.19388539e+03,\n", " 8.38816980e+03, 8.58245420e+03, 8.77673861e+03, 8.99236271e+03,\n", " 9.20798682e+03, 9.42361092e+03, 9.63923503e+03, 9.85485913e+03,\n", " 1.00704832e+04, 1.03087344e+04, 1.05469856e+04, 1.07852368e+04,\n", " 1.10234879e+04, 1.12617391e+04, 1.14999903e+04, 1.17673860e+04,\n", " 1.20347817e+04, 1.23021774e+04, 1.25695732e+04, 1.28369689e+04,\n", " 1.31043646e+04, 1.33717603e+04, 1.36391561e+04, 1.39065518e+04,\n", " 1.41739475e+04, 1.44000000e+04])\n", " t_events: None\n", " y: array([[0.00000000e+00, 0.00000000e+00, 1.19224476e-13, ...,\n", " 2.27072093e-01, 1.98648405e-01, 1.77412864e-01],\n", " [0.00000000e+00, 0.00000000e+00, 3.03909825e-13, ...,\n", " 1.54477792e+00, 1.35140976e+00, 1.20694535e+00],\n", " [0.00000000e+00, 0.00000000e+00, 1.96363686e-12, ...,\n", " 1.58519348e+00, 1.38676803e+00, 1.23851862e+00],\n", " ...,\n", " [0.00000000e+00, 2.91933365e-08, 1.16661187e-07, ...,\n", " 5.65008407e+00, 4.94283188e+00, 4.41444553e+00],\n", " [0.00000000e+00, 3.91109887e-01, 7.82207482e-01, ...,\n", " 5.63491242e+00, 4.92955788e+00, 4.40259621e+00],\n", " [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,\n", " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00]])\n", " y_events: None" ] }, "execution_count": 89, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol" ] }, { "cell_type": "code", "execution_count": 58, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0.00000000e+000, 9.46314230e+003, 6.45365940e+003, 5.67275261e+003,\n", " 6.26418952e-007, 1.79899429e-007, 1.79899429e-007, 3.22904947e-009,\n", " 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000,\n", " 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000,\n", " 0.00000000e+000, 5.56433140e+115, 0.00000000e+000, 0.00000000e+000,\n", " inf, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000,\n", " 0.00000000e+000, 0.00000000e+000, 4.49748572e-008, 4.49748572e-008,\n", " 7.19597716e-007, 4.49748572e-008, 4.49748572e-008, 0.00000000e+000,\n", " 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000,\n", " 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000,\n", " 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000,\n", " 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000,\n", " 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000,\n", " 0.00000000e+000, 0.00000000e+000, 2.15076749e+017, 0.00000000e+000,\n", " 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000,\n", " 3.11391673e-008, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000,\n", " 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000,\n", " 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000,\n", " 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000,\n", " 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000,\n", " 0.00000000e+000, inf, 0.00000000e+000, 0.00000000e+000,\n", " 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000,\n", " 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000,\n", " 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000,\n", " 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000,\n", " inf, 7.42698527e-313, 6.93839637e-310, 0.00000000e+000,\n", " 2.70839593e-316, 3.92400703e-013, 6.36598738e-313, 1.04267605e-017,\n", " 2.70844099e-316, 2.54639495e-313, 4.03179201e-313, 1.69759664e-313,\n", " 6.15378780e-313, 2.70852044e-316, 1.31563739e-312, 2.70883190e-316,\n", " 1.18575755e-322, 2.70855087e-316, 7.00258612e-313, 2.70857854e-316,\n", " nan, 5.09278990e-313, 0.00000000e+000, 0.00000000e+000,\n", " 2.70865008e-316, 2.70866984e-316, 1.33685735e-312, 0.00000000e+000,\n", " 6.93840202e-310, 1.69759664e-313, 1.18831764e-312, 1.18575755e-322,\n", " 4.94065646e-324, 8.06358401e-313, 6.93839888e-310, 9.54898107e-313,\n", " 0.00000000e+000, 0.00000000e+000, 1.61271680e-312, 6.93839890e-310,\n", " 2.70938723e-316, 1.18575755e-322, 3.18299369e-313, 2.70888368e-316,\n", " 6.93840200e-310, 2.54639496e-313, 0.00000000e+000, 0.00000000e+000,\n", " 4.94065646e-324, 0.00000000e+000, 4.94065646e-324, 6.93832955e-310,\n", " 4.29464389e-317, 2.54639495e-313, 6.93832870e-310, 4.29472294e-317,\n", " 4.94065646e-324, 6.93832948e-310, 4.29476246e-317, 8.29771096e+039,\n", " 4.94065646e-324, 1.48219694e-323, 1.18575755e-322, 2.54639496e-313,\n", " 2.70915561e-316, 2.70916510e-316, 4.94065646e-324, 0.00000000e+000,\n", " 2.70919553e-316, 7.85138443e-313, 9.88131292e-324, 4.94065646e-324,\n", " 1.08221785e-312, 6.93839637e-310, 6.79038654e-313, 6.57818696e-313,\n", " 4.94065646e-324, 4.94065646e-324, 1.48219694e-323, 1.48219694e-323,\n", " 1.48219694e-323, 4.24399159e-313, 1.90979621e-312, 0.00000000e+000,\n", " 4.94065646e-324, 6.93839886e-310, 7.42698528e-313, 1.08694442e-322,\n", " 0.00000000e+000, 4.94065646e-324, 8.39911598e-323, 0.00000000e+000,\n", " 9.88131292e-323, 6.93840203e-310, 2.70952043e-316, 2.54639496e-313,\n", " 0.00000000e+000, 4.45619117e-313, 0.00000000e+000, 0.00000000e+000,\n", " 1.90979824e-313])" ] }, "execution_count": 58, "metadata": {}, "output_type": "execute_result" } ], "source": [ "fse=s1[:,:,0]\n", "numpy.argwhere(numpy.isnan(fse))\n", "fse[:,0]" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[3. 2.]\n", " [2. 4.]\n", " [2. 6.]]\n", "[[3. 2.]\n", " [2. 4.]\n", " [2. 6.]]\n" ] } ], "source": [ "M=numpy.ones((3,2,2))\n", "M[0,0,1]=2\n", "M[1,1,0]=3\n", "M[2,1,1]=5\n", "v=numpy.ones(2)\n", "q=M.dot(v)\n", "q1=q.ravel()\n", "q2=numpy.reshape(q1,q.shape)\n", "print(q)\n", "print(q2)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.10" } }, "nbformat": 4, "nbformat_minor": 4 }