######################################################
#Galaktikate filamendi telje ümber paiknemise analüüs.

## Andmed on kättesaadavad järgnevatest kohtadest:
## Andmestike sees esimeses reas võib leiduda märk #, mille peab eemaldama

## 1
# andmestik a3 (andmestiku nimi on sdss_fcat_r05fin_gals.txt) on kättesaadav veebilehel:
# http://www.aai.ee/~elmo/sdss_fcat_r05fin_gals.txt.zip

## 2
# andmestikud a2 (andmestiku nimi on cat_table_a2.txt) ja a33 (andmestiku nimi on cat_table_a3.txt",) on kättesaadavd veebilehelt:
# http://www.aai.ee/~elmo/sdss-filaments/
# Ascii tables

## 3
# andmestik a34 (andmestiku nimi on dr8_18mar11_gal_v2.txt) on kättesaadav veebilehelt:
# http://www.aai.ee/~elmo/dr8gal.zip

## 1. Andmete filtreerimine, sidumine ja konstrueerimine

a2=read.table(".../cat_table_a2.txt",na.strings=" ",header=T,sep="")
dim(a2)
	#275599     14
	#Galaktikate koordinaadid, kaugus otspunktist:
a3=read.table(".../sdss_fcat_r05fin_gals.txt",na.strings=" ",header=T,sep="")
dim(a3)
	#576493     26
a33=read.table(".../cat_table_a3.txt",na.strings=" ",header=T,sep="")
	#Galaktikaparvede jaoks vajalik tunnus rank
a34=read.table(".../dr8_18mar11_gal_v2.txt",na.strings=" ",header=T,sep="")

length(unique(a2$id))
max(a2$npts)
dist_fil_end=a3[,c(19)]
koord=a3[,c(1,2,3)]
rank=a34$rank
a3=cbind(a33,dist_fil_end,koord,rank)

## 1.1 Vahemikku kuuluvate filamentide leidmine:

aid=a2$id[which(a2$dist>=100 & a2$dist<=250)]
aidd=unique(aid)
a2d=subset(a2, id %in% aidd)
length(unique(a2d$id))
#8568
dim(a2d)
#153960  14
a3u=subset(a3, fil_id!=0)

## 1.2 Loome uue tunnuse, mis annab filamendi siseselt puntkidele id-d

filpt=matrix(0,153960,1)
z=matrix(0,15421,1)
for(i in (1:153960)){
	z[a2d$id[i],1]=z[a2d$id[i],1]+1
	filpt[i,1]=z[a2d$id[i],1]
}
head(filpt)
filpt[1:15]

## 1.3 Järjestan iga id siseselt vastavalt kaugusele paremini ümber

a2b=cbind(a2d,filpt)
a2b[1:100,c(2,3,15)]
a2c=a2b[order(a2b$id), ]

## 1.4 Valime välja andmestikust a3u need galaktikad, mis asuvad meid huvitavas piirkonnas

filpn=a2c$idpts
	#Vektor, mis sisaldab filamentide ID koode, mis asuvad meile huvipakkuvas piirkonnas
a2cid=unique(a2c$id)
a3c=subset(a3u, fil_id %in% a2cid) 
dim(a3c)
#198861 23
a3c=a3c[order(a3c$fil_id,a3c$fil_idpts),]
	#Tükeldan andmestiku vastavalt sellele, millise punktikese ümber üldse on galaktikaid
filidpts=unique(a3c$fil_idpts)
	#võtan need filamendi punktid, mille juurde kuulub galaktika
a2d=subset(a2c, idpts %in% filidpts)
#80724    15

## 1.5 Lisame andmestikku tunnuse filpt

	#Kõikide nende filamendipunktide id, mille juurde kuulub galaktika
galidpts=a3c$fil_idpts 
x=as.data.frame(table(galidpts))
kordused=x$Freq
length(kordused)
length(a2d$filpt)
	#Nüüd kordame vektoris filpt olevaid arve vastavalt arvudele vektoris
filpt=a2d$filpt

punkt=rep(a2d$filpt,times=kordused)

a3d=cbind(a3c,punkt)
dim(a3d)

	
## 1.6 Valime välja npts väärtused, mille korral fil_id ja id langevad kokku

npts=a2d$npts[match(a3d$fil_id, a2d$id)]
a3e=cbind(a3d,npts)
length(a3e$fil_id)
#198861

## 1.7.1 Piiritlen ülevalt toru, galaktika võib olla teljest kuni 0.5 Mpc kaugusel.

a3f=subset(a3e,a3e$fil_dist<=0.5) #66785
dim(a3f)
#66785   25

## 1.7.2 Piiritlen toru lõpp-punktidest Galaktikad, mis on viimasest punktist maksimaalselt 0.25 Mpc kaugusel

a3f=subset(a3f,a3f$dist_fil_end<=0.25)
dim(a3f)
#64805    23
length(unique(a3f$fil_id))
#8545 nii palju erinevaid filamente

		
## 1.7 Võtame galaktikad, mis asuvad vähmalt 30 Mpc filamentidel valitud ümbruses

a3Mpc30=subset(a3f,a3f$npts>=60)
dim(a3Mpc30)
	#5193   23
max(a3Mpc30$npts)
filid30=unique(a3Mpc30$fil_id)
length(unique(a3Mpc30$fil_id))
	#[1] 165 Meil on 165 filamenti vaatluse all. Muidu on meil filamente 8568

## 1.8 Loon uue andmestiku, kus on filamendid jupitatud, esimesest ja viimasest punktist 

## 1.8.1 Pööran iga filamendi siseselt järjestuse ümber

max(a3Mpc30$npts)
a3Mpc30_2=a3Mpc30[order(a3Mpc30$fil_id,-a3Mpc30$fil_idpts),]
head(a3Mpc30_2)
uuspunkt=a3Mpc30_2$npts-a3Mpc30_2$punkt+1
a3Mpc30_2=cbind(a3Mpc30_2,uuspunkt)

## 1.8.2 Löön andmestikud kaheks: vaatan ühelt ja teiselt poolt filamenti
	#andmed1=uuspunkt<=60 ja andmed2=punkt<=60

andmed1=a3Mpc30[which(a3Mpc30$punkt<=60),]
head(andmed1)
dim(andmed1)
	#4403 20
andmed2=a3Mpc30_2[which(a3Mpc30_2$uuspunkt<=60),]
head(andmed2)
dim(andmed2)
	#4378   21
	#paneme punktide andmed kokku:
galfilamendil=c(andmed1$punkt,andmed2$uuspunkt)
length(galfilamendil)
# 8781

#Filamendi teljest 0.5 Mpc toru sees asuvate galaktikate jaotus mööda teljepunkte
hist(galfilamendil,breaks=c(seq(from=0,to=60,by=1)),cex.axis=1.3,cex.lab=1.3,main="",xlab="Filamendi teljepunkt",ylab="Galaktikate arv")
	
	#Kui palju asetseb täpselt:
galfiltabel=as.data.frame(table(galfilamendil))
galfiltabel_2=galfiltabel[order(-galfiltabel$Freq),]

## 2. Galaktikate projekteerimine filamendi teljele 

## 2.1. Galaktika projektsiooni kaugus lähimast punktist

andmed3=subset(andmed1,andmed1$dist_fil_end==0)
dim(andmed3)
#4372  30
andmed4=subset(andmed2,andmed2$dist_fil_end==0)
dim(andmed4)
#4355  31
	#Galaktika kaugus filamendi punktist
a2d_1=a2d[match(andmed3$fil_idpts, a2d$idpts),]
kaugus1=matrix(0,nrow(andmed3),1)
for (i in 1:(nrow(andmed3))){
	kaugus1[i]=sqrt((a2d_1$x[i]-andmed3$x[i])**2+(a2d_1$y[i]-andmed3$y[i])**2+(a2d_1$z[i]-andmed3$z[i])**2)
	}
galtelg1=sqrt(kaugus1**2-andmed3$fil_dist**2)
	#asendan kõik Nan väärtused nullidega.
sum(is.na(as.matrix(galtelg1)))
galtelg1=replace(galtelg1,is.na(galtelg1),0) 
andmed3=cbind(andmed3,galtelg1)

a2d_2=a2d[match(andmed4$fil_idpts, a2d$idpts),]
kaugus2=matrix(0,nrow(andmed4),1)
for (i in 1:(nrow(andmed4))){
	kaugus2[i]=sqrt((a2d_2$x[i]-andmed4$x[i])**2+(a2d_2$y[i]-andmed4$y[i])**2+(a2d_2$z[i]-andmed4$z[i])**2)
	}
galtelg2=sqrt(kaugus2**2-andmed4$fil_dist**2)
sum(is.na(as.matrix(galtelg2)))
galtelg2=replace(galtelg2,is.na(galtelg2),0) 
andmed4=cbind(andmed4,galtelg2)
	#Tunnused galtelg1 ja galtelg2 on nüüd galaktika kaugus mööda telge lähimast puntkist.

## 2.2 Galaktika kaugus lähima punkti naabritest  
	#Loon uued andmeveerud: parem, vasem ---need on vaadeldava punkti paremal ja
	#vasakul pool asuva punkti koordinaadid

#ANDMED 3

parem1=matrix(0,nrow=nrow(andmed3),ncol=3)
vasem1=matrix(0,nrow(andmed3),3)
parem1[,1]=andmed3$fil_id
vasem1[,1]=andmed3$fil_id
parem1[,3]=andmed3$npts
vasem1[,3]=andmed3$npts
for(i in (1:(nrow(andmed3)))){
	parem1[i,2]=andmed3$punkt[i]+1
	vasem1[i,2]=andmed3$punkt[i]-1
	}
for(i in (1:nrow(parem1))){
	if (parem1[i,2]>parem1[i,3]){
		parem1[i,2]=parem1[i,3]
		}
	}
for(i in (1:nrow(vasem1))){
	if(vasem1[i,2]==0){
		vasem1[i,2]=1
		}
	}
parem1=as.data.frame(parem1)
colnames(parem1)=c("id","filpt","npts")
vasem1=as.data.frame(vasem1)
colnames(vasem1)=c("id","filpt","npts")

install.packages("sqldf")
library(sqldf)
a2b_p1=sqldf("SELECT * FROM a2b,parem1 WHERE a2b.id==parem1.id and a2b.filpt==parem1.filpt")
a2b_v1=sqldf("SELECT * FROM a2b,vasem1 WHERE a2b.id==vasem1.id and a2b.filpt==vasem1.filpt")
sum(parem1$filpt-a2b_p1$filpt)

	#Leiame galaktikate kaugused vastavatest filamendi punktidest
	#Neid kauguseid võrreldes saame öelda, kummale poole lähimat punkti see galaktika jääb!
kaugus1_v=matrix(0,nrow(andmed3),1)
for (i in (1:nrow(andmed3))){
	kaugus1_v[i]=sqrt((andmed3$x[i]-a2b_v1$x[i])**2+(andmed3$y[i]-a2b_v1$y[i])**2+(andmed3$z[i]-a2b_v1$z[i])**2)
	}
