<< Statistics`LinearRegression`

<< Graphics`Graphics`

<< Graphics`Colors`

<< Graphics`MultipleListPlot`

<< Statistics`MultiDescriptiveStatistics`

<< LinearAlgebra`MatrixManipulation`

getRows[list_List, text_Real, column_] := 
        list[[First[Transpose[Position[list[[All, column]], text]]], 
            All]];

BlockDiagonalMatrix[mat : {_?MatrixQ .. }] := 
        Module[{colLengths, offsets, padfunc}, 
          colLengths = Last /@ Dimensions /@ mat; 
          offsets = Most[FoldList[Plus, 0, colLengths]]; 
          padfunc[x_, 
              y_] := PadRight[#, Tr[colLengths], 0, 
                  offsets[[y[[1]]]]] & /@ x; 
          Flatten[MapIndexed[padfunc, mat], 1]];

    
ProductMatrix[a_List, 
        b_List] := 
Transpose[
          AppendColumns[
            Table[a, {Length[b]}]]]×
AppendColumns[
          Table[b, {Length[a]}]]

Distance[One_List,  Two_List,  
          norm_]  := 
Module[{lowerBound,  upperBound,  
            diff,  pointSet}, 
diff[
              x_]  :=  Interpolation[One, InterpolationOrder -> 1][
                x] - Interpolation[Two, InterpolationOrder -> 1][
                x]; 
lowerBound  =  
            Max[{Min[One[[All, 1]]], 
                Min[Two[[All, 1]]]}]; 
upperBound  =  
            Min[{Max[One[[All, 1]]], 
                Max[Two[[All, 1]]]}]; 
pointSet = 
            Select[Select[
                Union[One[[All, 1]], 
                  Two[[All, 1]]], # ≥ 
                    lowerBound &], # ≤ 
                  upperBound &]; 
Norm[Map[diff, pointSet], 
              norm]/If[norm >  1, 1, Length[pointSet]]];

ConnT  =  
        Table[{0.5  i, 0.0}, {i, 0, 30}];

    LyreT  =  
        Table[{0.5  i, 0.0}, {i, 0, 30}];

    

    ConnF1 = 
        Join[ {{0.0, 0.0}, {0.5, 0.0}}, 
          Table[{0.5  +  0.5  i,  0.01  i}, {i, 1, 
              28}]];

    LyreF1 = 
        Join[ Table[{0.1 + 0.5  i, 0.14  - 0.005  i}, {i, 0, 
              28}],  {{14.6, 0.0}}];

DisplayTogether[
ListPlot[Union[ConnF1, LyreF1],  
        PlotJoined  -> True,  PlotStyle -> {Blue}], 
      ListPlot[Union[ConnT, LyreT],  PlotJoined  -> True,  
        PlotStyle -> {Red}],  
      PlotRange -> {-0.2, 0.4}]

DisplayTogether[
ListPlot[ConnF1, 
        PlotJoined  -> True,  
        PlotStyle -> {Green}],  
ListPlot[LyreF1,  
        PlotJoined  -> True,  PlotStyle -> {Blue}], 
      ListPlot[Union[ConnT, LyreT],  PlotJoined  -> True,  
        PlotStyle -> {Red}],  
      PlotRange -> {-0.1, 0.3}]

