function varargout = TMDesign(varargin)
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @TMDesign_OpeningFcn, ...
                   'gui_OutputFcn',  @TMDesign_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
function TMDesign_OpeningFcn(hObject, eventdata, handles, varargin)
handles.output = hObject;
Initial(hObject, eventdata, handles) ;
guidata(hObject, handles);
function varargout = TMDesign_OutputFcn(~, ~, handles) 
varargout{1} = handles.output;

%%  Input GroundMotion and Run OpenSEES Panel
function Initial(~, ~, handles) 
global MainDir minM minK minC maxM maxK maxC RBF_value Obj_Func Neq GMNames
clc;
RBF_value=1;
MainDir=pwd;
cd(MainDir);
GMData=readtable('GM_List.txt');
set(handles.GMList,'Data',GMData{:,2:3});
set(handles.GMList,'RowName',GMData{:,1});
minM=str2double(get(handles.Mmin,'String'));
minK=str2double(get(handles.Kmin,'String'));
minC=str2double(get(handles.Cmin,'String'));
maxM=str2double(get(handles.Mmax,'String'));
maxK=str2double(get(handles.Kmax,'String'));
maxC=str2double(get(handles.Cmax,'String'));
Obj_Func=get(handles.Max_Dr,'Value');
GMNames=get(handles.GMList,'RowName');
Neq=length(GMNames);
function ReadGM_Callback(hObject, eventdata, handles)%#ok
global MainDir
[GMName,GMPath]=uigetfile('*.txt','Select the GroundMotion .txt file');
cd(MainDir);
GMfile=[GMPath,GMName];
GMData=readtable(GMfile);
set(handles.GMList,'Data',GMData{:,2:3});
set(handles.GMList,'RowName',GMData{:,1});
function GMList_CellSelectionCallback(hObject, eventdata, handles)%#ok
global Dt NStep GMName MainDir Index GMNames
Index=eventdata.Indices(1);
GM_Data=get(handles.GMList,'Data');
Dt=GM_Data(Index,1);
NStep=GM_Data(Index,2);
GMNames=get(handles.GMList,'RowName');
GMName=GMNames{Index};
set(handles.GM_Name,'String',GMName);
x=Dt:Dt:NStep*Dt;
cd(MainDir);
GM_Path=['GroundM\',GMName,'.txt'];
Acc=readtable(GM_Path);
y=Acc.Var1;
axes(handles.axes_GM);
plot(x',y,'color','red')
function noTMD_Callback(hObject, eventdata, handles)%#ok
set(handles.MTMD,'Enable','Off');
set(handles.KTMD,'Enable','Off');
set(handles.CTMD,'Enable','Off');
function withTMD_Callback(hObject, eventdata, handles)%#ok
set(handles.MTMD,'Enable','on');
set(handles.KTMD,'Enable','on');
set(handles.CTMD,'Enable','on');
function Create_Input
global Dt NStep GMName MainDir MTMD KTMD CTMD
m1=['set GM ',GMName];              
m2=['set NStep ',num2str(NStep)];           
m3=['set Dt ',num2str(Dt)];              
m4=['set MTMD ',MTMD,'; # kg.s2/m'];
m5=['set KTMD ',KTMD,'; # kg/m  '];
m6=['set CTMD ',CTMD,'; # kg.s/m'];
m=char(m1,m2,m3,m4,m5,m6);
cd(MainDir);dlmwrite('Input.tcl',m,'');
function Run_OPS_Callback(hObject, eventdata, handles)%#ok
global MTMD KTMD CTMD MainDir GMName
TMD_Status=get(handles.withTMD,'Value');
MTMD=get(handles.MTMD,'String');
KTMD=get(handles.KTMD,'String');
CTMD=get(handles.CTMD,'String');
Create_Input;
Drift=zeros(1,9);
if TMD_Status==1
    cd(MainDir);
    dos('opensees.exe Run-withTMD.tcl');
    pathoutput=[MainDir,'\output\',GMName,'\M',MTMD,'-K',KTMD,'-C',CTMD,'\EnvDrift.out'];
    Answer=dlmread(pathoutput);
    for i=1:9
        Drift(i)=Answer(3,i*2);
    end
    maxDrift=max(Drift);
    set(handles.Driftmax,'String',num2str(maxDrift));
else
    cd(MainDir);
    dos('opensees.exe Run-noTMD.tcl'); 
    pathoutput=[MainDir,'\output\',GMName,'\no-TMD\EnvDrift.out'];
    Answer=dlmread(pathoutput);
    for i=1:9
        Drift(i)=Answer(3,i*2);
    end
    maxDrift=max(Drift);
    set(handles.Driftmax,'String',num2str(maxDrift));    
end

%%  Genetic Algorithm Panel
function Max_Dr_Callback(~, ~, ~)%#ok
global Obj_Func
Obj_Func=1;
function Mean_Dr_Callback(~, ~, ~)%#ok
global Obj_Func
Obj_Func=0;
function z=Drift(x)
global  Obj_Func Index Neq ANN01 ANN02 ANN03 ANN04 ANN05 ANN06 ANN07
    if Obj_Func==1
        switch Index
            case 1
                z=sim(ANN01,x');
            case 2
                z=sim(ANN02,x');
            case 3
                z=sim(ANN03,x');
            case 4
                z=sim(ANN04,x');
            case 5
                z=sim(ANN05,x');
            case 6
                z=sim(ANN06,x');
            case 7
                z=sim(ANN07,x');                
        end
    else
        z=(sim(ANN01,x')+sim(ANN02,x')+sim(ANN03,x')+sim(ANN04,x')+sim(ANN05,x')+sim(ANN06,x')+sim(ANN07,x'))/Neq;
    end
function RunGA_Callback(hObject, eventdata, handles)%#ok
global minM minK minC maxM maxK maxC Obj_Func MainDir GMName GMNames Neq
% Define Hadaf and Paramets
Hadaf=@Drift;       % Objective Function
N=3;                % Num variables

if Obj_Func==1
    DrnoTMD=zeros(1,9);
    pathoutput=[MainDir,'\output\',GMName,'\no-TMD\EnvDrift.out'];
    Answer=dlmread(pathoutput);
    for i=1:9
        DrnoTMD(i)=Answer(3,i*2);
    end
    maxDr_noTMD=max(DrnoTMD);
    set(handles.DriftnoTMD,'String',num2str(maxDr_noTMD));
else
    maxDr_noTMD=zeros(1,Neq);
    for eq=1:Neq
        DrnoTMD=zeros(1,9);
        pathoutput=[MainDir,'\output\',GMNames{eq},'\no-TMD\EnvDrift.out'];
        Answer=dlmread(pathoutput);
        for i=1:9
            DrnoTMD(i)=Answer(3,i*2);
        end
        maxDr_noTMD(eq)=max(DrnoTMD);
    end
    set(handles.DriftnoTMD,'String',num2str(mean(maxDr_noTMD)));    
end

maxIt=str2double(get(handles.maxIt,'String'));      % Maximun Iteration
pC=str2double(get(handles.Cross,'String'));         % Crossover Percent
pM=str2double(get(handles.Mutate,'String'));        % Mutate percent
NP=str2double(get(handles.NPopulation,'String'));   %Number of Population

nC=2*round(pC*NP/2);
nM=round(pM*NP);

minN=[minM minK minC];
maxN=[maxM maxK maxC];

NSize=[1 N];            % Number of Variable Size Matrix

% Initialization
Afrad.Pos=[];
Afrad.Val=[];
pop=repmat(Afrad,NP,1);
for i=1:NP
    pop(i).Pos=unifrnd(minN,maxN,NSize);
    pop(i).Val=Hadaf(pop(i).Pos);
end
% Sort Population
Vals=[pop.Val];
[Vals,So]=sort(Vals);
pop=pop(So);
WorstVal=Vals(end);

% GA Loop
SP=10;
axes(handles.axes_GA);
for it=1:maxIt
    Vals=[pop.Val];
    P=exp(-SP*Vals/WorstVal);
    P=P/sum(P);
    
% Create Crossover Population
pop2=repmat(Afrad,nC/2,2);
for Cross=1:nC/2
    % Selection Parents
    C=cumsum(P);
    i1=find(rand<=C,1);
    i2=find(rand<=C,1);
    
%     i1=randi([1 NP]);
%     i2=randi([1 NP]);
    
    p1=pop(i1);
    p2=pop(i2);
    
    % Crossover Function
    alpha=rand;
    x1=p1.Pos;
    x2=p2.Pos;
    ch1.Pos=alpha.*x1+(1-alpha).*x2;
    ch2.Pos=alpha.*x2+(1-alpha).*x1;

    ch1.Val=Hadaf(ch1.Pos);
    ch2.Val=Hadaf(ch2.Pos); 
    
    pop2(Cross,1)=ch1;
    pop2(Cross,2)=ch2;
    
end
pop2=pop2(:);

% Create Mutation Population
 pop3=repmat(Afrad,nM,1);
 for Mu=1:nM
    i=randi([1 NP]);
    p=pop(i);
    x=p.Pos;
    j=randi([1 N]);
    
    r=unifrnd(-1,1);
    y1=x-minN;
    y2=maxN-x;
    if r<0
        x(j)=x(j)+r*y1(j);
    else
        x(j)=x(j)+r*y2(j);
    end
    q.Pos=x;
    q.Val=Hadaf(q.Pos);
    pop3(Mu)=q;
 end


% Sumation Population

    % Merge
    pop=[pop
        pop2
        pop3];%#ok
    
    % Sort
    Vals=[pop.Val];
    [Vals, so]=sort(Vals);
    pop=pop(so);
    
    % Delete Extra
    pop=pop(1:NP);
    BestVal(it)=Vals(1);%#ok
    
    % Update WorstVal
    WorstVal=max(WorstVal,Vals(end));
    
% Results
BestSol=pop(1);
pause(0.0001);
semilogy(BestVal,'r','LineWidth',2);
hold off;
xlabel('Generation (Iteration)');

MKC=BestSol.Pos;
set(handles.Mopt,'String',num2str(round(MKC(1))));
set(handles.Kopt,'String',num2str(round(MKC(2))));
set(handles.Copt,'String',num2str(round(MKC(3))));

maxDr_wTMD=Vals(1);
Opt=100*(maxDr_noTMD-maxDr_wTMD)/maxDr_noTMD;
set(handles.DriftwTMD,'String',num2str(maxDr_wTMD));
if Opt<0
    set(handles.Optimization,'ForegroundColor',[0.6 0 0]);
else
    set(handles.Optimization,'ForegroundColor',[0 0.6 0]);    
end
set(handles.Optimization,'String',num2str(Opt));

end

%%  Artifitial Neural Network Panel
function RBFtype_Callback(hObject, eventdata, handles)%#ok
global RBF_value
RBF_value=1;
set(handles.Epoch_txt,'Enable','off');
set(handles.Epoch_num,'Enable','off');

set(handles.Valid_txt,'Enable','off');
set(handles.Valid_num,'Enable','off');

set(handles.Spread_txt,'Enable','on');
set(handles.Spread_num,'Enable','on');

set(handles.TrFunc_txt,'Enable','off');
set(handles.TrFunc_list,'Enable','off');
function MLPtype_Callback(hObject, eventdata, handles)%#ok
global RBF_value
RBF_value=0;
set(handles.Epoch_txt,'Enable','on');
set(handles.Epoch_num,'Enable','on');

set(handles.Valid_txt,'Enable','on');
set(handles.Valid_num,'Enable','on');

set(handles.Spread_txt,'Enable','off');
set(handles.Spread_num,'Enable','off');

set(handles.TrFunc_txt,'Enable','on');
set(handles.TrFunc_list,'Enable','on');
function Training_Callback(hObject, eventdata, handles)%#ok
global RBF_value minM maxM minK maxK minC maxC MainDir GMName net Index ...
    ANN01 ANN02 ANN03 ANN04 ANN05 ANN06 ANN07
    linM=str2double(get(handles.linM,'String'));
    linK=str2double(get(handles.linK,'String'));
    linC=str2double(get(handles.linC,'String'));
    
    Counter=0;
% Create Input and Target
    Input=zeros(3,linM*linK*linC);
    Target=zeros(1,linM*linK*linC);
    
    gamM=(maxM-minM)/(linM-1);
    gamK=(maxK-minK)/(linK-1);
    gamC=(maxC-minC)/(linC-1);
    
    for Mt=minM:gamM:maxM
        for Kt=minK:gamK:maxK
            for Ct=minC:gamC:maxC
                Dr=zeros(1,9);
                Counter=Counter+1;
                Input(1,Counter)=Mt;
                Input(2,Counter)=Kt;
                Input(3,Counter)=Ct; 
   pathoutput=[MainDir,'\output\',GMName,'\M',num2str(Mt),'-K',num2str(Kt),'-C',num2str(Ct),'\Envdrift.out'];
                Answer=dlmread(pathoutput);
                for i=1:9
                    Dr(i)=Answer(3,i*2);
                end
                Target(1,Counter)=max(Dr);
            end
        end
    end
        MaxN=str2double(get(handles.Neuron_num,'String'));
if RBF_value==1
    Goal=0;%#ok
    Spread=str2double(get(handles.Spread_num,'String'));
    
%     net=newrb(Input,Target, Goal,Spread,MaxN);
    net=newrbe(Input,Target,Spread);

    txtTitle='max Drift in SAC9 with RBF Network';
else
    numtrain=get(handles.TrFunc_list,'Value');
    switch numtrain
        case 1
            trainFcn='trainlm';
        case 2
            trainFcn='trainbr';            
        case 3
            trainFcn='trainscg';            
    end
            
    net=fitnet(MaxN,trainFcn);
    
    net.TrainParam.Epochs=str2double(get(handles.Epoch_num,'String'));
    net.TrainParam.Max_Fail=str2double(get(handles.Valid_num,'String'));
    
    net=train(net,Input,Target);
    txtTitle='max Drift in SAC9 with MLP Network';
end
    Output=net(Input);
    performance=perform(net,Target,Output);%#ok
   figure
   plot(Output,'b');
   hold on;
   plot(Target,'r:');
   legend('Output','Target');
   title(txtTitle);
   
   switch Index
       case 1
           ANN01=net;
       case 2
           ANN02=net; 
       case 3
           ANN03=net;
       case 4
           ANN04=net;
       case 5
           ANN05=net;
       case 6
           ANN06=net;
       case 7
           ANN07=net;           
   end
   disp([GMName,'  Trained']);
   set(handles.Training,'Value',0);
function TrainAll_Callback(hObject, eventdata, handles)%#ok
global Neq Index GMNames GMName
for i=1:Neq
    Index=i;
    GMName=GMNames{i};
    Training_Callback(hObject, eventdata, handles)
end
set(handles.TrainAll,'value',0);