kaugus1_p=matrix(0,nrow(andmed3),1)
for (i in (1:nrow(andmed3))){
	kaugus1_p[i]=sqrt((andmed3$x[i]-a2b_p1$x[i])**2+(andmed3$y[i]-a2b_p1$y[i])**2+(andmed3$z[i]-a2b_p1$z[i])**2)
	}

#ANDMED 4

parem2=matrix(0,nrow=nrow(andmed4),ncol=3)
vasem2=matrix(0,nrow(andmed4),3)
parem2[,1]=andmed4$fil_id
vasem2[,1]=andmed4$fil_id
parem2[,3]=andmed4$npts
vasem2[,3]=andmed4$npts
for(i in (1:(nrow(andmed4)))){
	parem2[i,2]=andmed4$punkt[i]+1
	vasem2[i,2]=andmed4$punkt[i]-1
	}
for(i in (1:nrow(parem2))){
	if (parem2[i,2]>parem2[i,3]){
		parem2[i,2]=parem2[i,3]
		}
	}
for(i in (1:nrow(vasem2))){
	if(vasem2[i,2]==0){
		vasem2[i,2]=1
		}
	}
parem2=as.data.frame(parem2)
colnames(parem2)=c("id","filpt","npts")
vasem2=as.data.frame(vasem2)
colnames(vasem2)=c("id","filpt","npts")

	#Kuna me oleme vaadeldava andmestiku ümber pööranud, vaatama filamenti lõpust algusesse
a2b_p2=sqldf("SELECT * FROM a2b,parem2 WHERE a2b.id==parem2.id and a2b.filpt==parem2.filpt")
a2b_p2=a2b_p2[order(a2b_p2$id,-a2b_p2$filpt),]
a2b_v2=sqldf("SELECT * FROM a2b,vasem2 WHERE a2b.id==vasem2.id and a2b.filpt==vasem2.filpt")
a2b_v2=a2b_v2[order(a2b_v2$id,-a2b_v2$filpt),]

	#Leiame galaktikate kaugused vastavatest filamendi punktidest
kaugus2_v=matrix(0,nrow(andmed4),1)
for (i in (1:nrow(andmed4))){
	kaugus2_v[i]=sqrt((andmed4$x[i]-a2b_v2$x[i])**2+(andmed4$y[i]-a2b_v2$y[i])**2+(andmed4$z[i]-a2b_v2$z[i])**2)
	}
kaugus2_p=matrix(0,nrow(andmed4),1)
for (i in (1:nrow(andmed4))){
	kaugus2_p[i]=sqrt((andmed4$x[i]-a2b_p2$x[i])**2+(andmed4$y[i]-a2b_p2$y[i])**2+(andmed4$z[i]-a2b_p2$z[i])**2)
	}

## 2.3 Galaktika asukoht lähima filamendi punkti suhtes: tunnus KOOD

	#ANDMED 3 korral:
	#Kui ta on vasempoolsele punktile lähemal siis tuleb koodiks -1 
	#Kui ta on parempoolsele punktile lähemal tuleb koodiks +1
kood1=matrix(0,nrow(andmed3),1)
for(i in (1:nrow(andmed3))){
	if (kaugus1_v[i]<kaugus1_p[i]){
		kood1[i]=-1
		}
	else if (kaugus1_v[i]>kaugus1_p[i]){
		kood1[i]=1
		}
		else kood1[i]=NA
	}
sum(is.na(as.matrix(kood1)))
#0 mitte ükski galaktika ei asu mõlemast punktist sama kaugel
kood2=matrix(0,nrow(andmed4),1)
for(i in (1:nrow(andmed4))){
	if (kaugus2_v[i]<kaugus2_p[i]){
		kood2[i]=-1
		}
	else if (kaugus2_v[i]>kaugus2_p[i]){
		kood2[i]=1
		}
		else kood2[i]=NA
	}
sum(is.na(as.matrix(kood2)))
#0 mitte ükski galaktika ei asu mõlemast punktist sama kaugel
andmed3=cbind(andmed3,kood1)
andmed4=cbind(andmed4,kood2)

## 2.4 Galaktika teine lähim punkt: selle järgi saab õigema järjestuse

kumbpool=matrix(0,nrow(andmed3),1)
for(i in (1:nrow(andmed3))){
	if (kaugus1_v[i]<kaugus1_p[i]){
		kumbpool[i]=a2b_v1$filpt[i]
		}
	else if (kaugus1_v[i]>kaugus1_p[i]){
		kumbpool[i]=a2b_p1$filpt[i]
		}
		else kumbpool[i]=NA
	}
andmed3=cbind(andmed3,kumbpool)
	#Kui galaktika asub esimese või viimase punkti juures
for(i in (1:nrow(andmed3))){
	if(andmed3$kumbpool[i]==andmed3$punkt[i]){
		andmed3$kood1[i]=0
		}
}
	#Sorteerime lähima teise punkti alusel
andmed3=andmed3[order(andmed3$fil_id,andmed3$kumbpool),]

kumbpool_2=matrix(0,nrow(andmed4),1)
for(i in (1:nrow(andmed4))){
	if (kaugus2_v[i]<kaugus2_p[i]){
		kumbpool_2[i]=a2b_v2$filpt[i]
		}
	else if (kaugus2_v[i]>kaugus2_p[i]){
		kumbpool_2[i]=a2b_p2$filpt[i]
		}
		else kumbpool_2[i]=NA
	}
kumbpool=andmed4$npts-kumbpool_2+1

andmed4=cbind(andmed4,kumbpool)
for(i in (1:nrow(andmed4))){
	if(andmed4$kumbpool[i]==andmed4$uuspunkt[i]){
		andmed4$kood2[i]=0
	}
}
andmed4=andmed4[order(andmed4$fil_id,andmed4$kumbpool),]	

## 3. Galaktikate kaugused otspunktist: tunnus GALKAUGUSED

galkaugused=matrix(0,nrow(andmed3),1)
for (i in (1:nrow(andmed3))){
	if (andmed3$punkt[i]==1){
		galkaugused[i]=andmed3$galtelg1[i]
	}
	else{
	galkaugused[i]=(andmed3$punkt[i]-1)*0.5+andmed3$kood1[i]*andmed3$galtelg1[i]
	}
}
andmed3=cbind(andmed3,galkaugused)

galkaugused=matrix(0,nrow(andmed4),1)
for (i in (1:nrow(andmed4))){
	if (andmed4$uuspunkt[i]==1){
		galkaugused[i]=andmed4$galtelg2[i]
	}
	else{
	galkaugused[i]=(andmed4$uuspunkt[i]-1)*0.5+andmed4$kood2[i]*andmed4$galtelg2[i]
	}
}
andmed4=cbind(andmed4,galkaugused)

	#KONTROLL
max(andmed3$galkaugused)
#[1] 29.76334
min(andmed3$galkaugused)
#[1] 0.001786505
max(andmed4$galkaugused)
#[1] 29.75707

## 4. Naabergalaktikate vahelised kaugused filamendi kaupa
## Lahutan kahe galaktika kauguse otspunktist selleks, et teada saada nende omavahelist kaugust

galgalkaug=matrix(0,nrow(andmed4),1)
for (i in (1:nrow(andmed4))){
	galgalkaug[i]=abs(andmed4$galkaugused[i]-andmed4$galkaugused[i+1])
	}
andmed4=cbind(andmed4,galgalkaug)
galgalkaug1=matrix(0,nrow(andmed3),1)
for (i in (1:nrow(andmed3))){
	galgalkaug1[i]=abs(andmed3$galkaugused[i]-andmed3$galkaugused[i+1])
	}
andmed3=cbind(andmed3,galgalkaug1)
	#Valed elemendid veergudes galgalkaug1 on need, mis näitavad filamendi viimase elemendi kaugust.
	#Sest ta arvutab erinevuse esimese filamendi viimasest hoopis teise filamendi esimese punkti kauguste vahel


## 5. Teljele projekteeritud naabergalaktikate kauguste vahelise korrelatsiooni leidmine

## 5.1. Galaktika kaugus paremast ja vaskust naabrist

#andmed1#
ngalgalkaug1=matrix(0,length(andmed3$galgalkaug1)-1,1)
	for(i in (1:length(andmed3$galgalkaug1)-1)){
		ngalgalkaug1[i]=andmed3$galgalkaug1[i+1]
	}
ngalgalkaug1[4372]=0
kaugusednihkes1=cbind(andmed3$fil_id,andmed3$punkt,andmed3$galkaugused,andmed3$galgalkaug1,ngalgalkaug1,andmed3$rank,andmed3$nrich)
colnames(kaugusednihkes1)=c("fil_id","punkt","galkaugused","galgalkaug1","ngalgalkaug1","rank","nrich")
kaugusednihkes1=kaugusednihkes1[-4372,]
kaugusednihkes1=as.data.frame(kaugusednihkes1)

#andmed2#
ngalgalkaug2=matrix(0,length(andmed4$galgalkaug)-1,1)
	for(i in (1:length(andmed4$galgalkaug)-1)){
		ngalgalkaug2[i]=andmed4$galgalkaug[i+1]
	}
ngalgalkaug2[4355]=0
kaugusednihkes2=cbind(andmed4$fil_id,andmed4$uuspunkt,andmed4$galkaugused,andmed4$galgalkaug,ngalgalkaug2,andmed4$rank,andmed4$nrich)
colnames(kaugusednihkes2)=c("fil_id","punkt","galkaugused","galgalkaug","ngalgalkaug2","rank","nrich")
kaugusednihkes2=kaugusednihkes2[-4355,]
kaugusednihkes2=as.data.frame(kaugusednihkes2)

## 5.2. Ebasobivate paaride eemaldamine (nihutamisel tekkiv viga kahe filamendi vahel)

	#andmed3#
kn1=kaugusednihkes1[!duplicated(kaugusednihkes1$fil_id, fromLast=TRUE),]
kaugusednihkes1=kaugusednihkes1[setdiff(rownames(kaugusednihkes1),rownames(kn1)),]

	#andmed4#
kn2=kaugusednihkes2[!duplicated(kaugusednihkes2$fil_id, fromLast=TRUE),]
kaugusednihkes2=kaugusednihkes2[setdiff(rownames(kaugusednihkes2),rownames(kn2)),]

## 5.3 Andmestike ühendamine (paremalt ja vasakult loetud filamendid)

## 5.3.1 Selleks, et need andmestikud kokku panna, peame looma filamentidele uued nimed 
## andmestikus kaugusednihkes2

a=length(unique(kaugusednihkes1$fil_id))+length(unique(kaugusednihkes2$fil_id))
	#a=330
b=unique(kaugusednihkes1$fil_id) #maksimaalne on 7141
	#genereerime arvud 7200st 7200+165
c=c(7200:7364)
x=as.data.frame(table(kaugusednihkes2$fil_id))
kordused=x$Freq
length(kordused)
fil_id=rep(c,times=kordused)
	#Võtan andmed 2 abil tehtud andmestikust fil_id tunnuse välja, sest see kattub andmed 1 omadega
kaugusednihkes2_2=kaugusednihkes2[,-1]
	#Asendan uue fil_id koodiga
kaugusednihkes2_2=cbind(fil_id,kaugusednihkes2_2)


## 5.3.2 Paneme andmestikud kokku

