% Delivered to the Government in Accordance with Contract FA8650-10-C-7080
% 
%  (c) 2011 United Technologies Corp.
% Developed with the sponsorship of the Defense Advanced Research Projects Agency (DARPA).
%  
% Permission is hereby granted, free of charge, to any person obtaining a copy of this data, including any software or models in source or binary form, as well as any drawings, specifications, and documentation (collectively "the Data"), to deal in the Data without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Data, and to permit persons to whom the Data is furnished to do so, subject to the following conditions:
%  
% The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Data.
%  
% THE DATA IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS, SPONSORS, DEVELOPERS, CONTRIBUTORS, OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE DATA OR THE USE OR OTHER DEALINGS IN THE DATA.
% 


function demo_dse1(dse_iter,dir_name,obj_error,arch_num)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%IBM copyrighted material (2011)
%IBM Confidential
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fclose all;
format long;
L_d=dse_iter;
eps_dse = obj_error; %in percentages
arch=arch_num;
source_name = dir_name;
destination_name = 'reference_arch';
copyfile(source_name, destination_name);
%D = dir('destination_name');
%D.name;
dirpath=strcat(destination_name,'\',dir_name);
for i=1:arch
    bpr=0;
    TSFC=0;
    ma_core=0;
    rc10=0;
    ind_arch = num2str(i,10);
    arch_name_read=strcat(dirpath,'\Arch_',ind_arch,'.mod');
    arch_name_write=strcat(dirpath,'\Arch_dse',ind_arch,'.mod');
    fid_wr = fopen('find_point.txt', 'w+');
    fprintf(fid_wr, '%s  \n', 'reset;');
    model_name=strcat('model', 32,'"',arch_name_read,'";'); %32 is ASCII space - to ovveride Matlab's feature.
    data_name=strcat('data', 32,'"',dirpath,'\Arch_',ind_arch,'.dat";');
    fprintf(fid_wr, '%s  \n', model_name);
    fprintf(fid_wr, '%s  \n', data_name);
    %fprintf(fid_wr, '%s  \n', 'option solver bonmin;');
    fprintf(fid_wr, '%s  \n', 'option solver ipopt;');
    fprintf(fid_wr, '%s  \n', 'solve;');
    fprintf(fid_wr, '%s  \n', 'display TSFC, bpr, ma_core, rc10, objective;');
    fclose(fid_wr);
    system('ampl find_point.txt > fp_results.txt');
    
    fid_rd = fopen('fp_results.txt'); %read results file
    tline_rd = fgetl(fid_rd);
    while ischar(tline_rd)
        tstr_rd=char(tline_rd);
        parsed_line=regexp(tstr_rd,'\s+', 'split');
    
        if regexp(parsed_line{1},'TSFC')
            TSFC(1)=str2num(parsed_line{3});
        end
   
        if regexp(parsed_line{1},'ma_core')
            ma_core(1)=str2num(parsed_line{3});
        end
    
        if regexp(parsed_line{1},'rc10')
            rc10(1)=str2num(parsed_line{3});
        end
   
        if regexp(parsed_line{1},'bpr')
            bpr(1)=str2num(parsed_line{3});
        end
        tline_rd = fgetl(fid_rd);
    end
    fclose(fid_rd);
    
%%%%%%%%%%%%%%%%%%%% DSE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    for k=1:L_d
        cost_len = length(bpr);
        eps_dse_str=num2str((TSFC(1)*eps_dse)/100,10);
        fid = fopen(arch_name_read);
        fid_wr = fopen(arch_name_write, 'w+');
        tstr_new = 'minimize objective: mfuel/Thto;';
        tstr_next = 'maximize objective: gamma_alg;';
        tline = fgetl(fid);
        first_line=tline;
        fprintf(fid_wr, '%s  \n', first_line);


    while ischar(tline)
        tstr=char(tline);
        m=regexp(tstr,tstr_new);
        t=regexp(tstr,tstr_next);
        
        if m == 1
            TSFC_1 = num2str(TSFC(1),10);
            printout_model_TSFC = strcat('subject to TSFC_1 : TSFC <= (1+',eps_dse_str,')*',TSFC_1,';');
            printout_model_gamma = strcat('var gamma_alg;');            
            fprintf(fid_wr, '%s  \n', printout_model_TSFC);  
            fprintf(fid_wr, '%s  \n', printout_model_gamma);
            bpr_j = num2str(bpr(cost_len),10);
            ma_core_j = num2str(ma_core(cost_len),10);
            rc10_j = num2str(rc10(cost_len),10);
            ind_i = num2str(cost_len,10);
            printout_model_square_sum_3 = strcat('((',bpr_j,'-bpr)/22)^2 + ((',ma_core_j,'-ma_core)/30)^2 +((',rc10_j,'-rc10)/28)^2');
            printout_gamma_3 = strcat('subject to gamma_',ind_i,' : gamma_alg <= sqrt(',printout_model_square_sum_3,');');
            fprintf(fid_wr, '%s  \n', printout_gamma_3);
            fprintf(fid_wr, '%s  \n', 'maximize objective: gamma_alg;');
            tline = fgetl(fid);
            continue;
        end
        if t == 1
                bpr_j = num2str(bpr(cost_len),10);
                ma_core_j = num2str(ma_core(cost_len),10);
                rc10_j = num2str(rc10(cost_len),10);
                ind_i = num2str(cost_len,10);
                printout_model_square_sum_3 = strcat('((',bpr_j,'-bpr)/22)^2 + ((',ma_core_j,'-ma_core)/30)^2 +((',rc10_j,'-rc10)/28)^2'); %%ma_core=85, range_bpr=28, range_rc=28
                printout_gamma_3 = strcat('subject to gamma_',ind_i,' : gamma_alg <= sqrt(',printout_model_square_sum_3,');');
                fprintf(fid_wr, '%s  \n', printout_gamma_3);
                fprintf(fid_wr, '%s  \n', 'maximize objective: gamma_alg;');
        elseif regexp(tstr,char(first_line));
        else
            fprintf(fid_wr, '%s  \n', tline);
        end
        tline = fgetl(fid);
    end

    fclose(fid_wr);
    fclose(fid);
    fclose all;
    %if i == 9 & k == 1
    %    pause(0.25);
    %end
    if i == 12 & k == 1
        pause(1);
    end
    movefile(arch_name_write, arch_name_read,'f');
    system('ampl find_point.txt > fp_results.txt');
    fid_rd = fopen('fp_results.txt'); %read results file
    tline_rd = fgetl(fid_rd);
    while ischar(tline_rd)
        tstr_rd=char(tline_rd);
        parsed_line=regexp(tstr_rd,'\s+', 'split');
    
        if regexp(parsed_line{1},'TSFC')
            TSFC(cost_len + 1)=str2num(parsed_line{3});
        end
   
        if regexp(parsed_line{1},'ma_core')
            ma_core(cost_len + 1)=str2num(parsed_line{3});
        end
    
        if regexp(parsed_line{1},'rc10')
            rc10(cost_len + 1)=str2num(parsed_line{3});
        end
   
        if regexp(parsed_line{1},'bpr')
            bpr(cost_len + 1)=str2num(parsed_line{3});
        end
        tline_rd = fgetl(fid_rd);
    end
    fclose(fid_rd);

    %bpr
    %rc10
    %ma_core
    subplot(4,3,i),
    %figure(i);
    if TSFC(1) < 0.65
        stem3(ma_core(1),bpr(1),rc10(1),'filled','ro-'); axis([35 130 0 30 0 30]); grid on; hold on; drawnow;
    elseif TSFC(1) < 0.82
        stem3(ma_core(1),bpr(1),rc10(1),'ro-'); axis([35 130 0 30 0 30]); grid on; hold on; drawnow;
    else
        stem3(ma_core(1),bpr(1),rc10(1),'ro:'); axis([35 130 0 30 0 30]); grid on; hold on; drawnow;
    end
    
    if TSFC(cost_len+1) < 0.65
        stem3(ma_core(cost_len+1),bpr(cost_len+1),rc10(cost_len+1),'filled','o-'); axis([35 130 0 30 0 30]); grid on; hold on; drawnow; %axis([40 100 2 30 2 30]);
    elseif TSFC(cost_len+1) < 0.82
        stem3(ma_core(cost_len+1),bpr(cost_len+1),rc10(cost_len+1),'o-'); axis([35 130 0 30 0 30]); grid on; hold on; drawnow; %axis([40 100 2 30 2 30]);
    else
        stem3(ma_core(cost_len+1),bpr(cost_len+1),rc10(cost_len+1),'o:'); axis([35 130 0 30 0 30]); grid on; hold on; drawnow; %axis([40 100 2 30 2 30]);
    end
    
    xlabel('CORE MASS FLOW');
    ylabel('BPR');
    zlabel('RC');
    title_arch_num=num2str(i,10);
    optimal=num2str(TSFC(1),10);
    arch_name=strcat('Arch(',title_arch_num,'); Opt. Fuel Consumption:', 32,optimal);
    title(arch_name);
    if i==2
        text(200,120,1,'Performance(TSFC): Best-Solid,Filled; Average-Solid,Empty; Worst-Dotted,Empty','fontsize',14);
    end   

    end
end