analyze[errorMole_,  errorLaser_] :=
 CompoundExpression[

  
errorBundle = 0.0; 
trackerConn  =  -1. ;trackerLyre  =  16. ;
PointErrorConn[y_, 
                error_] := 
((errorBundle^2) + 
                  (errorMole^2) + (((errorLaser (y  -  trackerConn)))^2)  + (error^2));
PointErrorLyre[y_, 
                error_] := 
((errorBundle^2) + 
                  (errorMole^2) + (((errorLaser (trackerLyre  -  y)))^2) + (error^2));
PointWeightConn[{y_, 
                  error_}] :=  
{y, 
                1. /((errorBundle^2) + 
                      (errorMole^2) + (((errorLaser (y  -  trackerConn)))^2)  + (error^2))};
PointWeightLyre[{y_, 
                  error_}] := 
{y, 
                1. /((errorBundle^2) + 
                      (errorMole^2) + (((errorLaser (trackerLyre  -  y)))^2) + (error^2))};

  
sawToothAvg  =  
              Distance[ConnF1,  LyreF1, 1];
sawToothMax =  
              Distance[ConnF1,  LyreF1, ∞];
distConnMax =  
              Distance[ConnF1, ConnT, ∞];
distLyreMax =  
              Distance[LyreF1, LyreT, ∞];

  
fitOneConn =Regress[ConnF1,Table[x^a, {a, 0, 1}],x,
  RegressionReport->{BestFit,BestFitParameters,AdjustedRSquared}];
fitOneLyre =Regress[LyreF1,Table[x^a, {a, 0, 1}],x,
  RegressionReport->{BestFit,BestFitParameters,AdjustedRSquared}];

  
analyzeAperture[ConnF1_List,  LyreF1_List, fitConn_,  
                fitLyre_, 
distConn_,  
                distLyre_] := 
Module[{ygrekConn, 
                  ygrekLyre, zetConn,  zetLyre, slope, shift, ref, 
                  shiftConn, shiftLyre, addConn, corrConn, addLyre, 
                  corrLyre}, 
ygrekConn  =  
                  Min[ConnF1[[All, 1]]] + 
                    distConn; 
ygrekLyre  =  
                  Max[LyreF1[[All, 1]]] - 
                    distLyre; 
zetConn  =  BestFit /. 
                      fitConn /. 
                    x ->  
                      ygrekConn; 
zetLyre  =  BestFit /. 
                      fitLyre /. 
                    x ->  
                      ygrekLyre; 
slope  =  (zetLyre - 
                        zetConn)/(ygrekLyre  -  
                        ygrekConn); 
shift  =  
                  zetConn - slope ygrekConn; 
ref  =  
                  slope x  +  
                    shift; 

shiftConn =  (BestFit /. fitConn) -  
                    ref; 
shiftLyre =  (BestFit /. 
                        fitLyre) -  
                    ref; 

addConn[{abs_,  ord_}]  :=  If[
                    abs  ≤  ygrekConn,  {abs, ord},  {abs,  
                      ord -  shiftConn /. 
                        x -> abs}]; 
corrConn = 
                  Map[addConn, 
                    ConnF1]; 
addLyre[{abs_,  
                      ord_}]  :=  
                  If[abs  ≥  ygrekLyre,  {abs, 
                      ord},  {abs,  
                      ord -  shiftLyre /. 
                        x -> abs}]; 
corrLyre = 
                  Map[addLyre, LyreF1]; 
Distance[
                  corrConn,  corrLyre,  ∞]
];

  
i =  0; 
results  =  {};
\(Do[
              CompoundExpression[
	
                i = i + 1; 
	
                j = 0; 
Do[
                  CompoundExpression[
	
                    j = j + 1; 
	
                    resultOne = 
                      analyzeAperture[
ConnF1, \ LyreF1, 
                        fitOneConn, \ 
                        fitOneLyre, 
0.2\  + \ \((i - 
                                1)\)\ 0.2, 
                        0.2 + \ \((j - 
                                1)\)\ 0.2]; 
results = \ 
                      Join[results, 
{Flatten[{0.2 + \ \((i - 1)\)\ 0.2, \ 0.2 + \ \((j - 1)\)\ 0.2, \ 
                              resultOne}]}]]
, {10}]], \ {10}]\);

lastCorrectPointOneConn = getRows[
                  results, Min[results[[All, {3}]]], 3][[1, 1]];
lastCorrectPointOneLyre = getRows[
                  results, Min[results[[All, {3}]]], 3][[1, 2]];

  

  
ygrekOneConn  =  
              Min[ConnF1[[All, 1]]] + lastCorrectPointOneConn;
ygrekOneLyre  =  
              Max[LyreF1[[All, 1]]] - lastCorrectPointOneLyre;
zetOneConn  =  BestFit /. 
                  fitOneConn /. x ->  ygrekOneConn;
zetOneLyre  =  BestFit /. 
                  fitOneLyre /. x ->  ygrekOneLyre;
slopeOne  =  (zetOneLyre - 
                    zetOneConn)/(ygrekOneLyre  -  ygrekOneConn);
shiftOne  =  
              zetOneConn - slopeOne ygrekOneConn;
refOne  =  
              slopeOne x  +  shiftOne;

 (* 
                original data with fits and reference line displayed *) 
orig :=DisplayTogether[
Plot[BestFit /. fitOneConn, {x,  -1,  16}, 
                    PlotStyle -> {Green}],
Plot[
                    BestFit /. fitOneLyre, {x,  -1,  16}, 
                    PlotStyle -> {Blue}],
Plot[
                    refOne,  {x,  -1,  16},  
                    PlotStyle -> {Black}],
ListPlot[ConnF1, 
                    PlotJoined  -> True,  PlotStyle -> {Green}], 
ListPlot[LyreF1,  
                    PlotJoined  -> True,  PlotStyle -> {Blue}],ListPlot[Union[ConnT, LyreT],  
                    PlotJoined  -> True,  PlotStyle -> {Red}], PlotRange -> {-0.8, 0.8},
PlotLabel->
\(ToString[
                        magnet]\  <> \ " A1, max saw tooth "\  <> \ \ ToString[SetAccuracy[sawToothA1Max, 
                          4]]\  <> 
"
 max dist conn "\  <> \ ToString[SetAccuracy[distA1ConnMax, 
                          4]] <> 
", lyre "\  <> \ 
                      ToString[SetAccuracy[distA1LyreMax, 4]]\)];

  