colnames(kaugusednihkes1)=c("fil_id","punkt","galkaugused","galgalkaug","ngalgalkaug","rank","nrich")
colnames(kaugusednihkes2_2)=c("fil_id","punkt","galkaugused","galgalkaug","ngalgalkaug","rank","nrich")

kaugusednihkes=rbind(kaugusednihkes1,kaugusednihkes2_2)

## 5.4 Andmete puhastamine 
## (Iga filamendi viimases reas on väär infot viimase galaktika paiknemise kohta)

kn5=kaugusednihkes[!duplicated(kaugusednihkes$fil_id, fromLast=TRUE),]
kaugusednihkes_m=kaugusednihkes[setdiff(rownames(kaugusednihkes),rownames(kn5)),]
	#Kui mitmel galaktikal on mõlemad naabrid
dim(kaugusednihkes_m)
	#8065 4

## 5.5 Korrelatsioonide leidmine 

install.packages("data.table")
library(data.table)
kaugusednihkes=as.data.table(kaugusednihkes_m)
setkey(kaugusednihkes,fil_id)
korrelatsioonid=kaugusednihkes[,list(corr = cor(galgalkaug,ngalgalkaug,method = 'pearson')), by =fil_id]
uusid=c(1:330)
korrelatsioonid=cbind(korrelatsioonid,uusid)
	#ega kuskil probleeme ei tekkinud:
sum(is.na(as.matrix(korrelatsioonid$corr)))
	#0 ei ole
#Filamendile projekteeritud naabergalaktikate vaheliste kauguste vaheline korrelatsioon
plot(korrelatsioonid$uusid,korrelatsioonid$corr,cex.axis=1.3,cex.lab=1.5,main="",xlab="Filament",ylab="Korrelatsioon")

#Joonis naabergalaktikate kaugustest vaadeldaval filamendil
plot(kaugusednihkes$galgalkaug,kaugusednihkes$ngalgalkaug,cex.axis=1.3,cex.lab=1.5,main="",xlab="Galaktikate gal_(i) ja gal_(i+1) vaheline kaugus",ylab="Galaktikate gal_(i+1) ja gal_(i+2) vaheline kaugus")

min(korrelatsioonid$corr)
[1] -0.5372933
max(korrelatsioonid$corr)
[1] 0.3726377
mean(korrelatsioonid$corr)
#-0.1522318

## 6. Paaris korrelatsioonifunktsiooni leidmine

galaktikaots_3=cbind(andmed3$fil_id,andmed3$galkaugused)
	#uued koodid filamentidele
a_2=length(unique(andmed4$fil_id))
	#genereerime arvud 7200st 7200+165
c_2=c(7200:7364)
x_2=as.data.frame(table(andmed4$fil_id))
kordused_2=x_2$Freq
length(kordused_2)
fil_id=rep(c_2,times=kordused_2)
galaktikaots_4=cbind(fil_id,andmed4$galkaugused)
galaktikaots=rbind(galaktikaots_3,galaktikaots_4)
colnames(galaktikaots)=c("fil_id","galkaugused")
galaktikaots=as.data.frame(galaktikaots)
dim(galaktikaots)
#8726   2

## 6.1 Leian iga filamendi jaoks kõikide galaktikate kaugused teistest galaktikatest

suurandmestik=as.matrix(0,nrow=1,ncol=1)
	suurvektor=c(0)
	indikaatorid=unique(galaktikaots$fil_id)				
	for (i in (1:length(indikaatorid))){
		B1=subset(galaktikaots,fil_id %in% indikaatorid[i])
		b1=B1$galkaugused
			andmed=as.matrix(c(rep(b1,length(b1))))
			andmed=matrix(andmed,nrow=length(b1),ncol=length(b1))
			lahutaja=as.matrix(c(rep(b1,each=length(b1))))
			lahutaja=matrix(lahutaja,nrow=length(b1),ncol=length(b1))
			suurandmestik=abs(andmed-lahutaja)
			diag(suurandmestik)=NA
			#eemaldan diagonaalielemendid:
			A1=matrix(0,nrow=nrow(suurandmestik)-1,ncol=1)
			for (j in (1:ncol(suurandmestik))){
				veerg=suurandmestik[,j]
				a1=as.matrix(veerg[which(!is.na(veerg))])
				A1=cbind(A1,a1)}
			A1=A1[,-1]
			vektor1=as.vector(A1)
			#eemaldan kõik galaktika kaugused iseendast
			#vektor1=vektor1[!(vektor1==0)]
			suurvektor=c(suurvektor,vektor1)
	}
	#eemaldan selle 0 algusest
suurvektor=suurvektor[-1]

#Galaktikate kaugused teineteisest 330 erineva filamendi pealt
hist(suurvektor,cex.lab=1.3,cex.axis=1.5,breaks=60,main="",xlab="Kaugus r (h^{-1} Mpc)",ylab="Galaktikate paaride arv")

	#Vaatame, kus on täpselt maksimumid:
as.data.frame(table(cut(suurvektor,breaks=c(seq(from=0,to=30,by=0.5)))))

	#Andmestikule suurvektor uued ID-d iga filamendi eraldamiseks
x_3=as.data.frame(table(galaktikaots$fil_id))
korda=x_3$Freq*(x_3$Freq-1)
uusid=c(1:330)
fil_id=rep(uusid,times=korda)
	#Kontroll
length(suurvektor)
	#250346
length(fil_id)
	#250346
galaktjaotus=cbind(fil_id,suurvektor)
galaktjaotus=as.data.frame(galaktjaotus)

## 6.2 Leiame filamentide galaktikate vaheliste kauguste jaotust
## Moodustame maatriksi DD:60 x 330

	#Leian DD iga filamendi siseselt:
DD=matrix(0,nrow=60,ncol=330)
for(i in (1:length(uusid))){
	B1=subset(galaktjaotus,fil_id %in% uusid[i])
	b1=B1$suurvektor
	DD[,i]=as.data.frame(table(cut(b1,breaks=c(seq(from=0,to=30,by=0.5)))))$Freq
}

## 6.3 Ühtlase jaotusega juhuslike suuruste genereerimine
## Genereerin igale filamendile vastavusse 10 korda rohkem ühtlase jaotusega asetsevaid punkte lõigul [0,30]

x_3=as.data.frame(table(galaktikaots$fil_id))
kordused_3=x_3$Freq
sum(kordused_3)
#8727
y=c(0)
for(i in (1:length(kordused_3))){
	x=runif(kordused_3[i]*10,0,30)
	y=c(y,x)}
y=y[-1]
length(y)
#87270

## 6.4 Leiame galaktikate kaugused juhuslikest punktidest	

suurandmestik_bgal=as.matrix(0,nrow=1,ncol=1)
galaktuhtlasest=c(0)				
x=array(0,1)						
for(i in (1:length(kordused_3))){			
	B1=subset(galaktikaots,fil_id %in% indikaatorid[i])	
	bgal=B1$galkaugused							
	x=kordused_3[i]*10
	lahutaja_bgal=as.matrix(c(rep(bgal,each=x)))
	lahutaja_bgal=matrix(lahutaja_bgal,nrow=x,ncol=kordused_3[i])
	z=y[(sum(kordused_3[0:(i-1)]*10)+1):(sum(kordused_3[0:(i)])*10)]
	andmedj=as.matrix(c(rep(z,kordused_3[i])))
	andmedj=matrix(andmedj,nrow=length(z),ncol=kordused_3[i])
	suurandmestik_bgal=abs(andmedj-lahutaja_bgal)
	vektor1_bgal=as.vector(suurandmestik_bgal)
	galaktuhtlasest=c(galaktuhtlasest,vektor1_bgal)
}
galaktuhtlasest=galaktuhtlasest[-1]
#Galaktika kaugused juhuslikult paiknevatest punktidest 30 ühikulisel sirgel
hist(galaktuhtlasest,cex.lab=1.4,cex.axis=1.5,breaks=60,main="",xlab="Kaugus r (h^{-1} Mpc)",ylab="Galaktikate ja punktide paaride arv")

sum(kordused_3*10*kordused_3)
#[1] 2590730 
length(galaktuhtlasest)
#[1] 2590730

## 6.5 Moodustame maatriksi DR:60 x 330, mis sisaldab iga filamendi galaktikate ja juhuslike punktide 
## vaheliste kauguste jaotust.

korda_uhtl=kordused_3*kordused_3*10
fil_id=rep(uusid,times=korda_uhtl)
galuhtljaotus=as.data.frame(cbind(fil_id,galaktuhtlasest))
DR=matrix(0,nrow=60,ncol=330)
for(i in (1:length(uusid))){
	B1=subset(galuhtljaotus,fil_id %in% uusid[i])
	b1=B1$galaktuhtlasest
	DR[,i]=as.data.frame(table(cut(b1,breaks=c(seq(from=0,to=30,by=0.5)))))$Freq
}

## 6.6 Leiame DD ja DR suhte iga filamendi jaoks

KSI=matrix(0,nrow=60,ncol=330)
for(i in (1:ncol(DD))){
	KSI[,i]=(DD[,i]/DR[,i])*10
}

## 6.8 Leiame paaris-korrelatsioonifunktsiooni

ksi=matrix(0,nrow=60,ncol=1)
for(i in (1:nrow(KSI))){
	ksi[i,1]=sum(KSI[i,])/330
}
	#Vaatame, kust algavad NAN:
ksi
	#algavad 50.vahest
ksi=ksi[1:50]
telgx=as.data.frame(table(cut(c(0:30),breaks=c(seq(from=0,to=30,by=0.5)))))$Var
telgx=telgx[1:50] #sest 45 elemendil on väärtused

	#Vaatame järjestust:
ksitelg=cbind(telgx,ksi)
ksitelg=ksitelg[order(-ksi),]

	#Joonistame: Paaris-korrelatsioonifunktsioon
plot(telgx,ksi,cex.lab=1.3,col="white",ylab="g(r, r+dr)",xlim=c(0,51),main="",xlab="Galaktikate vaheline kauguse vahemik r, r+dr")
lines(ksi,col="green")
lines(x=c(0,50),y=c(1,1),lty=b,col="sky blue")
min(ksi)
#[1] 0.6661287
mean(ksi)
#[1] 0.9490232
cbind(telgx,ksi)

## 7. Jackknife veahinnangu leidmine

	#1#
#jätan välja i-nda filamendi
jackksi=matrix(0,nrow=60,ncol=330)
JackKSI=matrix(0,nrow=60,ncol=329)
for(i in (1:ncol(DD))){
	DD_sub=DD[,-i]
	DR_sub=DR[,-i]
	for(j in (1:(ncol(DD_sub)))){
		JackKSI[,j]=(DD_sub[,j]/DR_sub[,j])*10
	}
	for(k in (1:(nrow(JackKSI)))){
		jackksi[k,i]=sum(JackKSI[k,])/329
	}
}
#jackksi algab alates 50ndast reast
jackksi=jackksi[1:50,]

	#2#
pseudoksi=matrix(0,50,330)
ksi=as.matrix(ksi,nrow=50,ncol=1)
for(i in (1:ncol(jackksi))){
	pseudoksi[,i]=330*ksi-329*jackksi[,i]
}

	#3#
