R/emgen.R

In mixAR: Mixture Autoregressive Models

## Do not edit this file manually.
## It has been automatically generated from *.org sources.

mixgenMstep <- function(y, tau, model, index, fix = NULL # 2012-11-02 est_shift=TRUE
                        , comp_sigma = FALSE
                        , method = "BBsolve", maxit = 100, trace = FALSE, lessverbose = TRUE,
                        ... ){
    verbose <- lessverbose && interactive()

    n <- length(y)

    g     <- .nmix(model)
    order <- model@order
    sigma <- model@scale
    p     <- order
    pm    <- max(order)
                                                         # 2012-11-02 tova beshe vav fit_mixAR
                                                  # ako ostavya taka, ste raboti bez izmenenie
    shift_ar_params <- unlist(c(model@shift, model@arcoef@a))

    ## 2021-03-31 new
    ##    TODO: this could be merged better with the code for the other cases
    if(identical(fix, "white.noise")) {
        ## TODO: the calling function should check that the last component has order pmax !!!
        white_noise <- TRUE
        fix <- NULL     # for now; could allow fixed params later!

        constr_params <- numeric(g + sum(model@order))
        constr_indx   <- c(g, g + sum(model@order[1:(g-1)]) + 1:model@order[g])
        unconstr_indx <- seq_along(constr_params)[-constr_indx]
    }else
        white_noise <- FALSE

    dont_g <- .nmix(model)
    dont_ar <- sum(model@order)
    dontfix_shift  <- if(identical(fix,"shift")) rep(FALSE, dont_g) else rep(TRUE, dont_g)
    dontfix_scale  <- rep(TRUE, dont_g)
    dontfix_arcoef <- rep(TRUE, dont_ar)
    if(is.list(fix)){ # should contain named logical vectors with TRUE for fixed elements
               # dontfix_prob   <- if(!is.null(fix$prob))   !fix$prob   else rep(TRUE, dont_g)
        if(!is.null(fix$shift))  dontfix_shift  <- !fix$shift
        if(!is.null(fix$scale))  dontfix_scale  <- !fix$scale
        if(!is.null(fix$arcoef)) dontfix_arcoef <- !fix$arcoef
    }
    dontfix_shift_ar <- c(dontfix_shift, dontfix_arcoef)

    est_shift  <- all(dontfix_shift_ar)                        # compatibility
    est_ar_all <- all(dontfix_arcoef) && all(!dontfix_shift)

    armaxit <- maxit # todo: make an argument for this?

    dist   <- get_edist(model)               # dist <- model@dist
    Fscore <- lapply( dist, function(x) x$Fscore )     # Fscore <- noise_dist(model, "Fscore")
    fpdf   <- lapply( dist, function(x) x$pdf    )     # fpdf   <- noise_dist(model, "pdf"   )

                       # 2012-10-26 TODO: slagam vremenno na TRUE, za da testvam novite raboti
                       #            todo: otchitay i novite raboti!
           # 2012-10-31 estdist_flag <- TRUE  # any( sapply(dist, function(x) x$any_param()) )
           # 2021-03-31 was: estdist_flag <- with(environment(dist[[1]]$pdf), any(param_flags))
           #     it seems that I have changed the way this is retrieved at the time      
#browser()
    ## 20211-06-26 TODO: this seems to assume that any_param() for the 1st component gives an
    ##                   overall result. Shouldn't it be any(sapply(dist, function(x) x$any_param)) ?
    ##
    ## For now leaving as is but this will not work correctly if the first component has no
    ## params to estimate but some of the others do.
    estdist_flag <- dist[[1]]$any_param()

    prob <- tau2probhat(tau)  # TODO: new mixing weights, but model@prob is updated after the other parameters,
                              #    i.e., prob is not used in their updates but rather model@prob.
                              #       TODO: why? (seems ok!)

    lcond <-

        if(white_noise){  # 2021-03-31 new - modified from the 'else' part below
            function(par){
                    # shift_ar_params[dontfix_shift_ar] <- par
                    # model@shift <- shift_ar_params[1:g]
                    # model@arcoef <- rag_modify(model@arcoef, shift_ar_params[-(1:g)])

                cur_probs <- model@prob  # TODO: Dali tova ne tryava da e 'prob'?  
                                         #       may ne, 'prob' e veroyatno initial value
                shift_ar_params[unconstr_indx] <- par
                model@shift <- c(shift_ar_params[1:(g-1)],
                                 - sum(cur_probs[-g] * shift_ar_params[1:(g-1)]) / cur_probs[g])

                    # tmp_m <- matrix(shift_ar_params[-(1:g)], nrow = g - 1, ncol = pm)
                tmp_m <- matrix(0, nrow = g - 1, ncol = pm)
                tmp_v <- shift_ar_params[-(1:g)] # only ar param
                for(i in 1:(g-1)){
                    r_ind <- 1:model@order[i]
		    tmp_m[i, r_ind] <- tmp_v[r_ind]
                    tmp_v <- tmp_v[-r_ind]
		}

                constr_ar_param <- numeric(pm)
                for(k in 1:(g-1)) {
                    for(i in 1:pm){
                        constr_ar_param[i] <-
                            constr_ar_param[i] + cur_probs[k] * tmp_m[k, i]
                    }
                }
                shift_ar_params[constr_indx[-1]] <- - constr_ar_param / cur_probs[g]
                
                model@arcoef <- rag_modify(model@arcoef, shift_ar_params[-(1:g)])

                stdetk <- mix_ek(model, y, index, scale=TRUE)    # standardised resid.
                sc <-  tau * (log %of% ((fpdf %of% stdetk)/ sigma))
                                        # tuk nyama nuzhda da vzemam podmnozhestvo ponezhe se
                                        # vrasta samo edna stoynost!
                sum(sc@m) + sum( (tau * log(model@prob))@m )
            }

        } else if(est_shift){
            function(par){
                model@shift <- par[1:g]
                model@arcoef <- rag_modify(model@arcoef, par[-(1:g)])

                stdetk <- mix_ek(model, y, index, scale=TRUE)    # standardised resid.
                sc <-  tau * (log %of% ((fpdf %of% stdetk)/ sigma))

                res <- sum(sc@m) + sum( (tau * log(model@prob))@m )
                res
            }

        }else if(est_ar_all){
            function(par){
                model@arcoef <- rag_modify(model@arcoef, par)

                stdetk <- mix_ek(model, y, index, scale=TRUE)    # standardised resid.

                log_etk <- log %of% ((fpdf %of% stdetk))
                            if(any(log_etk@m==-Inf)){
                                bad <- log_etk@m== -Inf
                                log_etk@m[bad] <- - 750  # todo: needs better thought!
                                                         # Arbitrary constant!
                            }                            # (but see R transcript above)

                sc <- tau * log_etk

                                        # if(any(is.infinite(sc))){
                                        #     bad <- sc==-Inf
                                        #     sc[bad] <- ifelse(tau[bad]<1e-10, 0, - 750 )
                                        # }

                res <- sum(sc@m) + sum( (tau * log(model@prob/sigma))@m )

                            if(any(is.infinite(res))){
                                ## 2020-06-12 was: cat("infinite values in res!\n")
                                ##                 browser()
                                if(verbose)
                                    message("infinite values in res!")
                            }

                res[res==-Inf] <- - 10000   # see the R snipped above

                res <- sum(res)

                            if(res == -Inf){
                                ## 2020-06-12 was: cat("res is: ", res, "\n")
                                ##                  browser()
                                if(verbose)
                                    message("res is: ", res, "\n")
                            }
                res
            }

        }else{
            function(par){
                shift_ar_params[dontfix_shift_ar] <- par
                model@shift <- shift_ar_params[1:g]
                model@arcoef <- rag_modify(model@arcoef, shift_ar_params[-(1:g)])

                stdetk <- mix_ek(model, y, index, scale=TRUE)    # standardised resid.
                sc <-  tau * (log %of% ((fpdf %of% stdetk)/ sigma))
                                        # tuk nyama nuzhda da vzemam podmnozhestvo ponezhe se
                                        # vrasta samo edna stoynost!
                sum(sc@m) + sum( (tau * log(model@prob))@m )
            }
        }

    eqns <- matrix(NA, nrow = g, ncol = 1 + pm)

    feqns <-

        if(est_shift){ # estimate all shift and ar parameters (ie all of them are not fixed)
            function(par){
                model@shift <- par[1:g]
                model@arcoef <- rag_modify(model@arcoef, par[-(1:g)])

                stdetk <- mix_ek(model, y, index, scale=TRUE)    # standardised resid.

                sc <-  (tau * (Fscore %of% stdetk)) / sigma
                                                   # TODO: tova tyrabva da se vidi vnimatelno!
                                  # 2012-10-31 smenyam znaka; todo: no dali e pravilno?
                eqns[,1] <- - colSums(sc@m)        # separate formula for intercepts
                for(i in seq_len(pm)){  # 1:pm
                    eqns[,i+1] <- - inner( y[(pm+1 -i):(n-i)], sc)
                }
                res <- eqns[,1]
                for(i in 1:g){
                    res <- c(res, eqns[i,1+seq_len(p[i])])
                }
                                        # cat("feqns\n")

                res #  -res   ne mozhe tuk da smenyam znaka, ponezhe smenya i na interceptite!
                    #         2012-10-31 todo: No zasto da ne go smenya za tyach? !!! check!!!
            }

        }else if(est_ar_all){ # all shift parameters fixed and all ar parameters not fixed.
            function(par){
                model@arcoef <- rag_modify(model@arcoef, par)
                stdetk <- mix_ek(model, y, index, scale=TRUE)    # standardised resid.

                                        # cat("feqns 2\n")

                        # todo: need methods for unitary "-"
                        #       currently the sign is changed just before returning the value.
                                           # sc <- - (tau * (Fscore[[1]] %of% stdetk)) / sigma
                sc <-  (tau * (Fscore %of% stdetk)) / sigma
                                        # cat("feqns 3\n")

                for(i in seq_len(pm)){   # 1:pm
                    eqns[,i] <- - inner( y[(pm+1 -i):(n-i)], sc)
                }
                                        # cat("feqns 4\n")
                res <- numeric(0)
                for(i in 1:g){
                    res <- c(res, eqns[i,seq_len(p[i])])   #   [i,1:p[i]]
                }
                                        # cat("feqns\n")
                res
            }

        }else{ # estimate subset of the shift-ar parameters
            function(par){
                shift_ar_params[dontfix_shift_ar] <- par
                model@shift <- shift_ar_params[1:g]
                model@arcoef <- rag_modify(model@arcoef, shift_ar_params[-(1:g)])

                stdetk <- mix_ek(model, y, index, scale=TRUE)    # standardised resid.
                sc <-  (tau * (Fscore %of% stdetk)) / sigma
                                                   # TODO: tova tyrabva da se vidi vnimatelno!
                                  # 2012-10-31 smenyam znaka; todo: no dali e pravilno?
                eqns[,1] <- - colSums(sc@m)        # separate formula for intercepts
                for(i in seq_len(pm)){  # 1:pm
                    eqns[,i+1] <- - inner( y[(pm+1 -i):(n-i)], sc)
                }
                res <- eqns[,1]
                for(i in 1:g){
                    res <- c(res, eqns[i,1+seq_len(p[i])])
                }

                res[dontfix_shift_ar] # return the required subset of the equations.
            }
        }

                                     #     2012-11-02 smen yam s dolnoto zaradi promenite
                                     # if(est_shift)
                                     #     param <- c( model@shift, ragged2vec(model@arcoef) )
                                     # else
                                     #     param <- ragged2vec(model@arcoef)
    if(white_noise){
        tmp_v <- ragged2vec(model@arcoef)
        param <- c( model@shift[-g], tmp_v[1:(length(tmp_v) - model@order[g])] )
    }else if(est_shift)
        param <- c( model@shift, ragged2vec(model@arcoef) )
    else if(est_ar_all)
        param <- ragged2vec(model@arcoef)
    else{ # subset of shift-ar parameters
        param <- c( model@shift, ragged2vec(model@arcoef) )
        param <- param[dontfix_shift_ar]
    }

                                           # print(param) # cat(feqns(param), "\n", sep = " ")
    if(verbose){
        cat("Estimating")
        cat(if(lessverbose) "." else " AR parameters...")
    }
    oldparam <- param
    param <- switch(method,
                        # BBsolve   = {print("Ouch!"); BBsolve(par=param, fn = feqns, ...) },
                    BBsolve   = {BBsolve(par=param, fn = feqns, quiet = TRUE,
                                         control = list(maxit=armaxit, trace = trace), ...) },
                                               # requires monotonicity (set M=1 in "control"?)
                    BBoptim   = BBoptim(par=param, fn = lcond,
                                        control=list(maximize=TRUE, trace = trace), ...),
                    BBoptimgr = BBoptim(par=param, fn = lcond, gr = feqns,
                                        control=list(maximize=TRUE, trace = trace), ...),
                    BBspg     = spg(par=param, fn = lcond, gr = feqns,
                                    control=list(maximize=TRUE, trace = trace), ...),
                    stop("no method specified")
                    )
    if(verbose)
        cat(if(lessverbose) param$convergence else c("convergence =", param$convergence, "fn.reduction =", param$fn.reduction, "\n"))
    if(param$convergence > 0  && !white_noise){             # TODO: tova e krapka!
        oldlik <- cond_loglik(model, y)
        tmpmodel <- model
        if(est_shift){
            tmpmodel@shift <- param$par[1:g]
            tmpmodel@arcoef <- rag_modify(model@arcoef, param$par[-(1:g)])
        }else if(est_ar_all){
            tmpmodel@arcoef <- rag_modify(model@arcoef, param$par)
        }else{
            shift_ar_params[dontfix_shift_ar] <- param$par

            tmpmodel@shift <- shift_ar_params[1:g]
            tmpmodel@arcoef <- rag_modify(tmpmodel@arcoef, shift_ar_params[-(1:g)])
        }

        tmplik <- cond_loglik(tmpmodel, y)
                                        # todo: more work needed, maybe enforce nondecrease?
        if(tmplik < oldlik){  # not only lack of convergence but also worse  likelihood,
                              # try again enforcing monotonicity (M=1)
            if(verbose)
                cat("tmplik:", tmplik, "is smaller than the old lik,", oldlik, "\n")

            tmpparam <- BBsolve(par=oldparam, fn = feqns, quiet = TRUE,
                                control = list(maxit=max(armaxit,100),
                                trace = trace, M = 1), ...)
            if(verbose)
                cat("\n\tconvergence =", tmpparam$convergence, "fn.reduction =", 
                    tmpparam$fn.reduction, "\n")
            param <- tmpparam
        }
    }
                             # if(est_shift){
                             #     model@shift <- param$par[1:g]
                             #     model@arcoef <- rag_modify(model@arcoef, param$par[-(1:g)])
                             # }else{
                             #     model@arcoef <- rag_modify(model@arcoef, param$par)
                             # }
    if(white_noise){  # 2021-03-31 new
        shift_ar_params[unconstr_indx] <- param$par

        model@shift <- shift_ar_params[1:g]
        model@arcoef <- rag_modify(model@arcoef, shift_ar_params[-(1:g)])
    }else if(est_shift){
        model@shift <- param$par[1:g]
        model@arcoef <- rag_modify(model@arcoef, param$par[-(1:g)])
    }else if(est_ar_all){
        model@arcoef <- rag_modify(model@arcoef, param$par)
    }else{
        shift_ar_params[dontfix_shift_ar] <- param$par

        model@shift <- shift_ar_params[1:g]
        model@arcoef <- rag_modify(model@arcoef, shift_ar_params[-(1:g)])
    }

    curlik <- cond_loglik(model, y)

    tmpmodel2 <- model
    if(white_noise){
        ## update mixing weights
        ui <- rbind(rep(-1, g-1), diag(g-1))
        ci <- c(-1, numeric(g-1))

        
        f_wn_pi <- function(par) {
            tmpmodel2@prob[1:(g-1)] <- par
            tmpmodel2@prob[g] <- 1 - sum(par)

            tmpmodel2@shift[g] <- - sum(par * tmpmodel2@shift[1:(g-1)]) / (1 - sum(par))

            constr_ar_param <- - par %*% tmpmodel2@arcoef[1:(g-1), ] / (1 - sum(par))
            tmpmodel2@arcoef[g, ] <- as.numeric(constr_ar_param)
            
            - cond_loglik(tmpmodel2, y)    # 'minus' since constrOptim minimises
        }

        new_probs <- constrOptim(tmpmodel2@prob[1:(g-1)], f_wn_pi, grad = NULL, ui = ui, ci = ci)
        
        tmpmodel2@prob[1:(g-1)] <- new_probs$par
        tmpmodel2@prob[g] <- 1 - sum(new_probs$par)

        tmpmodel2@shift[g] <- - sum(new_probs$par * tmpmodel2@shift[1:(g-1)]) / (1 - sum(new_probs$par))

        tmpmodel2@arcoef[g, ] <- 
            as.numeric( - new_probs$par %*% tmpmodel2@arcoef[1:(g-1), ] / (1 - sum(new_probs$par)) )
    }else{
        tmpmodel2@prob <- prob
    }

    problik <- cond_loglik(tmpmodel2, y)
    if(!white_noise && problik < curlik){
        if(verbose) cat("\tBoshwarning: using the estimated prob reduces the likelihood.\n")
    }

        # 2021-04-04 was: model@prob <- prob  # update the mixture probabilitties
    model <- tmpmodel2

    if(verbose) cat(if(lessverbose) "." else "   scale parameters...")
    etk <- mix_ek(model, y, index)

    if(any(dontfix_scale)){  # update the scale parameters
        sigmahat <- em_est_sigma(tau, etk, Fscore, model@scale,
                                 dontfix = dontfix_scale,
                                 compwise = comp_sigma)
                                                               # was: tauetk2sigmahat(tau,etk)
        model@scale <- sigmahat
    }

        

    if(estdist_flag){
        if(verbose)
            cat(if(lessverbose) "." else  "   noise parameter(s)...")

        nu <- noise_params(model)
                                        # print("nu is: ")
                                        # print(nu)
        parscore  <- lapply( dist, function(x) x$Parscore )
        wrklogpdf <- lapply( dist, function(x) x$logpdf )

                                        # cat("Entering em_est_dist\n")
        nu <- em_est_dist(tau, etk, parscore, sigma, nu, wrklogpdf)

        model <- set_noise_params(model, nu)
    }
    if(verbose)
        cat("\n")

    model
}

                                                              # 2011-11-03 name was em_sigma
                                                              # 2012-10-23 prerabotvam osnovno
