text:
##
#Binary DIRECTED Networks
##

struct  ErgmDirBin1 <: Ergm #Bin stands for (Binary) Adjacency matrix
    "ERGM for directed networks with in and out degrees as statistics and
    possibility to have groups of nodes associated with a single pair of in out
    parameters. "
    obs:: Array{<:Real,1} #[InDegs;OutDegs] degrees of each node (even if groups are present)
    Par::Array{<:Real,1} # One parameter per group of nodes  UNRESTRICTED VERSION  [InPar;OutPar]
    groupsInds::Array{<:Real,1} #The index of the groups each node belongs to
    # two parameters per node (In and out)
 end
fooErgmDirBin1 =  ErgmDirBin1(ones(6),zeros(6),ones(3))
ErgmDirBin1(deg:: Array{<:Real,1}) =  ErgmDirBin1(deg,zeros(length(deg[:])) , Int.(1:round(Int,length(deg[:])/2)))

function identify(Model::ErgmDirBin1,parIO::Array{<:Real,1};idType ="equalIOsums" )
    "Given a vector of parameters, return the transformed vector that verifies
    an identification condition. Do not do anything if the model is identified."
    #set the first of the in parameters equal to one (Restricted version)
    N,parI,parO = splitVec(parIO)

    if idType == "equalIOsums"
        Δ = sum(parI[isfinite.(parI)]) - sum(parO[isfinite.(parO)])
        shift = Δ/(2N)
    elseif idType == "firstZero"
        shift = parI[1]
    end
    parIO = [ parI .- shift ;  parO .+ shift ]
    return parIO
  end

function DegSeq2graphDegSeq(Model::ErgmDirBin1,SeqIO::Array{<:Real,1})
    #From a degree sequence find another one that is certanly graphicable
    #The matrix associated with the starting degree sequence is very hubs oriented

    # the largest node catches as many links as it needs to fill all its strubs
    N, SeqI,SeqO = splitVec(SeqIO)
    sum(SeqI) == sum(SeqO) ? () : error()
    inds = Vector(1:N)
    StubsI_inds = Int.(sortrows([SeqI inds],rev=true))#decreasing sort
    StubsO_inds = Int.(sortrows([SeqO inds],rev=true))

    mat = zeros(Int64,N,N)
    for i in 1:N
       j = 1
       while (StubsI_inds[i,1]>0)&(j<N)
           ind_i = StubsI_inds[i,2]
           ind_j = StubsO_inds[j,2]

           if (mat[ind_i,ind_j] == 0)& (ind_i!=ind_j) &(StubsO_inds[j,1]>0)

               mat[ind_i,ind_j] = 1
               StubsI_inds[i,1]  -= 1
               StubsO_inds[j,1]  -= 1
           end
           j += 1
       end
    end
    ((sum(mat[Array{Bool}(eye(mat))])) != 0)  ? error() : ()
    graphSeq= [sumSq(mat,2); sumSq(mat,1)]
    return graphSeq
 end
function zeroDegParFun(Model::ErgmDirBin1,N::Int;degsIO=zeros(Int,10),
            parIO=zeros(Float64,10), method="AVGSPACING" ,smallProb = 1e-1)
    #define a number big enough to play the role of Inf for
    # purposes of sampling N(N-1) bernoulli rvs with prob 1/(1+exp(  bigNumb))
    # prob of sampling degree > 0 (N-1)/(1+exp(bigNumb)) < 1e-6
    if uppercase(method) == "INF"
        # zero degrees have -Inf fitnesses
        zeroDegPar = - Inf #log((N-1)/smallProb - 1)
    elseif uppercase(method) == "SMALL"
        # such that the probability of observing a link is very small
        zeroDegPar = -10
    elseif uppercase(method) == "AVGSPACING"
        if (sum(degsIO)==0)|(sum(parIO) ==0)
            error()
        end
        ~,degI ,degO = splitVec(degsIO)
        ~,parI ,parO = splitVec(parIO)
        imin = findmin(degI[degI.!=0])[2]
        omin = findmin(degO[degO.!=0])[2]
        # find a proxy of the one degree difference between fitnesses
        diffDegs = degI[degI.!=0] .- degI[degI.!=0][imin]
        avgParStepI = mean(((parI[degI.!=0] .- parI[degI.!=0][imin])[diffDegs.!=0]) ./ diffDegs[diffDegs.!=0])
        !isfinite(avgParStepI) ? println(((parI .- parI[imin])[diffDegs.!=0]) ./ diffDegs[diffDegs.!=0])error() : ()
        diffDegs = degO[degO.!=0] .- degO[degO.!=0][omin]
        avgParStepO = mean((parO[degO.!=0] .- parO[degO.!=0][omin])[diffDegs.!=0] ./ diffDegs[diffDegs.!=0])
        !isfinite(avgParStepO) ? println(((parI .- parI[imin])[diffDegs.!=0]) ./ diffDegs[diffDegs.!=0])error() : ()
        #subtract a multiple of this average step to fitnesses of smallest node
        zeroDegParI = parI[degI.!=0][imin] .- avgParStepI .* degI[degI.!=0][imin]
        zeroDegParO = parO[degO.!=0][omin] .- avgParStepO .* degO[degO.!=0][omin]
    end
    return zeroDegParI,zeroDegParO
 end
function linSpacedPar(Model::ErgmDirBin1,Nnodes::Int;Ngroups = Nnodes,deltaN::Int=3,graphConstr = true,degMax = Nnodes-2)
    "for a given static model return a widely spaced sequence of groups parameters"
    N=Nnodes
    #generate syntetic degree distributions, equally spaced, and estimate the model
    tmpdegsIO = round.([ Vector(range(deltaN,stop=degMax,length=N));Vector(range(deltaN,stop=degMax,length=N))])
    if graphConstr
        tmpdegsIO  =   DegSeq2graphDegSeq(fooErgmDirBin1,tmpdegsIO)
    end

    if Ngroups>N#if more groups than nodes are required throw an error
        error()
    else
        groupsInds = distributeAinVecN(Vector(1:Ngroups),Nnodes)
    end
    out,it,estMod = estimate(ErgmDirBin1(tmpdegsIO))
    return out,tmpdegsIO
 end
function samplSingMatCan(Model::ErgmDirBin1,expMat::Array{<:Real,2})
    """
    given the expected matrix sample one ranUtilities matrix from the corresponding pdf
    """
    N1,N2 = size(expMat)
    N1!=N2 ? error() : N=N1
    out = zeros(Int8,N,N)
    #display((maximum(expMat),minimum(expMat)))
    for c = 1:N
        for r=1:N
            if r!=c
                out[r,c] = rand(Bernoulli(expMat[r,c]))
            end
        end
    end
    return out
 end