keskjack=matrix(0,nrow=50,ncol=1)
for(i in (1:nrow(jackksi))){
	keskjack[i,]=sum(pseudoksi[i,])/330
}

	#4#
#suurus s**2
sruut=matrix(0,nrow=50,ncol=1)
for(i in (1:nrow(pseudoksi))){
	sruut[i,]=1/329*((sum(pseudoksi[i,]**2))-1/330*(sum(pseudoksi[i,]))**2)
}

	#5#
#s_tärn**2=s**2/k
starn=sqrt(sruut)/sqrt(330)

	#6#
#Usaldusvahemikud
#ylem=ksi+qt(0.025,330-1)*starn
#alum=ksi-qt(0.025,330-1)*starn
ylem=keskjack+qt(0.025,330-1)*starn
alum=keskjack-qt(0.025,330-1)*starn

#Joonistame:
plot(telgx,ksi,cex.axis=1.3,cex.lab=1.3,col="white",ylab="g(r, r+dr)",xlim=c(0,51),main="",xlab="Galaktikate vaheline kauguse vahemik (r, r+dr] h^{-1} Mpc")
lines(ksi,col="green")
lines(ylem,col="black")
lines(alum,col="black")
legend(22,2,box.col="white",border="white",c("g(r,r+dr)","Jackknife vahemikhinnangud"), cex=1.3, fill=c("green","black"))
lines(x=c(0,50),y=c(1,1),lty=b,col="sky blue")

##########################################################
#Galaktikaparvede filamendi telje ümber paiknemise analüüs.


##1. Andmete filtreerimine, sidumine ja konstrueerimine
	# Andmestikes andmed1, andmed2 on galaktikad, mis asuvad vaadeldavast filamendi teljest
	# 0.5 Mpc kaugusel ja lõpp-punktidest on galaktika maks 0.25 Mpc kaugusel. 
	# Andmestikud luuakse koodis Galaktikad.R.

	#Nüüd vaatame nende andmestike korral valikusse jäänud galaktikaparvi
	#rank==1 annab meile peagalaktikad
galgrupp1=andmed1[which(andmed1$rank==1 & andmed1$nrich>1),]
dim(galgrupp1)
# 969  25

galgrupp2=andmed2[which(andmed2$rank==1 & andmed2$nrich>1),]
dim(galgrupp2)
#970  26

parvfilamendil=c(galgrupp1$punkt,galgrupp2$uuspunk)
#Filamendi teljest 0.5 Mpc toru sees asuvate galaktikaparvede jaotus mööda teljepunk
hist(parvfilamendil,breaks=c(seq(from=0,to=60,by=1)),cex.lab=1.3,cex.axis=1.3,cex.main=1.5,main="",xlab="Filamendi teljepunkt",ylab="Galaktikaparvede arv")
temp_parvfil=as.data.frame(table(parvfilamendil))
temp_parvfil=temp_parvfil[order(ohoo$Freq),]

## 2. Galaktikaparvede projekteerimine filamendi teljele 

## 2.1. Galaktikaparvede projektsiooni kaugus lähimast punktist

galgrupp3=subset(galgrupp1,galgrupp1$dist_fil_end==0)
galgrupp4=subset(galgrupp2,galgrupp2$dist_fil_end==0)
dim(galgrupp3)
#[1] 962  25
dim(galgrupp4)
#[1] 969  26

a2d_1=a2d[match(galgrupp3$fil_idpts, a2d$idpts),]
kaugus1=matrix(0,nrow(galgrupp3),1)
for (i in 1:(nrow(galgrupp3))){
	kaugus1[i]=sqrt((a2d_1$x[i]-galgrupp3$x[i])**2+(a2d_1$y[i]-galgrupp3$y[i])**2+(a2d_1$z[i]-galgrupp3$z[i])**2)
}
parvtelg1=sqrt(kaugus1**2-galgrupp3$fil_dist**2)

sum(is.na(as.matrix(parvtelg1)))
parvtelg1=replace(parvtelg1,is.na(parvtelg1),0) 
galgrupp3=cbind(galgrupp3,parvtelg1)

a2d_2=a2d[match(galgrupp4$fil_idpts, a2d$idpts),]
kaugus2=matrix(0,nrow(galgrupp4),1)
for (i in 1:(nrow(galgrupp4))){
	kaugus2[i]=sqrt((a2d_2$x[i]-galgrupp4$x[i])**2+(a2d_2$y[i]-galgrupp4$y[i])**2+(a2d_2$z[i]-galgrupp4$z[i])**2)
}
parvtelg2=sqrt(kaugus2**2-galgrupp4$fil_dist**2)

sum(is.na(as.matrix(parvtelg2)))
parvtelg2=replace(parvtelg2,is.na(parvtelg2),0) 
galgrupp4=cbind(galgrupp4,parvtelg2)

## 2.2 Galaktikaparvede kaugus lähima punkti naabritest  

#GALGRUPP 3

parem1=matrix(0,nrow=nrow(galgrupp3),ncol=3)
vasem1=matrix(0,nrow(galgrupp3),3)
parem1[,1]=galgrupp3$fil_id
vasem1[,1]=galgrupp3$fil_id
parem1[,3]=galgrupp3$npts
vasem1[,3]=galgrupp3$npts
for(i in (1:(nrow(galgrupp3)))){
	parem1[i,2]=galgrupp3$punkt[i]+1
	vasem1[i,2]=galgrupp3$punkt[i]-1
}
for(i in (1:nrow(parem1))){
	if (parem1[i,2]>parem1[i,3]){
		parem1[i,2]=parem1[i,3]
	}
}
for(i in (1:nrow(vasem1))){
	if(vasem1[i,2]==0){
		vasem1[i,2]=1
	}
}
parem1=as.data.frame(parem1)
colnames(parem1)=c("id","filpt","npts")
vasem1=as.data.frame(vasem1)
colnames(vasem1)=c("id","filpt","npts")

#Koodis Galaktikad.R sisse loetud pakett:
#install.packages("sqldf")
#library(sqldf)
a2b_p1=sqldf("SELECT * FROM a2b,parem1 WHERE a2b.id==parem1.id and a2b.filpt==parem1.filpt")
a2b_v1=sqldf("SELECT * FROM a2b,vasem1 WHERE a2b.id==vasem1.id and a2b.filpt==vasem1.filpt")
sum(parem1$filpt-a2b_p1$filpt)

	#Leiame galaktikaparvede kaugused vastavatest filamendi punktidest
kaugus1_v=matrix(0,nrow(galgrupp3),1)
for (i in (1:nrow(galgrupp3))){
	kaugus1_v[i]=sqrt((galgrupp3$x[i]-a2b_v1$x[i])**2+(galgrupp3$y[i]-a2b_v1$y[i])**2+(galgrupp3$z[i]-a2b_v1$z[i])**2)
}
kaugus1_p=matrix(0,nrow(galgrupp3),1)
for (i in (1:nrow(galgrupp3))){
	kaugus1_p[i]=sqrt((galgrupp3$x[i]-a2b_p1$x[i])**2+(galgrupp3$y[i]-a2b_p1$y[i])**2+(galgrupp3$z[i]-a2b_p1$z[i])**2)
}

#GALGRUPP 4

parem2=matrix(0,nrow=nrow(galgrupp4),ncol=3)
vasem2=matrix(0,nrow(galgrupp4),3)
parem2[,1]=galgrupp4$fil_id
vasem2[,1]=galgrupp4$fil_id
parem2[,3]=galgrupp4$npts
vasem2[,3]=galgrupp4$npts
for(i in (1:(nrow(galgrupp4)))){
	parem2[i,2]=galgrupp4$punkt[i]+1
	vasem2[i,2]=galgrupp4$punkt[i]-1
}
for(i in (1:nrow(parem2))){
	if (parem2[i,2]>parem2[i,3]){ #filamendi lõpppunkti korral
		parem2[i,2]=parem2[i,3]
	}
}
for(i in (1:nrow(vasem2))){
	if(vasem2[i,2]==0){ #filamendi alguspunkti korral
		vasem2[i,2]=1
	}
}
parem2=as.data.frame(parem2)
colnames(parem2)=c("id","filpt","npts")
vasem2=as.data.frame(vasem2)
colnames(vasem2)=c("id","filpt","npts")

	#Kuna me oleme vaadeldava andmestiku ümber pööranud, vaatama filamenti lõpust algusesse
#install.packages("sqldf")
#library(sqldf)
a2b_p2=sqldf("SELECT * FROM a2b,parem2 WHERE a2b.id==parem2.id and a2b.filpt==parem2.filpt")
a2b_p2=a2b_p2[order(a2b_p2$id,-a2b_p2$filpt),]
a2b_v2=sqldf("SELECT * FROM a2b,vasem2 WHERE a2b.id==vasem2.id and a2b.filpt==vasem2.filpt")
a2b_v2=a2b_v2[order(a2b_v2$id,-a2b_v2$filpt),]

	#Leiame galaktikaparvede kaugused vastavatest filamendi punktidest
kaugus2_v=matrix(0,nrow(galgrupp4),1)
for (i in (1:nrow(galgrupp4))){
	kaugus2_v[i]=sqrt((galgrupp4$x[i]-a2b_v2$x[i])**2+(galgrupp4$y[i]-a2b_v2$y[i])**2+(galgrupp4$z[i]-a2b_v2$z[i])**2)
}
kaugus2_p=matrix(0,nrow(galgrupp4),1)
for (i in (1:nrow(galgrupp4))){
	kaugus2_p[i]=sqrt((galgrupp4$x[i]-a2b_p2$x[i])**2+(galgrupp4$y[i]-a2b_p2$y[i])**2+(galgrupp4$z[i]-a2b_p2$z[i])**2)
}

## 2.3 Galaktikaparvede asukoht lähima filamendi punkti suhtes

	#GALGRUPP 3 korral:

kood1=matrix(0,nrow(galgrupp3),1)
for(i in (1:nrow(galgrupp3))){
	if (kaugus1_v[i]<kaugus1_p[i]){
		kood1[i]=-1
	}
	else if (kaugus1_v[i]>kaugus1_p[i]){
		kood1[i]=1
	}
		else kood1[i]=NA
}
sum(is.na(as.matrix(kood1)))
#1] 0 mitte ükski galaktikaparv ei asu mõlemast punktist sama kaugel
kood2=matrix(0,nrow(galgrupp4),1)
for(i in (1:nrow(galgrupp4))){
	if (kaugus2_v[i]<kaugus2_p[i]){
		kood2[i]=-1
	}
	else if (kaugus2_v[i]>kaugus2_p[i]){
		kood2[i]=1
	}
		else kood2[i]=NA
}
sum(is.na(as.matrix(kood2)))
#0 mitte ükski galaktikaparv ei asu mõlemast punktist sama kaugel
galgrupp3=cbind(galgrupp3,kood1)
galgrupp4=cbind(galgrupp4,kood2)

## 2.4 Galaktikaparve teine lähim punkt