em_est_sigma <- function(tau, etk, Fscore, sigma, dontfix = rep(TRUE, length(sigma)),
                         compwise = FALSE){
                              # cat("input:       ", sigma, "\n")
                              # cat("colSums tau: ", colSums(tau@m), "\n")
    tauetk <- tau * etk   # no need to compute it each time inside f
    mtauetk <- tauetk@m
    metk <- etk@m
    tausums <- colSums(tau@m)
                            # todo: check if Fscore is guaranteed to be a list of g components
    getFscore <- function(k){                  # min() here is in case Fscore is of length one
                     if(is.function(Fscore)) Fscore else Fscore[[min(length(Fscore),k)]]
                 }
                          # todo: transformations for 'par' to allow, e.g., common parameters.
                          #     no za da mozhe da se izpolzva  tova tryabva i obst mechanizam.
    if(all(dontfix)){
        f <- function(par){     # multi-variate
            wrk <-  tauetk * (Fscore %of% (etk/par))
            wrk <- - colSums(wrk@m) / tausums
            wrk - par    # equation is wrk - par = 0
        }

        fk <- function(par){    # component-wise                          # par here is scalar
            wrk <-  mtauetk[,k_cur] * sapply(metk[,k_cur]/par, getFscore(k_cur))
            wrk <- - sum(wrk) / tausums[k_cur]
            wrk - par  #  equation is wrk - par = 0
        }
    }else{
        si <- sigma
        f <- function(par){     # multi-variate
            si[dontfix] <- par
            wrk <-  tauetk * (Fscore %of% (etk/si))
            wrk <- - colSums(wrk@m) / tausums
            wrk[dontfix] - par    # equation is wrk - par = 0
        }
    }

    ## (2021-04-03) TODO: temporary patch !!! Find more permanent solution!
    sum_f_sq <- function(par){
        sum(f(par)^2)
    }

                                         # the equations for sigma can be solved independently
                                         # for each component.  need to change `f()' above
                                         # though...  a lazy solution would be (wrk-par)[k]
    if(compwise){
        newsigma <- sigma
        for(k_cur in seq_along(sigma)){
              sigmak_cur <- BBsolve(par=sigma[k_cur], fn = fk, quiet = TRUE,
                                    control=list(trace = FALSE))
            newsigma[k_cur] <- sigmak_cur$par
        }
                  # print(cbind( oldsigma=sigma, newsigma=newsigma, diffsigma=newsigma-sigma))
        sigma <- newsigma
    }else if(all(dontfix)){
        ## (2021-04-03) TODO: temporary patch !!! Find more permanent solution!
        cat("sigma = ", sigma, "\n")
        if(any(sigma < 0.01))
            print("inside!\n")
        sigma[sigma < 0.01] <- 0.01

        ## (2021-04-03) TODO: temporary patch !!! Find more permanent solution!
            # newsigma <- BBsolve(par=sigma, fn = f, quiet = TRUE,
            #                             control=list(trace = FALSE))
        newsigma <- try( BBsolve(par=sigma, fn = f, quiet = TRUE,
                                 control=list(trace = FALSE)) )
        if(inherits(newsigma, "try-error") || any( newsigma$par < 0.01 ))
            newsigma <- optim(par=sigma, fn = sum_f_sq, method = "L-BFGS-B", lower = 0.01)

        newsigma <- newsigma$par
                  # print(cbind( oldsigma=sigma, newsigma=newsigma, diffsigma=newsigma-sigma))
        sigma <- newsigma
    }else{
        ## (2021-04-03) TODO: temporary patch !!! Find more permanent solution!
        if(any(sigma[dontfix] < 0.01))
            sigma[dontfix][sigma[dontfix] < 0.01] <- 0.01
        
        newsigma <- BBsolve(par=sigma[dontfix], fn = f, quiet = TRUE,
                            control=list(trace = FALSE))

        newsigma <- newsigma$par
                  # print(cbind( oldsigma=sigma, newsigma=newsigma, diffsigma=newsigma-sigma))
        sigma[dontfix] <- newsigma
    }
                                      # todo:  krapka, uravnenieto se reshava v sigma i -sigma
    if(any(sigma<0)){                 # kogato Fscore e nechetna function.
        sigma <- abs(sigma)

        ## 2010-06-12 was: print("Changed sign of some sigma[k] to make them positive")
    }
               # cat("newsigma: ", sigma, "\n")
                                          # if(any(is.nan(wrk))) wrk[is.nan(wrk)] <- 0.0000001
    sigma
}

