%time	day	f1	f2	f3	w1
X=[1	7	0	4	9	049
2	8	5	5	8	558
3	9	1	8	6	186
4	10	7	4	3	743
5	11	0	1	5	015
6	12	1	2	7	127
7	13	4	9	8	498
8	14	1	8	9	189
9	15	7	4	6	746
10	16	7	1	3	713
11	17	6	4	2	642
12	18	3	2	2	322
13	19	8	2	0	820
14	20	2	3	9	239
15	21	8	2	1	821
16	22	8	6	6	866
17	23	9	8	7	987
18	24	5	5	7	557
19	25	7	1	4	714
20	26	9	4	9	949
21	27	2	6	7	267
22	28 	5	6	4	564
23	29 	1	8	7	187
24	30 	6	2	1	621
25	31 	0	2	6	026
26	1 	1	3	3	133
27	2 	4	4	6	446
28	3 	0	2	4	024
29	4	6	6	3	663
30	5	5	9	6	596
31	6 	8	7	3	873
32	7 	7	2	0	720
33	8 	9	6	9	969
34	9	6	2	2	622
35	10	6	7	1	671
36	11	4	8	9	489
37	12	7	9	2	792
38	13	0	5	5	055
39	14	0	2	9	029
40	15	9	8	0	980
41	16	2	9	3	293
42	17	9	1	9	919
43	18	4	1	7	417
44	19	4	9	5	495
45	20	3	7	7	377
46	21	9	5	2	952
47	22	5	5	6	556
48	23	1	2	4	124
49	24	7	6	7	767
50	25	0	5	6	056
51	26	8	7	0	870
52	27	6	8	1	681
53	28 	4	4	3	443
54	29	0	9	9	099
55	30	0	2	5	025
56	1	8	3	2	832
57	2	0	2	5	025
58	3	2	0	6	206
59	4	3	5	3	353
60	5	2	8	4	284
61	6	1	4	5	145
62	7	1	6	6	166
63	8	9	0	2	902
64	9	2	2	7	227
65	10	4	1	4	414
66	11	1	2	9	129
67	12	4	7	2	472
68	13	1	5	4	154
69	14	3	5	9	359
70	15	2	7	9	279
71	16	6	3	9	639
72	17	2	1	3	213
73	18	1	9	1	191
74	19	8	7	5	875
75	20	7	4	7	747
76	21	4	0	3	401
77	22	7	5	3	753
78	23	1	8	7	187
79	24	4	5	6	456
80	25	3	5	2	352
81	26	9	5	5	955
82	27	8	3	7	837
83	28	1	5	1	151
84	29	7	4	8	748
85 	30	2	5	0	250
86 	31	6	5	1	651
87	1	1	1	3	113
88	2	3	9	0	390
89	3	5	0	2	502
90	4	4	5	6	456
91	5	4	4	8	448
92	6	7	4	0	740
93	7	2	1	5	215
94	8	6	0	5	605
95	9	6	1	2	612
96	10	7	7	9	779
97	11	9	3	9	939
98	12	3	1	6	316
99	13	7	6	3	763
100	14	6	1	1	611
101	15	4	6	5	465
102	16	3	1	4	314
103	17	6	8	9	689
104	18	7	3	3	733
105	19	8	1	5	815
106	20	2	0	1	201
107	21	7	5	7	757
108	22	2	5	9	259
109	23	9	3	7	937
110	24	4	0	9	409
111	25	4	0	4	404
112	26	5	4	9	549
113	27	6	7	1	671
114	28	3	2	3	323
115	29	9	9	0	990
116	30	1	8	6	186
117	1	4	6	1	461
118	2	5	2	1	521
119	3	1	1	5	115
120	4	4	7	5	475
121	5	5	6	7	567
122	6	4	4	4	444
123	7	3	6	4	364
124	8	6	1	5	615
125	9	0	8	6	086
126	10	1	3	6	136];
[m,n]=size(X);
a=1;
b=m;
%time=X(:,1:1);
%day=X(:,2:2);
%f1=X(:,3:3);	
%f2=X(:,4:4);	
%f3=X(:,5:5);	
%w1=X(:,6:6);
%
%
% Locate and Generate Double and Triples
%
%
doubleyn=[];
ddigit=[];
 for i=1:length(X(a:b,:));
     if X(i:i,3:3)==X(i:i,4:4);
         doubleyn(i)=1;
         ddigit(i)=X(i:i,3:3);
     elseif X(i:i,3:3)==X(i:i,5:5);
         doubleyn(i)=1;
         ddigit(i)=X(i:i,3:3);
     elseif X(i:i,4:4)==X(i:i,5:5);
            doubleyn(i)=1;
            ddigit(i)=X(i:i,4:4);
        else
         doubleyn(i)=0;
         ddigit(i)=nan;
     end;
