# -*- coding: utf-8 -*-# noinspection PyUnresolvedReferencesr""".. _ap-mp:============Matrix proxy============With the fully discrete weak formulation ``wf``, we can bring it into its algebraic proxy by calling its method``mp``, standing for *matrix proxy*,>>> mp = wf.mp()which is an instance of :class:`MatrixProxy`, .. autoclass:: MatrixProxy :members:Similarly, its ``pr`` method can illustrate it properly,>>> mp.pr()<Figure size ..... _ap-ap:========================Algebraic representation========================Depend on ``mp`` is linear or nonlinear, an algebraic system can be producedthrough either method ``ls`` or ``nls`` of ``mp``,see :meth:`MatrixProxy.ls` and :meth:`MatrixProxy.nls`.Method ``ls`` gives an instance of :class:`MatrixProxyLinearSystem`, i.e., .. autoclass:: MatrixProxyLinearSystem :members:And method ``nls`` leads to an instance of :class:`MatrixProxyNoneLinearSystem`, namely, .. autoclass:: MatrixProxyNoneLinearSystem :members:In this case, ``mp`` is a linear system. Thus, we should call ``ls`` method of it,>>> ls = mp.ls()>>> ls.pr()<Figure size ...Eventually, a fully discrete abstract linear system is obtained. We can send it a particular implementationwhich will *objectivize* it, for example by making matrices 2-dimensional arrays and making the vectors1-dimensional arrays. These implementations will be introduced in the following section."""importmatplotlib.pyplotaspltimportmatplotlibplt.rcParams.update({"text.usetex":True,"font.family":"DejaVu Sans","text.latex.preamble":r"\usepackage{amsmath, amssymb}",})matplotlib.use('TkAgg')fromsrc.form.othersimport_find_formfromsrc.configimport_root_form_ap_vec_settingfromtools.frozenimportFrozenfromsrc.algebra.arrayimportAbstractArrayfromsrc.wf.term.apimportTermLinearAlgebraicProxyfromsrc.wf.term.apimportTermNonLinearOperatorAlgebraicProxyfromsrc.algebra.linear_systemimportBlockMatrix,BlockColVector,LinearSystemfromsrc.algebra.nonlinear_systemimportNonLinearSystemfromsrc.wf.mp.linear_systemimportMatrixProxyLinearSystemfromsrc.wf.mp.nonlinear_systemimportMatrixProxyNoneLinearSystem
[docs]classMatrixProxy(Frozen):""""""def__init__(self,wf):""""""ap=wf.ap()# make an algebraic proxy in real time.ifap._fully_resolved:assertap.linearityin('linear','nonlinear'), \
f"ap linearity must be linear or nonlinear when it is fully resolved."else:raiseException(f"there is(are) term(s) in the wf not-resolved as 'algebraic proxy', check 'ap' of "f"the 'wf' to see which term(s) is(are) not resolved!.")self._wf=wfself._ap=apself._num_equations=len(ap._term_dict)self._unknowns=ap.unknownsself._test_vectors=ap.test_vectorsself.___total_indexing_length___=Noneself._lbv=BlockColVector(self._num_equations)# left block vector partself._rbv=BlockColVector(self._num_equations)# right block vector partself._left_nonlinear_terms=[[]for_inrange(self._num_equations)]self._left_nonlinear_signs=[[]for_inrange(self._num_equations)]self._right_nonlinear_terms=[[]for_inrange(self._num_equations)]self._right_nonlinear_signs=[[]for_inrange(self._num_equations)]self._linearity='linear'foriinap._term_dict:terms=ap._term_dict[i]signs=ap._sign_dict[i]forjinrange(2):lr_terms=terms[j]lr_signs=signs[j]k=0forterm,signinzip(lr_terms,lr_signs):linearity,term=self._test_vector_remover(i,term)iflinearity=='linear':assertself._ap._linearity_dict[i][j][k]=='linear',f'Must be this case.'ifj==0:self._lbv._add(i,term,sign)else:self._rbv._add(i,term,sign)eliflinearity=='nonlinear':assertself._ap._linearity_dict[i][j][k]=='nonlinear',f'Must be this case.'ifself._linearity=='linear':self._linearity='nonlinear'else:passifj==0:# noinspection PyTypeCheckerself._left_nonlinear_terms[i].append(term)self._left_nonlinear_signs[i].append(sign)else:# noinspection PyTypeCheckerself._right_nonlinear_terms[i].append(term)self._right_nonlinear_signs[i].append(sign)else:raiseNotImplementedError()k+=1self._l_mvs=list()# left matrix@vector sectionsself._r_mvs=list()# right matrix@vector sectionsself.parse(self._unknowns)self._bc=wf._bcself._freeze()@propertydeflinearity(self):""""""assertself._linearityin('linear','nonlinear'),f"Linearity must be among ('linear', 'nonlinear')."returnself._linearitydef_test_vector_remover(self,i,term):""""""test_vector=self._test_vectors[i]ifterm.__class__isTermLinearAlgebraicProxy:aa=term._abstract_arrayfactor=aa._factorcomponents=aa._componentstransposes=aa._transposesassertcomponents[0]==test_vector,f"cannot remove test vector {test_vector} from term {term}."asserttransposes[0]isTrue,f"cannot remove test vector {test_vector} from term {term}."new_aa=AbstractArray(factor=factor,components=components[1:],transposes=transposes[1:],)return'linear',new_aaelifterm.__class__isTermNonLinearOperatorAlgebraicProxy:tf_pure_lin_repr=test_vector._pure_lin_repr[:-len(_root_form_ap_vec_setting['lin'])]tf=_find_form(tf_pure_lin_repr)ifterm._tfisNone:term.set_test_form(tf)else:assertterm._tfistf,f"double check the test form."return'nonlinear',termelse:raiseNotImplementedError(term.__class__)defparse(self,targets):""""""targets=list(targets)lbv=self._lbvrbv=self._rbvbb=BlockColVector(self._num_equations)fori,targetinenumerate(targets):iftarget.__class__.__name__=='Form':target=target.ap()targets[i]=targetbb._add(i,target,'+')forlor,bvinenumerate((lbv,rbv)):bm=BlockMatrix((self._num_equations,len(targets)))remaining_bv=BlockColVector(self._num_equations)fori,entryinenumerate(bv._entries):forj,terminenumerate(entry):sign=bv._signs[i][j]ifterm.__class__isAbstractArray:t=term._components[-1]trans=term._transposes[-1]iftintargetsandnottrans:# found a correct term to be put int matrix blockk=targets.index(t)components=term._components[:-1]transposes=term._transposes[:-1]factor=term._factorbm_aa=AbstractArray(factor=factor,components=components,transposes=transposes,)bm._add(i,k,bm_aa,sign)else:remaining_bv._add(i,term,sign)else:raiseNotImplementedError()iflor==0:ifbm._is_empty():passelse:self._l_mvs.append((bm,bb))self._lbv=remaining_bvelse:ifbm._is_empty():passelse:self._r_mvs.append((bm,bb))self._rbv=remaining_bvself.___total_indexing_length___=Nonedef_total_indexing_length(self):""""""ifself.___total_indexing_length___isNone:a=len(self._l_mvs)c=len(self._r_mvs)ifself._lbv._is_empty():b=0else:b=1ifself._rbv._is_empty():d=0else:d=1self.___total_indexing_length___=(a,b,c,d),a+b+c+dreturnself.___total_indexing_length___def__getitem__(self,index):"""To retrieve a linear term: do 'a-b' or 'a-b,c', where a, b, c are str(integer)-s. when index = 'a-b' `a` refer to the `ath` block. `b` refer to the `b`th entry of the ColVec of the block. when index = 'a-b,c' `a` refer to the `ath` block. `b,c` refer to the `b,c`th entry of the BlockMatrix of the block. So when `a`th block is a Matrix @ Vector, we can use either 'a-b' (vector block) or 'a-b,c' (matrix block). But when `a`th block is a ColVec, we can only use 'a-b'. To retrieve a nonlinear term, to be continued. """assertisinstance(index,str),f"pls put index in string."assertindex.count('-')==1,f"linear term index={index} is illegal."block_num,local_index=index.split('-')assertblock_num.isnumeric(),f"linear term index={index} is illegal."block_num=int(block_num)abcd,total_length=self._total_indexing_length()a,b,c,d=abcdassert0<=block_num<total_length,f"linear term index={index} is illegal; it beyonds the length."ifblock_num<a:# retrieving a term in left l_mvsblock=self._l_mvs[block_num]elifblock_num<a+b:# retrieving a term in left remaining vectorblock=self._lbvelifblock_num<a+b+c:# retrieving a term in right l_mvsblock_num-=a+bblock=self._r_mvs[block_num]else:# retrieving a term in right remaining vectorblock=self._rbvindices=eval('['+local_index+']')ifisinstance(block,tuple):iflen(indices)==2:block=block[0]eliflen(indices)==1:block=block[1]else:raiseException(f"linear term index={index} is illegal.")else:passtry:returnblock(*indices)except(IndexError,TypeError):raiseException(f"linear term index={index} is illegal.")def_pr_text(self):symbolic=''plus=''variant=0forbm_bbinself._l_mvs:bm,bb=bm_bbassertnotbm._is_empty()ifvariant==0:_v_plus=''else:_v_plus='+'symbolic+=_v_plussymbolic+=bm._pr_text()symbolic+=bb._pr_text()variant+=1plus='+'ifself._lbv._is_empty():passelse:symbolic+=plus+self._lbv._pr_text()nonlinear_text=''ifself.linearity=='nonlinear':nonlinear_terms=self._left_nonlinear_termsnonlinear_signs=self._left_nonlinear_signsnum_terms=0fortermsinnonlinear_terms:num_terms+=len(terms)ifnum_terms==0:passelse:# there are nonlinear terms on the left hand sideforterms,signsinzip(nonlinear_terms,nonlinear_signs):iflen(terms)==0:nonlinear_text+=r'\boldsymbol{0}'else:fori,terminenumerate(terms):sign=signs[i]ifi==0:ifsign=='-':passelse:sign=''else:passnonlinear_text+=sign+term._sym_reprnonlinear_text+=r'\\'nonlinear_text=nonlinear_text[:-len(r'\\')]nonlinear_text=r" + \begin{bmatrix}"+nonlinear_text+r"\end{bmatrix}"symbolic+=nonlinear_text+'='plus=''variant=0forbm_bbinself._r_mvs:bm,bb=bm_bbassertnotbm._is_empty()ifvariant==0:_v_plus=''else:_v_plus='+'symbolic+=_v_plussymbolic+=bm._pr_text()symbolic+=bb._pr_text()variant+=1plus='+'ifself._rbv._is_empty():passelse:symbolic+=plus+self._rbv._pr_text()nonlinear_text=''ifself.linearity=='nonlinear':nonlinear_terms=self._right_nonlinear_termsnonlinear_signs=self._right_nonlinear_signsnum_terms=0fortermsinnonlinear_terms:num_terms+=len(terms)ifnum_terms==0:passelse:# there are nonlinear terms on the left hand sideforterms,signsinzip(nonlinear_terms,nonlinear_signs):iflen(terms)==0:nonlinear_text+=r'\boldsymbol{0}'else:fori,terminenumerate(terms):sign=signs[i]ifi==0:ifsign=='-':passelse:sign=''else:passnonlinear_text+=sign+term._sym_reprnonlinear_text+=r'\\'nonlinear_text=nonlinear_text[:-len(r'\\')]nonlinear_text=r" + \begin{bmatrix}"+nonlinear_text+r"\end{bmatrix}"symbolic+=nonlinear_textreturnsymbolicdef_mp_seek_text(self):"""seek text"""seek_text=self._wf._mesh.manifold._manifold_text()seek_text+=r'seek $\left('form_sr_list=list()space_sr_list=list()foruninself._ap.unknowns:form_sr_list.append(rf' {un._sym_repr}')space_sr_list.append(rf"{un._shape_text()}")seek_text+=','.join(form_sr_list)seek_text+=r'\right) \in 'seek_text+=r'\times '.join(space_sr_list)seek_text+='$, such that\n'returnseek_text
[docs]defpr(self,figsize=(12,8)):"""Print the representation, a figure, of this weak formulation. Parameters ---------- figsize : tuple, optional The figure size. It has no effect when the figure is over-sized. A tight configuration will be applied when it is the case. The default value is ``(12, 8)``. """fromsrc.configimportRANK,MASTER_RANKifRANK!=MASTER_RANK:returnelse:passseek_text=self._mp_seek_text()symbolic=r"$"+self._pr_text()+r"$"ifself._bcisNoneorlen(self._bc)==0:bc_text=''else:bc_text=self._bc._bc_text()fig=plt.figure(figsize=figsize)plt.axis((0,1,0,1))plt.axis('off')plt.text(0.05,0.5,seek_text+symbolic+bc_text,ha='left',va='center',size=15)fromsrc.configimport_setting,_pr_cacheif_setting['pr_cache']:_pr_cache(fig,filename='matrixProxy')else:plt.tight_layout()plt.show(block=_setting['block'])returnfig
def_parse_ls(self):""""""assertself._lbv._is_empty(),f"Format is illegal, must be like Ax=b, do '.pr()' to check!"assertlen(self._l_mvs)==1,f"Format is illegal, must be like Ax=b, do '.pr()' to check!"A,x=self._l_mvs[0]b=BlockColVector(self._rbv._shape)forMvinself._r_mvs:M,v=Mvforiinrange(M._shape[0]):forjinrange(M._shape[1]):vj,sj=v(j)Mij,Sij=M(i,j)assertlen(vj)==len(sj)==1,f"Format is illegal, must be like Ax=b, do pr() to check!"vj,sj=vj[0],sj[0]formij,sijinzip(Mij,Sij):ifsj==sij:sign='+'else:sign='-'mij_at_vj=mij@vjb._add(i,mij_at_vj,sign)foriinrange(self._rbv._shape):bi,si=self._rbv(i)forbj,sjinzip(bi,si):b._add(i,bj,sj)ls=LinearSystem(A,x,b)returnls
[docs]defls(self):"""Convert self to an abstract linear system. Returns ------- ls : :class:`MatrixProxyLinearSystem` The linear system instance. """assertself.linearity=='linear',f"'mp' is not linear, try to use '.nls'."ls=self._parse_ls()returnMatrixProxyLinearSystem(self,ls,self.bc)
[docs]defnls(self):"""Convert self to an abstract nonlinear system. Returns ------- nls : :class:`MatrixProxyNoneLinearSystem` The nonlinear system instance. """assertself.linearity=='nonlinear',f"'mp' is linear, just use '.ls'."ls=self._parse_ls()mp_ls=MatrixProxyLinearSystem(self,ls,self.bc)len_right_nonlinear_terms=0fortermsinself._right_nonlinear_terms:len_right_nonlinear_terms+=len(terms)assertlen_right_nonlinear_terms==0, \
f"To initialize a nonlinear system, move all nonlinear terms to the left hand side first."nls=NonLinearSystem(ls,(self._left_nonlinear_terms,self._left_nonlinear_signs),)returnMatrixProxyNoneLinearSystem(self,mp_ls,nls)