diffPlots:=DisplayTogether[Plot[Interpolation[ConnF1, 
                          InterpolationOrder -> 1][
                        x] - Interpolation[LyreF1, 
                          InterpolationOrder -> 1][x], {x, 1, 14.0},  
                    PlotStyle -> {Red}], ImageSize -> 300,PlotLabel->
ToString[
                        magnet] <>  ", diff, A1 (red) and A2 (blue)", Background->Floral];

  
shiftOneConn =  (BestFit /. 
                    fitOneConn) -  refOne;
shiftOneLyre =  (BestFit /. 
                    fitOneLyre) -  refOne;

  
addOneConn[{abs_,  ord_}]  :=  
              If[abs  ≤  ygrekOneConn,  {abs, ord},  {abs,  
                  ord -  shiftOneConn /. x -> abs}];
corrOneConn = Map[addOneConn, ConnF1];
addOneLyre[{abs_,  ord_}]  :=  
              If[abs  ≥  ygrekOneLyre,  {abs, 
                  ord},  {abs,  ord -  shiftOneLyre /. x -> abs}];
corrOneLyre = Map[addOneLyre, LyreF1];

  
corrsawToothA1Max  =  
              Distance[corrOneConn,  corrOneLyre,  ∞];

  
corr :=DisplayTogether[
Plot[refOne,  {x,  -1,  16},  
                    PlotStyle -> {Black}],
ListPlot[corrOneConn, 
                    PlotJoined  -> True,  PlotStyle -> {Green}], 
ListPlot[corrOneLyre,  
                    PlotJoined  -> True,  PlotStyle -> {Blue}], PlotRange -> {-0.001, 0.01},
PlotLabel->
"Plied curves + reference line" ,
 Background->Floral, ImageSize ->  210];

  
modelFunction = a y  +  b;

 (* 
                caclulate shifts for the error covariation matrix *) 
addOneConn[{abs_,  ord_}]  :=  
              If[abs  ≤  ygrekOneConn,  {abs, ord},  {abs, 
                  shiftOneConn /. x -> abs}];
shiftConn = Map[addOneConn, ConnF1];
addOneLyre[{abs_,  ord_}]  :=  
              If[abs  ≥  ygrekOneLyre,  {abs, 
                  ord},  {abs,  shiftOneLyre /. x -> abs}];