end;
XD=[X(a:b,:),doubleyn',ddigit'];
%
%
% end  Locate and Generate Double and Triples
%
%
%
%  Begin validation to generate and compare confidence sets
%
%


validationsize=0;
sumtsuccess=0;
rootsuccess=0;
sumrate=0;
rootrate=0;
for iv=b-validationsize-1:b-1;


%
%
% Sum and Root calculations
%
%
sumt= sum(X(a:iv+1,3:5),2);
cumsumt=cumsum(sumt);
numstr=num2str(sumt);
ind=find(numstr(:,1)==' ');
numstr(ind,1)='0';
roott=str2num(numstr(:,1))+str2num(numstr(:,2));
ind=find(roott==10);
roott(ind)=1;
%
%
%
% Root total Time since last drawn 
% Root total Average time between draws
%
%
rootconfidenceset=[];
roottsld=[];
 for i=1:10;
     ind=find(roott==i-1);
     if length(ind)>1 
     mindiff= nanmin(sort(diff(ind)));
     maxdiff= nanmax(sort(diff(ind)));
     avediff= nanmean(sort(diff(ind)));
     meddiff = nanmedian(sort(diff(ind)));
     %
     %5 percent trimmed mean
     %
     trimdiff=trimmean(sort(diff(ind)),0);
     roottsld(i:i,1:1)=mindiff;
     roottsld(i:i,2:2)=maxdiff;
     roottsld(i:i,3:3)=avediff;
     roottsld(i:i,4:4)=trimdiff;
     %
     %time since last drawn
     %
     roottsld(i:i,5:5)=meddiff;
     roottsld(i:i,6:6)=(length(roott))-ind(length(ind));
     roottsld(i:i,7:7)=i-1;
     if roottsld(i:i,6:6) > roottsld(i:i,5:5) 
	     rootconfidenceset=[rootconfidenceset,i-1];
     end;
     if roottsld(i:i,6:6) > roottsld(i:i,4:4) 
	     rootconfidenceset=[rootconfidenceset,i-1];
     end;
     if roottsld(i:i,6:6) > roottsld(i:i,3:3) 
	     rootconfidenceset=[rootconfidenceset,i-1];
     end;
     else
     roottsld(i:i,1:1)=nan;
     roottsld(i:i,2:2)=nan;
     roottsld(i:i,3:3)=nan;
     roottsld(i:i,4:4)=nan;
     roottsld(i:i,5:5)=nan;
     roottsld(i:i,6:6)=nan;
     roottsld(i:i,7:7)=i-1;
     end; 
     end;
     roottsld
     sumtconfidenceset=[];
     sumttsld=[];
 for i=1:28;
     ind=find(sumt==i-1);
     if length(ind)>1 
     mindiff= nanmin(sort(diff(ind)));
     maxdiff= nanmax(sort(diff(ind)));
     avediff= nanmean(sort(diff(ind)));
     meddiff = nanmedian(sort(diff(ind)));
     %
     %5 percent trimmed mean
     %
     trimdiff=trimmean(sort(diff(ind)),0);
     sumttsld(i:i,1:1)=mindiff;
     sumttsld(i:i,2:2)=maxdiff;
     sumttsld(i:i,3:3)=avediff;
     sumttsld(i:i,3:3)=avediff;
     sumttsld(i:i,4:4)=trimdiff;     
     sumttsld(i:i,5:5)=meddiff;
     %
     %time since last drawn
     %
     sumttsld(i:i,6:6)=(length(sumt))-ind(length(ind));
     sumttsld(i:i,7:7)=i-1;
     if sumttsld(i:i,6:6) > sumttsld(i:i,5:5) 
	     sumtconfidenceset=[sumtconfidenceset,i-1];
     end;
     if sumttsld(i:i,6:6) > sumttsld(i:i,4:4) 
	     sumtconfidenceset=[sumtconfidenceset,i-1];
     end;
     if sumttsld(i:i,6:6) > sumttsld(i:i,3:3) 
	     sumtconfidenceset=[sumtconfidenceset,i-1];
     end;
     else
     sumttsld(i:i,1:1)=nan;
     sumttsld(i:i,2:2)=nan;
     sumttsld(i:i,3:3)=nan;
     sumttsld(i:i,4:4)=nan;
     %time since last drawn
     sumttsld(i:i,5:5)=nan; 
     sumttsld(i:i,6:6)=nan; 
     sumttsld(i:i,7:7)=i-1;  
     end;
     end;
     sumttsld
     nextpick=iv+2
     rootconfidenceset=unique(rootconfidenceset)
     rootsuccess=ismember(roott(iv),rootconfidenceset)+rootsuccess
     sumtconfidenceset=unique(sumtconfidenceset)
     sumtsuccess=ismember(sumt(iv),sumtconfidenceset)+sumtsuccess
   
