% Initialization for Adaptive DO solvers
% This is basically code that is common amongst the Adaptive DO solvers for
% automatic initialization. That is, initialization that is not handled through
% the setup script is handled through this script.

if ~exist('ui','var') && (app.S(1)>0 || length(app.S)>1)
    %If app.S=[0 #] we set app.S=1 and start with a very small variance
    
    %Choose the number of modes for the size of the basis from which new
    %unstable modes will be chosen. 
    NUMMODESINBASIS=round(min(50,max(max(app.S)^2/2,30)));
    d0=fileparts([pwd,filesep]);
    d0=d0(1:strfind(d0,'trunk')+4);
    
    if length(app.S)>1 %adaptive
        if ~isfield(app,'PlotScript')
            PlotScript = 'DOAdaptGenericPlots';
        end
        %Need 1 extra mode to have a properly normalized basis
        [ui, vi, var] = Stoch_Initial_Perturb(XU,YU,XV,YV,Nbcs,...
            length(bcsDP)+length(bcsOP),NodeP,NUMMODESINBASIS+1);
        if isfield(app,'var')
            var = var*app.var(1);
        else
            var = var*1e-6;
        end
        %XU and YU etc. contain points that may be masked out, so we have
        %to remove these. Also remove the extra mode.
        ui = ui([[1:Nbcs]';Nbcs+idsu],1:NUMMODESINBASIS);
        vi = vi([[1:Nbcs]';Nbcs+idsv],1:NUMMODESINBASIS);
        existbasis = length(dir([d0,'/Save/DOAdaptSave/Basis/orth_basis*']));
        save([d0,sprintf('/Save/DOAdaptSave/Basis/orth_basis%0.3d.mat',existbasis)],'ui','vi');
        close all
    else %non-adaptive
        if ~isfield(app,'PlotScript')
            PlotScript = 'DOAdaptGenericPlots';
        end
        %Need to have one extra mode to have a proper orthonormalization
        [ui, vi, var] = Stoch_Initial_Perturb(XU,YU,XV,YV,Nbcs,...
            length(bcsDP)+length(bcsOP),NodeP,app.S + 1);
        var
        if isfield(app,'var')
            var = var*app.var(1);
        else
            var = var*1e-6;
        end
        close all;
        %XU and YU etc. contain points that may be masked out, so we have
        %to remove these.
        ui = ui([[1:Nbcs]';Nbcs+idsu],:);
        vi = vi([[1:Nbcs]';Nbcs+idsv],:);
    end

    %If the number of Monte-Carlo samples was not set by the setup script, 
    %arbitrarily choose 10,000
    if ~isfield(app,'MC'),app.MC=10000;end
    
    if (app.S(1)==0) 
        %For app.S = [0 #] case we need to set small variance, and initialize 
        %arrays properly
        app.S(1)=1;
        var(1)=1e-6;
        fui = zeros(length(uid),app.S(1));
        fvi = zeros(length(vid),app.S(1));
        fuij= zeros(length(uid),app.S(1)^2);
        fvij= zeros(length(vid),app.S(1)^2);
    end
    j=1;
    randn('state',0);
    YYt=zeros(app.MC,app.S(1)+1);
    for i=1:app.S+1
        YYt(:,j)=normrnd(0,sqrt(var(i)),app.MC,1);
        j=j+1;
    end
    
    ui = ui(:,1:app.S(1)+1);
    vi = vi(:,1:app.S(1)+1);
    Pi=zeros(length(P),app.S(1)+1);
    %Ortho-Normalization of basis fields to avoid round-off errors
    [YYt ui vi Pi] = DOorthnorm(app, dx, dy, YYt, ui, uid, vi, vid, Pi, pid);
    
    %Only select the first S modes
    ui = ui(:,1:app.S(1));
    vi = vi(:,1:app.S(1));
    Pi=zeros(length(P),app.S(1));
    YYt = YYt(:,1:app.S(1));
end