shiftLyre = Map[addOneLyre, LyreF1];

 (* 
                diagonal matrices for the error covariation matrix *) 
CMConn[thisData_List,  lastCorrectPoint_,  shifts_]:=Module[{Cov, i,  j, rowI, uncertainty},
Cov = {};i = 0;
Do[CompoundExpression[
	i  =  i +  1;
	j  =  0;rowI = {};
	Do[CompoundExpression[
	j  =  j  +  1;
	\(uncertainty\  = \ 
                                      If[i == j, 
	
                                        If[thisData[\([i, 
                                        1]\)]\  ≤ \ 
                                        lastCorrectPoint, \ 
		PointErrorConn[thisData[\([i, 1]\)], 
                                        0], 
		PointErrorConn[thisData[\([i, 1]\)], \ shifts[\([i, 
                                        2]\)]]], 
	\((\ 
                                        shifts[\([i, 
                                        2]\)])\) \((\ \ shifts[\([j, 
                                        2]\)])\)]\);
	rowI = Append[rowI, uncertainty];],{Length[thisData]}];
	Cov = Append[Cov, rowI];],
{Length[thisData]}];
Return [Cov];];
CMLyre[thisData_List,  lastCorrectPoint_,  shifts_]:=Module[{Cov, i,  j, rowI, uncertainty},
Cov = {};i = 0;
Do[CompoundExpression[
	i  =  i +  1;
	j  =  0;rowI = {};
	Do[CompoundExpression[
	j  =  j  +  1;
	\(uncertainty\  = \ 
                                      If[i == j, 
	
                                        If[thisData[\([i, 
                                        1]\)]\  ≥ \ 
                                        lastCorrectPoint, \ 
		PointErrorLyre[thisData[\([i, 1]\)], 
                                        0], 
		PointErrorLyre[thisData[\([i, 1]\)], \ shifts[\([i, 
                                        2]\)]]], 
	\((\ 
                                        shifts[\([i, 
                                        2]\)])\) \((\ \ shifts[\([j, 
                                        2]\)])\)]\);
	rowI = Append[rowI, uncertainty];],{Length[thisData]}];
	Cov = Append[Cov, rowI];],
{Length[thisData]}];
Return [Cov];];
CMOneConn  =  
              CMConn[ConnF1,  ygrekOneConn,  shiftConn];
CMOneLyre  =  
              CMLyre[LyreF1,  ygrekOneLyre,  shiftLyre];
  

corrCoeff  =  1.0;
(*CorrelationMatrix[AppendRows[shiftConn[[All, {2}]], 
                      shiftLyre[[All, {2}]]]][[1, 2]];*)
CMConnLyre  =  
              corrCoeff ProductMatrix[shiftConn[[All, 2]], 
                  shiftLyre[[All, 2]]];
CMLyreConn  =  
              corrCoeff ProductMatrix[shiftLyre[[All, 2]], 
                  shiftConn[[All, 2]]];
CM =  
              BlockMatrix[{{CMOneConn, CMConnLyre}, {CMLyreConn, 
                    CMOneLyre}}];
WM = Inverse[CM];

  
SubtractModel[{abscissa_,  
                  ordinata_}]  :=   ordinata  -  modelFunction /. 
                y -> abscissa;
residuals = Partition[Map[
                  SubtractModel, 
Join[ConnF1, LyreF1]], 
                1];
chi  =  
              First[First[Transpose[residuals] . WM . residuals]];
midPlane = 
              FindMinimum[chi, {a, -0.001}, {b, -0.001}];
birgeRatio = (Sqrt[midPlane[[1]]/(Length[residuals] - 2)]);

  
parameters = parameters  =  {a /. 
                    midPlane[[2]], b /. midPlane[[2]]};
yShift = 0; zetShift = parameters[[1]];
cosAngle = (Sqrt[1. /(1 + parameters[[2]]^2)]);
sinAngle  = 
              parameters[[2]]/(Sqrt[(1 + parameters[[2]]^2)]);
