Add LoadPath "./interval".

Require Import ssreflect ssrbool ssrfun 
  eqtype ssrnat seq fintype ssralg  matrix bigops.
Require Import Reals.
Require Import Fourier.
Require Import Rstruct min_max real_matrix.



(*
- instantiates a canonical matrix norm: the maximum norm
*)


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

Open Scope R_scope.
Open Scope ring_scope.

Delimit Scope R_scope with Re.
Delimit Scope ring_scope with Ri.


Import GRing.Theory.

Section matNormMaxDef.

Variable m n: nat.

Definition sum_abs (A: 'M_(m, n)) i:R := \sum_j (Rabs (A i j)).
Definition norm_m (A: 'M_(m, n)):R:= \big[Rmax/R0]_i (sum_abs A i).

End matNormMaxDef.

Section matNormMaxProp.



Lemma sum_abs_pos: forall n p (a:'M_(n, p)) i, R0 <= sum_abs a i.
Proof.
move => n p a i; rewrite /sum_abs; elim index_enum.
by rewrite big_nil; right.
move => x s IH; rewrite big_cons.
by apply Rplus_le_le_0_compat; [apply Rabs_pos|].
Qed.

Lemma eq_sum_abs:
forall n p (a:'M_(n, p)) i, sum_abs a i = Rabs (sum_abs a i).
Proof.
move => *; rewrite Rabs_right; [done| ].
apply Rle_ge; apply sum_abs_pos.
Qed.

Lemma pres_ineq1:
forall (n p: nat) (f r:R)( v: 'cV_p),(forall j,  Rabs (v j ord0)<= r )-> ( 0 <= f  )->
\sum_i (Rabs (v i ord0)* f )<= \sum_(i< p)( r*f ).
Proof.
move => n p f v r Hv Hf.
elim: (index_enum ) => //=.
by rewrite !big_nil; right.
move => x s IH ; rewrite !big_cons.
apply Rplus_le_compat; last by done.
apply Rmult_le_compat_r; [by apply Hf|by apply Hv].
Qed.

Lemma pres_ineq:
forall  p (f: 'I_p->R)( v: 'cV_p) r,
  (forall j,  Rabs (v j ord0)<= r )-> (forall x, 0 <= f x )->
  \sum_i (Rabs (v i ord0)* f i)<= \sum_i( r*f i).
Proof.
move => m f v r Hv Hf.
elim: (index_enum ) => //=.
by rewrite !big_nil; right.
move => x s IH; rewrite !big_cons .
apply Rplus_le_compat; last by done.
apply Rmult_le_compat_r; [by apply Hf|by apply Hv].
Qed.


Lemma desc_norm_m:forall n p (a:'M_(n , p)) r, norm_m a = Rmax (sum_abs a r)
    (\big[Rmax/R0]_(k<n|k!=r) (sum_abs a k)).
Proof.
move=>n p a r; rewrite /norm_m.
have H:= sum_abs_pos a.
by apply sxbigD1 with (P:=(fun _  :'I_n=> true)) (j:=r) (F:= (fun t => sum_abs a t)).
Qed.

Lemma normM_ex:
forall n p (M: 'M_(n, p)), (0 < n)%Dn -> exists k, norm_m M = sum_abs M k.
Proof.
move => n p M Hm.
rewrite /norm_m.
apply (@max_exists n Hm).
apply sum_abs_pos.
Qed.

Lemma normM_pos:
forall n p (m: 'M_(n, p)), 0 <= norm_m m.
Proof.
move => n p M; rewrite /norm_m.
elim index_enum.
by rewrite big_nil; right.
move=> t s IH.
rewrite big_cons.
by apply Rmax_ge; [apply sum_abs_pos|] .
Qed.


Lemma normM_max:
forall n p (m: 'M_(n, p)) i, sum_abs m i <= norm_m m. 
Proof.
move => n p a i.
rewrite (desc_norm_m a i); apply RmaxLess1.
Qed.


Lemma ineq_norm_m:
forall n p q (a: 'M_(n, p))  (b: 'M_(p, q)),
norm_m (a *m b) <= norm_m a * norm_m b.
Proof.
move =>n p q a b.
rewrite {1}/norm_m /mulmx /sum_abs.
have H1: \big[Rmax/0]_i (\sum_j Rabs
                             ((\matrix_(i0, k) (\sum_j0 a i0 j0 * b j0 k))
                                i j)) = 
\big[Rmax/0]_i (\sum_j Rabs (\sum_s a i s*b s j)).
apply eq_bigr => i _.
apply eq_bigr => j _.
by rewrite mxE.
rewrite H1.
apply Rle_trans with (\big[Rmax/0]_i (\sum_j (\sum_s (Rabs (a i s * b s j))))).
apply (Rmax_le_f_le ) => i.
apply (sum_le_f_le ) => j.
apply abs_sum_le_sum_abs.
have H2: \big[Rmax/0]_i (\sum_j \sum_s Rabs (a i s * b s j))  = 
               \big[Rmax/0]_i (\sum_j \sum_s Rabs (a i s )*Rabs( b s j)) .
apply eq_bigr =>i _; apply eq_bigr => j _;apply eq_bigr => s _; apply Rabs_mult.
rewrite H2.
have H3: \big[Rmax/0]_i (\sum_j \sum_s Rabs (a i s) * Rabs (b s j)) =
\big[Rmax/0]_i (\sum_s Rabs (a i s)*\sum_j  Rabs (b s j)) .
have H4: \big[Rmax/0]_i (\sum_j \sum_s Rabs (a i s) * Rabs (b s j)) =
\big[Rmax/0]_i (\sum_s \sum_j Rabs (a i s) * Rabs (b s j))  .
apply eq_bigr => i _.
apply: exchange_big.
rewrite H4.
apply eq_bigr => i _; apply eq_bigr => s _.
rewrite -big_distrr //=.
rewrite H3.
apply Rle_trans with (\big[Rmax/0]_i (\sum_s Rabs (a i s) *norm_m b)).
apply (Rmax_le_f_le ) => i.
apply (sum_le_f_le ) => s.
apply Rmult_le_compat_l; first by apply Rabs_pos.
fold (sum_abs b s).
 apply normM_max.
have H5: \big[Rmax/0]_i (\sum_s Rabs (a i s) * norm_m b)  = 
(\big[Rmax/0]_i (\sum_s Rabs (a i s) ))* norm_m b.
rewrite ->Rmult_comm.
have H6 : \big[Rmax/0]_i (\sum_s Rabs (a i s) * norm_m b) =
\big[Rmax/0]_i (\sum_s norm_m b *Rabs (a i s) ) .
apply eq_bigr => i _;
apply eq_bigr => s _; by rewrite /(@GRing.mul) /= Rmult_comm.
rewrite H6.
rewrite distr_mult_max.
apply eq_bigr => i _.
by rewrite -big_distrr.
apply normM_pos.
rewrite H5.
have H7: \big[Rmax/0]_i (\sum_s Rabs (a i s)) = norm_m a.
by rewrite /norm_m /sum_abs.
rewrite H7.
apply Rle_refl.
Qed.


Lemma normM_plus_le:
forall n p (a b: 'M_(n, p)), norm_m (a + b) <= norm_m a + norm_m b.
Proof.
move =>n p a b.
rewrite {1}/norm_m /addmx /sum_abs.
have H1:
\big[Rmax/0]_i (\sum_j Rabs ((\matrix_(i0, j0) (a i0 j0 + b i0 j0)) i j))
= \big[Rmax/0]_i (\sum_j Rabs (a i j + b i j)).
apply eq_bigr => i _.
apply eq_bigr => j _.
by rewrite mxE.
rewrite H1.
apply Rle_trans with (\big[Rmax/0]_i (\sum_j (Rabs (a i j) +Rabs( b i j)))).
apply (Rmax_le_f_le )=> i.
apply (sum_le_f_le ) => s.
apply Rabs_triang.
have H2: \big[Rmax/0]_i (\sum_j (Rabs (a i j) + Rabs (b i j))) =
\big[Rmax/0]_i (\sum_j Rabs (a i j) +\sum_j Rabs (b i j)) .
apply eq_bigr => i _.
rewrite big_split //=.
rewrite H2.
apply Rle_trans with (\big[Rmax/0]_i (\sum_j Rabs (a i j)) +
\big[Rmax/0]_i ( \sum_j Rabs (b i j)) ).
apply (distr_plus_max ).
move => i; fold (sum_abs b i); apply sum_abs_pos.
move => i;  fold (sum_abs a i); apply sum_abs_pos.
rewrite /norm_m /sum_abs.
apply Rle_refl.
Qed.

Lemma hom_pos_normM:
forall n p (m: 'M_(n,p)) a, norm_m (a *m: m) = Rabs a * norm_m m.
Proof.
move =>n p v a.
rewrite /norm_m.
elim: index_enum.
by rewrite !big_nil mulr0.
move=> t s IH; rewrite !big_cons IH .
rewrite <- RmaxRmult; [| apply Rabs_pos].
f_equal.
replace ((Rabs a * sum_abs v t)%Re) with ((Rabs a * sum_abs v t)%Ri);
  last by toR.
rewrite /sum_abs big_distrr /=.
apply: eq_bigr => i _; rewrite mxE Rabs_mult; toR.
Qed.


Lemma zeroM_0:
forall n p, @norm_m n p 0 = 0.
Proof.
move=> n p; rewrite /norm_m; elim: index_enum.
by rewrite big_nil.
move=> t s IH; rewrite !big_cons IH Rmax_left; last by apply sum_abs_pos.
by rewrite  /sum_abs; apply: big1 => * ; rewrite mxE Rabs_R0.
Qed.


Lemma canonical1:
 forall p q (A:'M_(p, q)) i j, Rabs (A i j) <= norm_m A.
Proof.
move =>p q a i j.
rewrite /norm_m.
apply Rle_trans with (sum_abs a i);[|apply max_max; apply sum_abs_pos].
rewrite /sum_abs.
rewrite (bigD1 j ) //=.
rewrite -{1}(Rplus_0_r (Rabs (a i j))).
apply Rplus_le_compat_l.
elim: index_enum.
by rewrite big_nil; right.
move=> x s IH; rewrite !big_cons.
case C: (x!=j); last by done.
rewrite -{1}(Rplus_0_r R0).
apply Rplus_le_compat;[apply Rabs_pos|done].
Qed.


Lemma canonical2:
forall p q (A B:'M_(p, q)), (forall i j, Rabs (A i j) <= Rabs (B i j)) 
                                         ->norm_m A <= norm_m B.
Proof.
move=> p q a b H.
rewrite /norm_m; apply (Rmax_le_f_le ) => i.
by rewrite /sum_abs; apply (sum_le_f_le ) => j.
Qed.


End matNormMaxProp.

Canonical Structure the_norm :=
@MNorm  norm_m normM_pos  
 hom_pos_normM  normM_plus_le
 ineq_norm_m
canonical1 canonical2.

Section sqrMatNormMax.

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

Lemma Hp: (0 < p)%Dn.
by done. Qed.


Lemma norm_unitM:
|| (1 : 'M_p) : the_norm |  = 1.
Proof.
elim: (normM_ex  (1: 'M_p)); last by exact Hp.
move=> k Hk.
rewrite /= Hk /sum_abs /GRing.one  (bigD1 k) //= !mxE eq_refl Rabs_R1.
have H: \sum_(i < p | i != k) Rabs (1%:M k i) = 0.
elim: index_enum.
by rewrite big_nil.
move => x s IH; rewrite big_cons.
case X:(x!=k).
rewrite IH mxE .
have H: (k==x)=false.
apply/eqP.
move/eqP:X=>X.
rewrite /not in X |-*=> H.
have H1: x = k by apply sym_eq.
by apply X.
 rewrite H Rabs_R0.
 rewrite GRing.add0r  //= .
apply IH.
by rewrite H addr0. 
Qed.


Lemma  matr_inv_norm:
forall a: 'M_(p), norm_m  (1 - a) < 1 ->  \det a <>0.
Proof.
move=> a Ha.
apply (@cor1_th5 p' _ norm_unitM).
by rewrite //=.
Qed.


Lemma norm_inv_ser:
forall a: 'M_(p),
norm_m a <1->
norm_m (( 1 - a)^-1) <= 1/(1 - norm_m a).
Proof.
move=> a H.
apply (@norm_inv_ser0 p' _  norm_unitM a H ).
exact norm_unitM.
Qed.



End sqrMatNormMax.


Implicit Arguments norm_m [m n].