3樓: Originally posted by Nkxz at 2015-08-05 09:33:37
不給出程序，無法給你確切的建議。數(shù)據(jù)量不大的話幾層for循環(huán)都無所謂，只要你自己不被弄糊涂就行。如果大數(shù)據(jù)的話最好還是優(yōu)化下，matlab對大數(shù)據(jù)處理加上n個for循環(huán)，估計計算慢死了。

%function run_DFA_Du
clear
clc
regions={'Argentina';'Australia';'Benin';'Brazil';'Canada';'China';'European';...
    'GHCND';'India_resolution';'Japan';'Korea';'Mongolia';'New_Zealand';'Saudi Arabia';'SouthAfrica';'USA_resolution'};
region={'Arg';'Aus';'Ben';'Bra';'Can';'Chi';'Eur';'Ghc';'Ind';'Jap';'Kor';'Mon';'New';'Sau';'Sou';'Usa'};
lenregions=length(regions);
indir='F:\Process\Du_Haibo\Global station climate data\data_mat\';
datalist=dir([indir,'*.mat']);
datanum=length(datalist);
if lenregions~=datanum
    disp('Errors in the number of data and regions!');
    %break;
    return
end
for iregion=1:lenregions
    disp('請不要關(guān)Matlab�。�！程序運(yùn)行中。。請不要關(guān)機(jī)！��！已設(shè)置自動關(guān)機(jī)');
    %for i=1:40
    %    zprintf('_Red','請不要關(guān)Matlab�。。〕绦蜻\(yùn)行中。。。');
    %    disp('請不要關(guān)機(jī)�。�！已設(shè)置自動關(guān)機(jī)')
    %end
    outdir=['F:\Process\Du_Haibo\Global climate extremes\R-DFA\',regions{iregion},'\'];
    if ~isdir(outdir);mkdir(outdir);end
    outtfile=[outdir,region{iregion},'_'];
    DATA=load([indir,datalist(iregion).name]);
    da=DATA.data;
    clear DATA;
    lensites=length(da)-3;
    latlon=da{lensites+2};station=da{lensites+3};
    ymd0=da{lensites+1};
    sy=1961;ey=2010;leny=ey-sy+1;ymd=make_ymd(sy,ey);
    rowymd=ismember(ymd0,ymd,'rows');
    hightmax=1;hightmin=2;lowtmax=3;lowtmin=4;prcp=5;
    miss=-99.9;%missing data
    order=1;qth=2;
    %dfahightmax=cell(lensites,1);dfahightmin=cell(lensites,1);dfalowtmax=cell(lensites,1);dfalowtmin=cell(lensites,1);dfaprcp=cell(lensites,1);
    cols=[hightmax;hightmin;lowtmax;lowtmin;prcp];%here control the variables for calculating
    lencols=length(cols);
    misspercentile=99;%syear=ymd(1,1);eyear=ymd(length(ymd(:,1)),1);
    valid_site=zeros(lensites,lencols);
    num=zeros(lensites,lencols);
    subname={'hightmax','hightmin','lowtmax','lowtmin','prcp'};
    missratio=0.2;%allow the max missing ratio of 0.2 for each year
    wait=waitbar(0,'0%','Name',['Schedule for ',regions{iregion}, '...']);
    nlen=floor(lensites/50);
    if lensites<=50;ju=lensites;
    else ju=50:50:50*nlen;
        if rem(lensites,50)~=0;ju0=50:50:50*nlen;ju=[ju0 lensites];end
    end
    dfaindex=cell(lensites,lencols);
    for isite=1:lensites
        if ~isempty(find(ju==isite, 1));disp('請不要關(guān)機(jī)！�。∫言O(shè)置自動關(guān)機(jī)');end
        tic
        x=da{isite};x=round(x*10)/10;
        x(isnan(x))=miss;
        x=x(rowymd,:);
        for icol=1:lencols
            if icol==prcp;col=1;
            elseif icol==hightmax;col=2;
            elseif icol==hightmin;col=3;
            elseif icol==lowtmax;col=2;
            else  col=3;
            end%column of variables: prcp,tmax,or tmin
            if length(find(x(:,col)==miss))/length(x(:,1))>0.5
                num(isite,icol)=leny;
            else
                num(isite,icol)=missingyearjudge(ymd,x,misspercentile,miss,icol,col,sy,ey,missratio);
            end
            yymd=ymd;
            if num(isite,icol)<missratio*leny%if the ratio of missing years are less than the missratio
                valid_site(isite,icol)=1;
                datain=x(:,col);
                yymd(datain==miss,:)=[];datain(datain==miss)=[];
                if icol==prcp
                    yymd=yymd(:,2);
                    unid=unique(yymd,'rows');
                    lenid=length(unid);
                    for id=1:lenid
                        ed=datain(ismember(yymd,unid(id,1),'rows'));
                        ed(ed==0)=[];
                        if isempty(ed)
                            unid(id,3)=0;
                        else
                            unid(id,3)=mean(ed);
                        end
                    end
                    [~,meanrow]=ismember(yymd,unid(:,1),'rows');
                else
                    yymd=yymd(:,2:3);
                    unid=unique(yymd,'rows');
                    lenid=length(unid);
                    for id=1:lenid
                        unid(id,3)=mean(datain(ismember(yymd,unid(id,1:2),'rows')));
                    end
                    [~,meanrow]=ismember(yymd,unid(:,1:2),'rows');
                end
                meand=unid(meanrow,3);
                y1=cumsum(datain-meand)';
                uni=unique(datain);
                if icol==hightmax || icol==hightmin || icol==prcp
                    uni=flipud(uni);
                end
                if icol==prcp
                    thre=perthre_for_dfa([datain datain datain],miss,50,icol,col);
                    differe=abs(uni-thre);
                else
                    thre=perthre_for_dfa([datain datain datain],miss,1,icol,col);
                    differe=abs(uni-thre);
                end
                row=find(differe==min(differe));
                uni=uni(1:row(1));
                lenuni=length(uni);
                dfa=zeros(lenuni,1);
               
                warning('off','all');
                N=length(datain);
                S=1:0.25:4;S=4.^S;S=S';S=floor(S);S(S>N/4)=[];
                len1=length(S);
                fs=cell(len1,1);
                fn=cell(len1,1);
                if icol<lencols;diffv=cell(len1,1);end
                n=fix(N./S); % The number of segmention for each segmention length
                for ifn=1:len1
                    s=S(ifn);
                    n1=floor(N/s);
                    fitcoef1=zeros(n1,order+1);
                    Yn0=zeros(n(ifn),1);
                    for jfn=1:n(ifn)
                        xfit=((jfn-1)*s+1):jfn*s;
                        yfit=y1(((jfn-1)*s+1):jfn*s);
                        fn{ifn,1}(jfn,:)=xfit;
                        %fitcoef1(jfn,:)=polyfit(xfit,yfit,order);
                        %Yn0(jfn)=mean((polyval(fitcoef1(jfn,:),xfit)-yfit).^2);
                        xfit=xfit(:);yfit=yfit(:);
                        V=zeros(s,order+1);
                        V(:,order+1)=ones(s,1);
                        for ij1 = order:-1:1
                            V(:,ij1) = xfit.*V(:,ij1+1);
                        end
                        fitcoef1 = V\yfit;
                        Yn0(jfn)=mean((fitcoef1(1)*xfit+fitcoef1(2)-yfit).^2);
                        if icol>lencols;diffv{ifn,1}(jfn,:)=(fitcoef1(1)*xfit+fitcoef1(2)-yfit).^2;end
                    end
                    fs{ifn}=Yn0;
                end
                if icol>lencols
                    parfor iuni=1:lenuni
                        ivalue=uni(iuni);
                        y1;fn;fs;N;S;diffv;
                        warning('off','all');
                        len1=length(S);
                        F=zeros(len1,1);
                        if icol==lowtmax || icol==lowtmin;
                            r1=find(datain<=ivalue);
                        else
                            r1=find(datain>=ivalue);
                        end
                        for ifn=1:len1
                            s=S(ifn);
                            Yn1=fs{ifn};
                            [row,~]=find(ismember(fn{ifn,1},r1)==1);
                            unir=unique(row);
                            lenr=length(unir);
                            for ir=1:lenr
                                jfn=unir(ir);
                                rowj1=((jfn-1)*s+1):jfn*s;
                                if icol==lowtmax || icol==lowtmin
                                    subr=find(datain(rowj1)>=ivalue);
                                else
                                    subr=find(datain(rowj1)<=ivalue);%tmax
                                end
                                if isempty(subr)
                                    Yn1(jfn)=0;
                                else
                                    subdiff=diffv{ifn,1}(jfn,:);
                                    Yn1(jfn)=mean(subdiff(subr));
                                end
                            end
                            F(ifn)=mean(Yn1(Yn1~=0).^(qth/2)).^(1/qth);
                        end
                        v=zeros(len1,2);
                        v(:,2)=ones(len1,1);
                        v(:,1)=log2(S).*v(:,2);
                        p=v\log2(F);
                        dfa(iuni,1)=p(1);
                    end
                else
                    parfor iuni=1:lenuni
                        ivalue=uni(iuni);
                        y1;fn;fs;N;S;%diffv;
                        warning('off','all');
                        len1=length(S);
                        F=zeros(len1,1);
                        if icol==lowtmax || icol==lowtmin;
                            r1=find(datain<=ivalue);r2=find(datain>ivalue);
                        else
                            r1=find(datain>=ivalue);r2=find(datain<ivalue);
                        end
                        for ifn=1:len1
                            s=S(ifn);
                            Yn1=fs{ifn};
                            [row,~]=find(ismember(fn{ifn,1},r1)==1);
                            unir=unique(row);
                            lenr=length(unir);
                            for ir=1:lenr
                                jfn=unir(ir);
                                rowj1=((jfn-1)*s+1):jfn*s;
                                r01= ismember(rowj1,r2)==1;xfit1=rowj1(r01);
                                lenx=length(xfit1);
                                if lenx>1
                                    yfit1=y1(xfit1);
                                    %fitcoef=polyfit(xfit1,yfit1,order);
                                    %Yn1(jfn)=mean((polyval(fitcoef,xfit1)-yfit1).^2);
                                    xfit1=xfit1(:);yfit1=yfit1(:);
                                    V1=zeros(lenx,order+1);
                                    V1(:,order+1)=ones(lenx,1);
                                    for ij1 = order:-1:1
                                        V1(:,ij1) = xfit1.*V1(:,ij1+1);
                                    end
                                    fitcoef = V1\yfit1;
                                    Yn1(jfn)=mean((fitcoef(1)*xfit1+fitcoef(2)-yfit1).^2);
                                else
                                    Yn1(jfn)=0;
                                end
                            end
                            F(ifn)=mean(Yn1(Yn1~=0).^(qth/2)).^(1/qth);
                        end
                        v=zeros(len1,2);
                        v(:,2)=ones(len1,1);
                        v(:,1)=log2(S).*v(:,2);
                        p=v\log2(F);
                        dfa(iuni,1)=p(1);
                    end
                end
                dfainx=[dfa uni];
                dfainx(isnan(dfainx(:,1)),:)=[];dfainx(dfainx(:,1)==0,:)=[];
                if icol==prcp;dfainx(dfainx(:,2)==0,:)=[];end
                %plot(dfainx(:,2),dfainx(:,1),'.-');
                %print('-dtiff','-r300',['prcp',num2str(station(U)),'.tiff']);
                %close(figure(1));
                dfaindex{isite,icol}=dfainx;
            end
        end%icol
        toc
        waitbar(isite/lensites,wait,[num2str(isite/lensites*100,'%.2f\n'),'%']);
        if ~isempty(find(ju==isite, 1))
            for icol=lencols:lencols
                if icol==prcp;col=1;
                elseif icol==hightmax;col=2;
                elseif icol==hightmin;col=3;
                elseif icol==lowtmax;col=2;
                else  col=3;
                end%
                if ~isempty(find(ismember(cols,icol)==1, 1))
                    row1=find(valid_site(:,icol)==1);lenrow1=length(row1);
                    if isempty(row1);data=[];
                    else
                        data=dfaindex(row1,icol);data{lenrow1+1}=latlon(row1,:);data{lenrow1+2}=station(row1);
                    end
                    save([outtfile,subname{icol},'.mat'],'data');
                end
            end
        end
    end
    clear da;
    delete(wait)
end
%!shutdown -s

%%%Du Haibo, 8/7/15%%%

4樓: Originally posted by windoi at 2015-08-06 15:52:26
很多地方都可以用自定義函數(shù)或者矩陣形式計算的吧，我之前也遇到過，有的地方可以直接用矩陣形式把for循環(huán)代替，快了不少哇。