kumbpool=matrix(0,nrow(galgrupp3),1)
for(i in (1:nrow(galgrupp3))){
	if (kaugus1_v[i]<kaugus1_p[i]){
		kumbpool[i]=a2b_v1$filpt[i]
	}
	else if (kaugus1_v[i]>kaugus1_p[i]){
		kumbpool[i]=a2b_p1$filpt[i]
	}
		else kumbpool[i]=NA
}
galgrupp3=cbind(galgrupp3,kumbpool)
	#Kui galaktikaparv asub esimese või viimase punkti juures
for(i in (1:nrow(galgrupp3))){
	if(galgrupp3$kumbpool[i]==galgrupp3$punkt[i]){
		galgrupp3$kood1[i]=0
	}
}
	#Sorteerime lähima teise punkti alusel
galgrupp3=galgrupp3[order(galgrupp3$fil_id,galgrupp3$kumbpool),]

kumbpool_2=matrix(0,nrow(galgrupp4),1)
for(i in (1:nrow(galgrupp4))){
	if (kaugus2_v[i]<kaugus2_p[i]){
		kumbpool_2[i]=a2b_v2$filpt[i]
	}
	else if (kaugus2_v[i]>kaugus2_p[i]){
		kumbpool_2[i]=a2b_p2$filpt[i]
	}
		else kumbpool_2[i]=NA
}
kumbpool=galgrupp4$npts-kumbpool_2+1

galgrupp4=cbind(galgrupp4,kumbpool)
for(i in (1:nrow(galgrupp4))){
	if(galgrupp4$kumbpool[i]==galgrupp4$uuspunkt[i]){
		galgrupp4$kood2[i]=0
	}
}
galgrupp4=galgrupp4[order(galgrupp4$fil_id,galgrupp4$kumbpool),] 

## 3. Galaktikaparvede kaugused otspunktist

parvkaugused=matrix(0,nrow(galgrupp3),1)
for (i in (1:nrow(galgrupp3))){
	if (galgrupp3$punkt[i]==1){
		parvkaugused[i]=galgrupp3$parvtelg1[i]
	}
	else{
	parvkaugused[i]=(galgrupp3$punkt[i]-1)*0.5+galgrupp3$kood1[i]*galgrupp3$parvtelg1[i]
	}
}
galgrupp3=cbind(galgrupp3,parvkaugused)

parvkaugused=matrix(0,nrow(galgrupp4),1)
for (i in (1:nrow(galgrupp4))){
	if (galgrupp4$uuspunkt[i]==1){
		parvkaugused[i]=galgrupp4$parvtelg2[i]
	}
	else{
	parvkaugused[i]=(galgrupp4$uuspunkt[i]-1)*0.5+galgrupp4$kood2[i]*galgrupp4$parvtelg2[i]
	}
}
galgrupp4=cbind(galgrupp4,parvkaugused)

	#KONTROLL
max(galgrupp3$parvkaugused)
#[1] 29.73
min(galgrupp3$parvkaugused)
#[1]0.1036733
max(galgrupp4$parvkaugused)
#[1] 29.736

## 4. Naabergalaktikaparvede vahelised kaugused filamendi kaupa

parvparvkaug=matrix(0,nrow(galgrupp4),1)
for (i in (1:nrow(galgrupp4))){
	parvparvkaug[i]=abs(galgrupp4$parvkaugused[i]-galgrupp4$parvkaugused[i+1])
}
galgrupp4=cbind(galgrupp4,parvparvkaug)
parvparvkaug1=matrix(0,nrow(galgrupp3),1)
for (i in (1:nrow(galgrupp3))){
	parvparvkaug1[i]=abs(galgrupp3$parvkaugused[i]-galgrupp3$parvkaugused[i+1])
}
galgrupp3=cbind(galgrupp3,parvparvkaug1)
	#Valed elemendid veergudes parvparvkaug1 on need, mis näitavad filamendi viimase elemendi kaugust,
	#sest ta arvutab erinevuse esimese filamendi viimasest hoopis teise filamendi esimese punkti kauguste vahel

## 5. Teljele projekteeritud naabergalaktikaparvede kauguste vahelise korrelatsiooni leidmine

## 5.1. Galaktikaparvede kaugus paremast ja vaskust naabrist

#galgrupp3
nparvparvkaug1=matrix(0,length(galgrupp3$parvparvkaug1)-1,1)
	for(i in (1:length(galgrupp3$parvparvkaug1)-1)){
		nparvparvkaug1[i]=galgrupp3$parvparvkaug1[i+1]
}
nparvparvkaug1[962]=0
kaugusednihkes1_p=cbind(galgrupp3$fil_id,galgrupp3$punkt,galgrupp3$parvkaugused,galgrupp3$parvparvkaug1,nparvparvkaug1,galgrupp3$rank,galgrupp3$nrich)
colnames(kaugusednihkes1_p)=c("fil_id","punkt","parvkaugused","parvparvkaug1","nparvparvkaug1","rank","nrich")
kaugusednihkes1_p=kaugusednihkes1_p[-962,]
kaugusednihkes1_p=as.data.frame(kaugusednihkes1_p)

#galgrupp4
nparvparvkaug2=matrix(0,length(galgrupp4$parvparvkaug)-1,1)
	for(i in (1:length(galgrupp4$parvparvkaug)-1)){
		nparvparvkaug2[i]=galgrupp4$parvparvkaug[i+1]
}
nparvparvkaug2[969]=0
kaugusednihkes2_p=cbind(galgrupp4$fil_id,galgrupp4$uuspunkt,galgrupp4$parvkaugused,galgrupp4$parvparvkaug,nparvparvkaug2,galgrupp4$rank,galgrupp4$nrich)
colnames(kaugusednihkes2_p)=c("fil_id","punkt","parvkaugused","parvparvkaug","nparvparvkaug2","rank","nrich")
kaugusednihkes2_p=kaugusednihkes2_p[-969,]
kaugusednihkes2_p=as.data.frame(kaugusednihkes2_p)

## 5.2. Ebasobivate paaride eemaldamine (nihutamisel tekkiv viga kahe filamendi vahel)

	#galgrupp3
kn1_p=kaugusednihkes1_p[!duplicated(kaugusednihkes1_p$fil_id, fromLast=TRUE),]
kaugusednihkes1_p=kaugusednihkes1_p[setdiff(rownames(kaugusednihkes1_p),rownames(kn1_p)),]
#796   8
	#galgrupp4#
kn2_p=kaugusednihkes2_p[!duplicated(kaugusednihkes2_p$fil_id, fromLast=TRUE),]
kaugusednihkes2_p=kaugusednihkes2_p[setdiff(rownames(kaugusednihkes2_p),rownames(kn2_p)),]
#803	7
									

## 5.3 Andmestike ühendamine (paremalt ja vasakult loetud filamendid)

## 5.3.1 Loome kaugusednihkes2 tarvis uued filamentide koodid

a=length(unique(kaugusednihkes1_p$fil_id))+length(unique(kaugusednihkes2_p$fil_id))
b=unique(kaugusednihkes1_p$fil_id) #maksimaalne on 7141
c=c(7200:7364)
x=as.data.frame(table(kaugusednihkes2_p$fil_id))
kordused=x$Freq
length(kordused)
fil_id=rep(c,times=kordused)
kaugusednihkes2_2p=kaugusednihkes2_p[,-1]
kaugusednihkes2_2p=cbind(fil_id,kaugusednihkes2_2p)


## 5.3.2 Paneme andmestikud kokku

colnames(kaugusednihkes1_p)=c("fil_id","punkt","parvkaugused","parvparvkaug","nparvparvkaug","rank","nrich")
colnames(kaugusednihkes2_2p)=c("fil_id","punkt","parvkaugused","parvparvkaug","nparvparvkaug","rank","nrich")
kaugusednihkesparved=rbind(kaugusednihkes1_p,kaugusednihkes2_2p)
												
## 5.4 Andmete puhastamine 
## (Iga filamendi viimases reas on väär infot viimase galaktikaparve paiknemise kohta)

kn5=kaugusednihkesparved[!duplicated(kaugusednihkesparved$fil_id, fromLast=TRUE),]
kaugusednihkes_mp=kaugusednihkesparved[setdiff(rownames(kaugusednihkesparved),rownames(kn5)),]
	#Kui mitmel parvel on mõlemad naabrid
dim(kaugusednihkes_mp)
	#1269 

## 5.5 Korrelatsioonide leidmine

kaugusednihkesparved=as.data.frame(kaugusednihkes_mp)
indikaatorid=unique(kaugusednihkesparved$fil_id)

	#PEARSON
parved_korr2=matrix(0,nrow=length(indikaatorid))
for(i in (1:length(indikaatorid))){
	D1=subset(kaugusednihkesparved,fil_id %in% indikaatorid[i])
	d1=D1$parvparvkaug
	d2=D1$nparvparvkaug
	parved_korr2[i,]=cor(d1,d2,method="pearson")
}
uusid=c(1:330)

#Leiame keskmise korrelatsiooni:
parved_korr3=subset(parved_korr2,!is.na(parved_korr2))
mean(parved_korr3)
#-0.2313763

## 5.5.1 Uurime naaberitega galaktikaparvede hulkasid filamentidel

#Mitmel filamendil Pearsoni korrelatsiooni ei arvutatud?
sum(is.na(parved_korr2))
#17
#Maksimaalne ja minimaalne naaberparvede arv filamendil
sagedus_naaber=as.data.frame(table(kaugusednihkesparved$fil_id))
sagedus_naaber=sagedus_naaber[order(sagedus_naaber$Freq),]
sagedus_naaber[1:20,]
sagedus_naaber[300:330,]

	#Filamendile projekteeritud naabergalaktikaparvede vaheliste kauguste vaheline korrelatsioon
plot(uusid,parved_korr2,main="",cex.lab=1.3,cex.axis=1.5,xlab="Filament",ylab="Korrelatsioon")

	#Filamendi teljele projekteeritud galaktikaparvede kaugus naabergalaktikatest
plot(kaugusednihkesparved$parvparvkaug,kaugusednihkesparved$nparvparvkaug,cex.axis=1.3,cex.lab=1.5,main="",xlab="Galaktikaparve parv_(i) ja parv_(i+1) vaheline kaugus",ylab="Galaktikaparve parv_(i+1) ja parv_(i+2) vaheline kaugus")
	#Kui paljudel parvedel on mõlemad naabrid?
length(kaugusednihkesparved$parvparvkaug)
	#1269
	#Kui paljud omavad nii kaugel asuvaid naaberparvesid
parv_sag=as.data.frame(table(cut(kaugusednihkesparved$parvparvkaug,breaks=c(seq(from=0,to=30,by=1)))))
sum(parv_sag$Freq[1:11])
#1240
sum(parv_sag$Freq[11:30])
#42

## 6. Paaris korrelatsioonifunktsiooni leidmine
											
parvots_3=cbind(galgrupp3$fil_id,galgrupp3$parvkaugused)
	#uued koodid filamentidele
a_2=length(unique(galgrupp4$fil_id))
	#genereerime arvud 7200st 7200+165
c_2=c(7200:7364)
x_2=as.data.frame(table(galgrupp4$fil_id))
kordused_2=x_2$Freq
length(kordused_2)
fil_id=rep(c_2,times=kordused_2)
parvots_4=cbind(fil_id,galgrupp4$parvkaugused)
parvots=rbind(parvots_3,parvots_4)
colnames(parvots)=c("fil_id","parvkaugused")
parvots=as.data.frame(parvots)

