function[mas_newsost,mas_sredn,r_min,r_max,r_avg,s,errors]=prognoz_maxr(y_urez,dt,rmin_usr,rmax_usr,s,nprognoz,tbeg_progn,k,tbeg,tend,tbeg_sost,tend_sost);
% Функция прогнозирования, с использованием цепи максимального порядка.
% Пытаемся найти все повторения предыстории (не менее указанного количества
% повторений) и по их исходах прогнозируется вероятность следующего
% состояния.
%y_urez - ряд для прогноза ,dt - дискретизация ряда, которую нужно
%проделать,r - порядок цепи,s - количество уровней квантования,n -
%,bif_flag - флаг бифуркаций. Если 0 - выбирается меньшая, 1 - большая.
%tbeg - место начала прогноза
%
% errors - количество ситуаций, когда НЕВОЗМОЖНО найти следующий (правка
% 20081124
%Глобальная переменная.
%Блок дискретизации приращений. Из урезанного временного ряда получаем
%массив приращений с шагом дх.
r=rmin_usr;
N1=length(y_urez);
global type_prirash;
if isempty(type_prirash)
    type_prirash=0;
end
if type_prirash==0
    dy=y_urez(dt+1:N1)-y_urez(1:N1-dt);% абсолютные приращения.
else
    dy=(y_urez(dt+1:N1)-y_urez(1:N1-dt))./y_urez(dt+1:N1);% Относительные приращения
end
global type_s;
global s_const;
global statdt1;
statdt=getstdotklon(dy);
if isempty(statdt1)
    statdt1=statdt;
end
otnoshenie_sigma=statdt/statdt1;
lastdsigma=statdt;
if type_s==2
    s=round(s_const*otnoshenie_sigma);
    if (s<3)
        s=3;
    end
    global porog_smax;
    if (s>porog_smax)
        s=porog_smax;
    end
    global logs_skleyka;
    global hlogfile;
    if (isempty(hlogfile))
        hlogfile=1;
    end
    if (~isempty(logs_skleyka))
        fprintf(hlogfile,'dt=%d\tdispersiya=%d\tS=%d\n',dt,statdt,s);
    end
end

 
% Блок преобразования приращений в номера состояний, подсчет вероятностей
% переходов
%raspr=raspred_veroyatn(dy,10);
[ysost,sredn_sost,stdotklon_group,stdotklon_all,matrix_kol]=dy2sost(dy,s,r,dt,tbeg_sost,tend_sost);

% блок прогнозирования. 
mas_sredn=sredn_sost;
tbeg_progn=tbeg_progn-dt; % Состояние, соотв. точке tbeg_progn в процессе нахожд приращений, смещается на дт.
if (tbeg_progn>length(ysost))
    tbeg_progn=length(ysost);
end

%Готовим начальные состояния
kkk=ceil(tbeg_progn/dt)-1;
sost4prognoz=ysost(tbeg_progn-kkk*dt:dt:tbeg_progn);
minlen=floor(length(sost4prognoz))-1;
mas_newsost=[];
%global bif_flag; 
%bif_flag=1;
h=waitbar(0,strcat('Calculating dt=',num2str(dt)));
r_min=99999;
r_max=0;
r_avg=0;
sumr=0;
if (tend>length(ysost))
    tend=length(ysost);
end
ysost_yrez=ysost(tbeg:tend);
errors=0;
global hash_table_findcursost;% хеш-матрица для поиска значений уже найденных
%pack;
%hash_table_findcursost=sparse(4294967294,10000); % очистка хеш-матрицы
hash_table_findcursost=sparse(4294967294/2,10000); % очистка хеш-матрицы /2 по требованию матлаба 2009
global hash_results_findcursost;
hash_results_findcursost=[];
global hash_findnextver;
hash_findnextver=sparse(4294967294/2,42949672/2);
for i=1:nprognoz
    beg=length(sost4prognoz)-minlen;
    lastsost=sost4prognoz(beg:beg+minlen);
    [vectorver,rmax,place]=findnextver(ysost_yrez,rmin_usr, rmax_usr,dt,lastsost,k);
    nextsost=vector2sost(vectorver,dt,i);
    mas_newsost=[mas_newsost,nextsost];
    sost4prognoz=[sost4prognoz,nextsost];
    % Подсчет статистики по р для исследования параметров с и к.
    if (length(place)<k)% не нашлось данного элемента
        errors=errors+1;
    end
    if (rmax>r_max)
        r_max=rmax;
    end
    if (rmax<r_min)
        r_min=rmax;
    end
    sumr=sumr+rmax;
    %fprintf(1,'%d\n',nextsost);
    waitbar(i/nprognoz,h);
    %i
    %nprognoz
end
r_avg=sumr/nprognoz;
close(h);


function newsost=vector2sost(vector_ver,dt,i);
% Функция определения наиболее вероятного следующего состояния по вектору
% вероятностей перехода в другие состояния.
%verthreshold=max(matrix_ver(numsost_prev,1:numsost));
global sigma;
if isempty(sigma)
    sigma=0.01;
end
if sigma >1
    sigma=1;
end
if sigma <0
    sugma=0;
end
%sigma=0.01;
verthreshold=max(vector_ver)-sigma;
newsostsum=sum(vector_ver);
newsostrnd=round(random('unif',1,newsostsum));
summ=0;
newsost2=0;
numsost=length(vector_ver);
    sostmax=[];
    for jj=1:numsost
        if (vector_ver(jj)>=verthreshold)
            sostmax=[sostmax,jj];
      
        end
        summ=summ+vector_ver(jj);
        if ((summ>newsostrnd)&&(newsost2==numsost))
            newsost2=jj;
        end
    end
nmaxsost=length(sostmax);
        
%inogda weroyatnosti ravny 0. W etom sluchae my ne znaem, kuda poyti.
%Predlagaem 
    %global masstat;
    %r=length(masstat)+1;
sumthreshold=2; % количество появления события, которое считаем очень редким и вынуждены рассмотреть цепи меньших порядков.
mid=round(length(sostmax)/2); 
if ((mid>0)&&(mid<=length(sostmax)))
    newsost=sostmax(mid);%naibolee weroyatnoe - posredine
else
    -99999
end
% Ситуация нескольких наиболее вероятных состояний. Выбор одного или двух
% (в случае бифуркации).

if (nmaxsost>1)
    sostmax;
    global bif_filename;
    if isempty(bif_filename)
        bif_filename='biftimes.txt';
    end
    
    % Алгоритм: 1) найти наибольшие кластеры (группы соседних точек с
    % наибольшим количеством)
    maxclusters=getclusters(sostmax,numsost);
    sost_sredn=round(numsost/2);
    blizhnsost=numsost;
    bestcenters=[];
    for kk=1:length(maxclusters)
        if (abs(maxclusters(kk).center-sost_sredn)<blizhnsost)
            bestcenters=maxclusters(kk).center;
            blizhnsost=abs(maxclusters(kk).center-sost_sredn);
        else
            if (maxclusters(kk).center-sost_sredn==blizhnsost)
                bestcenters=[bestcenters,maxclusters(kk).center];
            end
        end
    end
    newsost=bestcenters;
    global bif_flag; %Выбираем одну из бифуркаций, в зависимости от глобальных настроек (значения глобальной переменной bif_flag)
    if (length(newsost)>1)
        if isempty(bif_flag) %не задана переменная для выбора вида бифуркации.
            newsost=newsost(1);
        else
            newsost=newsost(2);
        end
    end
    global bif_filename;
    if isempty(bif_filename)
        bif_filename='biftimes.txt';
    end
    fp=fopen(bif_filename,'a');
    fprintf(fp,'dt=%d,t=%d,probabilities=(',dt,i);
    for ii=1:length(vector_ver)
        fprintf(fp,'%d,',vector_ver(ii));
    end
    fprintf(fp,'), newsost=%d\n',newsost);
    fclose(fp);
    % 2) Из них выбрать тот, который ближе к центру.
    % 3) Если всего 2 кластера и оба на одинаковом расстоянии от центра -
    % делается бифуркация (выбирается 2 варианта)
    % 4) при выборе состояния из кластера выбирается точка середины
    % кластера.
