(* ::Package:: *)

BeginPackage["laporta2`"];

processLaportaFirst::usage="";
processLaportaSecond::usage="";
processLaportaThird::usage="";
writeLaportaSolution::usage="";
readLaportaSolution::usage="";
clearLaportaSolution::usage="";
setLaportaLinkFile::usage="";
unsetLaportaLinkFile::usage="";
getNecessaryId::usage="";
getNecessaryIdandPath::usage="";
getSolvedVariable::usage="";

$loopCount::usage="";
$numberOfSolvedVariable::usage="";
$stageLaporta::usage="";

$laportaCache::usage="";



Begin["`private`"];


<<"jTool`"

Unprotect["laporta2`*"];
Unprotect["laporta2`private`*"];

$laportaCache=False;


solveX[equation_,variable_]:=
	Module[{coefficient,divisor},
		divisor=Coefficient[equation,variable];
		$simplifySol[(-equation/.variable->0)/divisor]
	];

positionList[expr_,pattern_]:=
	Block[{i=1,repRule,tmp},
		repRule=Append[#:>{#,i++}&/@pattern,_:>(i++;Unevaluated[Sequence[]])]//Dispatch;
		tmp=List@@(Replace[expr,repRule,{1}]);
		Union[Flatten[Cases[tmp,{#,a_}:>a]]]&/@pattern
	];

origVarToRemappedVar[in_]:=
	Module[{data,mappedVariable,mapping,mapping2,mappingVariable},
		data=HoldComplete@@in;
		mappedVariable=List@@Intersection[$varMap,data];
		mapping=MapThread[Rule,{mappedVariable,$var[#[[1]]]&/@positionList[$varMap,mappedVariable]}];
		mappingVariable=SortBy[Complement[in,mappedVariable],$sortFunc];
		If[mappingVariable=!={},
			mapping2=MapThread[Rule,{mappingVariable,Array[$var,Length[mappingVariable],Length[$varMap]+1]}];
			$varMap=Join[$varMap,HoldComplete@@mappingVariable];
			$solMap=Join[$solMap,ConstantArray[0,Length[mappingVariable]]];
			mapping=Join[mapping,mapping2]
		];
		mapping
	];

idLess[a_,b_]:=
	Module[{tmp1,tmp2,ml},
		tmp1=getVariable[a,$var][[1;;,1]]//Reverse;
		tmp2=getVariable[b,$var][[1;;,1]]//Reverse;
		ml=Length/@{tmp1,tmp2};
		tmp1=PadRight[tmp1,ml];
		tmp2=PadRight[tmp2,ml];
		OrderedQ[{tmp1,tmp2}]
	];

readData[equation_]:=
	Module[{tmpEquation,tmpVariable,repRule,repRule2},
		$stageLaporta="stage 1";
		$loopCount=0;
		$numberOfSolvedVariable=0;
		tmpEquation=equation/.True->0/.(lhs__==rhs__):>lhs-rhs;
		tmpEquation=$simplifySol/@tmpEquation;
		$stageLaporta="stage 2";
		tmpVariable=SortBy[getVariable[equation,$origVar],$sortFunc];
		If[$sorting,
			repRule2=MapThread[Rule,{tmpVariable,Array[$var,Length[tmpVariable]]}];
			tmpEquation=tmpEquation/.Dispatch[repRule2];
			If[$vlStream=!=None,
				tmpEquation={tmpEquation,Range[Length[tmpEquation]]}//Transpose;
				tmpEquation=SortBy[tmpEquation,idLess];
				Write[$vlStream,{"r",tmpEquation[[1;;,2]]}];
				tmpEquation=tmpEquation[[1;;,1]],
				(*else*)
				tmpEquation=SortBy[tmpEquation,idLess];
			];
			repRule2=Reverse/@repRule2;		
			repRule=origVarToRemappedVar[tmpVariable];
			repRule=repRule2/.Dispatch[repRule],
			(*else*)
			repRule=origVarToRemappedVar[tmpVariable];
			If[$vlStream=!=None,Write[$vlStream,{"r",Range[Length[tmpEquation]]}]];
		];
		tmpEquation=tmpEquation/.Dispatch[repRule];
		Scan[If[MemberQ[Expand[#],
			HoldPattern[$var[_]^_|$var[_]$var[_]],{0,\[Infinity]}],unsetLaportaLinkFile[];Abort[]]&,tmpEquation];
		MapThread[Set,{$equation/@Range[1,Length[tmpEquation]],tmpEquation}];
		Length[tmpEquation]
	];

forwardElimination[equNumber_]:=
	Module[{currentEquation,currentVariable,currentSolution,usedSol,needSolution,repRule},
		$stageLaporta="stage 3";
		Do[
			$loopCount=i;
			currentEquation=$equation[i];
			$equation[i]=.;
			currentVariable=getVariable[currentEquation,$var];
			needSolution={$solMap[[currentVariable[[1;;,1]]]],currentVariable}//Transpose;
			needSolution=DeleteCases[needSolution,{0,_}];
			If[needSolution=!={},
				usedSol=needSolution[[1;;,1]];
				backSubstitution[needSolution[[1;;,2]]];
				If[$compress,
					repRule=#[[2]]->Uncompress[$solution[#[[1]]]]&/@needSolution//Dispatch,
					repRule=#[[2]]->$solution[#[[1]]]&/@needSolution//Dispatch;
				];
				currentEquation=$simplifySol[currentEquation/.repRule];
				currentVariable=getVariable[currentEquation,$var],
				(*else*)
				usedSol={};
			];
			If[!currentVariable==={},
				currentVariable=SortBy[currentVariable,$sortFunc[$origVar@@$varMap[[#[[1]]]]]&];
				currentSolution=solveX[currentEquation,currentVariable[[-1]]];
				$solMap[[currentVariable[[-1,1]]]]=++$solIndex;
				$solution[$solIndex]=If[$compress,Compress[currentSolution],currentSolution];
				If[$vlStream=!=None,Write[$vlStream,{"f",$solIndex,i,usedSol}]];
				$numberOfSolvedVariable=$solIndex,
				If[currentEquation=!=0,Print[currentEquation];clearLaportaSolution[];unsetLaportaLinkFile[]Abort[]];
			];
			If[$laportaCache===False,ClearSystemCache[]];
			If[Divisible[i,2500],Share[]],
			{i,equNumber}
		];
	];

backSubstitution[variable_]:=
	Module[{processingVariable,processedVariable,nextVariable,needSolution,
					processingSolution,currentSolution,currentVariable,repRule,j},
		processingVariable=variable[[1;;,1]];
		processedVariable={};
		processingSolution=Reap[
			While[processingVariable=!={},
				processedVariable=Join[processedVariable,processingVariable];
				needSolution=$solMap[[processingVariable]];
				needSolution=DeleteCases[needSolution,0];
				Sow[needSolution];
				nextVariable=Reap[
					If[$compress,
						Scan[Sow[getVariable[Uncompress[$solution[#]],$var][[1;;,1]]]&,needSolution],
						Scan[Sow[getVariable[$solution[#],$var][[1;;,1]]]&,needSolution];
					];
				][[2]]//Flatten//Union;
				processingVariable=Complement[nextVariable,processedVariable];
			];
		][[2]]//Flatten//Union//Sort//Reverse;
		j=1;
		Scan[
			Function[i,
				currentSolution=If[$compress,Uncompress[$solution[i]],$solution[i]];
				currentVariable=getVariable[currentSolution,$var];
				needSolution={$solMap[[currentVariable[[1;;,1]]]],currentVariable}//Transpose;
				needSolution=DeleteCases[needSolution,{0,_}];
				If[needSolution=!={},
					If[$vlStream=!=None,Write[$vlStream,{"b",i,needSolution[[1;;,1]]}]];
					If[$compress,
						repRule=#[[2]]->Uncompress[$solution[#[[1]]]]&/@needSolution//Dispatch,
						repRule=#[[2]]->$solution[#[[1]]]&/@needSolution//Dispatch;
					];
					currentSolution=$simplifySol[currentSolution/.repRule];
					$solution[i]=If[$compress,Compress[currentSolution],currentSolution];
				];
				If[$laportaCache===False,ClearSystemCache[]];
				If[Divisible[j++,2500],Share[]];
			],
			processingSolution
		];
	];


Options[processLaportaFirst]=
	{"sortFunc"->priorityOrder2,"compress"->False,"numerical"->False,"delta"->0,"equationSorting"->True};
processLaportaFirst[equation_List,variable_,OptionsPattern[]]:=
	Module[{delta,compress,numerical,equLength,filename},
		$origVar=variable;
		$numerical=(OptionValue["numerical"]===True);
		$delta=OptionValue["delta"];
		$compress=(OptionValue["compress"]===True);
		$sortFunc=OptionValue["sortFunc"];
		$solIndex=0;
		$sorting=(OptionValue["equationSorting"]===True);
		clearLaportaSolution[];
		If[$vlStream=!=None,
			filename=Close[$vlStream];
			DeleteFile[filename];
			$vlStream=OpenWrite[filename];
		];
		processLaportaSecond[equation];
	];

processLaportaSecond[equation_List]:=
	Module[{equLength},
		If[$numerical&&$delta===0,$simplifySol=Expand];
		If[$numerical&&$delta=!=0,$simplifySol=Chop[Expand[#],$delta]&];
		If[!$numerical&&$delta===0,$simplifySol=Collect[#,$origVar[__],Together]&];
		If[!$numerical&&$delta=!=0,$simplifySol=Chop[Collect[#,$origVar[__],Together],$delta]&];
		equLength=readData[equation];
		If[!$numerical&&$delta===0,$simplifySol=Collect[#,$var[__],Together]&];
		If[!$numerical&&$delta=!=0,$simplifySol=Chop[Collect[#,$var[__],Together],$delta]&];
		forwardElimination[equLength];
	];

processLaportaThird[variable_List]:=
	Module[{currentVariable,needSolution,result,repRule},
		$stageLaporta="stage 4";
		currentVariable=origVarToRemappedVar[variable][[1;;,2]];
		needSolution={$solMap[[currentVariable[[1;;,1]]]],currentVariable}//Transpose;
		needSolution=DeleteCases[needSolution,{0,_}];
		If[needSolution==={},$stageLaporta="Done.";Return[{}]];
		backSubstitution[needSolution[[1;;,2]]];
		If[$compress,
			result={#[[2]],Uncompress[$solution[#[[1]]]]}&/@needSolution,
			result={#[[2]],$solution[#[[1]]]}&/@needSolution;
		];
		currentVariable=getVariable[result,$var];
		repRule=Thread[currentVariable->List@@$varMap[[currentVariable[[1;;,1]]]]];
		result=result/.Dispatch[repRule];
		$stageLaporta="Done.";
		result
	];


clearLaportaSolution[]:=(
	$solIndex=0;
	$varMap=HoldComplete[];
	$solMap={};
	Clear[$equation,$solution];
);
clearLaportaSolution[];

clearLaportaSolution[var_List]:=
	Module[{varPos},
		varPos=positionList[$varMap,var]//Flatten;
		Quiet[Scan[($solution[#]=.)&,$solMap[[varPos]]]];
		$solMap[[varPos]]=0;
	];

writeLaportaSolution[filename_]:=
	Module[{stream},
		stream=OpenWrite[filename];
		If[stream===$Failed,Return[$Failed]];
		Write[stream,#]&/@{$origVar,$numerical,$delta,$compress,$sortFunc,$solIndex};
		Write[stream,$varMap];
		Write[stream,$solMap];
		Do[
			Write[stream,$solution[i]],
			{i,$solIndex}
		];
		Close[stream]
	];

readLaportaSolution[filename_]:=
	Module[{stream},
		stream=OpenRead[filename];
		If[stream===$Failed,Return[$Failed]];
		{$origVar,$numerical,$delta,$compress,$sortFunc,$solIndex}=Table[Read[stream],{6}];
		$varMap=Read[stream];
		$solMap=Read[stream];
		Clear[$solution];
		Do[
			$solution[i]=Read[stream],
			{i,$solIndex}
		];
		If[!$numerical&&$delta===0,$simplifySol=Collect[#,$var[__],Together]&];
		If[!$numerical&&$delta=!=0,$simplifySol=Chop[Collect[#,$var[__],Together],$delta]&];
		Close[stream]
	];

If[$vlStream=!=None,Close[$vlStream]//Quiet];
$vlStream=None;

setLaportaLinkFile[filename_]:=
	Module[{tmp},
		tmp=OpenAppend[filename];
		If[tmp===$Failed,Return[$Failed]];
		$vlStream=tmp;
		filename
	];

unsetLaportaLinkFile[]:=
	Module[{tmp},
		If[$vlStream===None,Return[$Failed]];
		tmp=Close[$vlStream];
		$vlStream=None;
		tmp
	];


take2[f_,sol_]:=
	Module[{tmp,old,currentSol,processingSol},
		processingSol={{sol,\[Infinity]}};
		While[processingSol=!={},
			currentSol=processingSol[[1]];
			processingSol=Rest[processingSol];
			tmp=f[currentSol[[1]]];
			Sow[currentSol[[1]],1];
			old=Select[tmp,#[[1]]<currentSol[[2]]&];
			f[currentSol[[1]]]=Complement[f[currentSol[[1]]],old];
			Scan[
				Function[v,
					If[v[[2]]==="f",
						Sow[v[[{1,4}]],2];
						Sow[v[[5]],1],
						(*else*)
						Sow[v[[4]],1];
					];
				],
				old
			];
			Scan[
				Function[v,
					processingSol=Join[{#,v[[1]]}&/@v[[-1]],processingSol];
				],
				old
			];
		];
	];

takeIdList[stream_,sol_]:=
	Module[{prevPos,sline,f,count,tmp,curId,newId,curSol,newSol,idMap,result},
		prevPos=0;
		sline=Reap[
			While[
				tmp=Read[stream];
				tmp=!=EndOfFile,
				If[tmp[[1]]==="r",Sow[prevPos]];
				prevPos=StreamPosition[stream];
			];
		][[2]]//Flatten//Union;
		curSol=sol;
		result={};
		Do[
			SetStreamPosition[stream,sline[[i]]];
			idMap=Read[stream][[2]];
			count=0;
			Clear[f];
			f[_]={};
			While[
				tmp=Read[stream];
				(tmp=!=EndOfFile)&&(tmp[[1]]=!="r"),
				AppendTo[f[tmp[[2]]],{++count,Sequence@@tmp}];
			];
			{newSol,newId}=Reap[take2[f,#]&/@curSol,{1,2}][[2]];
			curSol=newSol//Flatten//Union;
			If[newId=!={},
				newId=SortBy[newId[[1]],First];
				newId=newId[[1;;,2]]//DeleteDuplicates;
				PrependTo[result,idMap[[newId]]],
				(*else*)
				PrependTo[result,{}];
			],
			{i,Length[sline],1,-1}
		];
		{result,curSol}
	];

getNecessaryId[variable_,linkfile_,solutionfile_]:=getNecessaryIdandPath[variable,linkfile,solutionfile][[1]];

getNecessaryIdandPath[variable_,linkfile_,solutionfile_]:=
	Module[{streamL,streamS,varMap,solMap,mappedVariable,repRule,var,sol,id,res},
		streamL=OpenRead[linkfile];
		streamS=OpenRead[solutionfile];
		If[streamL===$Failed||streamS===$Failed,Return[$Failed]];
		Do[Read[streamS],{6}];
		varMap=List@@Read[streamS];
(*
HoldComplete\:306e\:307e\:307e\:4e00\:5ea6\:3082\:8a55\:4fa1\:3092\:3057\:306a\:3044\:5834\:5408\:3001\:4f8b\:3048\:3070 7/2 \:304c Times[7,Power[2,-1]] \:3068\:8868\:73fe\:3055\:308c\:308b\:3002
\:4e00\:65b9\:3001\:4e00\:5ea6\:3067\:3082\:8a55\:4fa1\:3059\:308b\:3068 Ratioal[7,2] \:3068\:8868\:3055\:308c\:308b\:3002 "variable" \:306b\:5206\:6570\:304c\:542b\:307e\:308c\:308b\:5834\:5408\:3001varMap\:3092
\:6700\:4f4e\:4e00\:5ea6\:8a55\:4fa1\:3059\:308b\:5fc5\:8981\:304c\:3042\:308b\:3002
*)
		mappedVariable=List@@Intersection[varMap,variable];
		var=positionList[varMap,mappedVariable]//Flatten//Union;
		solMap=Read[streamS];
		sol=DeleteCases[solMap[[var]],0];
		res=takeIdList[streamL,sol];
		res[[2]]=List@@varMap[[Flatten[positionList[solMap,res[[2]]]]]];
		Close[streamS];
		Close[streamL];
		res
	];

getSolvedVariable[]:=
	Module[{solvedVar},
		solvedVar=Position[$solMap,Except[0]]//Flatten//Rest;
		List@@$varMap[[solvedVar]]
	];


Protect["laporta2`*"];
Protect["laporta2`private`*"];
Unprotect["laporta2`$*"];
Unprotect["laporta2`private`$*"];


End[];

EndPackage[];
