R/stheoreme.R

Defines functions utild2group utild1group utild2filt utild1filt utild2clean utild1clean utild2bin utild1bin plot.pvalign print.pvalign pvalign.default pvalign pvaligner plot.d2spec print.d2spec d2spec.default d2spec d2spectrum d2nat.d1nat plot.d1spec print.d1spec d1spec.default d1spec d1spectrum plot.d1nat print.d1nat d1nat.default d1nat d1native d1char.d1nat cxds.stheorem crit.stheorem print.stheorem stheorem.default stheorem stheoremme

Documented in crit.stheorem cxds.stheorem d1char.d1nat d1nat d1nat.default d1native d1spec d1spec.default d1spectrum d2nat.d1nat d2spec d2spec.default d2spectrum plot.d1nat plot.d1spec plot.d2spec plot.pvalign print.d1nat print.d1spec print.d2spec print.pvalign print.stheorem pvalign pvalign.default pvaligner stheorem stheorem.default stheoremme utild1bin utild1clean utild1filt utild1group utild2bin utild2clean utild2filt utild2group

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, 
                       method='rad', meansub=TRUE) {
  
  # 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
    ####
    
    # radial method grouping
  } else {
    method='rad'
    quadrantos0 <- array(0,c(pw_size0,(2*pw_size0+1)))
    for (i in 1:pw_size0) {
      for (j in 1:(2*pw_size0+1)) {
        quadrantos0[i,j]=(i^2+(j-pw_size0-1)^2) 
      }
    }
    #
    quadrantos1 <- array(1,c(pw_size1,(2*pw_size1+1)))
    for (i in 1:pw_size1) {
      for (j in 1:(2*pw_size1+1)) {
        quadrantos1[i,j]=(i^2+(j-pw_size1-1)^2) 
      }
    }
    
    #4 arrays to define with radial space
    radius0 <- c(0:pw_size0)
    radius1 <- c(0:pw_size1)
    freq0 <- 0.5*c(0:(pw_size0-1))/pw_size0
    freq1 <- 0.5*c(0:(pw_size1-1))/pw_size1
    
    #radial power grouped by radii
    pwR0<-rep(0,pw_size0)
    for (k in 1:pw_size0) {
      pwR0[k] = sum(pw0[which(quadrantos0<((radius0[k+1])^2) 
                              & quadrantos0>=((radius0[k])^2))])
    }
    #
    pwR1<-rep(0,pw_size1)
    for (k in 1:pw_size1) {
      pwR1[k] = sum(pw1[which(quadrantos1<((radius1[k+1])^2) 
                              & quadrantos1>=((radius1[k])^2))])
    }
    
    # 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, 
                  method='rad', meansub=TRUE) UseMethod("d2spec")
d2spec.default <- function(d2arr0, d2arr1, band=c(0,0), 
                           brks=0, method='rad', meansub=TRUE) {
  
  # 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")
  } else {cat("radial freqs ($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"}
  if (x$method=='rad') {xlab <- "Radial_frequency, 1/pixel"}
  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))
  }
  if (method == 'iaddsd') {
    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
  radius <- ceiling(outsize)
  if (outsize==0) {sigma=1}
  kern <- d1gausscoeffs(sigma, radius)
  outlow0 <- arr0
  outlow1 <- arr1
  if (outsize==0) {radius=1}
  outlow0[(radius+1):(length(arr0)-radius)] = 
    convolve(arr0,kern,type='f')
  outlow1[(radius+1):(length(arr1)-radius)] = 
    convolve(arr1,kern,type='f')
  outlow0[1:(radius)] = outlow0[radius+1] #border values left
  outlow0[(length(arr0)-radius+1):length(arr0)] = outlow0[(length(arr0)-radius)] # also border 
  outlow1[1:(radius)] = outlow1[radius+1] # also border
  outlow1[(length(arr1)-radius+1):length(arr1)] = outlow1[(length(arr1)-radius)] # also border 
  
  
  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
  radius <- ceiling(outsize)
  if (outsize==0) {sigma=1}
  kern <- d2gausscoeffs(sigma, radius)
  if (outsize==0) {radius=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 (radius<=0) {stop('!radius is invalid (must be > 0)!')}
  if (radius>0.4*min(length(arr0),length(arr1))) {
    stop('! radius is too high!')
  }
  
  groupi=floor(2*radius+1)
  
  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!')
  }
  if (radius<=0) {stop('!radius is invalid (must be > 0)!')}
  if (radius>0.4*min(length(d2arr0[1,]),length(d2arr1[1,]),
                     length(d2arr0[,1]),length(d2arr1[,1]))) {
    stop('! radius is too high!')
  }
  
  groupi=floor(2*radius+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.