(*nP,nR,nS,nPR,nRS:dimensionsoftheP,R,S,PR,RSspaces,eMPS,eMP,eMR,eMS,eM:tensor,matricesandvectordefiningequalityconstraints,iMPS,iMP,iMR,iMS,iM:tensor,matricesandvectordefininginequalityconstraints,Y0:astatesatisfyingconstraints*){{nP,nR,nS,nPR,nRS},{eMPS,eMP,eMR,eMS,eM},{iMPS,iMP,iMR,iMS,iM}}=Import["input.mx"];Y0=Import["Y0.mx"];
qRS=SparseArray[IdentityMatrix[nRS]];qPR=SparseArray[IdentityMatrix[nPR]];projectRS[YY_]:=Module[{walP,walR,walS,iMaRS,iVaRS,eMaRS,eVaRS,valRS,solRS},(*projectionontheRSspace*){walP,walR,walS}=YY;iMaRS=Join[iMR,iMS+walP.iMPS,2];iVaRS=iM+iMP.walP;eMaRS=Join[eMR,eMS+walP.eMPS,2];eVaRS=eM+eMP.walP;valRS=-Join[walR,walS]RandomReal[{0.8,1.25},nRS];solRS=Check[QuadraticOptimization[{qRS,valRS},{iMaRS,iVaRS},{eMaRS,eVaRS}],Return[{}]];walR=Take[solRS,nR];walS=Drop[solRS,nR];{walP,walR,walS}];projectPR[YY_]:=Module[{walP,walR,walS,iMaPR,iVaPR,eMaPR,eVaPR,valPR,solPR},(*projectiononthePRspace*){walP,walR,walS}=YY;iMaPR=Join[iMP+Transpose[iMPS.walS],iMR,2];iVaPR=iM+iMS.walS;eMaPR=Join[eMP+Transpose[eMPS.walS],eMR,2];eVaPR=eM+eMS.walS;valPR=-Join[walP,walR]RandomReal[{0.8,1.25},nPR];solPR=Check[QuadraticOptimization[{qPR,valPR},{iMaPR,iVaPR},{eMaPR,eVaPR},Tolerance->10^(-3)],Return[{}]];walP=Take[solPR,nP];walR=Drop[solPR,nP];{walP,walR,walS}];(*pastT:previouslyfoundsolutions*)pastT={Y0};add[elem_]:=Module[{l},l=Length[pastT];pastT=If[l<20,AppendTo[pastT,elem],ReplacePart[pastT,RandomInteger[{1,l}]->elem]]];(*program:length[]>0meansthatasolutionhasbeenfound;seethecodeofprojectionmodules*)step[YT_]:=Module[{ZT,XT},Print[nit++];XT=projectPR[YT];ZT=If[Length[XT]>0,projectRS[XT],{}];If[Length[ZT]>0,add[ZT];ZT,RandomChoice[pastT]]];nit=0;Nest[step,Y0,200]