%
% End validation
%
     	end;
%
% %
%
%
%    
%Update Data matrix to include sum total and root total
XU=[XD,sumt,roott];
%looking at last 20
%
%
sumtsort=sort(sumt(length(sumt)-18:length(sumt)));
roottsort=sort(roott(length(roott)-18:length(roott)));
%one up one down
%Note this is reversed row-wise
oneuponedown=[X(m:m,6:6)+100,X(m:m,6:6)-100;
X(m:m,6:6)+10,X(m:m,6:6)-10;
X(m:m,6:6)+1,X(m:m,6:6)-1]'
sumtmat=[...
999 000 000 000 000 000 000 000 000 000 000 000 000 000 000 
899 000 000 000 000 000 000 000 000 000 000 000 000 000 000 
799 889	000 000 000 000 000 000 000 000 000 000 000 000 000 
699 789 888 000 000 000 000 000 000 000 000 000 000 000 000 
599 689 788 779 000 000 000 000 000 000 000 000 000 000 000 
499 589 679 688 778 000 000 000 000 000 000 000 000 000 000 
399 489 579 588 669 678 777 000 000 000 000 000 000 000 000 
299 389 479 488 569 578 668 677 000 000 000 000 000 000 000 
199 289 379 388 469 478 559 568 577 667 000 000 000 000 000 
099 189 278 288 369 378 459 468 477 558 567 666 000 000 000 
089 179 188 269 278 359 368 377 449 458 467 557 566 000 000 
079 088 169 178 259 268 277 349 358 367 448 457 466 556 000 
069 078 159 168 177 249 258 267 339 348 357 366 447 456 555 
059 068 077 149 158 167 239 248 257 266 338 347 356 446 455 
049 058 067 139 148 157 166 229 238 247 256 337 346 355 445 
039 048 057 066 129 138 147 156 228 237 246 255 336 345 444 
029 038 047 056 119 128 137 146 155 227 236 245 335 344 000 
019 028 037 046 055 118 127 136 145 226 235 244 334 000 000 
009 018 027 036 045 117 126 135 144 225 234 333 000 000 000 
008 017 026 035 044 116 125 134 224 233 000 000 000 000 000 
007 016 025 034 115 124 133 223 000 000 000 000 000 000 000 
006 015 024 033 114 123 222 000 000 000 000 000 000 000 000 
005 014 023 113 122 000 000 000 000 000 000 000 000 000 000 
004 013 022 112 000 000 000 000 000 000 000 000 000 000 000 
003 012 111 000 000 000 000 000 000 000 000 000 000 000 000 
002 011 000 000 000 000 000 000 000 000 000 000 000 000 000 
001 000 000 000 000 000 000 000 000 000 000 000 000 000 000 
000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 
];
s1=1:9;
aug=[s1';s1';s1'];
aug=10-aug;
aug=[aug;0];
roottmat=[sumtmat aug];
%missroot is the missing root
rootseq=(0:9)';
sumtseq=(0:27)';
missingsum=setdiff(sumtseq,sumtsort);
missingroot=setdiff(rootseq,roottsort);
%find(roottmat(:,16:16)==missingroot);
%sumtmat is layout to have (28 - row index)  as its
%missing sum total index. So if the missing total is 10
% the row index to look up in sumtmat is 18.
%misstsum=1;
%outlist=sumtmat(misstsum:misstsum,:)
%Saott's System
%
%Using the last K rows of the matrix X
XSS=X(a:b,:);
K=5;
SS=XSS(length(XSS)-(K-1):length(XSS),3:5);
SSout=[];
for i=2:1:K-1
%criss cross
row1=[SS(K-i:K-i,1:1),SS(K-1:K-1,2:2),SS(K:K,3:3)];
row2=[SS(K-i:K-i,3:3),SS(K-1:K-1,2:2),SS(K:K,1:1)];
%do the columns
row3=[SS(K-i:K-i,1:1);SS(K-1:K-1,1:1);SS(K:K,1:1)]';
row4=[SS(K-i:K-i,2:2);SS(K-1:K-1,2:2);SS(K:K,2:2)]';
row5=[SS(K-i:K-i,3:3);SS(K-1:K-1,3:3);SS(K:K,3:3)]';
%do the corners
row6=[SS(K-i:K-i,2:2),SS(K-i:K-i,3:3),SS(K-1:K-1,3:3)];
row7=[SS(K-1:K-1,3:3),SS(K:K,3:3),SS(K:K,2:2)];
row8=[SS(K:K,2:2),SS(K:K,1:1),SS(K-1:K-1,2:2)];
row9=[SS(K-1:K-1,2:2),SS(K-i:K-i,1:1),SS(K-i:K-i,2:2)];
SSout=[row1;row2;row3;row4;row5;row6;row7;row8;row9;SSout];
end
%Calculating sum total and root sum of the Saott numbers
sumSS= sum(SSout(:,1:3),2);
cumsumSS=cumsum(sumSS);
numstr=num2str(sumSS);
ind=find(numstr(:,1)==' ');
numstr(ind,1)='0';
ind=find(numstr(:,2)==' ');
numstr(ind,2)='0';
%
rootSS=str2num(numstr(:,1))+str2num(numstr(:,2));
ind=find(rootSS==10);
roott(ind)=1;
SSout=[SSout 100*SSout(:,1:1)+10*SSout(:,2:2)+SSout(:,3:3) rootSS];
%end
stringsum=num2str(sumtmat,'%0.3d %');
sumttab=tabulate(sumt);
roottab=tabulate(roott);
sumtdist=[1 1 2 3 4 5 7 8 10 12 13 14 15 15 15 15 14 13 12 10 8 7 5 4 3 2 1 1]';
sumtexpected=(sumtdist./sum(sumtdist))*sum([0;sumttab(:,2:2);0;0;0]);
chisquarestatistic=sum(([0;sumttab(:,2:2);0;0;0]-sumtexpected).^2./sumtexpected);
sumtpvalue=1-chi2cdf(chisquarestatistic,length(sumttab)-1);
roottdist=[1 24 24 25 24 24 25 24 24 25]';
roottexpected=(roottdist./sum(roottdist))*sum([0;roottab(:,2:2)]);
chisquarestatistic=sum(([0;roottab(:,2:2)]-roottexpected).^2./roottexpected);
roottpvalue=1-chi2cdf(chisquarestatistic,length(roottab)-1);
% tabulate(roott(1:length(roott)-2:length(sumt)));
tabulate(roott);
% tabulate(sumt(length(sumt)-18:length(sumt)));
