Add LoadPath "./interval".


Require Import Fourier Rdefinitions Raxioms RIneq Rbasic_fun.
Require Import  ssreflect ssrbool ssrnat eqtype ssrfun seq fintype finfun perm matrix
   bigops ssralg poly charpoly.
Require Import Rstruct min_max topology interval_R real_matrix.

(*
abstract the notion of eigenvalues and spectral radius,
admit results: Perron-Frobenius, properties of Rayleigh quotients
*)

Set Implicit Arguments.
Unset Strict Implicit.
Import Prenex Implicits.

Open Scope R_scope.
Delimit Scope R_scope with Re. 


Open Scope ring_scope.
Delimit Scope ring_scope with Ri.

Import GRing.Theory.

Section eigenDef.

(*need to talk about eigenvalues, spectral radius *)
Variable n': nat.
Notation n := n'.+1.

Variable T: comUnitRingType.
   
Variable A: matrix T n n.

Definition eigenv := root (char_poly A).
(* roots defined for unitRing, char_poly for comRing *)

End eigenDef.

Section detProp.

Variable n: nat.
(*Notation n := n'.+1.*)

Variable T: fieldType.
   

Lemma det0P_cv_n :
forall A: 'M[T]_n,  reflect (forall x: 'cV_n, A *m x = 0 -> x = 0)  (\det A != 0).
Proof.
move=> A.
have H0:= @det0P T n (trmx A).
case C: (\det A != 0); [apply: ReflectT |apply: ReflectF].
move=> x Hx.
replace  x  with ( (invmx A) *m (A *m x)).
by rewrite Hx mulmx0.
rewrite mulmxA mulVmx.
by rewrite mul1mx.
rewrite /unitmx -topredE /=.
rewrite unitfE.
apply: C.
move/negP: C.
move/negP => C.
rewrite -det_tr in C.
elim: (elimT H0 C) => x Hx1 Hx2 Hc.
have Hcx:= Hc (x^T).
have Htx: A *m x^T = 0.
by rewrite -(trmxK A) -trmx_mul Hx2 trmx_const /GRing.zero /=.
move/eqP: Hx1 => Hx1.
elim Hx1.
rewrite /GRing.zero /= -trmx_const -(trmxK x) (Hcx Htx); f_equal.
Qed.