## 6.1 Leian iga filamendi jaoks kõikide galaktikaparvede kaugused teistest galaktikaparvedest

suurandmestik_p=as.matrix(0,nrow=1,ncol=1)
	suurvektor_p=c(0)
	indikaatorid_p=unique(parvots$fil_id)				
	for (i in (1:length(indikaatorid_p))){
		B1=subset(parvots,fil_id %in% indikaatorid_p[i])
		b1=B1$parvkaugused
			andmed=as.matrix(c(rep(b1,length(b1))))
			andmed=matrix(andmed,nrow=length(b1),ncol=length(b1))
			lahutaja=as.matrix(c(rep(b1,each=length(b1))))
			lahutaja=matrix(lahutaja,nrow=length(b1),ncol=length(b1))
			suurandmestik_p=abs(andmed-lahutaja)
			diag(suurandmestik_p)=NA
			A1=matrix(0,nrow=nrow(suurandmestik_p)-1,ncol=1)
			for (j in (1:ncol(suurandmestik_p))){
				veerg=suurandmestik_p[,j]
				a1=as.matrix(veerg[which(!is.na(veerg))])
				A1=cbind(A1,a1)}
			A1=A1[,-1]
			vektor1_p=as.vector(A1)
			suurvektor_p=c(suurvektor_p,vektor1_p)
	}
suurvektor_p=suurvektor_p[-1]

#Galaktikaparvede kaugused teineteisest 330 erineva filamendi pealt
hist(suurvektor_p,cex.axis=1.3,cex.lab=1.5,breaks=30,main="",xlab="Kaugus r (h^{-1} Mpc)",ylab="Galaktikaparvede paaride arv")

	#Vaatame, kus on täpselt maksimumid:
suurvektorp_sag=as.data.frame(table(cut(suurvektor_p,breaks=c(seq(from=0,to=30,by=1)))))
suurvektorp_sag=suurvektorp_sag[order(-suurvektorp_sag$Freq),]

	#Andmestikule suurvektor_p uued ID-d iga filamendi eraldamiseks
uusid=c(1:330)
x_3=as.data.frame(table(parvots$fil_id))
korda=x_3$Freq*(x_3$Freq-1)
fil_id=rep(uusid,times=korda)
	#Kontroll:
length(suurvektor_p) #10186
length(fil_id) #10186

parvjaotus=cbind(fil_id,suurvektor_p)
parvjaotus=as.data.frame(parvjaotus)

## 6.2 Leiame filamentide galaktikaparvede vaheliste kauguste jaotust
## Moodustame maatriksi DD:30 x 330

	#Leian DD iga filamendi siseselt:
DD_p=matrix(0,nrow=30,ncol=330)
for(i in (1:length(uusid))){
	B1=subset(parvjaotus,fil_id %in% uusid[i])
	b1=B1$suurvektor_p
	DD_p[,i]=as.data.frame(table(cut(b1,breaks=c(seq(from=0,to=30,by=1)))))$Freq #panin sammuks 1
}

## 6.3 Ühtlase jaotusega juhuslike suuruste genereerimine
## Genereerin igale filamendile vastavusse 10 korda rohkem ühtlase jaotusega asetsevaid punkte lõigul [0,30]

x_p=as.data.frame(table(parvots$fil_id))
kordused_p=x_3$Freq
sum(kordused_p*10)
#[1] 19310

## 6.4 Leiame galaktikate kaugused juhuslikest punktidest	

suurandmestikj=as.matrix(0,nrow=1,ncol=1)
uhtlanesuur=c(0)
x=array(0,1)
for(i in (1:length(kordused_p))){
	x=runif(kordused_p[i]*10,0,30)
	andmedj=as.matrix(c(rep(x,length(x))))
	andmedj=matrix(andmedj,nrow=length(x),ncol=length(x))
	lahutajaj=as.matrix(c(rep(x,each=length(x))))
	lahutajaj=matrix(lahutajaj,nrow=length(x),ncol=length(x))
	suurandmestikj=abs(andmedj-lahutajaj)
	diag(suurandmestikj)=NA
	A1=matrix(0,nrow=nrow(suurandmestikj)-1,ncol=1)
		for (j in (1:ncol(suurandmestikj))){
			veerg=suurandmestikj[,j]
			a1=as.matrix(veerg[which(!is.na(veerg))])
			A1=cbind(A1,a1)}
		A1=A1[,-1]
	vektor1j=as.vector(A1)
	uhtlanesuur=c(uhtlanesuur,vektor1j)
}
uhtlanesuur=uhtlanesuur[-1]
#Ühtlase jaotusega juhuslike punktide vahelised kaugused 330 lõigul
hist(uhtlanesuur,cex.lab=1.5,cex.axis=1.3,breaks=30,main="",xlab="Kaugus r (h^{-1} Mpc)",ylab="Ühtlase jaotusega punktipaaride arv")

## 6.5 Moodustame maatriskis RR_p, mis sisaldab igal lõigul olevate juhuslike punktide kauguseid üksteisest
## RR_p:30 x 330

sum(kordused_p*10*(kordused_p*10-1))
#[1] 1192390
length(uhtlanesuur)
#1192390
kordused_uhtlane=kordused_p*10*(kordused_p*10-1)
fil_id=rep(uusid,times=kordused_uhtlane)
uhtljaotus=as.data.frame(cbind(fil_id,uhtlanesuur))

RR_p=matrix(0,nrow=30,ncol=330)
for(i in (1:length(uusid))){
	B1=subset(uhtljaotus,fil_id %in% uusid[i])
	b1=B1$uhtlanesuur
	RR_p[,i]=as.data.frame(table(cut(b1,breaks=c(seq(from=0,to=30,by=1)))))$Freq #panin sammuks 1
}

## 6.6 Leiame DD_p ja RR_p suhte iga filamendi jaoks

KSI_p=matrix(0,nrow=30,ncol=330)
for(i in (1:ncol(DD_p))){
	KSI_p[,i]=(DD_p[,i]/RR_p[,i])*100
}

## 6.8 Leiame paaris korrelatsioonifunktsiooni

ksi_p=matrix(0,nrow=30,ncol=1)
for(i in (1:nrow(KSI_p))){
	ksi_p[i,1]=sum(KSI_p[i,])/330
}

ksi_p 
	#alates 25ndast Mpc on jagamine nulliga
ksi_p=ksi_p[1:23]

telgx=as.data.frame(table(cut(c(0:30),breaks=c(seq(from=0,to=30,by=1)))))$Var
telgx=telgx[1:23]

	#Vaatame jäjestust:
telgksi_p=cbind(telgx,ksi_p)
telgksi_p=telgksi_p[order(-(ksi_p)),]

	#Joonistame: Paaris-korrelatsioonifunktsioon
plot(telgx,ksi_p,col="white",xlim=c(0,25),cex.lab=1.3,main="",xlab="Parvede vaheline kauguse vahemik r, r+dr",ylab="g(r,r+dr)")
lines(ksi_p,col="green")

## 7. Jackknife veahinnangu leidmine

	#1#
#jätan välja i-nda filamendi
jackksi_p=matrix(0,nrow=30,ncol=330)
JackKSI_p=matrix(0,nrow=30,ncol=329)
for(i in (1:ncol(DD_p))){
	DD_psub=DD_p[,-i]
	RR_psub=RR_p[,-i]
	for(j in (1:(ncol(DD_psub)))){
		JackKSI_p[,j]=(DD_psub[,j]/RR_psub[,j])*100
	}
	for(k in (1:(nrow(JackKSI_p)))){
		jackksi_p[k,i]=sum(JackKSI_p[k,])/329
	}
}

#jackksi algab alles 23 reast:
jackksi_p=jackksi_p[1:23,]

	#2#
pseudoksi_p=matrix(0,23,330)
for(i in (1:ncol(jackksi_p))){
	pseudoksi_p[,i]=330*ksi_p-329*jackksi_p[,i]
}

	#3#
keskjack_p=matrix(0,nrow=23,ncol=1)
for(i in (1:nrow(jackksi_p))){
	keskjack_p[i,]=sum(pseudoksi_p[i,])/330
} 					
	
	#4#
#suurus s**2
sruut_p=matrix(0,nrow=23,ncol=1)
for(i in (1:nrow(pseudoksi_p))){
	sruut_p[i,]=1/329*((sum(pseudoksi_p[i,]**2))-1/330*(sum(pseudoksi_p[i,]))**2)
}

	#5#
#Kuidas leida s_tärn?
#s_tärn**2=s**2/k
starn_p=sqrt(sruut_p)/sqrt(330)

	#6#
#Usaldusvahemikud
#ylem=ksi_p+qt(0.025,330-1)*starn_p
#alum=ksi_p-qt(0.025,330-1)*starn_p
ylem=keskjack_p+qt(0.025,330-1)*starn_p
alum=keskjack_p-qt(0.025,330-1)*starn_p

#Joonistame:
plot(telgx,ksi_p,cex.axis=1.3,cex.lab=1.3,col="white",xlim=c(0,25),cex.lab=1.3,main="",xlab="Parvede vaheline kauguse vahemik (r,r+dr] h^{-1} Mpc",ylab="g(r,r+dr)")
lines(ksi_p,col="green")
lines(ylem,col="black")
lines(alum,col="black")
legend(10,0.6,box.col="white",border="white",c("g(r,r+dr)","Jackknife vahemikhinnang"), cex=1.3, fill=c("green","black"))
lines(x=c(0,50),y=c(1,1),lty=b,col="sky blue")
	

#######################################################
## Filamentide otspunktide uurimine

a2=read.table(".../cat_table_a2.txt",na.strings=" ",header=T,sep="")
a3=read.table(".../sdss_fcat_r05fin_gals.txt",na.strings=" ",header=T,sep="")
a33=read.table(".../cat_table_a3.txt",na.strings=" ",header=T,sep="")
a34=read.table(".../dr8_18mar11_gal_v2.txt",na.strings=" ",header=T,sep="") 
rank=a34$rank

dist_fil_end=a3[,c(19)]
koord=a3[,c(1,2,3)]
a3=cbind(a33,dist_fil_end,koord,rank)

	## 1. Andmete filtreerimine, sidumine ja loomine
	
## 1.1 Leiame igale filamendile iga punkti jaoks tunnuse filpt (iga filamendi sisene punkti indentifikaator)

filpt=matrix(0,275599,1)
z=matrix(0,275599,1)
for(i in (1:275599)){
	z[a2$id[i],1]=z[a2$id[i],1]+1
	filpt[i,1]=z[a2$id[i],1]
}
a2x=cbind(a2,filpt)
a2x=a2x[order(a2x$id,a2x$idpts), ]

## 1.2 Korrastan galaktikate andmestikku

a3=a3[order(a3$fil_id,a3$fil_idpts),]
ilmafilament=subset(a3,fil_id==0)
ilmafilament$punkt=0
dim(ilmafilament)
#180248     24
a3=subset(a3,fil_id!=0)
dim(a3)
#396245     23
punkt=a2x$filpt[match(a3$fil_idpts, a2x$idpts)]
a3x=cbind(a3,punkt)
dim(a3x)
#396245     24

