(* ::Package:: *)

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)



(* :Name: PwM`Chap08Numerics` *)


(* :Book: Programming with Mathematica: An Introduction *)


(* :Author: Paul Wellin *)


(* :Summary: This package contains the functions defined in Chapter 8, Numerics. *)


(* :Context: PwM`Chap08Numerics` *)


(* :Package Version: 1.0 *)


(* :Copyright: Copyright 2012, Paul R. Wellin *)


(* :History: Version 1.0 by Paul Wellin, Dec 2012 *)


(* :Keywords: autocorrelation, condition number, correlogram, digit sum, eigenvector, root finding, lag plot, radius of gyration, sparse array, tridiagonal matrix *)


(* :Source: Chap08Numerics.nb *)


(* :Mathematica Version: 8.0.4 *)


Remove[RadiusOfGyrationTensor]


BeginPackage["PwM`Chap08Numerics`"]


AutoCorrelation::usage="AutoCorrelation[data, lag] returns the auto-correlation statistic of data with time lag, lag.";


ConditionNumber::usage="ConditionNumber[mat] computes the condition number of the matrix mat.";


Correlogram::usage="Correlogram[data, {lagmin, lagmax, incr}] displays a correlogram of data with time lags running from lagmin to lagmax in increments of incr.";


DigitSum::usage="DigitSum[n, b] returns the digit sum of the number n in base b.";


EigenvectorPlot::usage="EigenvectorPlot[data] displays an eigenvector plot of data.";


findRoot::usage="findRoot[fun, {var, init}] implements Newton's root-finding method to find a zero of the function fun[var] starting with an initial estimate of init.";


LagPlot::usage="LagPlot[data, lag] returns a lag plot of data.";


MiddleSquareGenerator::usage="MiddleSquareGenerator[n, seed] generates n + 1 random numbers starting with the seed and using von Neumann's middle square algorithm.";


RadiusOfGyrationTensor::usage="RadiusOfGyrationTensor[lis] computes the radius of gyration tensor of the numerical list lis.";


RandomSparseArray::usage="RandomSparseArray[n] returns an n \[Cross] n sparse array with random numbers between 0 and 1 on the diagonal and zeros everywhere else.";


TridiagonalMatrix::usage="TridiagonalMatrix[n, p, q] returns an n \[Cross] n tridiagonal matrix with element p on the diagonal and q on the off-diagonal.";


Options[Correlogram]=Join[{Coefficient->0.05},Options[ListPlot]];


Options[findRoot]={
MaxIterations:>$IterationLimit,
PrecisionGoal->Automatic,
WorkingPrecision->Automatic
};


Begin["`Private`"]


AutoCorrelation[data_,lag_:1]:=Correlation[Drop[data,lag],Drop[data,-lag]]


ConditionNumber[mat_?MatrixQ]:=Module[{sv=SingularValueList[mat]},
First[sv]/Last[sv]]


Correlogram[data_,{lagmin_,lagmax_,incr_:1},opts:OptionsPattern[]]:=
Module[{rh,corrs},
rh=OptionValue[Coefficient];
corrs=Table[{lag,AutoCorrelation[data,lag]},{lag,lagmin,lagmax,incr}];
ListPlot[corrs,
FilterRules[{opts},Options[ListPlot]],AspectRatio->.4,Frame->True,Axes->False,PlotRange->Automatic,FrameTicks->{{Automatic,False},{Automatic,False}},
Epilog->{Thin,Dashed,Line[{{0,rh},{(lagmax-lagmin+1)/incr,rh}}],Line[{{0,-rh},{(lagmax-lagmin+1)/incr,-rh}}]}]
]


DigitSum[n_,base_:10]:=Total[IntegerDigits[n,base]]


EigenvectorPlot[data:{{_,_}..},tensor_]:=Module[{T=tensor,cmx,cmy,x,\[Lambda]1,\[Lambda]2,v1x,v1y,v2x,v2y,l1,l2},
{\[Lambda]1,\[Lambda]2}=Eigenvalues[T];
{cmx,cmy}=Mean[data];
{{v1x,v1y},{v2x,v2y}}=Eigenvectors[T];
l1=v1y/v1x (x-cmx)+cmy;
l2=v2y/v2x (x-cmx)+cmy;
Show[{
ListLinePlot[data,PlotStyle->LightGray],
Plot[l1,{x,cmx-\[Lambda]1,cmx+\[Lambda]1},PlotStyle->{Gray,Thick}],
Plot[l2,{x,cmx-\[Lambda]2,cmx+\[Lambda]2},PlotStyle->{Gray,Thick}],
Graphics[{PointSize[Large],Red,Point[Mean[data]]}]
},AspectRatio->Automatic]]


findRoot[fun_,{var_,init_?NumericQ},opts:OptionsPattern[]]:=Module[{maxIterations,precisionGoal,workingPrecision,initx,df=Derivative[1][fun],result},
{maxIterations,precisionGoal,workingPrecision}=OptionValue[{MaxIterations,PrecisionGoal,WorkingPrecision}];
If[precisionGoal===Automatic,precisionGoal=Precision[init]];If[workingPrecision===Automatic,workingPrecision=precisionGoal+10];initx=SetPrecision[init,workingPrecision];result=SetPrecision[FixedPoint[#1-fun[#1]/df[#1]&,initx,maxIterations],precisionGoal];
{var->result}]


LagPlot[data_,lag_:1,opts:OptionsPattern[ListPlot]]:=ListPlot[Transpose[{Drop[data,lag],Drop[data,-lag]}],opts]


MiddleSquareGenerator[n_,seed_:1234567890]:=Module[{tmp={seed},s2,len,s=seed},
Do[
s2=IntegerDigits[s^2];
len=Length[s2];
s=FromDigits[If[len<20,PadLeft[s2,20,0],s2][[6;;15]]];
AppendTo[tmp,s],
{n}];
tmp
]


RadiusOfGyrationTensor[lis_]:=Module[{cmx,cmy,xcoords,ycoords,xy},
{cmx,cmy}=Mean[lis];
{xcoords,ycoords}=Transpose[lis];
xy=(xcoords-cmx).(ycoords-cmy)/Length[lis];
{{Mean[(xcoords-cmx)^2],xy},{xy,Mean[(ycoords-cmy)^2]}}]


RandomSparseArray[n_Integer]:=SparseArray[{Band[{1,1}]:>RandomReal[]},{n,n}]


TridiagonalMatrix[n_,p_,q_]:=SparseArray[{Band[{1,1}]->p,Band[{1,2}]->q,Band[{2,1}]->q},{n,n}]


End[]


EndPackage[]