function expMatrix2(Model::ErgmDirBin1, parNodesIO::Array{<:Real,1}  )
    "Given the vector of model parameters (groups), return the product matrix(often useful
    in likelihood computation) and the expected matrix"
    N,parI,parO = splitVec(parNodesIO)
    parMat =  Utilities.putZeroDiag_no_mut(exp.(parI  .+ parO') )
    expMat = parMat ./ (1 .+ parMat)
    return parMat,expMat
 end
expMatrix(Model::ErgmDirBin1, parNodesIO::Array{<:Real,1}  ) = expMatrix2(Model, parNodesIO )[2]# return only the expected matrix given the parameters
function expValStatsFromMat(Model::ErgmDirBin1,expMat::Array{<:Real,2}  )
    "Expected value of the statistic for a model, given the expected matrix"
    expDeg = [sumSq(expMat,2);sumSq(expMat,1)]
    expValStats = expDeg
    return expValStats
 end
expValStats(Model::ErgmDirBin1, parNodesIO::Array{<:Real,1} ) = expValStatsFromMat(Model, expMatrix(Model,parNodesIO )) # if only the paramters are given as input
function firstOrderCond(Model::ErgmDirBin1;degIO::Array{<:Real,1} = Model.obs,
                    degGroups::Array{<:Real,1} = zeros(4), parGroupsIO::Array{<:Real,1}=Model.Par, groupsInds::Array{<:Real,1} = Model.groupsInds )
    "Given the model, the degree sequence, the parameters and groups assignements,
    return the First order conditions. Gradient of the loglikelihood, or system
    of differences between expected and observed values of the statistics (either
    degrees(In and Out) of single nodes or of a group). "
    N,degI,degO = splitVec(degIO)
    groupsNames = unique(groupsInds)
    NG = length(groupsNames)
    NGcheck,parI,parO = splitVec(parGroupsIO)
    NG==NGcheck ?  () : error()
    if sum(degGroups)==0
        degGroups = zeros(2NG) ; for i =1:N degGroups[groupsInds[i]] += degI[i];degGroups[NG + groupsInds[i]] += degO[i]; end
    end
    parNodesIO = [parI[groupsInds];parO[groupsInds]]
    N,expDegI,expDegO = splitVec(expValStats(Model,parNodesIO))
    expGroupsDeg = zeros(2NG) ; for i =1:N expGroupsDeg[groupsInds[i]] += expDegI[i];expGroupsDeg[NG + groupsInds[i]] += expDegO[i]; end
    #println(groupsInds)
    #println(degGroups)
    #println(expGroupsDeg)
    foc = degGroups .- expGroupsDeg
    return foc
 end

function estimateSlow(Model::ErgmDirBin1; degIO::Array{<:Real,1} = Model.obs,
                    groupsInds = Int.(1:round(length(degIO)/2)), targetErr::Real=targetErrValStaticNets,
                    startVals =zeros(2))
    "Given model type, observations of the statistics and groups assignments,
    estimate the parameters."
    N,degI,degO = splitVec(degIO)
    groupsNames = unique(groupsInds)
    NG = length(groupsNames)
    tmpDic = countmap(groupsInds)
    #number of nodes per group
    NgroupsMemb = sortrows([[i for i in keys(tmpDic)] [  tmpDic[i] for i in keys(tmpDic) ] ])
    # agrregate degree of each group
    degGroups = zeros(2NG) ; for i =1:N degGroups[groupsInds[i]] += degI[i]; degGroups[NG + groupsInds[i]] += degO[i]; end
    ~,degGroupsI , degGroupsO = splitVec(degGroups)

    #sort the groups for degrees
    indsDegGroupsOrderI = Int.(sortrows([degGroupsI Vector(1:NG)], rev = true)[:,2])
    indsDegGroupsOrderO = Int.(sortrows([degGroupsO Vector(1:NG)], rev = true)[:,2])

    NzerI = sum(degGroupsI.==0)
    NzerO = sum(degGroupsO.==0)
    (sum(degGroupsI) - sum(degGroupsO) )< 10^Base.eps() ?  L =sum(degGroupsO) : error()
    LperLink  = L/(N*N-N)
    # The parameters inside this function are Restricted to be in [0,1] hence
    # they need to be converted (log(par)) in order to be used outside

    #uniform parameter as starting value
    unifstartval =  sqrt((LperLink/((1-LperLink))));
    if NG==1
        outParUN = [ 0 ; log(unifstartval)]
        outMod = ErgmDirBin1(Model.obs,outParUN,groupsInds)
        i = 0
    else
        parI = unifstartval*ones(NG)
        parO = unifstartval*ones(NG)
        #bigNumber = 1e3 ./ targetErr
        eps =  targetErr*1e-12 #10*Base.eps()

        i=0
        maxIt = 25
        relErr = 1
        tol_x = 10*Base.eps()
        sum(startVals)==0 ? () : (parI=startVals[1:N];parO=startVals[1+N:end];)
        function tmpFindZero(objFun::Function,brak::Array{<:Real,1},start::Real)
            #    wrapper find_zero options
            #println((brak[1],brak[2]))
            #println((objFun(brak[1]),objFun(brak[2])))
            #println()
            out = find_zero(objFun,(brak[1],brak[2]),Bisection())
            #println((out,objFun(out)))
            return out
        end
        while (i<maxIt)&(relErr>targetErr)
            parOMax = maximum(parO)
            parOMin = minimum(parO[parO.!=0])
            for nInd=1:NG
                ng = indsDegGroupsOrderI[nInd]
            #    println((1,n))
                #put to zero each component of the gradient1
                function objfunI_ng(parI_ng::Real)
                    parVecOthGroups = parI_ng.*parO[1:end .!=ng]
                    parO_ng = parO[ng]
                    expGroupDegI = NgroupsMemb[ng,2]*( sum( NgroupsMemb[1:end .!=ng ,2].*( parVecOthGroups ./ (1 + parVecOthGroups)) )    + (NgroupsMemb[ng,2]-1)* (parI_ng*parO_ng)/(1+(parI_ng*parO_ng)) )
                    return  degGroupsI[ng] - expGroupDegI

                end
                #find zero of given component

                if degGroupsI[ng] == NgroupsMemb[ng,2]*(N-1-NzerO)
                    #if a group has the maximum possible in degree set par to a
                    #very high value
                    # println((1,ng))
                     parI[ng] = (1-eps)/(eps* parOMin)
                    # objfunI_ngMinEps(parI_ng::Real) = (objfunI_ng(parI_ng::Real) - eps)
                    # tmpBnd = bigNumber #(degGroupsO[ng]/(NgroupsMemb[ng,2]*(N-1-NzerI)) )/(1-degGroupsO[ng]/(NgroupsMemb[ng,2]*(N-1-NzerI)) )
                    # brakets = [ 0,2* tmpBnd/parOMin]
                    # parI[ng] = tmpFindZero(objfunI_ngMinEps,brakets,parI[ng])
                elseif degGroupsI[ng]==0
                    #error("Remove zero col and rows")
                    parI[ng] = 0
                else
                    tmpBnd = (degGroupsI[ng]/(NgroupsMemb[ng,2]*(N-1-NzerO)) )/(1-degGroupsI[ng]/(NgroupsMemb[ng,2]*(N-1-NzerO)) )
                    #println((tmpBnd,degGroupsI[ng],NgroupsMemb[ng,2],NzerO))
                    brakets = [ 0,2* tmpBnd/parOMin]
                    parI[ng] = tmpFindZero(objfunI_ng,brakets,parI[ng])
                end
            end
            #parUn = log.([parI;parO])
            #err = firstOrderCond(Model;deg = deg,par = parUn , groupsInds = groupsInds)
            #println(err)
            parIMax = maximum(parI)
            parIMin = minimum(parI[parI.!=0])
            for nInd=1:NG
                ng = indsDegGroupsOrderO[nInd]
                #println((2,n))
                #put to zero each component of the gradient1
                function objfunO_ng(parO_ng::Real)
                    parVecOthGroups = parO_ng.*parI[1:end .!=ng]
                    parI_ng = parI[ng]
                    expGroupDegO = NgroupsMemb[ng,2]*( sum( NgroupsMemb[1:end .!=ng ,2].*( parVecOthGroups ./ (1 + parVecOthGroups)) )    + (NgroupsMemb[ng,2]-1)* (parI_ng*parO_ng)/(1+(parI_ng*parO_ng)) )
                    return  degGroupsO[ng] - expGroupDegO
                end
                #find zero of given component
                if degGroupsO[ng] == NgroupsMemb[ng,2]*(N-1-NzerI)
                    #if a group has the maximum possible out degree set par to a
                    #very high value
                    # println((2,ng))
                    parO[ng] = (1-eps)/(eps* parIMin)
                    # objfunO_ngMinEps(parO_ng::Real) = (objfunO_ng(parO_ng::Real) - eps)
                    # tmpBnd = bigNumber #(degGroupsO[ng]/(NgroupsMemb[ng,2]*(N-1-NzerI)) )/(1-degGroupsO[ng]/(NgroupsMemb[ng,2]*(N-1-NzerI)) )
                    # brakets = [ 0,2* tmpBnd/parIMin]
                    # parO[ng] = tmpFindZero(objfunO_ngMinEps,brakets,parO[ng])
                elseif degGroupsO[ng]==0
                    #error("Remove zero col and rows")
                    parO[ng] = 0
                else
                    tmpBnd = (degGroupsO[ng]/(NgroupsMemb[ng,2]*(N-1-NzerI)) )/(1-degGroupsO[ng]/(NgroupsMemb[ng,2]*(N-1-NzerI)) )
                    brakets = [ 0,2* tmpBnd/parIMin]
                    parO[ng] = tmpFindZero(objfunO_ng,brakets,parO[ng])
                end

            end
            parReIO = [parI;parO]
            parUnIO = -bigConstVal.*ones(parReIO)
            tmpInds = parReIO.>0
            #show(tmpInds)
            parUnIO[tmpInds] = log.(parReIO[tmpInds])
            err = firstOrderCond(Model;degIO = degIO,parGroupsIO = parUnIO , groupsInds = groupsInds)
        #    tmpMat = par.*par';expMat = putZeroDiag!(tmpMat ./ (1+tmpMat)) ;expDeg = sumSq(expMat,2)
            relErrVec = abs.(err ./ degGroups)
            #println(relErrVec)
            relErr =  maximum(relErrVec[degGroups.!=0])
            isnan(relErr) ? error() : ()
            # display(err')
             #println((mean(relErrVec[deg.!=0]),relErr))
            i+=1
        end
        #println([parI parO])
        #i<maxIt ? () : (println(i);error())
        parIO = [parI;parO]
        parUN = log.(parIO)
    end
    return parUN, i
 end
function estimate(model::ErgmDirBin1; degIO::Array{<:Real,1} = model.obs,
                    targetErr::Real=targetErrValStaticNets,identIter = false)
    ## Write and test a new way to estimate DirBOIn1 following Chatterjee, Diaconis and  Sly
    N,degI,degO = splitVec(degIO)
    ldegI = log.(degI);ldegO = log.(degO)
    (sum(degI) - sum(degO) )< 10^Base.eps() ?  L =sum(degO) : error()
    LperLink  = L/(N*N-N)
    unifstartval =  0.5*log((LperLink/((1-LperLink))))
    #inds = (ldegI.==0)
    #function that defines the iteration
    function phiFunc(vParUn::Array{<:Real,1})
        ~,vParReI,vParReO = splitVec(exp.(vParUn))
        matI = putZeroDiag(1 ./ (vParReI  .+ (1 ./ vParReO)' ))
        matO = putZeroDiag(1 ./ ( (1 ./ vParReI)  .+ vParReO'  ))
        outPhiI = (ldegI - log.(sumSq(matI,2)))
        outPhiO = (ldegO - log.(sumSq(matO,1)))
        #println(vParReI[ldegI.==0])
        #println(log.(sumSq(matI,2))[inds])
        #println(matI[inds,:][1,1:20])
        #println(matI[inds,:][2,1:20])
        return [outPhiI;outPhiO]
     end
    function checkFit(vParUn::Vector{<:Real})
        expMat = expMatrix(model,vParUn)
        nnzInds = degIO.!=0
        errsIO = [sumSq(expMat,2) - degI; sumSq(expMat,1) - degO][nnzInds]
        relErrsIO = abs.(errsIO) ./ (degIO[nnzInds])
        return prod(relErrsIO .< targetErr)
     end

    maxIt = 20000
    it=0
    it2 = 0
    vParUn = [ldegI;ldegO]# 0*unifstartval.*ones(2N)
    #starting from an identified set of par. and then identify at each step
    #amounts to doing an identified step. i.e. moving summing an identified
    #vector
    identIter ? vParUn =  identify(model,vParUn) : ()
    while (!checkFit(vParUn)) & (it<maxIt)
        vParUn = phiFunc(vParUn); it+=1
        identIter ? vParUn =  identify(model,vParUn) : ()
        #println(vParUn[degIO.==1][1:5])
    end
    #if did not obtain the required precision try with the slower method
    if it==maxIt
        outSlow = estimateSlow(model; degIO = degIO,targetErr = targetErr,
                            startVals =exp.(vParUn))

        vParUn = identify(model,outSlow[1])
        it2 = outSlow[2]
    end
    #zeroInds = degIO.==0
    #vParUn[zeroInds] = bigNegParFun(model,N)
    return vParUn,it,it2
 end
 #

function sampl(Model::ErgmDirBin1,Nsample::Int;  parGroupsIO::Array{<:Real,1}=Model.Par, groupsInds = Model.groupsInds )
    NGcheck,parI,parO = splitVec(parGroupsIO)
    parNodesIO = [parI[groupsInds];parO[groupsInds]]
    expMat = expMatrix(Model, parNodesIO)
    N = length(expMat[:,1])
    SampleMats =zeros(Int8,N,N,Nsample)
    for s=1:Nsample
        SampleMats[:,:,s] = samplSingMatCan(Model,expMat)
    end
    return SampleMats
 end

# Snapshot sequence Functions-----------------------------------------------------------------

struct SnapSeqNetDirBin1
    obsT:: Array{<:Real,2}  # In degrees and out Degrees for each t . The last
                            # dimension is always time. 2NxT
    Par::Array{<:Real,2} # In and Out Parameters [InW_T;outW_T] for each t
 end
SnapSeqNetDirBin1(obs::Array{<:Real,2}) = SnapSeqNetDirBin1(obs,zero(obs))
fooSnapSeqNetDirBin1 = SnapSeqNetDirBin1(zeros(2,3))

function estimate(Model::SnapSeqNetDirBin1; degsIO_T::Array{<:Real,2}=Model.obsT,
                    targetErr::Real=1e-4,identPost=false,identIter=false,
                    zeroDegFit="small")
   N2,T = size(degsIO_T);N = round(Int,N2/2)
   UnParT = zeros(2N,T)
   prog = 0
   for t = 1:T
       degs_t = degsIO_T[:,t]
      # println(degs_t)
       UnPar_t , Niter1,Niter2  = estimate(fooErgmDirBin1; degIO = degs_t ,
                                    targetErr =  targetErr,identIter=identIter)
       zeroParI,zeroParO = zeroDegParFun(fooErgmDirBin1,N;
                                        degsIO=degsIO_T[:,t],parIO = UnPar_t )
      # println((zeroParI,zeroParO))
       !isfinite(zeroParI) ? error() : ()
       !isfinite(zeroParO) ? error() : ()
       UnPar_t[findall(degs_t[1:N].==0)] .= zeroParI
       UnPar_t[N .+ findall(degs_t[1+N:2N].==0)] .= zeroParO
       identPost ? UnParT[:,t] = identify(fooErgmDirBin1,UnPar_t) :
                                    UnParT[:,t] = UnPar_t
       #round(t/T,1)>prog ?  (prog = round(t/T,1);println(prog)) : ()

        #println((t,Niter1,Niter2))
   end
   sum(.!isfinite.(UnParT))>0 ? error() : ()

 return UnParT
 end

function sampl(Model::SnapSeqNetDirBin1, vParUnIO_T::Array{<:Real,2}; Nsample::Int=1  )
    N2,T = size(vParUnIO_T);N = round(Int,N2/2)
    SampleMats =zeros(Int8,Nsample,N,N,T)
    for s=1:Nsample
        for t = 1:T
            SampleMats[s,:,:,t] = samplSingMatCan(fooErgmDirBin1,expMatrix(fooErgmDirBin1,vParUnIO_T[:,t]))
        end
    end
    Nsample == 1  ?  (return squeeze(SampleMats,1)) : (return SampleMats )
  end

 #---------------------------------Binary DIRECTED Networks with predetermined COVARIATES

 #
 # struct  ErgmDirBin1X0 <: ErgmBin #Bin stands for (Binary) Adjacency matrix
 #                                        # X for regressors and 0(lower level pox) for link specific parameters
 #    "ERGM with in and out node specific parameters and dependence on link
 #    specific regressors trough link specific parameters (0)  "
 #    Amats:: BitArray{3}   #Different realizations of the adjqcency matrix Needed!
 #                        #Because the covariates parameters are link specific. Time
 #                        #is the last dimension
 #    regressors0:: Array{<:Real,3}
 #    parDegs:: Array{<:Real,2} # One parameter per node  UNRESTRICTED VERSION  [InPar OutPar]
 #    parRegs:: Array{<:Real,2} #one parameter for each link associated with each regressor
 # end
 # fooErgmDirBin1X0 =  ErgmDirBin1X0(tmpArrayAdjmatrEqual(6,10),ones(6,6,10),zeros(6,2),ones(6,6))
 # ErgmDirBin1X0(Amats:: BitArray{3},regs0::Array{<:Real,3}) =  (N =length(Amats[:,1,1]); ErgmDirBin1X0(Amats,regs0, zeros(N ,2) , zeros(N,N) ))

function defineConstDegs(degsIO_T;thVarConst = 0.005)
        N2 = length(degsIO_T[:,1])
        isConst =  [StatsBase.trimvar(degsIO_T[i,:];prop = 0.001)  for i =1:N2].< thVarConst

        ~,isConI,isConO = splitVec(isConst)
        @show sum(isConst)
        return sum(isConst),isConI,isConO
    end
#Some useful wrap functions
function estManyAR(obs::Array{<:Real,2};p::Int=1,isCon::BitArray{1} = falses(length(obs[:,1])))
    #obs is a matrix with a time series for each row (time is last dimension)
    D,T = size(obs)
    parEst = zeros(D,p+1)
    arp = AReg.ARp(ones(p+1),.1,[])

    for d=1:D
        if !isCon[d]
            parEst[d,:] = AReg.fitARp(obs[d,:],arp).estimates
        else
            #if the ts is costant, use only the unconditional mean with zero persistence
            parEst[d,1] = mean(obs[d,:])
        end
    end
    return parEst
 end
function oneStepForAR1(obs_t::Float64,AR1Par::Vector{<:Real};steps::Int = 1)
    #given w and β AR1 parameters, compute the one step ahead forecast from obs_t
    if steps==1
        obstpSteps = AR1Par[1] + AR1Par[2] * obs_t
    else
        obstpSteps = AR1Par[1] * sum(AR1Par[2].^Vector(0:steps-1)) + AR1Par[2].^steps * obs_t
    end
    return obstpSteps
 end
#Estimate AR1 on the training sample
function forecastEvalAR1Net(allFit::Array{Float64,2}, obsNet_T::BitArray{3},Ttrain::Int;mat2rem=falses(obsNet_T[:,:,1]) )
    N2,T = size(allFit)
    Ttest = T-Ttrain
    trainFit = allFit[:,1:Ttrain]
    degsIO_T = [sumSq(obsNet_T,2);sumSq(obsNet_T,1)]
    Nconst,isConIn,isConOut = defineConstDegs(degsIO_T[:,1:Ttrain],thVarConst =0.005 )
    estAR1= estManyAR(trainFit;isCon = [isConIn;isConOut])
    ## One step ahead forecasts for each fitness for each time in Test sample
    foreFit = zeros(allFit)
    for n=1:N2 ,t=1:T-1
        #compute forecasting for all links, diagonal will be disregarded after
        foreFit[n,t+1] = oneStepForAR1(allFit[n,t],estAR1[n,:])
    end
    tp,fp,AUC,realVals,foreVals =    nowCastEvalFitNet(foreFit[:,Ttrain+1:end],obsNet_T[:,:,Ttrain+1:end];mat2rem=mat2rem)
    return tp,fp,AUC,realVals,foreVals
 end
function nowCastEvalFitNet(allFit::Array{Float64,2}, obsNet_T::BitArray{3};
            expMat_T :: Array{Float64,3}=zeros(10,10,50),mat2rem=falses(obsNet_T[:,:,1]), shift::Int=0,plotFlag = true )

    sum(expMat_T) == 0  ? inputExpMat =false : inputExpMat=true
    N2,T = size(allFit)
    inputExpMat  ? T = size(expMat_T)[3] : ()
    println((N2,T))
    #println(size(expMat_T))
    obsNet_Test_T= obsNet_T

    # se l'expected matrices sono date come input nel seguto non calcolarle
    #ricorda che se expMat_T vengono da multi steps ahead forecasting,non va usato
    # il temporal shift
    inputExpMat ? (shift!=0 ? error : ()) : ()
    ##use input fitnesses
    testFit = allFit
    #disregard diagonal elements and those that are in mat2rem
    noDiagIndnnc = putZeroDiag((.!mat2rem)    )
    @show Nlinksnnc =sum(noDiagIndnnc)
    if shift>=0
            println(Nlinksnnc*(T-shift))
            println(T)
        foreVals = 100 * ones(Nlinksnnc*(T-shift))
        realVals = falses( Nlinksnnc*(T-shift))
        lastInd=1
        for t=1:T-shift
            adjMat =  obsNet_Test_T[:,:,t+shift]
            if inputExpMat
                expMat =  expMat_T[:,:,t]
                #println(size(expMat))
            else
                expMat = expMatrix(fooErgmDirBin1,testFit[:,t])
            end
            foreVals[lastInd:lastInd+Nlinksnnc-1] = expMat[noDiagIndnnc]
            realVals[lastInd:lastInd+Nlinksnnc-1] = adjMat[noDiagIndnnc]

            lastInd += Nlinksnnc
        end
    else
        foreVals = 100 * ones(Nlinksnnc*(T+shift))
        realVals = falses( Nlinksnnc*(T+shift))
        lastInd=1
        for t=1:T+shift
            adjMat =  obsNet_Test_T[:,:,t]
            if inputExpMat
                expMat =  expMat_T[:,:,t]
            else
                expMat = expMatrix(fooErgmDirBin1,testFit[:,t-shift])
            end

            foreVals[lastInd:lastInd+Nlinksnnc-1] = expMat[noDiagIndnnc]
            realVals[lastInd:lastInd+Nlinksnnc-1] = adjMat[noDiagIndnnc]

            lastInd += Nlinksnnc
        end
    end
 # println(realVals)
    tp,fp,AUC = rocCurve(realVals,foreVals;plotFlag = plotFlag)
    return tp,fp,AUC,realVals,foreVals
 end

function forecastEvalAR1Net(      obsMat_allT::BitArray{3},Ttrain::Int;identPost = false,identIter=true,mat2rem = falses(obsMat_allT[:,:,1]) )
        #if only the matrices obs are given then estimate the one fit model for each t
        degsIO_T = [sumSq(obsMat_allT,2);sumSq(obsMat_allT,1)]
        snapModData = SnapSeqNetDirBin1(degsIO_T)
        allFit =  estimate(snapModData; identPost = identPost,identIter=identIter)
        return   forecastEvalAR1Net(allFit, obsMat_allT,Ttrain;mat2rem = mat2rem)

end

function dgpFitVarN(N::Int,degMin::Int,degMax::Int;exponent = 1)

    # exponent determines the convexity of the distributions of degrees between
    # deg min and deg max. increasing exponent --> decreasing density
    # if we let deg max grow with N than
    #generate a degree distribution
    tmp = (range(degMin,stop=degMax,length=N)).^exponent

    degIOUncMeans = repeat(round.(Int,Vector(tmp.*(degMax/tmp[end]) )),2)
    degIOUncMeans[degIOUncMeans.<degMin] = degMin
    meanFit ,~ = StaticNets.estimate(StaticNets.ErgmDirBin1(degIOUncMeans),degIO = degIOUncMeans)
    expDegs = StaticNets.expValStats(fooErgmDirBin1,meanFit)
    return meanFit, expDegs
end

function dgpDynamic(model::ErgmDirBin1,mode,N::Int,T::Int;NTV = 2,degIOUncMeans::Array{Float64,1} =zeros(10),
                    degb = [10, 40])
    if sum(degIOUncMeans)==0 # if the unconditional means of degrees are not Given then fix them
        degIOUncMeans = range(1,stop=10,length=N)
        error()
    else
        N = round(Int,length(degIOUncMeans)/2)

    end
    #data unconditional mean definisci ampiezza della variazione, che sarà ampiezza seno e varianza AR1

    ## generate the dynamical fitnesses

    #println((N,NTV))
    indsTV = round.(Int,range(2, stop=N, length=NTV))
    indLogTV = falses(N);indLogTV[indsTV] .= true; indLogTV = repeat(indLogTV,2)

    meanFit ,~ = StaticNets.estimate(StaticNets.ErgmDirBin1(degIOUncMeans),degIO = degIOUncMeans)
    #define the oscillations bounds for the fitnesses
    periodOsc = round(Int,T/2)

    # all fits start at the same point with different fases and amplitudes
    #define a variation
    ~,meanfitI,meanfitO = splitVec(meanFit)
    #shuffle ranUtilitiesly the fases of the sines but fix the ranUtilities seed for reproducibility

    #solve for x :deglb =  (N-1) * 1/(1+exp( - x - meanfitO))
    fitbI = - mean(meanfitO) .- log.( (N-1) ./ degb .-1)
    fitbO = - mean(meanfitI) .- log.( (N-1) ./ degb .-1)
    ampI = ones(NTV).*Vector(range(0.4,stop=1, length=NTV).*(diff(fitbI)))#(diff(fitbI))#
    ampO = ones(NTV).* Vector(range(0.4,stop=1,length=NTV).*(diff(fitbO))) # (diff(fitbO)) #
    #center the dynamical parameters at the middle of the bounds
    meanFit[indsTV] =  ampI ./ 2  .+ fitbI[1]
    meanFit[N.+indsTV] = ampO ./ 2  .+ fitbO[1]

    times = Array{Float64,1}(1:T)
    dynFit = repeat(meanFit,1,T)

    if uppercase(mode) ==  uppercase("sin")
        println(NTV)
        phases =  rand( MersenneTwister(1234),2NTV)*((2pi))
        for n=1:NTV
            dynFit[indsTV[n],:]   = meanFit[indsTV[n]]  .+ ampI[n]*sin.( times.*((2pi/periodOsc)) .+ phases[n] )
            dynFit[N+indsTV[n],:] = meanFit[N+indsTV[n]]  .+ ampO[n]*sin.( times.*((2pi/periodOsc)) .+ phases[n+NTV] )
        end
    elseif uppercase(mode) ==  uppercase("const")
        for n=1:NTV
            dynFit[indsTV[n],:]   = meanFit[indsTV[n]] .*ones(T)
            dynFit[N+indsTV[n],:] = meanFit[N+indsTV[n]]  .*ones(T)
        end
    elseif uppercase(mode) ==  uppercase("steps")
        Nsteps = 2
        phases =  rand( MersenneTwister(1234),1:T,2NTV)

        for n=1:NTV
            minPar = meanFit[indsTV[n] ] .- ampI[n]/2
            maxPar = meanFit[indsTV[n]]  .+ ampI[n]/2
            #steps between min and max values
            dynFit[indsTV[n],:]  =  randSteps(minPar,maxPar,Nsteps,T)
            maxPar = meanFit[indsTV[n]+N] + ampO[n]/2
            #steps between min and max values
            dynFit[indsTV[n]+N,:]  = randSteps(minPar,maxPar,Nsteps,T)
            #circyularly shift time index by the ranUtilities phase of each
            if false
                #circyularly shift time index by the ranUtilities phase of each
                dynFit[indsTV[n],:]  = circshift(dynFit[indsTV[n],:]  ,phases[n])
                dynFit[indsTV[n]+N,:]  = circshift(dynFit[indsTV[n]+N,:],phases[n])
            end
            minPar = meanFit[indsTV[n]+ N] - ampO[n]/2
        end

    elseif uppercase(mode) == uppercase("ar1")

        sigma = 0.01
         B = 0.999 #rand([0.7,0.99])

        for n=1:NTV
            scalType = "uniform"
            minPar = meanFit[indsTV[n] ] - ampI[n]/2
            maxPar = meanFit[indsTV[n]] + ampI[n]/2

            dynFit[indsTV[n],:]  = dgpAR(meanFit[indsTV[n]],B,sigma,T,minMax=[minPar,maxPar];scaling = scalType)
            #println((meanFit[indsTV[n]]  ,[minPar,maxPar]) )
            minPar = meanFit[indsTV[n]+ N] - ampO[n]/2
            maxPar = meanFit[indsTV[n]+N] + ampO[n]/2

            dynFit[indsTV[n]+N,:]   = dgpAR(meanFit[indsTV[n]+N],B,sigma,T,minMax=[minPar,maxPar];scaling = scalType)

        end

    end

    #

    for t=1:T dynFit[:,t] = identify(model,dynFit[:,t]) end
    return dynFit, indLogTV
end


text:
function sample(assembly, chr, pos, temp_len)
    start = max(1, pos - div(temp_len, 2))
    if !haskey(assembly, chr)
        error("The chosen assembly doesn't have $chr")
    end
    start = min(start, length(assembly[chr]) - temp_len)
    seq = uppercase(assembly[chr].sequence.seq[start:start+temp_len-1])
    return dna(seq), start
end

ext: jl


text:
# hartley space
# @manipulate for kb in 0:4, ke in 1:5, dk in 0.4:0.2:1, p in 0:0.005:1, shape in [:square, :circle]
#     plothartleysubspace(hartleysubspace(kbegin=kb,kend=ke,dk=dk,phase=p,shape=shape),floor(Int,ke/dk),dk,color=:coolwarm)
# end

hs = hartleysubspace(kbegin=0.2,kend=6.6,dk=0.2,addhalfcycle=true)
hgs = map(i -> begin
          ss = cas2sin(i...)
          grating(θ = ss.θ, sf = ss.f, phase = ss.phase, size = (5, 5), ppd = 30)
          end, hs)
# mean of gratings[0, 1] of a hartley subspace should be a uniform gray
minv,maxv = extrema(reduce((i,j)->i.+j,hgs)/length(hs))
@test minv≈maxv≈0.5

hs = hartleysubspace(kbegin=0.2,kend=6.6,dk=0.2,addhalfcycle=true)
hgs = map(i -> begin
          ss = cas2sin(i...)
          grating(θ = ss.θ, sf = ss.f, phase = ss.phase, size = (2, 2), ppd = 10)
          end, hs)
Xᵀ = mapreduce(vec,hcat,hgs)
X = Xᵀ'
a=Xᵀ*X
b=inv(Xᵀ*X)

ggs = [randn(30,30) for _ in 1:10000]
Xᵀ = mapreduce(vec,hcat,ggs)
X = Xᵀ'
a=Xᵀ*X
b=inv(Xᵀ*X)

heatmap(a,yflip=true,aspect_ratio=1,frame=:none,color=:fire)

## grating
ori=0.5π;sf=2;ppd=30
img = grating(θ=ori,sf=sf,ppd=ppd)
img .+= rand(size(img)...)/3
heatmap(img,yflip=true,aspect_ratio=:equal,frame=:none)
ps,f1,f2 = powerspectrum2(img,ppd,freqrange=[-6,6])
heatmap(f2,f1,ps,aspect_ratio=:equal,frame=:grid)
eori,esf = f1orisf(ps,f1,f2)
@test sf≈esf
@test ori≈eori

## image resize
img = grating(θ=0.25π,sf=6,ppd=60,size=(3,3))
heatmap(img,yflip=true,aspect_ratio=:equal,frame=:none)

imgra = imresize_antialiasing(img,(64,64))
heatmap(imgra,yflip=true,aspect_ratio=:equal,frame=:none)

## image mask and blend
@test alphablend(Gray(0.7),Gray(0.5)) == Gray(0.7)
@test alphablend(GrayA(0.7,0),Gray(0.5)) == Gray(0.5)
@test alphablend(RGBA(0.5,1,1,1),RGB(0.6,0.3,0.2)) == RGB(0.5,1,1)

@test alphablend(10,5,1) == 10
@test alphablend(fill(10,10,10),fill(5,10,10),0.5) == fill(7.5,10,10)

alphablend.(alphamask(Gray.(img)).y,Gray(0.5))
alphablend.(alphamask(Gray.(img),masktype="Gaussian",sigma=0.4).y,Gray(0.5))
alphablend.(alphamask(Gray.(img),masktype="DiskFade",sigma=15,radius=0.45).y,Gray(0.5))

## ImagerIO
readrawim_Mono8(joinpath(@__DIR__,"Mono8_100.Raw"),100,100)
readrawim_Mono12Packed(joinpath(@__DIR__,"Mono12Packed_100.Raw"),100,100)


text:
# # Freely decaying barotropic turbulence on a latitude-longitude strip

using Oceananigans
using Oceananigans.Grids

using Oceananigans.BoundaryConditions: fill_halo_regions!

using Oceananigans.Coriolis: HydrostaticSphericalCoriolis

using Oceananigans.Models.HydrostaticFreeSurfaceModels:
    HydrostaticFreeSurfaceModel,
    VerticalVorticityField,
    VectorInvariant,
    ExplicitFreeSurface,
    ImplicitFreeSurface

using Oceananigans.TurbulenceClosures

using Oceananigans.Utils: prettytime, hours, day, days, years, year
using Oceananigans.OutputWriters: JLD2OutputWriter, TimeInterval, IterationInterval
using Oceananigans.AbstractOperations: KernelFunctionOperation

using Statistics
using JLD2
using Printf
using CUDA

#####
##### Grid
#####

precompute = true

latitude = (-80, 80)
Δφ = latitude[2] - latitude[1]

resolution = 1/3 # degree
Nx = round(Int, 360 / resolution)
Ny = round(Int, Δφ / resolution)

# A spherical domain
@show grid = LatitudeLongitudeGrid(GPU(), size = (Nx, Ny, 1),
                                   longitude = (-180, 180),
                                   latitude = latitude,
                                   halo = (2, 2, 2),
                                   z = (-100, 0),
                                   precompute_metrics = precompute)

#####
##### Physics and model setup
#####

free_surface = ExplicitFreeSurface(gravitational_acceleration=1.0)

CUDA.allowscalar(true)

equatorial_Δx = grid.radius * deg2rad(mean(grid.Δλᶜᵃᵃ))
diffusive_time_scale = 120days

@show const νh₂ =        equatorial_Δx^2 / diffusive_time_scale
@show const νh₄ = 1e-5 * equatorial_Δx^4 / diffusive_time_scale

closure = ScalarBiharmonicDiffusivity(ν=νh₄)

coriolis = HydrostaticSphericalCoriolis()
Ω = coriolis.rotation_rate / 20
coriolis = HydrostaticSphericalCoriolis(rotation_rate=Ω)

model = HydrostaticFreeSurfaceModel(grid = grid,
                                    momentum_advection = VectorInvariant(),
                                    free_surface = free_surface,
                                    coriolis = coriolis,
                                    tracers = nothing,
                                    buoyancy = nothing,
                                    closure = closure)

#####
##### Initial condition: two streamfunction
#####

# Random noise
ψ★ = Field(Face, Face, Center, model.architecture, model.grid)
set!(ψ★, (x, y, z) -> rand())
fill_halo_regions!(ψ★)

# Zonal wind
step(x, d, c) = 1/2 * (1 + tanh((x - c) / d))
polar_mask(y) = step(y, -5, 60) * step(y, 5, -60)
zonal_ψ(y) = (cosd(4y)^3 + 0.5 * exp(-y^2 / 200)) * polar_mask(y)

ψ̄ = Field((Face, Face, Center), model.grid)
set!(ψ̄, (x, y, z) -> zonal_ψ(y))
fill_halo_regions!(ψ̄)

ψ_total = 40 * ψ★ + Ny * ψ̄

u, v, w = model.velocities
η = model.free_surface.η

if !isnothing(model.coriolis)
    using Oceananigans.Coriolis: fᶠᶠᵃ
    f = KernelFunctionOperation{Face, Face, Center}(fᶠᶠᵃ, model.grid, parameters=model.coriolis)
    g = model.free_surface.gravitational_acceleration
    η .= f * ψ_total / g
end

u .= - ∂y(ψ_total)
v .= + ∂x(ψ_total)

#####
##### Shenanigans for rescaling the velocity field to
#####   1. Have a magnitude (ish) that's a fixed fraction of
#####      the surface gravity wave speed;
#####   2. Zero volume mean on the curvilinear LatitudeLongitudeGrid.
#####

# Time-scale for gravity wave propagation across the smallest grid cell
g = model.free_surface.gravitational_acceleration
gravity_wave_speed = sqrt(g * grid.Lz) # hydrostatic (shallow water) gravity wave speed

minimum_Δx = grid.radius * cosd(maximum(abs, grid.φᵃᶜᵃ)) * deg2rad(maximum(abs, grid.Δλᶜᵃᵃ))
minimum_Δy = grid.radius * deg2rad(minimum(abs, grid.Δφᵃᶜᵃ))
wave_propagation_time_scale = min(minimum_Δx, minimum_Δy) / gravity_wave_speed

@info "Max speeds prior to rescaling:"
@show max_u = maximum(u)
@show max_v = maximum(v)
max_speed_ish = sqrt(max_u^2 + max_v^2)

target_speed = 0.5 * gravity_wave_speed
u .*= target_speed / max_speed_ish
v .*= target_speed / max_speed_ish

# Zero out mean motion
using Oceananigans.AbstractOperations: volume

u_cpu = XFaceField(CPU(), grid)
v_cpu = YFaceField(CPU(), grid)
set!(u_cpu, u)
set!(v_cpu, v)

@show max_u = maximum(u)
@show max_v = maximum(v)

u_dV = u_cpu * volume
u_reduced = AveragedField(u_dV, dims=(1, 2, 3))
compute!(u_reduced)
mean!(u_reduced, u_dV)
integrated_u = u_reduced[1, 1, 1]

v_dV = v_cpu * volume
v_reduced = AveragedField(v_dV, dims=(1, 2, 3))
mean!(v_reduced, v_dV)
integrated_v = v_reduced[1, 1, 1]

# Calculate total volume
u_cpu .= 1
v_cpu .= 1
compute!(u_reduced)
compute!(v_reduced)

u_volume = u_reduced[1, 1, 1]
v_volume = v_reduced[1, 1, 1]

@info "Max speeds prior zeroing out volume mean:"
@show maximum(u)
@show maximum(v)

u .-= integrated_u / u_volume
v .-= integrated_v / v_volume

@info "Initial max speeds:"
@show maximum(u)
@show maximum(v)

#####
##### Simulation setup
#####

ζ = VerticalVorticityField(model)
compute!(ζ)
Δt = 0.2 * minimum_Δx / gravity_wave_speed

@info """
    Maximum vertical vorticity: $(maximum(ζ))
    Inverse maximum vertical vorticity: $(prettytime(1/maximum(ζ)))
    Minimum wave propagation time scale: $(prettytime(wave_propagation_time_scale))
    Time step: $(prettytime(Δt))
"""

mutable struct Progress; interval_start_time::Float64; end

function (p::Progress)(sim)
    wall_time = (time_ns() - p.interval_start_time) * 1e-9

    compute!(ζ)

    @info @sprintf("Time: %s, iteration: %d, max(|ζ|): %.2e s⁻¹, wall time: %s",
                   prettytime(sim.model.clock.time),
                   sim.model.clock.iteration,
                   maximum(abs, ζ),
                   prettytime(wall_time))

    p.interval_start_time = time_ns()

    return nothing
end

simulation = Simulation(model,
                        Δt = Δt,
                        stop_time = 100days,
                        iteration_interval = 1000,
                        progress = Progress(time_ns()))

output_fields = merge(model.velocities, (η=model.free_surface.η, ζ=ζ))

output_prefix = "rotating_freely_decaying_barotropic_turbulence_Nx$(grid.Nx)_Ny$(grid.Ny)"

simulation.output_writers[:fields] = JLD2OutputWriter(model, (ζ = ζ,),
                                                      schedule = TimeInterval(10day),
                                                      filename = output_prefix,
                                                      overwrite_existing = true)

# Let's goo!
run!(simulation)


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Variogram{T,D}

A theoretical variogram model (e.g. Gaussian variogram)
with parameters of type `T` and distance of type `D`.
"""
abstract type Variogram{T,D} end

"""
    isstationary(γ)

Check if variogram `γ` possesses the 2nd-order stationary property.
"""
isstationary(γ::Variogram) = isstationary(typeof(γ))

"""
    range(γ)

Return the range of the variogram `γ` when defined.
"""
Base.range(γ::Variogram) = γ.range

"""
    sill(γ)

Return the sill of the variogram `γ` when defined.
"""
sill(γ::Variogram) = γ.sill

"""
    nugget(γ)

Return the nugget of the variogram `γ` when defined.
"""
nugget(γ::Variogram) = γ.nugget

"""
    distance(γ)

Return the distance of the variogram `γ`.
"""
distance(γ::Variogram) = γ.distance

"""
    γ(u, v)

Evaluate the variogram at points `u` and `v`.
"""
(γ::Variogram)(u::Point, v::Point) =
  γ(evaluate(distance(γ), coordinates(u), coordinates(v)))

"""
    γ(U, v)

Evaluate the variogram at geometry `U` and point `v`.
"""
function (γ::Variogram)(U::Geometry, v::Point)
  us = _sample(γ, U)
  mean(γ(u, v) for u in us)
end

"""
    γ(u, V)

Evaluate the variogram at point `u` and geometry `V`.
"""
(γ::Variogram)(u::Point, V::Geometry) = γ(V, u)

"""
    γ(U, V)

Evaluate the variogram at geometries `U` and `V`.
"""
function (γ::Variogram)(U::Geometry, V::Geometry)
  us = _sample(γ, U)
  vs = _sample(γ, V)
  mean(γ(u, v) for u in us, v in vs)
end

"""
    result_type(γ, u, v)

Return result type of γ(u, v).
"""
result_type(γ::Variogram, u, v) = typeof(γ(u, v))

#------------------
# IMPLEMENTATIONS
#------------------

"""
    GaussianVariogram(sill=s, range=r, nugget=n, distance=d)

A Gaussian variogram with sill `s`, range `r` and nugget `n`.
Optionally, use a custom distance `d`.
"""
@with_kw struct GaussianVariogram{T,D} <: Variogram{T,D}
  range::T  = 1.0
  sill::T   = 1.0
  nugget::T = 0.0
  distance::D = Euclidean()
end
(γ::GaussianVariogram)(h) = begin
  # add small eps to nugget
  # for numerical stability
  s = γ.sill
  r = γ.range
  n = γ.nugget + convert(typeof(s), 1e-6)
  (s - n) * (1 - exp(-3(h/r)^2)) + (h > 0) * n
end
isstationary(::Type{<:GaussianVariogram}) = true

"""
    ExponentialVariogram(sill=s, range=r, nugget=n, distance=d)

An exponential variogram with sill `s`, range `r` and nugget `n`.
Optionally, use a custom distance `d`.
"""
@with_kw struct ExponentialVariogram{T,D} <: Variogram{T,D}
  range::T  = 1.0
  sill::T   = 1.0
  nugget::T = 0.0
  distance::D = Euclidean()
end
(γ::ExponentialVariogram)(h) = (γ.sill - γ.nugget) * (1 - exp(-3(h/γ.range))) + (h > 0) * γ.nugget
isstationary(::Type{<:ExponentialVariogram}) = true

"""
    MaternVariogram(sill=s, range=r, nugget=n, order=ν, distance=d)

A Matérn variogram with sill `s`, range `r` and nugget `n`. The parameter
ν is the order of the Bessel function. Optionally, use a custom distance `d`.
"""
@with_kw struct MaternVariogram{T,D} <: Variogram{T,D}
  range::T  = 1.0
  sill::T   = 1.0
  nugget::T = 0.0
  order::T  = 1.0
  distance::D = Euclidean()
end
(γ::MaternVariogram)(h) = begin
  s = γ.sill
  r = γ.range
  n = γ.nugget
  ν = γ.order

  # shift lag by machine precision to
  # avoid explosion at the origin
  h2 = h + eps(typeof(h))
  h3 = sqrt(2.0ν)h2/r

  (s - n) * (1 - 2.0^(1 - ν)/gamma(ν) * h3^ν * besselk(ν, h3)) + (h > 0) * n
end
isstationary(::Type{<:MaternVariogram}) = true

"""
    SphericalVariogram(sill=s, range=r, nugget=n, distance=d)

A spherical variogram with sill `s`, range `r` and nugget `n`.
Optionally, use a custom distance `d`.
"""
@with_kw struct SphericalVariogram{T,D} <: Variogram{T,D}
  range::T  = 1.0
  sill::T   = 1.0
  nugget::T = 0.0
  distance::D = Euclidean()
end
(γ::SphericalVariogram)(h) = begin
  s = γ.sill
  r = γ.range
  n = γ.nugget

  (h < r) * (s - n) * (1.5(h/r) - 0.5(h/r)^3) + (h ≥ r) * (s - n) + (h > 0) * n
end
isstationary(::Type{<:SphericalVariogram}) = true

"""
    CubicVariogram(sill=s, range=r, nugget=n, distance=d)

A cubic variogram with sill `s`, range `r` and nugget `n`.
Optionally, use a custom distance `d`.
"""
@with_kw struct CubicVariogram{T,D} <: Variogram{T,D}
  range::T  = 1.0
  sill::T   = 1.0
  nugget::T = 0.0
  distance::D = Euclidean()
end
(γ::CubicVariogram)(h) = begin
  s = γ.sill
  r = γ.range
  n = γ.nugget

  (h < r) * (s - n) * (7*(h/r)^2 - (35/4)*(h/r)^3 + (7/2)*(h/r)^5 - (3/4)*(h/r)^7) +
  (h ≥ r) * (s - n) + (h > 0) * n
end
isstationary(::Type{<:CubicVariogram}) = true

"""
    PentasphericalVariogram(sill=s, range=r, nugget=n, distance=d)

A pentaspherical variogram with sill `s`, range `r` and nugget `n`.
Optionally, use a custom distance `d`.
"""
@with_kw struct PentasphericalVariogram{T,D} <: Variogram{T,D}
  range::T  = 1.0
  sill::T   = 1.0
  nugget::T = 0.0
  distance::D = Euclidean()
end
(γ::PentasphericalVariogram)(h) = begin
  s = γ.sill
  r = γ.range
  n = γ.nugget

  (h < r) * (s - n) * ((15/8)*(h/r) - (5/4)*(h/r)^3 + (3/8)*(h/r)^5) +
  (h ≥ r) * (s - n) + (h > 0) * n
end
isstationary(::Type{<:PentasphericalVariogram}) = true

"""
    PowerVariogram(scaling=s, exponent=a, nugget=n, distance=d)

A power variogram with scaling `s`, exponent `a` and nugget `n`.
Optionally, use a custom distance `d`.
"""
@with_kw struct PowerVariogram{T,D} <: Variogram{T,D}
  scaling::T  = 1.0
  nugget::T   = 0.0
  exponent::T = 1.0
  distance::D = Euclidean()
end
(γ::PowerVariogram)(h) = γ.scaling*h^γ.exponent + (h > 0) * γ.nugget
isstationary(::Type{<:PowerVariogram}) = false

"""
    SineHoleVariogram(sill=s, range=r, nugget=n, distance=d)

A sine hole variogram with sill `s`, range `r` and nugget `n`.
Optionally, use a custom distance `d`.
"""
@with_kw struct SineHoleVariogram{T,D} <: Variogram{T,D}
  range::T  = 1.0
  sill::T   = 1.0
  nugget::T = 0.0
  distance::D = Euclidean()
end
(γ::SineHoleVariogram)(h) = begin
  s = γ.sill
  r = γ.range
  n = γ.nugget

  # shift lag by machine precision to
  # avoid explosion at the origin
  h = h + eps(eltype(h))

  (s - n) * (1 - sin(π*h/r)/(π*h/r)) + (h > 0) * n
end
isstationary(::Type{<:SineHoleVariogram}) = true

"""
    NuggetEffect(n)
    NuggetEffect(nugget=n, distance=d)

A pure nugget effect variogram with nugget `n`.
Optionally use a custom distance `d`.
"""
@with_kw struct NuggetEffect{T,D} <: Variogram{T,D}
  nugget::T = 0.0
  distance::D = Euclidean()
end
NuggetEffect(n) = NuggetEffect(nugget=n)
(γ::NuggetEffect)(h) = (h > 0) * γ.nugget
Base.range(::NuggetEffect{T,D}) where {T,D} = zero(T)
sill(γ::NuggetEffect) = γ.nugget
isstationary(::Type{<:NuggetEffect}) = true


text:
# ==================================================
# Author: Hayk Aleksanyan
# ==================================================
module Sand

output = include("Output.jl") # supports output functionality (e.g. array to csv)

const N_size = 700;    # the radius of the lattice Z^2, the actual size becomes (2*N+1)x(2*N+1)
const dx = [1,0,-1,0]; # for a given (x,y) in Z^2, (x + dx, y + dy) for all (dx,dy) covers the neighborhood of (x,y)
const dy = [0,1,0,-1];

struct L_coord
  # represents a lattice coordinate
  x::Int
  y::Int
end

struct Graph_rep
  # used in representing a Graph, where each neigh array shows the indices of neighboring vertices to the given vertex
  neigh::Array{Int,1}
end

function Init_CheckArray()
  # returns an extended neighborhood of the origin of Z^2 lattice

  v=zeros(Int,(13,2));

  v[1,1]=0;     v[1, 2] = 0;
  v[2, 1] = 0;  v[2, 2] = 1;
  v[3, 1] = 0;  v[3, 2] = 2;
  v[4, 1] = 0;  v[4, 2] = -1;
  v[5, 1] = 0;  v[5, 2] = -2;
  v[6, 1] = 1;  v[6, 2] = 0;
  v[7, 1] = 2;  v[7, 2] = 0;
  v[8, 1] = -1; v[8, 2] = 0;
  v[9, 1] = -2; v[9, 2] = 0;
  v[10, 1] = 1;  v[10, 2] = 1;
  v[11, 1] = -1; v[11, 2] = 1;
  v[12, 1] = -1; v[12, 2] = -1;
  v[13, 1] = 1;  v[13, 2] = -1;

  return v;

end

function FindCoordinate(Z::Array{L_coord,1}, a::Int, b::Int)
   # in the given array Z of coordinates finds the (first) index of the tuple (a,b);
   # if no match, returns -1

  for i=1:length(Z)
    if (Z[i].x==a)&&(Z[i].y==b)
      return i;
    end
  end

  return -1;
end

function moveSand(T_mass)
    # the main function moving a sand of mass T_mass out of the origin of Z^2
    # implements Algorthim 2 of https://arxiv.org/abs/1607.01525 Discrete balayage and boundary sandpile
    # provides a precise compuation of the final distribution of the boundary sandpile

    Z_lat = zeros(2*N_size+1,2*N_size+1);     # models the integer lattice Z^2
    V_sites = falses(2*N_size+1,2*N_size+1);  # all sites which are visited by the sandpile process, are being marked here
    Odometer = zeros(2*N_size+1, 2*N_size+1); # stores the values of the odometer function

    w_Bound = L_coord[];  # the coordinates of boundary sites which need to move
    w_Int = L_coord[];    # the coordinates of interior sites which need to move

    v_check = Init_CheckArray();   #Array(Int,13,2);

    Norm_mass=T_mass^0.5;  # the capacity of the boundary

    Z_lat[N_size+1,N_size+1] = T_mass;  # the initial mass is stored at the origin
    V_sites[N_size+1,N_size+1] = true;

    push!(w_Bound,L_coord(N_size+1,N_size+1));  # the origin needs to topple

    i, j = 1, 0 ;

    # i1, ... j2  -> show the boundaries of the box which is visited by the sandpile process
    i1, i2, j1, j2 = N_size + 1, N_size + 1, N_size + 1, N_size + 1 ;

    iter_ = 0;   # number of total iterations

    while 0 == 0
        # we iterate until the process becomes stable
        # it is proved in our paper https://arxiv.org/abs/1607.01525
        # that the process stabilises after finitely many moves
        top_ = length(w_Bound);

        iter_ = iter_+1;
        println("\nIteration N:: ", iter_);
        println("The boundaries are:: ", (i2-i1+1),"x",(j2-j1+1), "\n");
        println("Moving the boundary...");

        t1 = time_ns();

        while (top_>0)
            x=w_Bound[top_].x;
            y=w_Bound[top_].y;

            # any mass at the (x,y) needs to be moved out
            dist_mass = Z_lat[x,y]/4;    # the mass is split evenly among lattice neighbors

            Z_lat[x,y + 1] = Z_lat[x,y + 1] + dist_mass; V_sites[x,y + 1] = true;
            Z_lat[x,y - 1] = Z_lat[x,y - 1] + dist_mass; V_sites[x,y - 1] = true;
            Z_lat[x - 1,y] = Z_lat[x - 1,y] + dist_mass; V_sites[x - 1,y] = true;
            Z_lat[x + 1,y] = Z_lat[x + 1,y] + dist_mass; V_sites[x + 1,y] = true;

            Odometer[x,y] = Odometer[x,y] + Z_lat[x,y];
            Z_lat[x,y] = 0;

            pop!(w_Bound);

            i1 = min(i1,x-1);
            i2 = max(i2,x+1);
            j1 = min(j1,y-1);
            j2 = max(j2,y+1);

            for i=1:size(v_check,1)
                a = x+v_check[i,1];
                b = y+v_check[i,2];

                if Z_lat[a,b]!=0
                    is_Interior = ((V_sites[a,b + 1]==true)&&(V_sites[a,b - 1]==true)&&(V_sites[a - 1,b]==true)&&(V_sites[a + 1,b]==true));

                    if (is_Interior==false) && (Z_lat[a,b]<=Norm_mass)
                        continue
                    end

                   if is_Interior == true
                       j = FindCoordinate(w_Bound,a,b);

                       if j > 0
                          splice!(w_Bound,j); #removes the element with index j
                      end
                  end

                 # case 2. if (a,b) needs to be walked out,
                 # but it's not in the boundary stack. we check, to see if it needs to be added
                 if is_Interior==false
                    j = FindCoordinate(w_Bound,a,b);
                    if j == -1
                       push!(w_Bound,L_coord(a,b));
                    end
                 end

               end # Z_lat is non zero
             end # end of going over v_check

             top_ = length(w_Bound);

        end # end of while

        t2 = time_ns();
        print( "time elapsed: " , (t2 - t1)/1.0e9, "\n" );

        # =  = = we now move the interior

        println("Moving interiors.....");
        println("..... Step1. Populating the interior....");
        t1 = time_ns();

        # we first populate the interior
        w_Int=L_coord[];

        for i=i1:i2
            for j=j1:j2
                is_Interior = ((V_sites[i,j + 1]==true)&&(V_sites[i,j - 1]==true)&&(V_sites[i - 1,j]==true)&&(V_sites[i + 1,j]==true));
                if is_Interior == true
                   k = FindCoordinate(w_Int,i,j);
                      if k==-1
                          push!(w_Int,L_coord(i,j));
                      end
                end
            end
        end

        t2 = time_ns();
        print( "time elapsed: " , (t2 - t1)/1.0e9, "\n" );
        println("..... Step2. Creating transformations....");

        t1 = time_ns();

        Transforms_ = createTransforms(w_Int); #Transforms_ is a pair of 2 matrices [T_matrix, Delta_matrix] in this order
        # T_matrix is the coeff  matrix of the mass transform for the interior walk
        # Delta_matrix is the coeff matrix of the mass distribution to boundary for the interior walk

        weights_ = zeros(Float64, length(w_Int),1 );
        for i = 1:length(w_Int)
          weights_[i]=Z_lat[w_Int[i].x, w_Int[i].y];
        end

        t2 = time_ns();
        print( "time elapsed: " , (t2 - t1)/1.0e9, "\n" );

        println("..... Step3. Computing distributed weigths with inversion operator....");
        t1 = time_ns();
        id_M = zeros(Float64, length(w_Int), length(w_Int) );
        for ii = 1:length(w_Int)
           id_M[ii,ii] = 1;
        end

        weights_=(Transforms_[2]*(inv( id_M  - Transforms_[1]   )) )*weights_;

        t2 = time_ns();
        print( "time elapsed: " , (t2 - t1)/1.0e9, "\n" );

        println("..... Step4. Distributing the weigths....");

        t1 = time_ns();

        for i = 1:length(w_Int)
            a = w_Int[i].x;
            b = w_Int[i].y;

            Z_lat[a+1,b] = Z_lat[a+1,b] + weights_[i];
            Z_lat[a-1,b] = Z_lat[a-1,b] + weights_[i];
            Z_lat[a,b+1] = Z_lat[a,b+1] + weights_[i];
            Z_lat[a,b-1] = Z_lat[a,b-1] + weights_[i];

            Odometer[a,b] = Odometer[a,b]+4*weights_[i];
         end

         for i=1:length(w_Int)
             a = w_Int[i].x;
             b = w_Int[i].y;
             Z_lat[a,b]=0;
         end

        t2 = time_ns();
        print( "time elapsed: " , (t2 - t1)/1.0e9,"\n" );

        # it is left to re-populate the boundary

        println("..... Step5. Populating the new boundary....");

        t1 = time_ns();

        for i=1:length(w_Int)
           #checking neighbors of interior points
            for k = 1:4
               a = w_Int[i].x + dx[k];
               b = w_Int[i].y + dy[k];

              is_Interior = ((V_sites[a,b + 1]==true)&&(V_sites[a,b - 1]==true)&&(V_sites[a - 1,b] == true)&&(V_sites[a + 1,b]==true));

              if (is_Interior == false)&&(Z_lat[a,b] > Norm_mass)
                  j=FindCoordinate(w_Bound,a,b);
                  if j == -1
                      push!(w_Bound,L_coord(a,b));
                  end
              end
            end
        end # for loop over interior

        t2 = time_ns();
        print( "time elapsed: " , (t2 - t1)/1.0e9, "\n" );

        w_Int = L_coord[];

        if length(w_Bound) == 0
           break;  # WE're done!!!
        end

    end #end of the main while

    Z_lat = output.printIntoFile(Z_lat,  0, string("BSand_Z_", T_mass) );
    Odometer = output.printIntoFile(Odometer, 1,string("BSand_OD_", T_mass)  );

    output.saveAsGrayImage(Z_lat, string("BSand_Z_", T_mass), 20, 2);
    output.saveAsGrayImage(Odometer, string("BSand_OD_", T_mass), 20, 2);

    # for the total elapsed time, it's better to use the @time macros on the main call

    return Z_lat, Odometer; # these are trimmed in output module

end # end of function moveSand

function createGraph(coord_List::Array{L_coord,1})
  # given a list of lattice sites, creates a graph based on the lattice neighborhood
  # e.g. (0,0) will be connected to all of (0,1), (1,0), (-1,0) , (0, -1)
  # the graph representation is of the form i -> array_of_integers,  where i is the index of the site
  # and the array_of_integers is the list of indices of the neighbors of i

  N = length(coord_List);

  res_= Graph_rep[]; # the result

  for i=1:N
    push!(res_, Graph_rep(Int[]));
  end

  for i=1:N
    for j=1:N

      if (i!=j) && (abs(coord_List[i].x-coord_List[j].x) + abs(coord_List[i].y-coord_List[j].y)==1 )
        push!(res_[i].neigh,j);
      end

    end
  end

  return res_;
end

function createTransforms(coord_List::Array{L_coord,1})
  # T_matrix=zeros(Float64,1,1);  # coeff of the mass transform for the interior walk
  # Delta_matrix=zeros(Float64,1,1); # coeff of the mass distribution to boundary for the interior walk
  # definitions of these matrices are explained in Section 3 of our paper, see Algorithm 2

  N=length(coord_List);
  Delta_matrix=zeros(Float64,N,N);

  if N==1
    Delta_matrix[1,1]=0.25;
  end

  G = createGraph(coord_List);

  T_matrix = zeros(Float64, N, N);
  for i = 1:N
    T_matrix[i,i] = 1; # make T_matrix as an identity matrix
  end

  #println("........... Step2b. Creating the matrices....");

  for i=1:N

    for p=1:N
      Delta_matrix[i,p]=Delta_matrix[i,p]+T_matrix[i,p]/4;
    end

    for p=1:N
      for j=1:length(G[i].neigh)
        k=G[i].neigh[j];
        T_matrix[k,p]=T_matrix[k,p]+T_matrix[i,p]/4;
      end
      T_matrix[i,p]=0;

    end

  end

  return T_matrix, Delta_matrix;
end

function moveSand_Standard(T_mass, max_iter_Count = 50000)
    # performes divisible sandpile process in a standard way
    # each time checking for lattice sites which need to move and adding those in the moving stack

    Z_lat = zeros(2*N_size+1,2*N_size+1);   #Array(Float64, 2*N_size+1,2*N_size+1); models the Z^2 lattice
    V_sites = falses(2*N_size+1,2*N_size+1);  # all sites which are visited by the sandpile process, are being marked here
    Odometer = zeros(2*N_size+1, 2*N_size+1);

    w_Int = L_coord[];  # the stack of lattice sites which need to move (topple)

    Norm_mass = T_mass^0.5;  # the maximal allowed mass on the boundary
    Z_lat[N_size+1,N_size+1] = T_mass; # the initial mass is set at the origin
    V_sites[N_size + 1 , N_size + 1 ] = true;
    push!(w_Int,L_coord(N_size+1,N_size+1));

    i, j = 1, 0;
    #typed_Zeros = typeof(Z_lat[1,1])(0);

    i1, i2, j1, j2 = N_size+1, N_size+1, N_size+1, N_size+1;

    iter_ = 0;
    ex_mass_max = 0;
    ex_mass_max_new = 0;
    ex_mass = 0;

    while iter_ < max_iter_Count

        iter_ = iter_ + 1;
        if (iter_%1000 == 0)
            println("Iteration N:: ", iter_);
            println("The boundaries are:: ", (i2-i1+1),"x",(j2-j1+1), "\n");
        end

        ex_mass_max = 0;

        for  i = 1:length(w_Int)
           x, y = w_Int[i].x, w_Int[i].y ;
           ex_mass_max = max(Z_lat[x,y] , ex_mass_max);
        end

      for tt=1:1500
          for i=1:length(w_Int)
               # the toppling
               x, y = w_Int[i].x, w_Int[i].y;
               if Z_lat[x,y] == 0.0
                  continue
               end

               dist_mass = Z_lat[x,y]/4;

               Z_lat[x,y + 1] += dist_mass;
               Z_lat[x,y - 1] += dist_mass;
               Z_lat[x - 1,y] += dist_mass;
               Z_lat[x + 1,y] += dist_mass;

               Odometer[x,y] += Z_lat[x,y];

               Z_lat[x,y] = 0 ;
          end

          ex_mass_max_new = 0 ;
          for i = 1:length(w_Int)
              x, y = w_Int[i].x, w_Int[i].y;
              ex_mass_max_new = max(Z_lat[x,y]  - 1, ex_mass_max_new);
          end
          if ex_mass_max_new < 0.5*ex_mass_max
             break;
          end
       end # interior cycles

        for i = 1:length(w_Int)
            x , y = w_Int[i].x, w_Int[i].y;

            i1, i2, j1, j2 = min(i1,x-1), max(i2,x+1), min(j1,y-1), max(j2,y+1);

            V_sites[x, y + 1] = true;
            V_sites[x, y - 1] = true;
            V_sites[x - 1, y] = true;
            V_sites[x + 1, y] = true;
        end

        w_Int = L_coord[];
        q = false;

        for i=i1:i2
            for j=j1:j2
                is_Interior = ((V_sites[i,j + 1]==true)&&(V_sites[i,j- 1]==true)&&(V_sites[i - 1,j] == true)&&(V_sites[i + 1,j]==true));
                if (is_Interior)
                    push!(w_Int, L_coord(i,j));
                    if Z_lat[i,j] > 0
                      q = true;
                    end
                else
                    if Z_lat[i,j] > Norm_mass
                        push!(w_Int, L_coord(i,j));
                        q = true;
                    end
                end
            end
        end

        if q == false
          break;  # We are done
        end

    end #end of the main while

  println("saving into CSV...")
  Z_lat = output.printIntoFile(Z_lat,  0, string("BSand_Z_", T_mass), true );
  Odometer = output.printIntoFile(Odometer, 1,string("BSand_OD_", T_mass) , true );

  println("saving as Gray images...")
  output.saveAsGrayImage(Z_lat, string("BSand_Z_", T_mass), 20, 0, true);
  output.saveAsGrayImage(Odometer, string("BSand_OD_", T_mass), 20, 0, true);

  return Z_lat, Odometer; # these are trimmed in output module

end # end of function moveSand_Standard

end # end of the module


text:
type MonteCarloTestSolution{S,T} <: AbstractMonteCarloSolution
  solution_data::S
  errors::Dict{Symbol,Vector{T}}
  error_means::Dict{Symbol,T}
  error_medians::Dict{Symbol,T}
  elapsedTime::Float64
end

type MonteCarloSolution{T} <: AbstractMonteCarloSolution
  solution_data::Vector{T}
  elapsedTime::Float64
end

function calculate_monte_errors(sim::AbstractMonteCarloSolution)
  solution_data = sim.solution_data
  errors = Dict{Symbol,Vector{eltype(solution_data[1].u[1])}}() #Should add type information
  error_means  = Dict{Symbol,eltype(solution_data[1].u[1])}()
  error_medians= Dict{Symbol,eltype(solution_data[1].u[1])}()
  for k in keys(solution_data[1].errors)
    errors[k] = [sol.errors[k] for sol in solution_data]
    error_means[k] = mean(errors[k])
    error_medians[k]=median(errors[k])
  end
  return MonteCarloTestSolution(solution_data,errors,error_means,error_medians,sim.elapsedTime)
end

Base.length(sim::AbstractMonteCarloSolution) = length(sim.solution_data)
Base.endof( sim::AbstractMonteCarloSolution) = length(sim)
Base.getindex(sim::AbstractMonteCarloSolution,i::Int) = sim.solution_data[i]
Base.getindex(sim::AbstractMonteCarloSolution,i::Int,I::Int...) = sim.solution_data[i][I...]
Base.size(sim::AbstractMonteCarloSolution) = (length(sim),)
Base.start(sim::AbstractMonteCarloSolution) = 1
function Base.next(sim::AbstractMonteCarloSolution,state)
  state += 1
  (sim[state],state)
end
Base.done(sim::AbstractMonteCarloSolution,state) = state >= length(sim)

@recipe function f(sim::AbstractMonteCarloSolution)

  for sol in sim
    @series begin
      legend := false
      sol
    end
  end
end


text:
# In this tutorial, we will learn
#  - How to
#
#

a = 1+1

# Another section

max_stars_repo_name: gridap/TutorialsForDevelopers

max_issues_repo_name: gridap/TutorialsForDevelopers

max_forks_repo_name: gridap/TutorialsForDevelopers


text:
# Julia is a high-level, high-performance, dynamic programming language. While it is a general-purpose language and can be used to write any application, many of its features are well suited for numerical analysis and computational science

# Println basically prints out `Hello, World!` to your console/terminal and puts the cursor to the NEXT line. So if you print something else, if will be
# Done on the next line
println("Hello, World!")

max_stars_repo_name: AaryanKhClasses/Hello-World

max_issues_repo_name: AaryanKhClasses/Hello-World

max_forks_repo_name: AaryanKhClasses/Hello-World


text:
using Test
using Libxc
using DFTK: xc_fallback!, energy_per_particle

@testset "Fallback LDA" begin
    for func_name in (:lda_x, :lda_c_vwn, :lda_c_pw)
        func = Libxc.Functional(func_name)

        # Create reference
        ρ = abs.(randn(100))
        Eref  = similar(ρ)
        Vref  = similar(ρ)
        V2ref = similar(ρ)
        Libxc.evaluate!(func, rho=ρ, zk=Eref, vrho=Vref, v2rho2=V2ref)

        @testset "$func_name" begin
            # Compute in DFTK in elevated precision
            ρbig = Array{BigFloat}(ρ)
            E    = similar(ρbig)
            V    = similar(ρbig)
            V2   = similar(ρbig)
            xc_fallback!(func, Val(:lda), ρbig; zk=E, vrho=V, v2rho2=V2)

            @test E  ≈ Eref  atol=5e-15
            @test V  ≈ Vref  atol=5e-15
            @test V2 ≈ V2ref atol=5e-15

            # TODO Check floating-point type consistency:
            # ε = energy_per_particle(Val(func_name), rand(Float32))
            # @test Float32 == eltype(ε)
        end
    end
end

@testset "Fallback GGA" begin
    gga_fallback = (
        :gga_x_pbe,     :gga_c_pbe,     :gga_x_pbe_r,
        :gga_x_pbe_sol, :gga_c_pbe_sol,
        :gga_x_pbefe,   :gga_c_pbefe,
        :gga_x_xpbe,    :gga_c_xpbe,
        :gga_x_pbe_mol, :gga_c_pbe_mol,
        :gga_x_apbe,    :gga_c_apbe,
    )

    for func_name in gga_fallback
        func = Libxc.Functional(func_name)

        # Create reference
        ρ = abs.(randn(100))
        σ = abs.(randn(100))
        Eref  = similar(ρ)
        Vρref = similar(ρ)
        Vσref = similar(ρ)
        Libxc.evaluate!(func, rho=ρ, sigma=σ, zk=Eref, vrho=Vρref, vsigma=Vσref)

        @testset "$func_name" begin
            # Compute in DFTK in elevated precision
            ρbig = Array{BigFloat}(ρ)
            σbig = Array{BigFloat}(σ)
            E    = similar(ρbig)
            Vρ   = similar(ρbig)
            Vσ   = similar(ρbig)
            xc_fallback!(func, Val(:gga), ρbig; sigma=σbig, zk=E, vrho=Vρ, vsigma=Vσ)

            @test E  ≈ Eref   atol=5e-15
            @test Vρ ≈ Vρref  atol=5e-15
            @test Vσ ≈ Vσref  atol=5e-15

            # TODO Check floating-point type consistency:
            # ε = energy_per_particle(Val(func_name), rand(Float32))
            # @test Float32 == eltype(ε)
        end
    end
end