end
% wariant # 2 сгенерировать следующее состояние случайно, но по
% вероятностям, указанных в векторе переходов.
if (newsost2>0)
    %newsost=newsost2;%Принимаем во внимание следуюшее по 2-й методике.
    %Если закоментировать предыдущую строку, будет считаться 1-я методика -
    %наиболее вероятное состояние.
end

function  maxclusters=getclusters(sostmax,nsost);
% находим все кластеры 
n=length(sostmax);
elements=sostmax(1);
ncluster=1;
clusters=[];
maxelems=0;
for i=2:n
    if (sostmax(i-1)+1==sostmax(i))%кластер продолжаетсz
        elements=[elements,sostmax(i)];
        ncluster=ncluster+1;
    else % кластер закончился
        str=[];
        str.n=ncluster;
        str.elements=elements;
        if (mod(ncluster,2)==0) %Если кол-во элементов в кластере четное - неизвестно какой из центральных брать.
           if (elements(ceil(ncluster/2))<round(nsost/2))
               str.center=elements(ceil(ncluster/2)+1);%слева от центра. Берем правую точку.
           else
               str.center=elements(ceil(ncluster/2));%справа от центра. Берем левую точку.  
           end 
        else % Если нечетное - точно есть центральное значение.
           str.center=elements(ceil(ncluster/2));
        end
        clusters=[clusters,str];
        if (ncluster>maxelems)
            maxelems=ncluster;
        end
        elements=sostmax(i);
        ncluster=1;
    end
