module Sim3
include("./lib/probs.jl")
include("./lib/pccf.jl")
include("./detector.jl")
include("./lib/roc.jl")
include("./lib/utils.jl")
#using Probs, Detector, Pccf

function calcTmax(mu, n)
    round(Int,(mu*(n+1)))
end

function calcWidth(mu)
    1 + round(Int, 1.33 * mu)
end

# 1
function changesLocations(n, mu, sigma)
    changes1::Array{Int,1}
    tMax = calcTmax(mu,n)
    width = calcWidth(mu)
    changes1 = Probs.generateChanges(n,mu,sigma,width)
end
# 2
function errorsLocations(n,mu,errRate)
    tMax = calcTmax(mu,n)
    avgNumOfErrors = errRate * n
    #round(Int, rand(1:tMax, avgNumOfErrors))
    rand(1:tMax, avgNumOfErrors)
end


## ;;;;;;;;;;;;;;;;;Pccf/v1.0/src/recurrency_modelling2/sim3.clj;;;;;;;
## (defn gen-vals-in-interval-uniform [n from-v to-v]
##   "Function to generate values within interval inclusive"
##   (let [rng (inc (- to-v from-v))]
##     (map #(+ % from-v) (repeatedly n #(rand-int rng)))))
##
## (defn gen-vals-in-interval-gauss [n from-v to-v]
##   (let [mu (/ (+ to-v from-v) 2.0) sigma (inc (- to-v from-v))]
##     (utils/rnorm n mu sigma)))
## ;; (gen-vals-in-interval 5 7 8)
## ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

function genValsInIntervalGauss(n,st,en)
    # used to simulate changes
    mu = (en+st)/2 #old: (n+st)/2
    sigma = en-st+1
    randn(n) * sigma + mu
end

function updVec(vec, ids::Array{Int,1}, newVals)
    v=copy(vec)
    for i in 1:length(ids)
        v[ids[i]]=newVals[i]
    end
    return v
end

function genSig(n,mu,sigma,errRate,tvf,tvt,evf,evt)
    # Generate signal. For changes values are from tvf-tvt; errors in rvf-evt
    tMax        = calcTmax(mu,n)
    changesLoc  = changesLocations(n,mu,sigma)
    errorsLoc   = errorsLocations(n,mu,errRate)
    nChps       = length(changesLoc)
    nErrors     = length(errorsLoc)
    valsChanges = genValsInIntervalGauss(nChps, tvf, tvt)
    valsErrors  = genValsInIntervalGauss(nErrors, evf, evt)
    vals        = zeros(tMax)
    changes1    = updVec(vals, errorsLoc, valsErrors)
    changes2    = updVec(changes1, changesLoc, valsChanges)
    Dict("changes"=>changesLoc, "sig"=>changes2)
end

function avgPairs(listOfTuples)
    f=[x[1] for x in listOfTuples]
    s=[x[2] for x in listOfTuples]
    (mean(f),mean(s))
end
# julia> avgPairs([(1,2),(5,7)])-> (3.0,4.5)

function simForH(nSim,nChps, mu,sigma,errRate,tvFrom, tvTo, evFrom,evTo,h)
    # Perform n simulations for a given threshold h
    performance=[]
    for i in 1:nSim
        dat=genSig(nChps,mu,sigma,errRate,tvFrom,tvTo,evFrom,evTo)
        newsignal=dat["sig"]
        tChps=dat["changes"]
        maxDelay=1
        detections=Detector.naiveDetector(newsignal,h)
        tprate=Roc.tpRate(detections,tChps,maxDelay)
        farate=Roc.faRate(detections,tChps,maxDelay,length(newsignal))
        push!(performance,(tprate,farate))
    end
    avgPairs(performance)
end

function simDynamicH(probsFile,nSim,nChps,mu,sigma,errRate,tvFrom,tvTo,
                     evFrom,evTo,delayInMu,hmin,hmax)
    # Calculate TP/FA rates for a given pair of hmin hmax
    maxDelay=1
    tMax=calcTmax(mu,nChps)
    hForProbsFunction=0.1
    T=round(Int, mu*nChps*2)
    probsNoConfirm=Pccf.pccf(T,mu,sigma)
    performance=[]
    for i in 1:nSim
        dat = genSig(nChps,mu,sigma,errRate,tvFrom,tvTo,evFrom,evTo)
        newsignal = dat["sig"]
        tChps = dat["changes"]
        subProbsNoConfirm = probsNoConfirm[1:length(newsignal)]
        probs = Pccf.pccfWithConfirmation(subProbsNoConfirm,delayInMu,tChps)
        dynamicHvec = Probs.dynamicH(probs,hForProbsFunction,hmin,hmax, tMax)

        detections = Detector.naiveDetectorDynamic(newsignal,dynamicHvec)

        tprate = Roc.tpRate(detections, tChps, maxDelay)
        farate = Roc.faRate(detections, tChps, maxDelay,length(newsignal))
        push!(performance, (tprate,farate))
    end
    avgPairs(performance)
end

function simNotAdaptive(nSim,nChps,mu,sigma,
                        errRate,tvFrom,tvTo,evFrom,
                        evTo,h1,h2)
    # For every threshold h from h1 to h2 perform n-sim simulations
    performance=[]
    h=h1
    while h<=h2
        avgPerf=simForH(nSim,nChps,mu,sigma,errRate,tvFrom,tvTo,evFrom,evTo, float(h))
        push!(performance, (h, avgPerf[1],avgPerf[2]))
        h=h+1
    end
    performance
end

function runSimNotAdaptive(nSim, mainpath, settings)
    s=settings
    notAdaptivePerf=simNotAdaptive(nSim,s["nchps"],s["mu"],s["sigma"],
                                   s["errRate"],s["tvFrom"],s["tvTo"],s["evFrom"],
                                   s["evTo"],20,50)

    write_items(mainpath*"NotadAptivePerf.txt", notAdaptivePerf, "Delay,TP,FA")
end

function runSimAdaptive(nSim,delay1,delay2,hmin,hmax,mainpath,settings)
    s=settings
    performance=[]
    for confirmDelay in delay1:delay2
        adaptivPerf=simDynamicH(mainpath*"recprobs.txt",nSim,s["nchps"],s["mu"],
                                s["sigma"],s["errRate"],s["tvFrom"],s["tvTo"],
                                s["evFrom"],s["evTo"],confirmDelay,hmin,hmax)
        push!(performance, (confirmDelay,adaptivPerf[1],adaptivPerf[2]))
    end
    outpath=mainpath*"DelayVsPerformance.txt"
    write_items(outpath, performance, "Delay,TP,FA")
    println("... runSimAdaptive(): results are written into $outpath")
end

function run(settings)
    mainpath="D:/_tmp/out/"
    stCommon=Dict("nsim"=>200)
    stAdaptive=Dict("delay1"=>1, "delay2"=>3,"hmin"=>27.0,"hmax"=>30.0)
    runSimNotAdaptive(stCommon["nsim"], mainpath, settings)
    runSimAdaptive(stCommon["nsim"],stAdaptive["delay1"],stAdaptive["delay2"],
                   stAdaptive["hmin"],stAdaptive["hmax"], mainpath, settings)
end

end