## 1.2.1 Valime ainult parvede peagalaktikad

a3x=subset(a3x,rank==1)
dim(a3x)
# 233261    24

## 1.2.2 Valime parved, mis asuvad meist 98 Mpc kuni 252 Mpc kaugusel

kaugus=matrix(0,nrow(a3x),1)
for(i in (1:nrow(a3x))){
	kaugus[i,1]=sqrt(a3x$x[i]**2+a3x$y[i]**2+a3x$z[i]**2)
}
length(kaugus)
# 233261
dim(a3x)
# 233261   24
a3y=cbind(a3x,kaugus)
	#Valime parved sobivalt:
a3z=a3y[which(kaugus>=98 & kaugus<=252),]
dim(a3z)
#109729  erinevat sobivat galatkikaparve


## 1.2.3 Valime parved, mis ei asu filamendi telgedel sobivalt:

ilmafilament2=subset(ilmafilament,rank==1)
kaugus=matrix(0,nrow(ilmafilament2),1)
for(i in (1:nrow(ilmafilament2))){
	kaugus[i,1]=sqrt(ilmafilament2$x[i]**2+ilmafilament2$y[i]**2+ilmafilament2$z[i]**2)
}
ilmafilament2=cbind(ilmafilament2,kaugus)
ilmafilament2=ilmafilament2[which(kaugus>=98 & kaugus<=252),]
dim(ilmafilament2)
#4457 25

## 1.3 Loome kõikide filamentide otspunktide andmestiku

koik_otspunktid=array(0)
ID_otskoik=unique(a2x$id)
for(i in (1:length(ID_otskoik))){
	IDD=ID_otskoik[i]
	vek=a2x[which(a2x$id==IDD),]
	al=vek[which(vek$filpt==1),]
	yl=vek[which(vek$filpt==vek$npts),]
	koik_otspunktid=rbind(koik_otspunktid,al,yl)
}
koik_otspunktid=koik_otspunktid[-1,]
dim(koik_otspunktid)
# 30842    15

## 1.4 Nüüd jätame andmestikust välja need filamendid, mille otspunktid ei asu
## 100 Mpc kuni 250 Mpc kaugusel
## Andmestik otspunktid on vaatluse all olevate filamentide andmestik

E1=koik_otspunktid[which(koik_otspunktid$dist<=250 & koik_otspunktid$dist>=100),]
	#id-de sagedustabel
sagedus=as.data.frame(table(E1$id))
	#Ainult need filamendid, mille mõlemad otspunktid sattusid vahemikku:
vajalikud=sagedus[which(sagedus$Freq==2),]

	#Valime sobivad otspunktid
otspunktid=subset(E1,id %in% vajalikud$Var1)
	#Kontroll
nrow(vajalikud)
#[1] 8144
nrow(otspunktid)
#[1] 16288
8144*2
#[1] 16288
max(otspunktid$dist)
# 249.9965
min(otspunktid$dist)
# 100.0289

			## 2. Vaatleme galaktikaparvi 2 Mpc raadiusega keras

## 2.1 Leiame mitu galaktikaparve (andmestikust a3z) asub filamendi otspunktist 2 Mpc kaugusel

## 2.1.1 Valime sobivaks teljest 2 Mpc kaugusel asuvad parved

a3q=a3z[which(a3z$fil_dist<=2),]
dim(a3q)
#56880    25

## 2.1.2 Leiame galaktikaparved, mis asuvad filamentide otspunktist 2 Mpc kaugusel

A=matrix(0,1,3)
for(i in (1:nrow(otspunktid))){
	ylemine=otspunktid$dist[i]+2
	alumine=otspunktid$dist[i]-2
	a3q_sub1=a3q[which(a3q$kaugus<=ylemine & a3q$kaugus>=alumine),]
	h=matrix(0,nrow(a3q_sub1),3)
	h[,1]=rep(otspunktid$id[i],times=nrow(a3q_sub1))
	h[,2]=a3q_sub1$id
	for(j in (1:nrow(a3q_sub1))){
		h[j,3]=sqrt((otspunktid$x[i]-a3q_sub1$x[j])**2+(otspunktid$y[i]-a3q_sub1$y[j])**2+(otspunktid$z[i]-a3q_sub1$z[j])**2)
	}
	a=h[which(h[,3]<=2),]
	A=rbind(A,a)
}
A=A[-1,]
colnames(A)=c("otspunkt","galaktika","h")

## 2.2 Võtame otspunktide lähedal asuvate galaktikate andmed
	
	#Andmestik A sisaldab filamente, mille otspunkti lähistel asub galaktika
dim(A)
# 42994     3
A=as.data.frame(A)
	#galaktikad_sub sisaldab galaktikaid, mis paiknevad filamentide otspunktide keras
galaktikad_sub=subset(a3q,id %in% A$galaktika)
dim(galaktikad_sub)
# 30470     25
max(galaktikad_sub$nrich)
#256
	#Mitu nii rikast parve filamentide otsas paikneb?
otspunktis_parv=as.data.frame(table(galaktikad_sub$nrich))

## 2.3 Võtame kokku parved, mis on rikkamad kui 15

for(i in (1:nrow(galaktikad_sub))){
	if(galaktikad_sub$nrich[i]>=15){
		galaktikad_sub$nrich[i]=15
	}
	else{galaktikad_sub$nrich[i]=galaktikad_sub$nrich[i]
	}		
}

## 2.4 Loeme üle mitu filamenti

IDD=unique(otspunktid$id)
length(IDD) 
#8144

## 2.5 Loeme üle mitu filamenti neist ei oma otspunktides galaktikaparvesid

IDDD=unique(A$otspunkt)
fil_ilma_id=IDD[!IDD %in% IDDD]
length(fil_ilma_id)
#31

		## 3. Vaatleme galaktikaparvi 1 Mpc raadiusega keras

## 3.1 Valime sobivaks teljest 2 Mpc kaugusel asuvad parved

a3w=a3z[which(a3z$fil_dist<=1),]
dim(a3w)
#36298     25

## 3.2 Leiame galaktikaparved, mis asuvad filamentide otspunktist 1 Mpc kaugusel

B=matrix(0,1,3)
for(i in (1:nrow(otspunktid))){
	ylemine=otspunktid$dist[i]+1
	alumine=otspunktid$dist[i]-1
	a3w_sub1=a3w[which(a3w$kaugus<=ylemine & a3w$kaugus>=alumine),]
	h=matrix(0,nrow(a3w_sub1),3)
	h[,1]=rep(otspunktid$id[i],times=nrow(a3w_sub1))
	h[,2]=a3w_sub1$id
	for(j in (1:nrow(a3w_sub1))){
		h[j,3]=sqrt((otspunktid$x[i]-a3w_sub1$x[j])**2+(otspunktid$y[i]-a3w_sub1$y[j])**2+(otspunktid$z[i]-a3w_sub1$z[j])**2)
	}
	b=h[which(h[,3]<=1),]
	B=rbind(B,b)
}
B=B[-1,]
colnames(B)=c("otspunkt","galaktika","h")

## 3.3 Võtame otspunktide lähedal asuvate galaktikaparvede andmed
	
	#Andmestik A sisaldab filamente, mille otspunkti lähistel asub galaktikaparv
dim(B)
# 11031    3
B=as.data.frame(B)
	#galaktikad_sub2 sisaldab galaktikaid, mis paiknevad filamentide otspunktide keras
galaktikad_sub2=subset(a3w,id %in% B$galaktika)
dim(galaktikad_sub2)
# 10335    25
max(galaktikad_sub2$nrich)
#169
	#Mitu nii rikas parve filamentide otsas paikneb?
otspunktis_parv2=as.data.frame(table(galaktikad_sub2$nrich))

	#Võtame alates rikkus=15 kokku
for(i in (1:nrow(galaktikad_sub2))){
	if(galaktikad_sub2$nrich[i]>=15){
		galaktikad_sub2$nrich[i]=15
	}
	else{galaktikad_sub2$nrich[i]=galaktikad_sub2$nrich[i]
	}		
}

## 3.4 Joonistame ühise pildi 2 Mpc ja 1 Mpc kaugusel asuvate galaktikaparvede kumulatiivsest jaotusest

## 3.4.1 Andmed 2 Mpc kaugusel asuvatest galaktikaparvedest

sagedus1=as.data.frame(table(galaktikad_sub$nrich))
parvedotstes=c(cumsum(sagedus1$Freq))
xtelg1=c(1:15)

## 3.4.2 Andmed 1 Mpc kaugusel asuvatest galaktikaparvedest

sagedus2=as.data.frame(table(galaktikad_sub2$nrich))
parvedotstes2=c(cumsum(sagedus2$Freq))

## 3.4.3 Joonistame

matplot(xtelg1, cbind(parvedotstes,parvedotstes2),xaxt="n",cex.axis=1.5,cex.lab=1.3,pch=c("o","o"), col=c("blue","green"),xlab="Galaktikaparve rikkus", ylab="Kumulatiivne galaktikaparvede arv")
axis(1, 1:15, c("1","2","3","4","5","6","7","8","9","10","11","12","13","14","15+"))
matlines(xtelg1,parvedotstes,col="blue")
matlines(xtelg1,parvedotstes2,col="green")
legend(7,20000,box.col="white",border="white",cex=1.5,c("2 h^(-1)Mpc kera sees","1 h^(-1) Mpc kera sees"),pch="oo",col=c("blue","green"))

## 3.5 Loeme mitu filamenti ei oma otspunktides galaktikaparvesid

IDDD2=unique(B$otspunkt)
length(IDDD2)
#6439
8144-6439
#1705	

## 4. Leian vaatluse all olevate filamentide otspunktide ümber asuvate teiste otspunktide jaotuse

## 4.1 Valin kõikide filamentide seast need otspunktid, mis asuvad vaatlejast
## 98 Mpc kuni 252 Mpc kaugusel

otspunktid_sub=koik_otspunktid[which(koik_otspunktid$dist>=98 & koik_otspunktid$dist<=252),]
dim(otspunktid_sub)
# 17141    15

## 4.2 Leian mitu otspunkti asub vaadeldava filamendi otspunkti 2 Mpc raadiusega keras.

E1=matrix(0,1,4)
for(i in (1:nrow(otspunktid))){
	ylemine=otspunktid$dist[i]+2
	alumine=otspunktid$dist[i]-2
	otspunktid_sub1=otspunktid_sub[which(otspunktid_sub$dist<=ylemine & otspunktid_sub$dist>=alumine),]
	r=otspunktid$id[i]
	otspunktid_sub1=otspunktid_sub1[drop(otspunktid_sub1$id != r),]
	h=matrix(0,nrow(otspunktid_sub1),4)
	h[,1]=rep(otspunktid$id[i],times=nrow(otspunktid_sub1))
	h[,2]=rep(otspunktid$filpt[i],times=nrow(otspunktid_sub1))
	h[,3]=otspunktid_sub1$id
	for(j in (1:nrow(otspunktid_sub1))){
		h[j,4]=sqrt((otspunktid$x[i]-otspunktid_sub1$x[j])**2+(otspunktid$y[i]-otspunktid_sub1$y[j])**2+(otspunktid$z[i]-otspunktid_sub1$z[j])**2)
	}
	e=h[which(h[,4]<=2),]
	E1=rbind(E1,e)
}
E1=E1[-1,]
colnames(E1)=c("otspunkt","punktikood","otspunktlahedal","h")
	#Kui palju erinevaid filamentide otspunkte kerasse sattus?