end
str=[];
str.n=ncluster;
str.elements=elements;
if (mod(ncluster,2)==0) %Если кол-во элементов в кластере четное - неизвестно какой из центральных брать.
    if (elements(ceil(ncluster/2))<round(nsost/2))
        str.center=elements(ceil(ncluster/2)+1);%слева от центра. Берем правую точку.
    else
        str.center=elements(ceil(ncluster/2));%справа от центра. Берем левую точку.  
    end 
else % Если нечетное - точно есть центральное значение.
    str.center=elements(ceil(ncluster/2));
end
clusters=[clusters,str];
if (ncluster>maxelems)
    maxelems=ncluster;
end
% Выбор максимальных кластеров
ncl=length(clusters);
maxclusters=[];
for i=1:ncl
    if (clusters(i).n==maxelems)
        maxclusters=[maxclusters,clusters(i)];
    end
end

%function centers=getcenters(maxclusters);
% возвращает центры кластеров

%% Ошибка! В случае, когда кластер с четным количеством эл-тов попадает
%% ровно на центр распределения. Необходимо брать 2 центральных варианта.
%% вариант 2 - рассмотрен в рукописях номер 13. 

function [vectorver,rmin,place]=findnextver(ysost,rmin_usr, rmax_usr,dt,sost4prognoz,minnum);
% Функция определения следующего состояния.
%Тестирование процесса выбора значения.
%vectorver=[0.3 0.5 0.1 0.6 0.6]; %тест выбора максимума Наче работает.
len=length(sost4prognoz);
s=max(ysost);
maxlength=length(ysost)/dt;
i=len-1;
place=[];
vector_col=0;
sposob=0; %способ приближения к нужному р. 0 - уменьшением в 2 р, 1 - уменьшение на 1 единицу
rmin=rmin_usr;
rmax=rmax_usr;
if rmax >len-1
    rmax=len-1;
end
if (len>0)
    mashistories=[];
    while (rmax-rmin>=0) % Исправление 20081126 с целью заставить р уменьшаться до 1, а не до 2. Убраноменьшеравно. Приводит к беск. циклу
        %for i=len-1:-1:0 %цикл для длины предыстории.
        %i
        i=floor((rmin+rmax)/2);% srednee
        cursost=sost4prognoz(len-i+1:len);
        mashistories{length(mashistories)+1}=cursost;
        hashnum=gethashnumber2(cursost, s);
        %if (hashnum~=0) % !!!! На время отключим хеширование. Ошибка!!!!
        %    rmin=hashnum;
        %    break;
        %end
        [vector_col,place]=findcursost(ysost,dt,cursost);
        %[exactk,newplaces]=getexactk(place,dt);
        if (rmax-rmin)<2
            break;
        end
        if (sum(vector_col)>=minnum)
            rmin=i;
        else 
            rmax=i;
        end
    end
    i=rmin;
    if (i>=len) 
        i=len-1;
    end
    sethashmatrix(mashistories,rmin,s); % Ошибка! Не все истории нужно фиксировать в хеше, а только те, которые больше rmin.
    cursost=sost4prognoz(len-i+1:len);
    [vector_col,place]=findcursost(ysost,dt,cursost);
    if (sum(vector_col)>0)
        vectorver=vector_col./(sum(vector_col));
    else
        vectorver=ones(1,max(ysost));
        place=[];
    end
else
    vectorver=ones(1,max(ysost));
    place=[];
end
rmax=i+1;