newY[{y_, 
                  z_}] := (y - yShift) cosAngle + (z - zetShift) 
                  sinAngle;
newZet[{y_, 
                  z_}] := -(y - yShift) sinAngle + (z - 
                      zetShift) cosAngle;

  
plyOneConn  =  
              Transpose[
{ Map[newY, corrOneConn],  
                  Map[newZet, corrOneConn]}];
plyOneLyre  =  
              Transpose[
{ Map[newY, corrOneLyre],  
                  Map[newZet, corrOneLyre]}];

plydistA1ConnMax =  
              Distance[plyOneConn, ConnT, ∞];
plydistA1LyreMax =  
              Distance[plyOneLyre, LyreT, ∞];

ply :=DisplayTogether[
ListPlot[plyOneConn, PlotJoined  -> True,  
                    PlotStyle -> {Green}], 
ListPlot[plyOneLyre,  
                    PlotJoined  -> True,  PlotStyle -> {Blue}],ListPlot[Union[ConnT, LyreT],  
                    PlotJoined  -> True,  PlotStyle -> {Red}], PlotRange -> {-0.8, 0.8},
PlotLabel->
ToString[
                        magnet]  <>  " A1" <> 
"\n \
max dist conn "  <>  
                      ToString[
                        SetAccuracy[plydistA1ConnMax, 
                          4]] <> 
", lyre "  <>  
                      ToString[SetAccuracy[plydistA1LyreMax, 4]]];

  
fitCorr=Regress[Join[corrOneConn,  corrOneLyre],{1, x},x,
Weights  -> 
Join[
                        1. /((errorLaser^2) ((corrOneConn[[All, 
                                        1]] -  
                                      trackerConn)^2)), 
1. /((errorLaser^2) ((trackerLyre - 
                                      corrOneLyre[[All, 1]])^2))],
RegressionReport->{BestFit, FitResiduals, SequentialSumOfSquares,PartialSumOfSquares,ParameterCITable,
SinglePredictionCITable,ParameterConfidenceRegion}];

  
fitGA:=DisplayTogether[Plot[BestFit /. fitCorr,{x, 1, 16}, PlotStyle->{Blue}],
ListPlot[ConnF1, PlotJoined  -> True, PlotStyle->{Red}],
ListPlot[LyreF1, PlotJoined  -> True, PlotStyle->{Red, Dashing[{0.01, 0.01}]}],
Plot[
                    parameters[[1]] x  +  
                      parameters[[2]], {x, -1, 16},  
                    PlotStyle -> {Green, AbsoluteThickness[1.5]}],
\(PlotLabel -> \ 
"Original Data with Corrected GA

weighted fit : "\  <> \ ToString[parameters[\([1]\)] x\  + \ parameters[\([2]\)]\ , 
                        StandardForm]\  <> " (green)" <> \ "

correction : "\  <> ToString[BestFit /. fitCorr\ , 
                        StandardForm] <> 
" (blue)"\) ,
   Background->Floral
]; 
];

analyze[0.07, 5.  (10^(-3))];
fitGA;  {lastCorrectPointOneConn, lastCorrectPointOneLyre, midPlane, 
      birgeRatio, 
      Sqrt[(Total[(FitResiduals /. fitCorr) (FitResiduals /. 
                    fitCorr)]/(Length[corrOneConn] + 
                Length[corrOneLyre]  -  2))]}

av  = Fit[Join[ConnF1, LyreF1], {1, x}, x];DisplayTogether[Plot[BestFit /. fitCorr,{x, 1, 16}, PlotStyle->{Blue}],
Plot[av,{x, 1, 16}, PlotStyle->{Orange}],
ListPlot[ConnF1, PlotJoined  -> True, PlotStyle->{Red}],
ListPlot[LyreF1, PlotJoined  -> True, PlotStyle->{Red, Dashing[{0.01, 0.01}]}],
Plot[
            parameters[[1]] x  +  parameters[[2]], {x, -1, 
              16},  PlotStyle -> {Green, AbsoluteThickness[1.5]}],
\(PlotLabel -> \ 
"Original Data with Corrected GA

weighted fit : "\  <> \ ToString[parameters[\([1]\)] x\  + \ parameters[\([2]\)]\ , 
                StandardForm]\  <> " (green)" <> \ "
correction : "\  <> ToString[BestFit /. fitCorr\ , 
                StandardForm] <> 