dim(E1)
#[1] 9808   3

## 4.3 Loon filamentide esimeste punktide andmestiku. Milles on tunnus, mis näitab filamendi esimese
## otspunkti ümber loodud kerasse sattunud teiste otspunktide arvu: FILFIL1

E1=as.data.frame(E1)
E2=subset(E1,punktikood==1)
dim(E2)
#4867    4
sagedus_otsp1=as.data.frame(table(E2$otspunkt))
#Lisan otspunkti koodi:
sagedus_otsp1$filpt=1

#Loon esimese punkti andmestiku:
filfil1=subset(otspunktid,filpt==1)
filfil1$punkte_otsas=sagedus_otsp1$Freq[match(filfil1$id,sagedus_otsp1$Var1)]
filfil1$punkte_otsas=replace(filfil1$punkte_otsas,is.na(filfil1$punkte_otsas),0)

dim(filfil1)

## 4.4 Loon filamentide teiste punktide andmestiku. Milles on tunnus, mis näitab filamendi teise
## otspunkti ümber loodud kerasse sattunud teiste otspunktide arvu: FILFIL2

E3=subset(E1,punktikood!=1)
dim(E3)
#4941    4
sagedus_otsp2=as.data.frame(table(E3$otspunkt))
#Lisan otspunkti koodi:
sagedus_otsp2$filpt=E3$punktikood[match(sagedus_otsp2$Var1,E3$otspunkt)]

#Loon teise punkti andmestiku:
filfil2=subset(otspunktid,filpt!=1)
filfil2$punkte_otsas=sagedus_otsp2$Freq[match(filfil2$id,sagedus_otsp2$Var1)]
filfil2$punkte_otsas=replace(filfil2$punkte_otsas,is.na(filfil2$punkte_otsas),0)

## 4.5 Loon ühise andmestiku FILFIL

filfil=rbind(filfil1,filfil2)
filfil=filfil[order(filfil$id),]

## 4.6 Joonistame

esimene=filfil1$punkte_otsas
teine=filfil2$punkte_otsas
#värvide pakett
install.packages("broman")
library(broman)


#MOSAIC
install.packages("vcd", dependencies=TRUE)
library(vcd)
crosstab <- table(esimene, teine)
crosstab <- crosstab[,] #Reorder or and subset rows
crosstab=addmargins(crosstab)

values <- c(crosstab)
rowvarcat <- c("0","1","2","3","4","5","Summa")
columnvarcat <- c("0","1","2","3","4","5","Summa")
names=c("Filamendi esimene otspunkt", "Filamendi teine otspunkt")
dims <- c(7,7) #columns then rows 

TABS <- structure( c(values), .Dim = as.integer(dims), .Dimnames = structure( list(rowvarcat,columnvarcat ),
.Names = c(names) ) , class = "table") 
mosaic(TABS,spacing=spacing_dimequal(1:2 / 1), pop=FALSE, main="", sub="",highlighting=c(names),highlighting_fill = c("sky blue"),labeling_args=list(gp_labels=(gpar(fontsize=14))))
labeling_cells(text=TABS , clip_cells=FALSE)(TABS )

## 5. Mitu filamenti läbib otspunkti ümber konstrueeritud kera

## 5.1 Vähendan kõikide punktide andmestikku (a2x)

a2y=a2x[which(a2x$dist<=252 & a2x$dist>=98),]
dim(a2y)
#[1] 152235     15

## 5.2 Leian kõik punktid, mis jäävad kerade sisse

U=matrix(0,1,4)
for(i in (1:nrow(otspunktid))){
	ylemine=otspunktid$dist[i]+2
	alumine=otspunktid$dist[i]-2
	filamentide_punktid=a2y[which(a2y$dist<=ylemine & a2y$dist>=alumine),]
	r=otspunktid$id[i]
	filamentide_punktid=filamentide_punktid[drop(filamentide_punktid$id != r),]
	h=matrix(0,nrow(filamentide_punktid),4)
	h[,1]=rep(otspunktid$id[i],times=nrow(filamentide_punktid))
	h[,4]=rep(otspunktid$filpt[i],times=nrow(filamentide_punktid))
	h[,2]=filamentide_punktid$id
	for(j in (1:nrow(filamentide_punktid))){
		h[j,3]=sqrt((otspunktid$x[i]-filamentide_punktid$x[j])**2+(otspunktid$y[i]-filamentide_punktid$y[j])**2+(otspunktid$z[i]-filamentide_punktid$z[j])**2)
	}
	f=h[which(h[,3]<=2),]
	U=rbind(U,f)
}
F=U[-1,]
colnames(F)=c("otspunkt","punktlahedal","h","filpt")
dim(F)
#[1] 39393     3

## 5.3 Eemaldan need punktid, mis on otspunktide korral juba arvesse võetud

G=as.data.frame(F)
indikaatorid=unique(G$otspunkt)
R=matrix(0,1,ncol(G))
colnames(R)=c("otspunkt","punktlahedal","h","filpt")
for(i in (1:length(indikaatorid))){
	ident=indikaatorid[i]
	temp=subset(G,otspunkt==ident)
	temp2=subset(E1,otspunkt==ident)
	r=temp[!(temp$punktlahedal %in% temp2$otspunktlahedal),]
	R=rbind(R,r)
}
R=R[-1,]
#Vaatame andmeid:
head(R)
dim(R)
#16703 	
#Andmestik, mis sisaldab kera läbivaid filamente on R. Igat filamenti on arvesse võetud iga punkt korda

indikaatorid=unique(R$otspunkt)
length(indikaatorid)
# 2657 filamendi otspunkti läbib teine filament

## 5.4 Loen iga filamendi, mis satus kerasse sisse ainult ühekordselt (andmestikus R oli neid nii palju kui
## vaadeldava filamendi punkte sinna sattus)

K=matrix(0,1,ncol(R))
colnames(K)=c("otspunkt","punktlahedal","h","filpt")
for(i in (1:length(indikaatorid))){
	ident=indikaatorid[i]
	temp=subset(R,otspunkt==ident)
	k=temp[!duplicated(temp$punktlahedal),]
	K=rbind(K,k)
}
K=K[-1,]

## 5.5 Vaatleme andmeid:

dim(K)
# 2969  erinevat filamenti läbib kera
length(unique(K$otspunkt)) #KONTROLL
# 2657

## 5.6 Loon filamentide esimeste punktide andmestiku. Milles on tunnus, mis näitab filamendi esimese
## otspunkti ümber loodud kera läbivate filamentide arv: FILFIL1

J1=as.data.frame(K)
J2=subset(J1,filpt==1)
dim(J2)
#1567 

sagedus_labiv1=as.data.frame(table(J2$otspunkt))
#Lisan otspunkti koodi:
sagedus_labiv1$filpt=1
max(sagedus_labiv1$Freq)
# 2
#Loon esimese punkti andmestiku:
filfil1$punkte_labib=sagedus_labiv1$Freq[match(filfil1$id,sagedus_labiv1$Var1)]
filfil1$punkte_labib=replace(filfil1$punkte_labib,is.na(filfil1$punkte_labib),0)

## 5.6 Loon filamentide teiste punktide andmestiku. Milles on tunnus, mis näitab filamendi teise
## otspunkti ümber loodud kera läbivate filamentide arvu: FILFIL2

J3=subset(J1,filpt!=1)
dim(J3)
#1402   4
sagedus_labiv2=as.data.frame(table(J3$otspunkt))
#Lisan otspunkti koodi:
sagedus_labiv2$filpt=J3$punktikood[match(sagedus_labiv2$Var1,J3$otspunkt)]

#Loon teise punkti andmestiku:
filfil2$punkte_labib=sagedus_labiv2$Freq[match(filfil2$id,sagedus_labiv2$Var1)]
filfil2$punkte_labib=replace(filfil2$punkte_labib,is.na(filfil2$punkte_labib),0)

algusesfil=filfil1$punkte_labib
lopusfil=filfil2$punkte_labib

## 5.7 Joonistame:

#MOSAIC

crosstab <- table(algusesfil, lopusfil)
crosstab <- crosstab[,] #Reorder or and subset rows
crosstab=addmargins(crosstab)

values <- c(crosstab)
rowvarcat <-c("0","1","2","Summa")
columnvarcat <-c("0","1","2","Summa")
names=c("Filamendi esimene otspunkt", "Filamendi teine otspunkt")
dims <- c(4,4) #columns then rows 

TABS <- structure( c(values), .Dim = as.integer(dims), .Dimnames = structure( list(rowvarcat,columnvarcat ),
.Names = c(names) ) , class = "table") 
mosaic(TABS,spacing=spacing_dimequal(1:2 / 1), pop=FALSE, main="", sub="",highlighting=c(names),highlighting_fill = c("yellow green"),labeling_args=list(gp_labels=(gpar(fontsize=14))))
labeling_cells(text=TABS , clip_cells=FALSE)(TABS )

	## Vaatleme kera ümber asuvaid otspunkte ja läbivaid filamente üheskoos

## 6.1 Meil on kaks andmestikku filfil1 ja filfil2, mis sisaldavad vajalikke tunnuseid

head(filfil1)
dim(filfil1)
#8144   17
max(filfil1$punkte_otsas)
# 5
max(filfil1$punkte_labib)
#2
head(filfil2)
dim(filfil2)
#8144   17

## 6.2 Vaatleme objektina kera:

#Kontrollime:
kontroll=filfil1$id-filfil2$id
max(kontroll)
#[1] 0
min(kontroll)
#[1] 0   Järjestus on sama
filfil4=cbind(filfil2,c(15422:23565))
osa1=filfil1[,c(1,16,17)]
osa2=filfil4[,c(18,16,17)]
colnames(osa2)=c("id","punkte_otsas","punkte_labib")
kerad=rbind(osa1,osa2)
dim(kerad)
# 16288     3

## 6.3 Joonistame:

crosstab <- table(kerad$punkte_otsas, kerad$punkte_labib)
crosstab <- crosstab[,] #Reorder or and subset rows
crosstab=addmargins(crosstab)

values <- c(crosstab)
rowvarcat <- c("0","1","2","3","4","5","Summa")
columnvarcat <- c("0","1","2","Summa")
names=c("Otspunkti ümber asuvate filamentide otspunktide arv", "Otspunkti ümbrust läbivate filamentide arv")
dims <- c(7,4) #columns then rows 

TABS <- structure( c(values), .Dim = as.integer(dims), .Dimnames = structure( list(rowvarcat,columnvarcat ),
.Names = c(names) ) , class = "table") 
mosaic(TABS,spacing=spacing_dimequal(1:2 / 1), pop=FALSE, main="", sub="",highlighting=c(names),highlighting_fill = c("salmon"),labeling_args=list(gp_labels=(gpar(fontsize=14))))
labeling_cells(text=TABS , clip_cells=FALSE)(TABS )