# R/stheoreme.R In stheoreme: Klimontovich's S-Theorem Algorithm Implementation and Data Preparation Tools

#### Documented in crit.stheoremcxds.stheoremd1char.d1natd1natd1nat.defaultd1natived1specd1spec.defaultd1spectrumd2nat.d1natd2specd2spec.defaultd2spectrumplot.d1natplot.d1specplot.d2specplot.pvalignprint.d1natprint.d1specprint.d2specprint.pvalignprint.stheorempvalignpvalign.defaultpvalignerstheoremstheorem.defaultstheoremmeutild1binutild1cleanutild1filtutild1grouputild2binutild2cleanutild2filtutild2group

```stheoremme <- function(distribution0, distribution1, criterionly=FALSE){

# Date: 13-02-2015

# Create and check some items
f0 <- distribution0/sum(distribution0)
f1 <- distribution1/sum(distribution1)
if (abs(length(f1)-length(f0))>0) {
stop('!Distributions have not equal lenghts!')
}
if (length(f0[which(f0<0)])>0) {
stop('!Distribution0 is not valid: probably contains negatives')
}
if (length(f1[which(f1<0)])>0) {
stop('!Distribution1 is not valid: probably contains negatives')
}
#if (alpha>=1 | alpha<=0) {stop('Alpha value must be 0<|alpha|<1!')}

#new distribution functions with zeros repaired
f0[which(f0<1e-15)]=1e-15
f0 = f0/sum(f0)
f1[which(f1<1e-15)]=1e-15
f1 = f1/sum(f1)

### func entropy calulation
dS <- function(f0,f1){

#new distribution functions with zeros repaired
f0[which(f0<1e-15)]=1e-15
f0 = f0/sum(f0)
f1[which(f1<1e-15)]=1e-15
f1 = f1/sum(f1)

# zero step of analysis
f0s <- f0^(1/1); f0s = f0s/sum(f0s)
dH0to1 <- sum(-f0s*log(f0))-sum(-f1*log(f0))
f1s <- f1^(1/1); f1s = f1s/sum(f1s)
dH1to0 <- sum(-f1s*log(f1))-sum(-f0*log(f1))
greenFlag <- TRUE
D <- 1001 #D is temperature

#check if s-criterion works at all
if (dH0to1*dH1to0 > 0) {
greenFlag = FALSE
dS <- 1001 #absurd high if criterion doesn't work
}
if (dH0to1<0 & greenFlag) {
Tm <- 1+c(0:200)/10 #{if Ham.mean diff is negative => D up(?)}
dH <- NaN
for (i in 1:200) {
D <- Tm[i]
f0s <- f0^(1/D); f0s = f0s/sum(f0s)
dH[i] = abs(sum(-f0s*log(f0))-sum(-f1*log(f0)))
}
D=Tm[which(dH==min(dH))]
f0s <- f0^(1/D); f0s = f0s/sum(f0s)
dS <- sum(-f1*log(f1))-sum(-f0s*log(f0s))
}
if (dH0to1>0 & greenFlag) {
Tm <- 1+c(0:200)/10
dH <- NaN
for (i in 1:200) {
D <- Tm[i]
f1s <- f1^(1/D); f1s = f1s/sum(f1s)
dH[i] = abs(sum(-f1s*log(f1))-sum(-f0*log(f1)))
}
D=Tm[which(dH==min(dH))]
f1s <- f1^(1/D); f1s = f1s/sum(f1s)
dS <- sum(-f1s*log(f1s))-sum(-f0*log(f0))
}
if (abs(dH0to1)==0 | abs(dH1to0)==0) {
dS <- 0
D=1
}

return(list(greenFlag=greenFlag,val=dS,D=D))
}
###

# P_value when if direct evolution is possible:
p_val <- 1

#search for the medium state if direct evolution is forbidden
ds <- dS(f0,f1)
ref0_index <- 50 #initial params
ref1_index <- 50
ref0_ds <- ds\$val
ref1_ds <- 0

if (ds\$greenFlag == FALSE) {
i <- c(1)
while (i<100) { #generate new functions in between
fbtw1 <- ((i)*f0+(100-i)*f1)/100; fbtw1=fbtw1/sum(fbtw1)
ds_forward <- dS(f0,fbtw1)
if (ds_forward\$greenFlag == FALSE) {
i=i+1
} else{
ref0_index <- i
ref0_ds <- ds_forward\$val
i=101
}
}
i <- c(1)
while (i<100) {#generate new functions in between
fbtw1 <- ((i)*f1+(100-i)*f0)/100; fbtw1=fbtw1/sum(fbtw1)
ds_back <- dS(fbtw1,f1)
if (ds_back\$greenFlag  == FALSE) {
i=i+1
} else {
ref1_index <- 100-i
ref1_ds <- ds_back\$val
i=101
}
}

if ((ref0_index-ref1_index) == 1) {
p_val = (max(abs(50-ref0_index),abs(50-ref1_index))/50)^2
#a bit speculative residual function
}

if ((ref0_index-ref1_index) < 1) {
#need to calculate fbtw1 again to match the indices to be at
#the same medium point
fbtw1 <- ((ref0_index)*f0+(100-ref0_index)*f1)/100; fbtw1=fbtw1/sum(fbtw1)
ds_back <- dS(fbtw1,f1)
ref1_ds <- ds_back\$val
p_val = (max(abs(50-ref0_index),abs(50-ref1_index))/50)^2
#a bit speculative residual function
}

if ((ref0_index-ref1_index) > 1) {
ref0_ds = 1001
p_val = 0
}
}

###

#outputs:
dS_02 <- ref0_ds;  if (dS_02==1001) {dS_02=1e-15} #protection patch
dS_21 <- ref1_ds;  if (dS_21==1001) {dS_21=1e-15} #protection patch
dH_val <- -sum(f1*log(f1))+sum(f0*log(f0))
dS_val <- dS_02+dS_21;  if (dS_val==1001) {dS_val=1e-15} #protection patch
dH_ext <- dH_val-dS_val

is_direct <- ds\$greenFlag

# Protection from abnormals when testing
if (p_val<0) {p_val=1e-15}
p_val = p_val
###

ref_indices <- c(ref0_index,ref1_index)

if (criterionly==FALSE) {
formula <- paste(as.character(dH_val),c('{0to1} ='),as.character(dS_val),
c('{0to1} +'),as.character(dH_ext),c('{0to1}'))
} else formula <- {'criterionly'}
#identity <- 1
#if (p_val < alpha) {identity = 0}

return(list(is_direct=is_direct, r2_val = p_val, dH_val=dH_val, dS_val=dS_val,
dH_ext=dH_ext, dS_02=dS_02, dS_21=dS_21, formula=formula,
ref_indices=ref_indices))

}

stheorem<- function(distribution0, distribution1, criterionly) UseMethod("stheorem")
stheorem.default <- function(distribution0, distribution1, criterionly=FALSE) {

distribution0 <- as.vector(distribution0)
distribution1 <- as.vector(distribution1)
est <- stheoremme(distribution0, distribution1, criterionly)
est\$is_direct <- as.logical(est\$is_direct)
est\$r2_val <- as.numeric(est\$r2_val)
est\$dH_val <- as.numeric(est\$dH_val)
est\$dS_val <- as.numeric(est\$dS_val)
est\$dH_ext <- as.numeric(est\$dH_ext)
est\$dS_02 <- as.numeric(est\$dS_02)
est\$dS_21 <- as.numeric(est\$dS_21)
est\$formula <- as.character(est\$formula)
est\$ref_indices <- as.vector(est\$ref_indices)

est\$call <- match.call()
class(est) <- "stheorem"
est
}

print.stheorem <- function(x,...){

if (x\$formula=='criterionly') {
cat("\n")
cat(" S-theorem convergence criterion\n")
cat("\nSystem evolution from state0 to state1 is thermodynamically")
if (x\$r2_val==1) {
cat(" allowed ")
cat("(R^2 = ")
cat(as.character(x\$r2_val)); cat(").\n")
}
if (x\$r2_val<0.00043) {
cat(" forbidden \n")
cat("(R^2 = 0).\n")
}
if (x\$r2_val<1 & x\$r2_val>=0.00043) {
cat(" possible through an indirect\n medium state2 ")
cat("(R^2 = ")
cat(as.character(x\$r2_val)); cat(").\n")
}
}

###
if (x\$formula!='criterionly') {
if (x\$r2_val<=1 & x\$r2_val>=0.00043) {
cat("\n S-theorem open system evolution model\n")
cat("\n")
cat("Overall entropy shift ")
cat("{H1-H0 = dS + dI}:\n")
cat(as.character(x\$formula))
}
if (x\$r2_val<1 & x\$r2_val>=0.00043) {
cat("\nwhere dS consists of two:\n")
cat(as.character(x\$dS_val),c('{0to1} ='),as.character(x\$dS_02),
c('{0to2} +'),as.character(x\$dS_21),c('{2to1}\n'))
}
if (x\$r2_val<0.00043) {
cat("\nSystem evolution is thermodynamically forbidden. No valid model\n")
}
}
cat("\n")
}

crit.stheorem <- function(distribution0, distribution1) {
return(stheorem(distribution0, distribution1, criterionly=TRUE))
}

cxds.stheorem <- function(distribution0, distribution1) {
return(stheorem(distribution0, distribution1, criterionly=FALSE))
}

d1char.d1nat <- function(farr0, farr1, reject=c("")) {

# Date: 23-02-2015

for (i in 1:length(reject)) {
farr0=farr0[which(farr0!=reject[i])]
farr1=farr1[which(farr1!=reject[i])]
}

cat("Separate objects with corresponding object numbers:\n")
print(t(levels(as.factor(c(farr0,farr1)))))

out <- d1nat(as.numeric(as.factor(c(farr0,farr1)))[1:length(farr0)],
as.numeric(as.factor(c(farr0,farr1)))[(length(farr0)+1):(length(farr0)+length(farr1))],
brks=max(as.numeric(as.factor(c(farr0,farr1)))))
return(out)
}

d1native <- function(sample0, sample1, band=c(0,0), brks=0){

# Date: 12-02-2015

# set the window & adjust data to it:
if (band[2]==band[1] & band[2]==0) {
band = c(min(sample1, sample0), max(sample1, sample0))
}
data_0 = sample0[which(sample0>=band[1] & sample0<=band[2])]
data_1  = sample1[which(sample1>=band[1] & sample1<=band[2])]
if (length(data_1)>1 & length(data_0)>1) {
full_range = (max(data_1, data_0)-min(data_1, data_0))
} else stop("Band size doesn't match!")

# set the breaks
if (brks == 0) {
num_breaks = floor(sqrt(length(data_0)))  #dance from zero!
break_size = (max(data_0)-min(data_0))/num_breaks
num_breaks = floor(full_range/break_size)
} else {num_breaks = brks}
break_size = full_range/num_breaks
the_breaks = min(data_1, data_0)+c(0:num_breaks)*break_size

#mid points
the_mids <- the_breaks[1]-break_size/2+c(1:num_breaks)*break_size

#distributions final - as quazi-histogram grouping
f0 <- rep(0,num_breaks)
for (i in 1:(num_breaks)) {
f0[i] = sum(which(data_0>the_breaks[i] & data_0<=the_breaks[i+1])/
which(data_0>the_breaks[i] & data_0<=the_breaks[i+1]))
}
f0[1]=f0[1]+sum(which(data_0==the_breaks[1])/
which(data_0==the_breaks[1])) #take also zeros missed previously
#
f1 <- rep(0,num_breaks)
for (i in 1:(num_breaks)) {
f1[i] = sum(which(data_1>the_breaks[i] & data_1<=the_breaks[i+1])/
which(data_1>the_breaks[i] & data_1<=the_breaks[i+1]))
}
f1[1]=f1[1]+sum(which(data_1==the_breaks[1])/
which(data_1==the_breaks[1])) #take also zeros missed previously
#
f0=f0/sum(f0)
f1=f1/sum(f1)

# absolute chaotic entropy calc
fneutral <- rep(1,length(f0))/length(f0)
entropy_max <- -sum(log(fneutral)*fneutral)

#statistic summary for new distributions
f0_n_expected = sum(f0*the_mids)
f0_n_stdev = sqrt(sum(f0*the_mids^2)-f0_n_expected^2)
f0_n_entropy = -sum(f0*log(f0+1e-15))
if (f0_n_entropy<0) {f0_n_entropy=0} #protection patch
f0_n_1stmod = min(the_mids[which(f0==max(f0))]) #use 1st if where is > 1 modes
ftemp<-f0; ftemp[which(the_mids==f0_n_1stmod)]=0 #temp without 1st mode
f0_n_2ndmod = min(the_mids[which(ftemp==max(ftemp))])
ftemp[which(the_mids==f0_n_2ndmod)]=0 #temp without 1st & 2nd modes
f0_n_3rdmod = min(the_mids[which(ftemp==max(ftemp))])
#
f1_n_expected = sum(f1*the_mids)
f1_n_stdev = sqrt(sum(f1*the_mids^2)-f1_n_expected^2)
f1_n_entropy = -sum(f1*log(f1+1e-15))
if (f1_n_entropy<0) {f1_n_entropy=0} #protection patch
f1_n_1stmod = min(the_mids[which(f1==max(f1))]) #use 1st if where is > 1 modes
ftemp<-f1; ftemp[which(the_mids==f1_n_1stmod)]=0 #temp without 1st mode
f1_n_2ndmod = min(the_mids[which(ftemp==max(ftemp))])
ftemp[which(the_mids==f1_n_2ndmod)]=0 #temp without 1st & 2nd modes
f1_n_3rdmod = min(the_mids[which(ftemp==max(ftemp))])

stat_summary <- cbind(expctd=c(f0=f0_n_expected,f1=f1_n_expected),
var=c(f0=f0_n_stdev^2,f1=f1_n_stdev^2),
fsum=c(f0=sum(f0),f1=sum(f1)),
xmin=c(f0=min(the_breaks),f1=min(the_breaks)),
xmax=c(f0=max(the_breaks),f1=max(the_breaks)),
n=c(f0=num_breaks,f1=num_breaks),
mod1=c(f0=f0_n_1stmod,f1=f1_n_1stmod),
mod2=c(f0=f0_n_2ndmod,f1=f1_n_2ndmod),
mod3=c(f0=f0_n_3rdmod,f1=f1_n_3rdmod),
H_val=c(f0=f0_n_entropy,f1=f1_n_entropy),
H_max=c(f0=entropy_max,f1=entropy_max))

essence='Native distributions'

#many happy returns!
return(list(f0=f0, f1=f1, midpoints=the_mids,
stat_summary=stat_summary, essence=essence))
}

d1nat<- function(sample0, sample1, band=c(0,0), brks=0) UseMethod("d1nat")
d1nat.default <- function(sample0, sample1, band=c(0,0), brks=0) {

sample0 <- as.vector(sample0)
sample1 <- as.vector(sample1)
band <- as.vector(band)
brks <- as.numeric(brks)
est <- d1native(sample0, sample1, band, brks)
est\$f0 <- as.vector(est\$f0)
est\$f1 <- as.vector(est\$f1)
est\$midpoints <- as.vector(est\$midpoints)
est\$stat_summary <- as.matrix(est\$stat_summary)
est\$essence <- as.character(est\$essence)
est\$call <- match.call()
class(est) <- "d1nat"
est
}

print.d1nat <- function(x,...){
cat("Call:\n")
print(x\$call)
cat("Two probability mass functions (\$f0,\$f1) have ")
cat("been generated at the common scale of values (\$midpoints)\n")
cat("Statistics summary:\n")
print(x\$stat_summary)
plot(x)
}

plot.d1nat <- function(x,...){

min_val <- x\$stat_summary[1,4]
max_val <- x\$stat_summary[1,5]
the_mids <- x\$midpoints
the_mids_1 <- the_mids+(max_val-min_val)/(x\$stat_summary[1,6])/7
plot(-1,-1, xlim = c(min_val,max_val),
ylim = c(0,(1.4*max(x\$f0,x\$f1))), ylab="Probability Mass",
xlab = "midpoints")
lines(the_mids, x\$f0, type='h', lwd = 3, col=1)
lines(the_mids_1, x\$f1, type='h', lwd = 3, col=8)
legend(mean(the_mids),(1.4*max(x\$f0,x\$f1)), c("\$f0 to sample0","\$f1 to sample1"),
lty=c(1,1), col=c(1,8))
title(x\$essence)
}

d1spectrum <- function(sample0, sample1, band=c(0,0), brks=0, meansub=TRUE) {

# Date: 08-02-2015

#mean subtraction
if (meansub==TRUE){
sample0 = sample0-mean(sample0)
sample1 = sample1-mean(sample1)
}

# power spectrum and freq spaces calc
pwR0 <- abs(fft(sample0))[1:(floor(length(sample0)/2))]
freq0 <- 0.5*c(1:length(pwR0))/length(pwR0)
pwR1 <- abs(fft(sample1))[1:(floor(length(sample1)/2))]
freq1 <- 0.5*c(1:length(pwR1))/length(pwR1)

# set the window & adjust data to it:
if (band[2]==band[1] & band[2]==0) {
band = c(0, 0.5) # freq band by default: full
}

#check the band
if (min(band[1],band[2])<0) {stop("Wrong band!")}
if (max(band[1],band[2])>0.5) {stop("Wrong band!")}

#choose powers and freqs from within band only
pwR0 = pwR0[which(freq0>=band[1] & freq0<=band[2])]
freq0 = freq0[which(freq0>=band[1] & freq0<=band[2])]
pwR1  = pwR1[which(freq1>=band[1] & freq1<=band[2])]
freq1  = freq1[which(freq1>=band[1] & freq1<=band[2])]
if (length(pwR0)>1 & length(pwR1)>1) {
full_range = band[2]-band[1]
} else stop("Band size doesn't match!")

# set the breaks
if (brks == 0) {
num_breaks = min(length(freq0),length(freq1)) # choose as min of 2 posibilities
break_size = full_range/num_breaks
num_breaks = floor(full_range/break_size)
} else {num_breaks = brks}
break_size = full_range/num_breaks
the_breaks = band[1]+c(0:num_breaks)*break_size

#midpoints
the_mids <- band[1]-break_size/2+c(1:num_breaks)*break_size

#final quazi-histogramming distributions from within the breaks
f0 <- rep(0,num_breaks)
for (i in 1:(num_breaks)) {
f0[i] = sum(pwR0[which(freq0>the_breaks[i] & freq0<=the_breaks[i+1])])
}
f0[1]=f0[1]+sum(pwR0[which(
freq0==the_breaks[1])]) #take also zeros missed previously
#
f1 <- rep(0,num_breaks)
for (i in 1:(num_breaks)) {
f1[i] = sum(pwR1[which(freq1>the_breaks[i] & freq1<=the_breaks[i+1])])
}
f1[1]=f1[1]+sum(pwR1[which(
freq1==the_breaks[1])]) #take also zeros missed previously
#

if (sum(f0)==0 | sum(f1)==0) {
stop('One of the samples is a constant. Analysis is not relevant!')
}

f0=f0/sum(f0)
f1=f1/sum(f1)

# absolute chaotic entropy calc
fneutral <- rep(1,length(f0))/length(f0)
entropy_max <- -sum(log(fneutral)*fneutral)

#statistic summary for new distributions
f0_n_expected = sum(f0*the_mids)
f0_n_stdev = sqrt(sum(f0*the_mids^2)-f0_n_expected^2)
f0_n_entropy = -sum(f0*log(f0+1e-15))
if (f0_n_entropy<0) {f0_n_entropy=0} #protection patch
f0_n_1stmod = min(the_mids[which(f0==max(f0))]) #use 1st if where is > 1 modes
ftemp<-f0; ftemp[which(the_mids==f0_n_1stmod)]=0 #temp without 1st mode
f0_n_2ndmod = min(the_mids[which(ftemp==max(ftemp))])
ftemp[which(the_mids==f0_n_2ndmod)]=0 #temp without 1st & 2nd modes
f0_n_3rdmod = min(the_mids[which(ftemp==max(ftemp))])
#
f1_n_expected = sum(f1*the_mids)
f1_n_stdev = sqrt(sum(f1*the_mids^2)-f1_n_expected^2)
f1_n_entropy = -sum(f1*log(f1+1e-15))
if (f1_n_entropy<0) {f1_n_entropy=0} #protection patch
f1_n_1stmod = min(the_mids[which(f1==max(f1))]) #use 1st if where is > 1 modes
ftemp<-f1; ftemp[which(the_mids==f1_n_1stmod)]=0 #temp without 1st mode
f1_n_2ndmod = min(the_mids[which(ftemp==max(ftemp))])
ftemp[which(the_mids==f1_n_2ndmod)]=0 #temp without 1st & 2nd modes
f1_n_3rdmod = min(the_mids[which(ftemp==max(ftemp))])

stat_summary <- cbind(expctd=c(f0=f0_n_expected,f1=f1_n_expected),
var=c(f0=f0_n_stdev^2,f1=f1_n_stdev^2),
fsum=c(f0=sum(f0),f1=sum(f1)),
xmin=c(f0=min(the_breaks),f1=min(the_breaks)),
xmax=c(f0=max(the_breaks),f1=max(the_breaks)),
n=c(f0=num_breaks,f1=num_breaks),
mod1=c(f0=f0_n_1stmod,f1=f1_n_1stmod),
mod2=c(f0=f0_n_2ndmod,f1=f1_n_2ndmod),
mod3=c(f0=f0_n_3rdmod,f1=f1_n_3rdmod),
H_val=c(f0=f0_n_entropy,f1=f1_n_entropy),
H_max=c(f0=entropy_max,f1=entropy_max))

#many happy returns!
return(list(f0=f0, f1=f1, freqs=the_mids,
stat_summary=stat_summary, essence='Spectral distributions'))
}

d1spec<- function(sample0, sample1, band=c(0,0), brks=0, meansub=TRUE) UseMethod("d1spec")
d1spec.default <- function(sample0, sample1, band=c(0,0), brks=0, meansub=TRUE) {

sample0 <- as.vector(sample0)
sample1 <- as.vector(sample1)
band <- as.vector(band)
brks <- as.numeric(brks)
meansub <- as.logical(meansub)
est <- d1spectrum(sample0, sample1, band, brks, meansub)
est\$f0 <- as.vector(est\$f0)
est\$f1 <- as.vector(est\$f1)
est\$freqs <- as.vector(est\$freqs)
est\$stat_summary <- as.matrix(est\$stat_summary)
est\$essence <- as.character(est\$essence)
est\$call <- match.call()
class(est) <- "d1spec"
est
}

print.d1spec <- function(x,...){
cat("Call:\n")
print(x\$call)
cat("Two spectral probability mass functions (\$f0,\$f1) have ")
cat("been generated at the common scale of values (\$freqs)\n")
cat("Statistics summary:\n")
print(x\$stat_summary)
plot(x)
}

plot.d1spec <- function(x,...){
min_val <- x\$stat_summary[1,4]
max_val <- x\$stat_summary[1,5]
the_mids <- x\$freqs
the_mids_1 <- the_mids+(max_val-min_val)/(x\$stat_summary[1,6])/7
plot(-1,-1, xlim = c(min_val,max_val),
ylim = c(0,(1.4*max(x\$f0,x\$f1))), ylab="Probability Mass",
xlab = "Frequency, 1/(sampling_interval)")
lines(the_mids, x\$f0, type='h', lwd = 3, col=1)
lines(the_mids_1, x\$f1, type='h', lwd = 3, col=8)
legend(mean(the_mids),(1.4*max(x\$f0,x\$f1)), c("\$f0 to sample0","\$f1 t0 sample1"),
lty=c(1,1), col=c(1,8))
title(x\$essence)
}

d2nat.d1nat <- function(d2arr0, d2arr1, band=c(0,0), brks=64, method='default'){

# Date: 08-02-2015

#3 methods: default image hist, cols and rows
d2arr0.as.d1 <- c(d2arr0) #default
d2arr1.as.d1 <- c(d2arr1) #default
if (method=='cols'){
d2arr0.as.d1 <- NaN
for (i in 1:length(d2arr0[1,])){
d2arr0.as.d1[i]=mean(d2arr0[,i])
}
d2arr1.as.d1 <- NaN
for (i in 1:length(d2arr1[1,])){
d2arr1.as.d1[i]=mean(d2arr1[,i])
}
}
if (method=='rows'){
d2arr0.as.d1 <- NaN
for (i in 1:length(d2arr0[,1])){
d2arr0.as.d1[i]=mean(d2arr0[i,])
}
d2arr1.as.d1 <- NaN
for (i in 1:length(d2arr1[,1])){
d2arr1.as.d1[i]=mean(d2arr1[i,])
}
}

#now simply apply d1nat
out <- d1nat(sample0=d2arr0.as.d1, sample1=d2arr1.as.d1, band=band, brks=brks)

return(out)
}

d2spectrum <- function(d2arr0, d2arr1, band=c(0,0), brks=0,

# Date: 09-02-2015

# check if arrays are squared
L0 <- c(length(d2arr0[1,]),length(d2arr0[,1]))
L1<- c(length(d2arr1[1,]),length(d2arr1[,1]))
if ((max(L0)-min(L0)) > 0) {cat('Caution! Array0 isnt square. Will be cut!\n')}
if ((max(L1)-min(L1)) > 0) {cat('Caution! Array1 isnt square. Will be cut!\n')}

# cut to the same size/square
L0 = min(L0)
L1 = min(L1)
d2arr0 <- d2arr0[1:L0,1:L0]
d2arr1 <- d2arr1[1:L1,1:L1]

#mean subtraction
if (meansub==TRUE){
d2arr0 = d2arr0-mean(d2arr0)
d2arr1 = d2arr1-mean(d2arr1)
}

#calc the size of future fft quarter space
pw_size0 <- floor(L0/2-1)
pw_size1 <- floor(L1/2-1)

#spectr_calc
pw0 <- abs(fft(d2arr0))
pleft <- pw0[1:(pw_size0),(pw_size0+1):(1)]
pright <- pw0[1:(pw_size0),L0:(L0-pw_size0+1)]
pw0 = cbind(pleft,pright)
pw0=pw0^2
#
pw1 <- abs(fft(d2arr1))
pleft <- pw1[1:(pw_size1),(pw_size1+1):(1)]
pright <- pw1[1:(pw_size1),L1:(L1-pw_size1+1)]
pw1 = cbind(pleft,pright)
pw1=pw1^2

#angular pixel gouping
if (method == 'ang' | method == 'ang90') {
quanglos0 <- array(0,c(pw_size0, (2*pw_size0+1)))
for (i in 1:pw_size0){
for (j in 1:(2*pw_size0+1)) {
quanglos0[i,j]=(180/pi)*atan((j-pw_size0-1)/(i-0.999999999999999))
}
}
#
quanglos1 <- array(0,c(pw_size1, (2*pw_size1+1)))
for (i in 1:pw_size1){
for (j in 1:(2*pw_size1+1)) {
quanglos1[i,j]=(180/pi)*atan((j-pw_size1-1)/(i-0.999999999999999))
}
}

#3 arrays to define with angular space
angle <- c(1:180)
freq0 <- c(0:180)
freq1 <- c(0:180)

#angular powers calc
pwR0<-rep(0,180); pwR1<-rep(0,180)
for (k in 1:180) {
pwR0[k]=sum(pw0[which(angle[k]==(91+as.integer(quanglos0)))])
pwR1[k]=sum(pw1[which(angle[k]==(91+as.integer(quanglos1)))])
}

# set the window & adjust data to it:
if (band[2]==band[1] & band[2]==0) {
band = c(0, 180)  #freq band by default: full
}
if (min(band[1],band[2])<0) {stop("Wrong band!")}
if (max(band[1],band[2])>180) {stop("Wrong band!")}
pwR0 = pwR0[which(freq0>=band[1] & freq0<=band[2])]
freq0 = freq0[which(freq0>=band[1] & freq0<=band[2])]
pwR1  = pwR1[which(freq1>=band[1] & freq1<=band[2])]
freq1  = freq1[which(freq1>=band[1] & freq1<=band[2])]
if (length(pwR0)>1 & length(pwR1)>1) {
full_range = band[2]-band[1]
} else stop("Band size doesn't match!")

# set the breaks:
if (brks == 0) {
num_breaks = full_range #it is 0:180 degrees
break_size = full_range/num_breaks
num_breaks = floor(full_range/break_size)
} else {num_breaks = brks}
break_size = full_range/num_breaks
the_breaks = band[1]+c(0:num_breaks)*break_size
####

} else {
for (i in 1:pw_size0) {
for (j in 1:(2*pw_size0+1)) {
}
}
#
for (i in 1:pw_size1) {
for (j in 1:(2*pw_size1+1)) {
}
}

#4 arrays to define with radial space
freq0 <- 0.5*c(0:(pw_size0-1))/pw_size0
freq1 <- 0.5*c(0:(pw_size1-1))/pw_size1

pwR0<-rep(0,pw_size0)
for (k in 1:pw_size0) {
}
#
pwR1<-rep(0,pw_size1)
for (k in 1:pw_size1) {
}

# set the window & adjust data to it:
if (band[2]==band[1] & band[2]==0) {
band = c(0, 0.5) #freq band by default: full
}
if (min(band[1],band[2])<0) {stop("Wrong band!")}
if (max(band[1],band[2])>0.5) {stop("Wrong band!")}
pwR0 = pwR0[which(freq0>=band[1] & freq0<=band[2])]
freq0 = freq0[which(freq0>=band[1] & freq0<=band[2])]
pwR1  = pwR1[which(freq1>=band[1] & freq1<=band[2])]
freq1  = freq1[which(freq1>=band[1] & freq1<=band[2])]
if (length(pwR0)>1 & length(pwR1)>1) {
full_range = band[2]-band[1]
} else stop("Band size doesn't match!")

# set the breaks:
if (brks == 0) {
num_breaks = min(length(pwR0),length(pwR1)) # dance from what is less size!
break_size = full_range/num_breaks
num_breaks = floor(full_range/break_size)
} else {num_breaks = brks}
break_size = full_range/num_breaks
the_breaks = band[1]+c(0:num_breaks)*break_size
}

#distribution midpoints
the_mids <- band[1]-break_size/2+c(1:num_breaks)*break_size

#quazi-histogram grouping to form final f0 & f1 distributions
f0 <- rep(0,num_breaks)
for (i in 1:(num_breaks)) {
f0[i] = sum(pwR0[which(freq0>=the_breaks[i] & freq0<the_breaks[i+1])])
}
#f0[num_breaks]=f0[num_breaks]+sum(pwR0[which(
#freq0==the_breaks[num_breaks+1])]) #take also zeros missed previously
f1 <- rep(0,num_breaks)
for (i in 1:(num_breaks)) {
f1[i] = sum(pwR1[which(freq1>=the_breaks[i] & freq1<the_breaks[i+1])])
}
#special case : angle is shifted to 90degrees
if (method=='ang90') {
the_mids = the_mids-90
for (i in 1:num_breaks) {
if (the_mids[i]<0) {the_mids[i]=180+the_mids[i]}
}
}

#
if (sum(f0)==0 | sum(f1)==0) {
stop('One of the samples is a constant. Analysis is not relevant!')
}

f0=f0/sum(f0)
f1=f1/sum(f1)

# absolute chaotic entropy calc
fneutral <- rep(1,length(f0))/length(f0)
entropy_max <- -sum(log(fneutral)*fneutral)

#statistic summary for new distributions
f0_n_expected = sum(f0*the_mids)
f0_n_stdev = sqrt(sum(f0*the_mids^2)-f0_n_expected^2)
f0_n_entropy = 0.0001*round(-sum(f0*log(f0+1e-15))*10000)
if (f0_n_entropy<0) {f0_n_entropy=0} #protection patch
f0_n_1stmod = min(the_mids[which(f0==max(f0))]) #use 1st if where is > 1 modes
ftemp<-f0; ftemp[which(the_mids==f0_n_1stmod)]=0 #temp without 1st mode
f0_n_2ndmod = min(the_mids[which(ftemp==max(ftemp))])
ftemp[which(the_mids==f0_n_2ndmod)]=0 #temp without 1st & 2nd modes
f0_n_3rdmod = min(the_mids[which(ftemp==max(ftemp))])
#
f1_n_expected = sum(f1*the_mids)
f1_n_stdev = sqrt(sum(f1*the_mids^2)-f1_n_expected^2)
f1_n_entropy = 0.0001*round(-sum(f1*log(f1+1e-15))*10000)
if (f1_n_entropy<0) {f1_n_entropy=0} #protection patch
f1_n_1stmod = min(the_mids[which(f1==max(f1))]) #use 1st if where is > 1 modes
ftemp<-f1; ftemp[which(the_mids==f1_n_1stmod)]=0 #temp without 1st mode
f1_n_2ndmod = min(the_mids[which(ftemp==max(ftemp))])
ftemp[which(the_mids==f1_n_2ndmod)]=0 #temp without 1st & 2nd modes
f1_n_3rdmod = min(the_mids[which(ftemp==max(ftemp))])

#
stat_summary <- cbind(expctd=c(f0=f0_n_expected,f1=f1_n_expected),
var=c(f0=f0_n_stdev^2,f1=f1_n_stdev^2),
fsum=c(f0=sum(f0),f1=sum(f1)),
xmin=c(f0=min(the_breaks),f1=min(the_breaks)),
xmax=c(f0=max(the_breaks),f1=max(the_breaks)),
n=c(f0=num_breaks,f1=num_breaks),
mod1=c(f0=f0_n_1stmod,f1=f1_n_1stmod),
mod2=c(f0=f0_n_2ndmod,f1=f1_n_2ndmod),
mod3=c(f0=f0_n_3rdmod,f1=f1_n_3rdmod),
H_val=c(f0=f0_n_entropy,f1=f1_n_entropy),
H_max=c(f0=entropy_max,f1=entropy_max))

#many happy returns!
return(list(f0=f0, f1=f1, efs=the_mids, method=method,
stat_summary=stat_summary, essence=c('Spectral distributions')))
}

d2spec<- function(d2arr0, d2arr1, band=c(0,0), brks=0,
d2spec.default <- function(d2arr0, d2arr1, band=c(0,0),

# check they are matrices
if (class(d2arr0)!='matrix' | class(d2arr1)!='matrix') {
stop('!one or both arrays are not of matrix class!')
}
d2arr0 <- as.matrix(d2arr0)
d2arr1 <- as.matrix(d2arr1)
band <- as.vector(band)
brks <- as.numeric(brks)
method <- as.character(method)
meansub <- as.logical(meansub)
est <- d2spectrum(d2arr0, d2arr1, band, brks, method, meansub)
est\$f0 <- as.vector(est\$f0)
est\$f1 <- as.vector(est\$f1)
est\$efs <- as.vector(est\$efs)
est\$method <- as.character(est\$method)
est\$stat_summary <- as.matrix(est\$stat_summary)
est\$essence <- as.character(est\$essence)
est\$call <- match.call()
class(est) <- "d2spec"
est
}

print.d2spec <- function(x,...){
cat("Call:\n")
print(x\$call)
cat("Two spectral probability mass functions (\$f0,\$f1) have ")
cat("been generated at the common scale of ")
if (x\$method=='ang90' | x\$method=='ang') {
cat("spectral angles (\$efs)\n")
cat("Statistics summary:\n")
print(x\$stat_summary)
plot(x)
}

plot.d2spec <- function(x,...){
min_val <- x\$stat_summary[1,4]
max_val <- x\$stat_summary[1,5]
the_mids <- x\$efs
the_mids_1 <- the_mids+(max_val-min_val)/(x\$stat_summary[1,6])/7

if (x\$method=='ang90') {xlab <- "Spectr_angles-90, degrees"}
if (x\$method=='ang') {xlab <- "Spectr_angles, degrees"}
plot(-1,-1, xlim = c(min_val,max_val),
ylim = c(0,(1.4*max(x\$f0,x\$f1))), ylab="Probability Mass",
xlab=xlab)
lines(the_mids, x\$f0, type='h', lwd = 3, col=1)
lines(the_mids_1, x\$f1, type='h', lwd = 3, col=8)
legend(mean(the_mids),(1.4*max(x\$f0,x\$f1)), c("\$f0 to d2arr0","\$f1 to d2arr1"),
lty=c(1,1), col=c(1,8))
title(x\$essence)
}

pvaligner <- function(distribution0, distribution1,
outcomes0=c(0), outcomes1=c(0), tolerance=0) {

# Date: 04-02-2015

f0 <- distribution0
f1 <- distribution1
# check the consistence of the things
if (min(f0)<0 | min(f1)<0) {
stop('!One of the distributions is invalid (<0?)')
}
if (tolerance == 0) {tolerance = 0.001} #default tolerance <-0.001
if (tolerance >= 1 | tolerance <= 0) {
stop('!Tolerance must be: 0<|tolerance|<1 !')
}

# set f0,f1 and new scale
if (length(outcomes0)<=1) {
x0 <- seq(1:length(f0))
} else {x0 <- outcomes0}
if (length(outcomes1)<=1) {
x1 <- seq(1:length(f1))
} else {x1 <- outcomes1}

#check other things
if ((abs(length(x0)-length(f0))>0) | (abs(length(x1)-length(f1))>0)) {
stop('!Lengths of distribution vector(s) dont match the length(s)
of vector(s) of outcomes')
}

#create equidistant integer scale of values
absolute_min <- min(x0,x1)
if (min(x0,x1) <= 0) {
x0=x0-absolute_min+1 #shift values all to positives
x1=x1-absolute_min+1 #shift values all to positives
}

#check if scales of value arrayas are equidistant and compatible
#using the tolerance value
if ((x0[1]>x0[2]) | (x1[1]>x1[2])) {
stop('Vector of outcomes is in decreasing order! Re-sort please')
}
dif_array1 <- NaN
for (i in 1:length(x0)-1) {dif_array1[i]=x0[i+1]-x0[i]}
ref_value1 <- (max(dif_array1)-min(dif_array1))/mean(dif_array1)
if (ref_value1>tolerance) {
stop('Scale of 1st vector of outcomes is not equidistant!')
}
dif_array2 <- NaN
for (i in 1:length(x1)-1) {dif_array2[i]=x1[i+1]-x1[i]}
ref_value2 <- (max(dif_array2)-min(dif_array2))/mean(dif_array2)
if (ref_value2>tolerance) {
stop('Scale of 2nd vector of outcomes is not equidistant!')
}
if (abs(mean(dif_array2)-mean(dif_array1))>tolerance*mean(dif_array1)) {
stop('Vectors of outcomes are not compatible!')
}

#create the common scale of values
xvector <- c(1:(1+as.integer(round((max(x0,x1)-min(x0,x1))/((x0[2]-x0[1]))))))

#new distributions at common scale
if (min(x0)<min(x1)) {
temp_array <- rep(0,length(xvector)-length(x0))
f0 <- c(f0,temp_array)
#f0[which(f0<1e-15)]=0
temp_array <- rep(0,length(xvector)-length(x1))
f1 <- c(temp_array,f1)

} else {
temp_array <- rep(0,length(xvector)-length(x1))
f1<- c(f1,temp_array)
temp_array <- rep(0,length(xvector)-length(x0))
f0 <- c(temp_array,f0)
}
#
f0=f0/sum(f0)
f1=f1/sum(f1)

#common x scale:
x0=x0-min(x0,x1)+1
the_mids=xvector+absolute_min-1

# absolute chaotic entropy calc
fneutral <- rep(1,length(f0))/length(f0)
entropy_max <- -sum(log(fneutral)*fneutral)

#statistic summary for new distributions
f0_n_expected = sum(f0*the_mids)
f0_n_var = (sum(f0*the_mids^2)-f0_n_expected^2)
f0_n_entropy = -sum(f0*log(f0+1e-15))
if (f0_n_entropy<0) {f0_n_entropy=0} #protection patch
f0_n_1stmod = min(the_mids[which(f0==max(f0))]) #use 1st if where is > 1 modes
ftemp<-f0; ftemp[which(the_mids==f0_n_1stmod)]=0 #temp without 1st mode
f0_n_2ndmod = min(the_mids[which(ftemp==max(ftemp))])
ftemp[which(the_mids==f0_n_2ndmod)]=0 #temp without 1st & 2nd modes
f0_n_3rdmod = min(the_mids[which(ftemp==max(ftemp))])
#
f1_n_expected = sum(f1*the_mids)
f1_n_var = (sum(f1*the_mids^2)-f1_n_expected^2)
f1_n_entropy = -sum(f1*log(f1+1e-15))
if (f1_n_entropy<0) {f1_n_entropy=0} #protection patch
f1_n_1stmod = min(the_mids[which(f1==max(f1))]) #use 1st if where is > 1 modes
ftemp<-f1; ftemp[which(the_mids==f1_n_1stmod)]=0 #temp without 1st mode
f1_n_2ndmod = min(the_mids[which(ftemp==max(ftemp))])
ftemp[which(the_mids==f1_n_2ndmod)]=0 #temp without 1st & 2nd modes
f1_n_3rdmod = min(the_mids[which(ftemp==max(ftemp))])

stat_summary <- cbind(expctd=c(f0=f0_n_expected,f1=f1_n_expected),
var=c(f0=f0_n_var,f1=f1_n_var),
fsum=c(f0=sum(f0),f1=sum(f1)),
xmin=c(f0=min(the_mids),f1=min(the_mids)),
xmax=c(f0=max(the_mids),f1=max(the_mids)),
n=c(f0=length(the_mids),f1=length(the_mids)),
mod1=c(f0=f0_n_1stmod,f1=f1_n_1stmod),
mod2=c(f0=f0_n_2ndmod,f1=f1_n_2ndmod),
mod3=c(f0=f0_n_3rdmod,f1=f1_n_3rdmod),
H_val=c(f0=f0_n_entropy,f1=f1_n_entropy),
H_max=c(f0=entropy_max,f1=entropy_max))

#many happy returns!
return(list(new_scale=the_mids, f0=f0, f1=f1, stat_summary=stat_summary))
}

pvalign<- function(distribution0, distribution1,
outcomes0=c(0), outcomes1=c(0), tolerance=0) UseMethod("pvalign")
pvalign.default <- function(distribution0, distribution1,
outcomes0=c(0), outcomes1=c(0), tolerance=0) {

distribution0 <- as.vector(distribution0)
distribution1 <- as.vector(distribution1)
outcomes0 <- as.vector(outcomes0)
outcomes1 <- as.vector(outcomes1)
tolerance <- as.numeric(tolerance)
est <- pvaligner(distribution0, distribution1,
outcomes0, outcomes1, tolerance)
est\$new_scale <- as.vector(est\$new_scale)
est\$f0 <- as.vector(est\$f0)
est\$f1 <- as.vector(est\$f1)
est\$stat_summary <- as.matrix(est\$stat_summary)
est\$call <- match.call()
class(est) <- "pvalign"
est
}

print.pvalign <- function(x,...){
cat("Call:\n")
print(x\$call)
cat("Two spectral probability mass functions (\$f0,\$f1) have ")
cat("been generated at the new common scale (\$new_scale)\n")
cat("Statistics summary:\n")
print(x\$stat_summary)
plot(x)
}

plot.pvalign <- function(x,...){

min_val <- x\$stat_summary[1,4]
max_val <- x\$stat_summary[1,5]
the_mids <- x\$new_scale
the_mids_1 <- the_mids+(max_val-min_val)/(x\$stat_summary[1,6])/7

plot(-1,-1, xlim = c(min_val,max_val),
ylim = c(0,(1.4*max(x\$f0,x\$f1))), ylab="Probability Mass",
xlab = "New scale of outcomes")
lines(the_mids, x\$f0, type='h', lwd = 3, col=1)
lines(the_mids_1, x\$f1, type='h', lwd = 3, col=8)
legend(mean(the_mids),(1.4*max(x\$f0,x\$f1)),
c("f0 to outcomes0","f1 to outcomes1"),
lty=c(1,1), col=c(1,8))
title("Re-scaled distributions")
}

utild1bin <- function(arr0, arr1, method='mean', trsh=0, inverted=FALSE, d2=FALSE){

# Date: 12-02-2015

thr0<-mean(c(arr0,arr1)) #default thresholds is mean for two
thr1<-thr0

if (method == 'mean') {thr1=mean(c(arr0,arr1));thr0=thr1} #default
if (method == 'mean0') {thr1=mean(c(arr0));thr0=thr1}
if (method == 'mean1') {thr1=mean(c(arr1));thr0=thr1}
if (method == 'med')  {thr1=median(c(arr0,arr1));thr0=thr1}
if (method == 'med0')  {thr1=median(c(arr0));thr0=thr1}
if (method == 'med1')  {thr1=median(c(arr1));thr0=thr1}
if (method == 'spec') {
if (length(trsh)==2) {thr0=trsh[1]; thr1=trsh[2]}
if (length(trsh)==1) {thr0=trsh; thr1=trsh}
if (length(trsh)>2) {stop('!more than 2 threshold levels introduced')}
}
if (method == 'imean') {thr0=mean(c(arr0)); thr1=mean(c(arr1))}
if (method == 'imed') {thr0=median(c(arr0)); thr1=median(c(arr1))}
if (method == 'i1stQ') {
thr0=summary(c(arr0))[2]
thr1=summary(c(arr1))[2]
}
if (method == 'i3rdQ') {
thr0=summary(c(arr0))[5]
thr1=summary(c(arr1))[5]
}
if (method == 'isubsd') {
thr0=mean(c(arr0))-sd(c(arr0))
thr1=mean(c(arr1))-sd(c(arr1))
}
thr0=mean(c(arr0))+sd(c(arr0))
thr1=mean(c(arr1))+sd(c(arr1))
}

if (d2==FALSE) {
b0 <- NaN
b1 <- NaN
if (class(arr0)=='matrix' | class(arr1)=='matrix') {
stop ('One of the arrays has > 1 dimensions')
}
} else {
b0 <- array(0,dim(arr0))
b1 <- array(0,dim(arr1))
}

if (inverted) {
b0[which(arr0<=thr0)]=TRUE
b0[which(arr0>thr0)]=FALSE
b1[which(arr1<=thr1)]=TRUE
b1[which(arr1>thr1)]=FALSE
} else {
b0[which(arr0<=thr0)]=FALSE
b0[which(arr0>thr0)]=TRUE
b1[which(arr1<=thr1)]=FALSE
b1[which(arr1>thr1)]=TRUE
}

#counts:
ntrue0 = sum(b0[which(b0==TRUE)])
ntrue1 = sum(b1[which(b1==TRUE)])
nfalse0 = length(b0)-ntrue0
nfalse1 = length(b1)-ntrue1

#
counts<-cbind(n_false=c(bin0=nfalse0,bin1=nfalse1),
n_true=c(bin0=ntrue0, bin1=ntrue1))

return(list(bin0=b0,bin1=b1,counts=counts))
}

########
utild2bin <- function(d2arr0, d2arr1, method='mean', trsh=0, inverted=FALSE){

if (class(d2arr0)!='matrix' | class(d2arr1)!='matrix') {
stop ('One of the arrays has < 2 dimensions')
}

return(utild1bin(d2arr0, d2arr1,
method=method, trsh=trsh, inverted=inverted, d2=TRUE))
}

utild1clean <- function(arr0, arr1, method='bothends', nsigma=3, d2=FALSE) {

# Date 15-02-2015

M0=mean(arr0); M1=mean(arr1)
S0=sd(arr0); S1=sd(arr1)

left0 <- M0-nsigma*S0; right0 <- M0+nsigma*S0
left1 <- M1-nsigma*S1; right1 <- M1+nsigma*S1

if (d2==FALSE) {
c0 <- arr0
c1 <- arr1
if (class(arr0)=='matrix' | class(arr1)=='matrix') {
stop ('One of the arrays has > 1 dimensions')
}
} else {
c0 <- arr0
c1 <- arr1
}

Ncount0 <- 0; Ncount1 <- 0

if (method=='left') {
Ncount0=length(c0[which(c0<left0)])
c0[which(arr0<left0)]=left0
Ncount1=length(c1[which(c1<left1)])
c1[which(arr1<left1)]=left1
}
if (method=='right') {
Ncount0=length(c0[which(c0>right0)])
c0[which(arr0>right0)]=right0
Ncount1=length(c1[which(c1>right1)])
c1[which(arr1>right1)]=right1
}
if (method!='left' & method!='right') {
Ncount0=length(c0[which(c0<left0)])+length(c0[which(c0>right0)])
c0[which(arr0<left0)]=left0
c0[which(arr0>right0)]=right0
Ncount1=length(c1[which(c1<left1)])+length(c1[which(c1>right1)])
c1[which(arr1<left1)]=left1
c1[which(arr1>right1)]=right1
}

cat(as.character(Ncount0)); cat(" outliers in array0 and ");
cat(as.character(Ncount1)); cat(" outliers in array1");
cat(" have been replaced by critical values\n"); cat("\n")

return(list(clean0=c0, clean1=c1))
}

utild2clean <- function(d2arr0, d2arr1, method='bothends', nsigma=3){

if (class(d2arr0)!='matrix' | class(d2arr1)!='matrix') {
stop ('One of the arrays has < 2 dimensions')
}

return(utild1clean(d2arr0, d2arr1, method=method, nsigma=nsigma, d2=TRUE))
}

utild1filt <- function(arr0, arr1, outsize=2, strong=1) {

# Date 23-02-2015

if (outsize<0) {stop('! bandpass frame is invalid (must be >= 0)')}
if (outsize>0.4*min(length(arr0),length(arr1))) {
stop('! outsize is too high')
}
if (strong<=0) {strong=1} #stupid protection
if (class(arr0)=="matrix" | class(arr1)=="matrix") {
stop('! one/both arrays have >1 dims!')
}

#gauss kernel calc func
d1gausscoeffs <- function(s, n) {
t <- seq(-n,n,1)
out <- exp(-(t^2/(2*s^2))); out=out/sum(out)
return (out)
}

#convolution for lowpass
sigma <- ceiling(outsize)*strong
if (outsize==0) {sigma=1}
outlow0 <- arr0
outlow1 <- arr1
convolve(arr0,kern,type='f')
convolve(arr1,kern,type='f')

out0<-outlow0
out1<-outlow1

return(list(filt0=out0, filt1=out1))
}

#####

utild2filt <- function(d2arr0, d2arr1, outsize=2, strong=1) {

# Date 23-02-2015

if (class(d2arr0)!="matrix" | class(d2arr1)!="matrix") {
stop('! one/both arrays have <2 dims!')
}
if (outsize<0) {stop('! bandpass frame is invalid (must be >= 0)')}
if (outsize>0.4*min(length(d2arr0[,1]),length(d2arr1[,1]),
length(d2arr0[1,]),length(d2arr1[1,]))) {
stop('! outsize is too high')
}
if (strong<=0) {strong=1} #stupid protection

#gauss kernel calc func
# 2D
d2gausscoeffs <- function(s, n) {
t<-array(0, c((2*n+1),(2*n+1)))
for (i in 1:(2*n+1)) {
for (j in 1:(2*n+1)) {
t[i,j] <- sqrt((i-n-1)^2+(j-n-1)^2)
}
}
out <- exp(-(t^2/(2*s^2))); out=out/sum(out)
return (out)
}

# convolve 2d function
conv2d <- function(arr2d, kernel2) {
gap <- (length(kernel2[1,])-1)/2
Wn <- length(arr2d[1,])-2*gap
Hn <- length(arr2d[,1])-2*gap

out <- arr2d
for (i in (1+gap):(Hn+gap)) {
for (j in (1+gap):(Wn+gap)) {
out[i,j]=0
for (k in -gap:gap) {
for (l in -gap:gap) {
out[i,j]=out[i,j]+arr2d[i+k,j+l]*kernel2[k+gap+1,l+gap+1]
}
}
}
}
return(out)
}

#convolution for lowpass
sigma <- ceiling(outsize)*strong
if (outsize==0) {sigma=1}
outlow0 = conv2d(d2arr0,kern)
outlow1 = conv2d(d2arr1,kern)

out0<-outlow0
out1<-outlow1

return(list(filt0=out0, filt1=out1))

}

utild1group <- function(arr0, arr1, radius=1, method='split1') {

#Date 24-02-15

}

if (class(arr0)=="matrix" | class(arr1)=="matrix") {
stop('! one/both arrays have >1 dims!')
}

if (method!='splitN') {

new_size0 <- floor(length(arr0)/groupi)
new_size1 <- floor(length(arr1)/groupi)
s0 <- c(1:new_size0)*0
s1 <- c(1:new_size1)*0

for (i in 1:new_size0) {
s0[i] = mean(arr0[((i-1)*groupi+1):(i*groupi)])
}
for (i in 1:new_size1) {
s1[i] = mean(arr1[((i-1)*groupi+1):(i*groupi)])
}

}

if (method=='splitN') {

new_size0 <- length(arr0)-groupi+1
new_size1 <- length(arr1)-groupi+1
s0 <- c(1:new_size0)*0
s1 <- c(1:new_size1)*0

for (i in 1:new_size0) {
s0temp <- c(1:groupi)*0
for (k in 1:groupi) {
#j = ((k-1)*(new_size0/radius)+i) # index of a new array
#shift = k+i-1
s0temp[k] = arr0[i+k-1]
}
s0[i] = mean(s0temp)
}

for (i in 1:new_size1) {
s1temp <- c(1:groupi)*0
for (k in 1:groupi) {
#j = ((k-1)*(new_size0/radius)+i) # index of a new array
#shift = k+i-1
s1temp[k] = arr1[i+k-1]
}
s1[i] = mean(s1temp)
}

}

out <- list(group0=s0, group1=s1)
return(out)

}

######
utild2group <- function(d2arr0, d2arr1, radius=1, method='split1') {

#Date 24-02-15

if (class(d2arr0)!="matrix" | class(d2arr1)!="matrix") {
stop('! one/both arrays have <2 dims!')
}
length(d2arr0[,1]),length(d2arr1[,1]))) {
}

if (method!='splitN') {

new_width0 <- floor(length(d2arr0[,1])/groupi)
new_width1 <- floor(length(d2arr1[,1])/groupi)
new_heigth0 <- floor(length(d2arr0[1,])/groupi)
new_heigth1 <- floor(length(d2arr1[1,])/groupi)

s0 <- array(0, c(new_width0, new_heigth0))
s1 <- array(0, c(new_width1, new_heigth1))

for (i in 1:new_width0) {
for (j in 1:new_heigth0) {
s0[i,j] = mean(d2arr0[((i-1)*groupi+1):(i*groupi),
((j-1)*groupi+1):(j*groupi)])
}
}
for (i in 1:new_width1) {
for (j in 1:new_heigth1) {
s1[i,j] = mean(d2arr1[((i-1)*groupi+1):(i*groupi),
((j-1)*groupi+1):(j*groupi)])
}
}

}

if (method=='splitN') {

new_width0 <- length(d2arr0[,1])-groupi+1
new_width1 <- length(d2arr1[,1])-groupi+1
new_heigth0 <- length(d2arr0[1,])-groupi+1
new_heigth1 <- length(d2arr1[1,])-groupi+1

s0 <- array(0, c(new_width0, new_heigth0))
s1 <- array(0, c(new_width1, new_heigth1))

for (i in 1:new_width0) {
for (j in 1:new_heigth0) {

s0temp <- array(0, c(groupi,groupi))

for (k in 1:groupi) {
for (m in 1:groupi) {
s0temp[k,m] = d2arr0[i+k-1,j+m-1]
}
}

s0[i,j] = mean(s0temp)

}
}

for (i in 1:new_width1) {
for (j in 1:new_heigth1) {

s1temp <- array(0, c(groupi,groupi))

for (k in 1:groupi) {
for (m in 1:groupi) {
s1temp[k,m] = d2arr1[i+k-1,j+m-1]
}
}
s1[i,j] = mean(s1temp)
}
}

}

out <- list(group0=s0, group1=s1)
return(out)

}
```

## Try the stheoreme package in your browser

Any scripts or data that you put into this service are public.

stheoreme documentation built on May 2, 2019, 9:33 a.m.