em_est_dist <- function(tau, etk, parscore, sigma, nu, logpdf){
                                        # cat("input:       ", nu, "\n")
                                        # cat("colSums tau: ", colSums(tau@m), "\n")
    lp <- logpdf

    xeval <- etk/sigma

    f <- function(nu){
        wrk <-  tau * (parscore %of% xeval)
        wrk <- colSums(wrk@m)
        wrk     # equation is wrk  = 0
    }
                      # origlpnu <- get("nu", environment(lp)) - this was when lp =logpdf[[1]]
    parenv <- parent.env(environment(lp[[1]]))

    if(exists("creator_fun", parent.env(parenv), inherits=FALSE))   # todo: krapka!
        parenv <- parent.env(parenv)
    else if(exists("creator_fun", parent.env(parent.env(parenv)), inherits=FALSE))
        parenv <- parent.env(parent.env(parenv))  # TODO: krapka! 2012-10-30 !!!
                                                  # tova stana s vavezhdaneto na ed_skeleton


    origlpnu <- get("nu", parenv)        # cat("origlpnu is: ", origlpnu, "\n")

    if(exists("set_non_fixed", parenv, inherits=FALSE)){  # todo: krapka 2012-10-26
        assign("nu", nu, parenv)  # todo: krapka, 'nu' currently comes as the
        nu <- evalq(get_non_fixed(), parenv)  #       whole nu, not as non-fixed only.

        f0 <- function(nu){
                    # evalq does not do here since we need nu to be evaluated before the call.
                    # evalq(set_non_fixed(nu), parenv) # assign("nu", nu, parenv)
            do.call("set_non_fixed", list(nu), envir = parenv)
            wrk <-  tau * (lp %of% xeval)

            if(any(is.nan(wrk@m)))         # todo: krapka, no: NaN appear from 0*(-Inf)
                wrk@m[is.nan(wrk@m)] <- 0  #     i.e. elem of tau = 0 and corresp value = -Inf
                                           # todo: tova sigurno ne e edinstvenata vazmozhnost.

            wrk <- colSums(wrk@m)
            sum(wrk)
        }
    }else
        f0 <- function(nu){
            assign("nu", nu, parenv)        # assign("nu", nu, environment(lp))
            wrk <-  tau * (lp %of% xeval)

            if(any(is.nan(wrk@m)))         # todo: krapka, no: NaN appear from 0*(-Inf)
                wrk@m[is.nan(wrk@m)] <- 0  #     i.e. elem of tau = 0 and corresp value = -Inf
                                           # todo: tova sigurno ne e edinstvenata vazmozhnost.
            wrk <- colSums(wrk@m)
            sum(wrk)
        }
                       # the equations for nu can be solved independently for each component.
                       # need to change `f()' above though...  a lazy solution would be wrk[k]
                               # todo: nu must be > 2 !!!
          # newnu <- BBsolve(par=nu, fn = f)

    newnu <- BBoptim(par=nu, fn = f0, lower = 2.2,
                     control = list(maximize=TRUE, trace = FALSE), quiet = TRUE)

    assign("nu", origlpnu, parenv)            # restore the original parameter(s) of lp

    newnu <- newnu$par
                    # print(cbind( oldnu=nu, newnu=newnu, diffnu=newnu-nu))

    if(exists("set_non_fixed", parenv, inherits=FALSE)){  # todo: krapka 2012-10-26
                                                         # vzh i v nachaloto na tazi funktsiya
        do.call("set_non_fixed", list(newnu), envir = parenv)
        nu <- evalq(get_nu(), parenv)
    }else
        nu <- newnu

    nu
}