function  [vector_col,place]=findcursost(ysost,dt,cursost);
% Находит все повторения данной предыстории в истории и возвращает вектор
% кол-ва выпадения для всех вариантов следующих событий. 
numsost=max(ysost);
history=length(cursost);
vector_col=zeros(1,numsost);
lasti=length(ysost)-(history+1)*dt;
hashnum=gethashnumber(cursost, numsost);
global hash_table_findcursost;
global hash_results_findcursost;
% хеш-таблица состоит из след. размерностей: 
% хеш-номера
% ер
% ответ
% При новом дт или ес должна очищаться.
if (hashnum<2147483647)
    if(hash_table_findcursost(hashnum,history)==0)
        place=[];
        for i=1:lasti;
            currrsost=ysost(i:dt:i+(history+1)*dt);
            if (sum(abs(cursost-currrsost(1:history)))/history<=0.04) %найденно нужное состояние
                nextsost=currrsost(history+1);
                if isempty(place)
                    place=[place,i];
                    vector_col(nextsost)=vector_col(nextsost)+1;
                else
                    if (i-place(length(place))>(dt/2))
                        place=[place,i];
                        vector_col(nextsost)=vector_col(nextsost)+1;
                    end
                end
            end
        end
        answer=[];
        answer.vector_col=vector_col;
        answer.place=place;
        hash_results_findcursost=[hash_results_findcursost,answer];%{vector_col,place};
        hash_table_findcursost(hashnum,history)=length(hash_results_findcursost);
    else
        answer_num=hash_table_findcursost(hashnum,history);
        answer=hash_results_findcursost(answer_num);
        vector_col=answer.vector_col;
        place=answer.place;
    end
else
    place=[];
    for i=1:lasti;
        currrsost=ysost(i:dt:i+(history+1)*dt);
        if (sum(abs(cursost-currrsost(1:history)))/history<=0.04) %найденно нужное состояние
            nextsost=currrsost(history+1);
            if isempty(place)
                place=[place,i];
                vector_col(nextsost)=vector_col(nextsost)+1;
            else
                if (i-place(length(place))>(dt/2))
                    vector_col(nextsost)=vector_col(nextsost)+1;
                    place=[place,i];
                end
            end
        end
    end
end



function num=gethashnumber(cursost, s);
% Возвращает хеш-значение для искомой последовательности состояний
n=length(cursost);
num=1;
for i=1:n
    % Кодируем с конца
    num=num*(s+1)+cursost(i);
end

function [rbest,masnum]=gethashnumber2(cursost, s);
% Находит хеш-значение для искомой последовательности состояний и
% возвращает лучшую длину р для данного к
n=length(cursost);
masnum=[];
num=1;
k=1;
for i=1:n
    % Кодируем с конца
    num=num*(s+1)+cursost(i);
    k=k+1;
    if ((s+1)^(k+1)>=2147483647)
        masnum=[masnum,num];
        num=1;
        k=1;
    end
end
masnum=[masnum,num];
global hash_findnextver;

if (length(masnum)>2)
    rbest=0;
else
    if length(masnum)>1
        if (masnum(1)<size(hash_findnextver,1))&&(masnum(2)<size(hash_findnextver,2))
            rbest=hash_findnextver(masnum(1),masnum(2));
        else
            rbest=0;
            -111;
        end
    else
        rbest=hash_findnextver(masnum(1),1);
    end
end

function sethashmatrix(mashistories,r,s);
% Добавляет в хеш-таблицу лучших р новые значения, соотв. массиву историй.
global hash_findnextver;
n_hist=length(mashistories);
for i=1:n_hist
    % Ошибка! Не все истории нужно фиксировать в хеше, а только те, которые больше rmin.
    if (length(mashistories{i})>=r)
        [rbest,masnum]=gethashnumber2(mashistories{i}, s);
        if (length(masnum)>2)%||(memory.MaxPossibleArrayBytes<100)
        else
            try
                if (length(masnum)>1)
                    hash_findnextver(masnum(1),masnum(2))=r;
                else
                    hash_findnextver(masnum(1),1)=r;
                end
            catch
                fprintf(1,'Low available memory! Number of records=%d!\n');
            end
        end
    end
end


function [exactk,newplaces]=getexactk(places,dt);
% Предполагается, что при сглаживании и поиске состояний, соседние по дт=1
% состояния будут одинаковые и ложно усиливают статистику. Предлагается
% выбрасывать близко идущие повторяющиеся последовательности состояний.
if (isempty(places))
    exactk=0;
    newplaces=[];
else
    newplaces=places(1);
    exactk=1;
    for i=2:length(places)
        if (newplaces(exactk)-places(i)<(dt/2))
            % Пропускаем одинаковые состояния.
        else
            newplaces=[newplaces,places(i)];
            exactk=exactk+1;
        end
    end
end