Lemma det0P_cv:
forall A: 'M[T]_n, reflect (exists2 v : 'cV_n, v != 0 & A *m v = 0) (\det A == 0).
Proof.
move=> A.
have H0:= @det0P T n (trmx A).
case C: (\det A == 0); [apply: ReflectT |apply: ReflectF].
rewrite -det_tr in C.
elim: (elimT H0 C) => x Hx1 Hx2 .
exists (x^T).
rewrite /GRing.zero /=  -trmx_const .
apply/eqP => Htx.
move/eqP: Hx1=> Hx1; elim: Hx1.
apply/matrixP => i j; rewrite mxE; move/matrixP: Htx => Htx.
by have:= Htx j i; rewrite /trmx !mxE.
by rewrite -(trmxK A) -trmx_mul Hx2 trmx_const /GRing.zero /=.
rewrite -det_tr in C.
move/eqP: C => Hdet  [x [Hex1 Hex2]].
have H: (exists2 v : 'rV_n, v != 0 & v *m A^T = 0).
exists (x^T).
apply/eqP => Htx.
move/eqP: Hex1=> Hx1; elim: Hx1.
 apply/matrixP => i j; rewrite mxE; move/matrixP: Htx => Htx.
by have:= Htx j i; rewrite /trmx !mxE.
by rewrite -trmx_mul Hex2 trmx_const /GRing.zero /=.
by elim: Hdet; apply/eqP; apply/H0.
Qed.

End detProp.


Section eigenProp.

Variable n': nat.
Notation n := n'.+1.

Variable T: fieldType.
   

Lemma lam_equ:
forall (A: 'M[T]_n) lam, lam \in eigenv A -> exists2 x: 'cV_n, x != 0 & A*m x = lam *m: x.
Proof.
move=> A lam Hlam.
rewrite /eigenv /root /char_poly in Hlam.
have H:= det0P_cv (lam%:M - A).
have Heq: (\det ('X%:M - matrixC A)).[lam] = \det (lam%:M - A).
rewrite /determinant horner_sum.
apply eq_bigr => s _ .
rewrite horner_mul .
congr (_ * _ ).
case (odd_perm s); rewrite /=.
by rewrite !expr1 horner_opp; f_equal; rewrite hornerC.
by rewrite !expr0 hornerC.
elim: index_enum.
by rewrite !big_nil hornerC.
move=> i t IH.
rewrite !big_cons -IH horner_mul.
congr (_ * _).
rewrite !mxE  horner_add horner_opp hornerC.
congr (_ + _).
by case (i == s i); rewrite /=; [rewrite !mulr1n hornerX| rewrite !mulr0n hornerC].
move: Hlam; rewrite -topredE /= => Hlam.
rewrite Heq in Hlam.
have Hr: (exists2 v : 'cV_n, v != 0 & (lam%:M - A) *m v = 0) by apply/H.
elim Hr => x [Hx1 Hx2].
exists x; first by done.
rewrite mulmx_subl  -scalemx1 -scalemxAl mul1mx in Hx2.
by move/eqP: Hx2; rewrite subr_eq0; move/eqP=> ->.
Qed.

End eigenProp.

Section eigenR.

Variable n': nat.
Notation n := n'.+1.

Variable eigens: matrix R n n -> seq R.
Hypothesis Heig: forall A lam, lam \in (eigenv A) <-> (lam \in (eigens A)).
(* all elements in eigens A, satisfy the predicate, i.e. are eigenvalues of A *)





(* spectral radius: we can define it with big[Rmax] on eigens *) 

Definition specr A := \big[Rmax/Rabs (last 0 (eigens A))]_(lam <- (eigens A)) Rabs lam.

Lemma lam_le_specr:
  forall A lam, lam \in eigens A -> Rabs lam <= specr A.
Proof.
move=> A lam Hlam.
rewrite /specr.
move: Hlam; elim (eigens A).
rewrite in_nil //.
move=> t s IH Hlam.
rewrite big_cons.
rewrite in_cons in Hlam.
move/orP: Hlam=> Hlam; elim Hlam; move/eqP => Hlam1.
rewrite Hlam1; apply RmaxLess1.
apply Rle_trans with ((\big[Rmax/Rabs (last 0 (t :: s))]_(j <- s) Rabs j)).
have Hs: s != [::].
apply/negP; move/eqP => H.
have H1 := in_nil lam.
move/eqP: Hlam1 => Hlam1.
rewrite H H1 in Hlam1.
discriminate Hlam1.
rewrite -(last_cons t Hs).
apply: IH; apply/eqP ; apply Hlam1.
apply RmaxLess2.
Qed.


(*Lemma lam_le_specr_val:
  forall A lam, lam \in eigenv A -> Rabs lam <= specr A.
Proof.
move=> A lam Hlam.
have H:= @lam_le_specr A lam.
have He:= Heig A.
move/allP : He => He.
have Hq := He lam.
apply H.



Hypothesis def_specr:
  forall A lam, eigenv A lam -> Rabs lam <= specr A.
*)


Definition lam_min A := \big[Rmin/ (last 0 (eigens A))]_(i <- eigens A) i.
Definition lam_max A :=  \big[Rmax/(last 0 (eigens A))]_(i <- eigens A) i .



Variable mnorm: can_norm.


Lemma lam_equ_norm:
forall (A: 'M[R]_n) (lam: R), eigenv A lam -> 
  exists2 x: 'cV_n,  || x: mnorm | = 1 & A *m x = lam *m: x .
Proof.
move=> A lam Hlam.
have Hl:= @lam_equ n' real_field A lam Hlam.
elim Hl => x [Hx1 Hx2].
have Hdifx: ||x : mnorm| != 0.
by apply/eqP; apply Rgt_not_eq; apply norm_pos_st.
have Hunit: GRing.unit (|| x : mnorm | : real_field).
by rewrite unitfE.
exists (1 / || x : mnorm | *m: x).
rewrite hom_pos_norm Rabs_right.
by rewrite -mulrA  mulVr; first by rewrite mul1r.
apply Rle_ge; apply Rmult_pos; first by toR; fourier.
by left; toR; rewrite /Rinvx Hdifx; apply Rinv_0_lt_compat; apply norm_pos_st.
by rewrite scalemxA (mulrC lam _) -(scalemxA _ lam) -Hx2  -scalemxAr.
Qed.


Lemma eigenv_le_norm:
forall (A: 'M[R]_n) (lam: R), eigenv A lam ->
Rabs lam <= || A: mnorm|.
Proof.
move=> A lam Hlam.
elim ( lam_equ_norm Hlam) => x [Hxnorm Hxeq].
rewrite -(mulr1 (Rabs lam)) -Hxnorm -hom_pos_norm -Hxeq.
apply Rle_trans  with (|| A : mnorm | * || x : mnorm |); 
  first by apply prod_ineq_norm.
by rewrite Hxnorm mulr1; right.
Qed.

Corollary specr_le_norm:
forall (A: 'M[R]_n), specr A <= || A : mnorm |.
Proof.
move=> A.
rewrite /specr.
case C: (eigens A != [::]).
have H1 : (forall i : real_eqType, i \in eigens A -> 0 <= Rabs i) .
move=>*; apply Rabs_pos.
elim: ( @max_exists_seq2 _ Rabs C) => k [Hkex Hkeq].
by rewrite Hkeq; apply eigenv_le_norm; rewrite Heig.
move/eqP: C => ->.
 rewrite big_nil /last Rabs_R0.
 apply pos_norm.
Qed.


Lemma pos_def_eigen_pos:
forall (A: 'M[R]_n) (lam: R), eigenv A lam -> 
 pos_def_gen A -> 0 < lam.
Proof.
move=> A lam Hlam HA.
rewrite /pos_def_gen in HA.
elim: (@lam_equ n' real_field A lam Hlam) => x Hx0 Heq.
have HAx:= HA x Hx0.
rewrite -mulmxA Heq -scalemxAr in HAx.
have HAx0 := HAx ord0 ord0.
move/vec_prod_pos: Hx0  => Hx0.
have Hxtx:= Hx0 ord0 ord0.
rewrite mxE mxE in HAx0, Hxtx.
apply Rmult_lt_reg_l with  ((x^T *m x) ord0 ord0); first by done.
by rewrite Rmult_0_r Rmult_comm.
Qed.


Theorem Perron_Frobenius:
forall A, Mle  0 A -> (forall lam, lam \in (eigenv A) <-> (lam \in (eigens A))) ->
  (eigenv A (specr A) /\ exists x: 'cV_n, 0 <=m: x /\ x !=0 /\ 
   A *m x =  (specr A) *m: x).
Admitted.


Corollary corPF:
forall A a, Mle 0 A -> (forall lam, lam \in (eigenv A) <-> (lam \in (eigens A))) ->
  (a <= (specr A) <-> exists u: 'cV_n,  0 <=m: u  /\  (a *m: u) <=m: A *m u).
Admitted.


Lemma Rayleigh_min:
forall (A:'M_n) (x: 'cV_n), x != 0 ->  
 (forall lam,  lam \in (eigenv (A^T*A)) <-> (lam \in (eigens (A^T*A))) ) -> 
  (lam_min (A^T * A))%:M <=m:  (x^T *m (A^T * A) *m x) / (x^T *m x).
Admitted.

Lemma Rayleigh_max:
forall (A:'M_n) (x: 'cV_n), x != 0 -> 
 (forall lam,  lam \in (eigenv (A^T*A)) <-> (lam \in (eigens (A^T*A))) ) -> 
   (x^T *m (A^T * A) *m x) / (x^T *m x) <=m: (lam_max (A^T * A))%:M  .
Admitted.


 
End eigenR.



Implicit Arguments eigenv[n' T].
Implicit Arguments specr[n'].