mixAR_cond_probs <- function(model, y, indx = NULL){
    if(is.null(indx))
        indx <- (max(model@order) + 1) : length(y)
    stdetk <- mix_ek(model, y, indx, scale=TRUE)  ## em_tau needs standardised resid.

                                       # em_tau_safe below takes care of this but is currently
                                       # commented out. Probably need better solution.
        # if(any(!is.finite(stdetk)))               # todo: krapka! appropriate if matching
        #     stdetk@m[!is.finite(stdetk)] <- Inf   #   p_k, scale_k are both (close to) zero.

    fpdf <- noise_dist(model, "pdf")  # note: this may change if noise dist. has parameters
                                      #       so, in general needs to be computed every time

         # em_tau_safe(stdetk, model@prob, model@scale, fpdf)
    em_tau(stdetk, model@prob, model@scale, fpdf) # todo: dali za slozha i argument indx?
}

mixARgenemFixedPoint <- function(y, model               # crit   # opts
                                 # ne, zasega chrez '...', fix = "" # 2012-11-02 new arg.
                                                  # 2012-11-01 comment out, they are for '...'
                                                  # , est_shift = TRUE
                                                  # , comp_sigma = FALSE
                                 , crit = 1e-14   # maybe set to a larger value?
                                 , maxniter = 200
                                 , minniter = 10
                                 , verbose = FALSE
                                 , ... ){
    verbose <- verbose && interactive()

    y <- as.numeric(y)

    oldmodel <- model

    p <- oldmodel@order
    pm <- max(oldmodel@order)
    n <- length(y)
    indx <- (pm+1):n

    nmp <- n - pm

    oldvallogf <- cond_loglik(oldmodel, y)
    newvallogf <- NA
    relchange <- crit + 1

    all_relchange <- numeric(0)  # for testing
    all_vallogf <- oldvallogf

    if(verbose)
        cat("niter: ", 0,  "\t", "vallogf: ", oldvallogf, "\n")

    niter <- 0
    while( niter <= minniter || niter <= maxniter && relchange > crit ){
        niter <- niter + 1
                        # 2012-10-23 zamenyam nyakolko komandi s edin call na mixAR_cond_probs
        newtau <- mixAR_cond_probs(oldmodel, y, indx)

                          # if(any(colSums(newtau@m) < 2)){  # todo: arbitrary threshold here!
                          #     print(newtau)
                          #     stop("One or more components are too improbable")
                          # }
        newmodel <- mixgenMstep(y, newtau, oldmodel, indx, ...)

        oldvallogf <- newvallogf
        newvallogf <- cond_loglik(newmodel, y)
        relchange <- abs((oldvallogf - newvallogf)/oldvallogf)

        all_vallogf   <- c(all_vallogf, newvallogf)
        all_relchange <- c(all_relchange, relchange)

        if(niter %% 10 == 1){  # todo: replace 10 with an argument.
            if(verbose){
                cat("niter: ", niter, "\t", "vallogf: ", newvallogf, "\n")
                show(newmodel)
            }
        }

        oldmodel <- newmodel
    }
    list(model=newmodel, vallogf=newvallogf, niter = niter
         , all_vallogf = all_vallogf
         , all_relchange = all_relchange)
}