" (blue)
simple average : "\  <> \ ToString[av, StandardForm]\) ,
ImageSize  ->  400, Background->Floral
];

ParameterCITable/.fitCorr

Show[Graphics[{RGBColor[0, 1, 0],Disk[{b /. midPlane[[2]], 
                  a /. midPlane[[2]]}, 
{0.00006, 
                  0.000004}]}],
Graphics[{RGBColor[0, 0, 1],  
            ParameterConfidenceRegion /. fitCorr}],Axes -> True,AxesLabel -> {"shift", "slope"},AxesOrigin -> {-0.002, 0. },
PlotRange ->  {{-0.002, 
              0.0042}, {-0.00012,  0.00002}},Background->Floral]

BestFit /. fitCorr

(a /. midPlane[[2]]) y  +  b /. midPlane[[2]]

DisplayTogether[
Plot[refOne,  {x,  -1,  16},  
            PlotStyle -> {Black}],
ListPlot[Union[ConnF1, LyreF1],  
            PlotJoined  -> True,  PlotStyle -> {Red}], ListPlot[Union[corrOneConn, corrOneLyre],  
            PlotJoined  -> True,  PlotStyle -> {Blue}], PlotRange -> {-0.005, 0.3},
PlotLabel->
"Plied curves + reference line" ,
 Background->Floral, ImageSize ->  210];

calculateWeightOneConn[{abs_,  ord_}]  :=  
            If[abs  ≤  
                ygrekOneConn,  
PointWeightConn[{abs, 
                  0}], 
PointWeightConn[ {abs, 
                  shiftOneConn /. x -> abs}]];
weightOneConn = 
            Map[calculateWeightOneConn, ConnF1];
calculateWeightOneLyre[{abs_,  
                ord_}]  :=  
            If[abs  ≥  
                ygrekOneLyre,  
PointWeightLyre[{abs, 
                  0}], 
PointWeightLyre[ {abs, 
                  shiftOneLyre /. x -> abs}]];
weightOneLyre = 
            Map[calculateWeightOneLyre, LyreF1];
DisplayTogether[
ListPlot[weightOneConn, PlotJoined -> True, PlotStyle->{Blue}],
ListPlot[weightOneLyre, PlotJoined -> True,PlotStyle->{Blue}], PlotRange -> {0, 300}]

DisplayTogether[
Plot[BestFit /. fitOneConn,{x,  -1,  16}, PlotStyle->{Orange, Dashing[{0.01, 0.01 }]}],
Plot[BestFit /. fitOneLyre,{x,  -1,  16},PlotStyle->{Orange, Dashing[{0.01, 0.01 }]}],
Plot[refOne,  {x,  -1,  16},  
            PlotStyle -> {Black}],
ListPlot[ConnF1, 
            PlotJoined  -> True,  PlotStyle -> {Green}], 
ListPlot[LyreF1,  
            PlotJoined  -> True,  PlotStyle -> {Blue}],ListPlot[Union[ConnT, LyreT],  PlotJoined  -> True,  
            PlotStyle -> {Red}], PlotRange -> {-0.01, 0.01},
PlotLabel ->  "Original data + \
reference line",Background->Floral,
ImageSize ->  210];

Min[results]

ListPlot3D[
      Partition[Flatten[results[[All, {3}]]], 
        10], 
PlotLabel  ->  "Distance between \
plied curves\nsearch by 10 steps of 20cm"]

getRows[results, Min[results], 3]

analyzeAperture[ConnF1, LyreF1, fitOneConn,  fitOneLyre, 0.5, 
      0.5]

DisplayTogether[ListPlot[shiftConn], ListPlot[shiftLyre]]