diff --git a/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/Demo_Test_DnCNN_DAG.m b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/Demo_Test_DnCNN_DAG.m
new file mode 100644
index 00000000..fe7d5bbe
--- /dev/null
+++ b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/Demo_Test_DnCNN_DAG.m
@@ -0,0 +1,120 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% @article{zhang2017beyond,
+%   title={Beyond a {Gaussian} denoiser: Residual learning of deep {CNN} for image denoising},
+%   author={Zhang, Kai and Zuo, Wangmeng and Chen, Yunjin and Meng, Deyu and Zhang, Lei},
+%   journal={IEEE Transactions on Image Processing},
+%   year={2017},
+%   volume={26}, 
+%   number={7}, 
+%   pages={3142-3155}, 
+% }
+
+% by Kai Zhang (1/2018)
+% cskaizhang@gmail.com
+% https://github.com/cszn
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% clear; clc;
+
+%% testing set
+addpath(fullfile('utilities'));
+
+folderModel = 'model';
+folderTest  = 'testsets';
+folderResult= 'results';
+imageSets   = {'BSD68','Set12'}; % testing datasets
+setTestCur  = imageSets{2};      % current testing dataset
+
+
+showresult  = 1;
+gpu         = 1;
+
+
+noiseSigma  = 25;
+
+% load model
+epoch       = 50;
+
+modelName   = 'DnCNN';
+
+% case one: for the model in 'data/model'
+%load(fullfile('data',folderModel,[modelName,'-epoch-',num2str(epoch),'.mat']));
+
+% case two: for the model in 'utilities'
+load(fullfile('utilities',[modelName,'-epoch-',num2str(epoch),'.mat']));
+
+
+
+
+net = dagnn.DagNN.loadobj(net) ;
+
+net.removeLayer('loss') ;
+out1 = net.getVarIndex('prediction') ;
+net.vars(net.getVarIndex('prediction')).precious = 1 ;
+
+net.mode = 'test';
+
+if gpu
+    net.move('gpu');
+end
+
+% read images
+ext         =  {'*.jpg','*.png','*.bmp'};
+filePaths   =  [];
+for i = 1 : length(ext)
+    filePaths = cat(1,filePaths, dir(fullfile(folderTest,setTestCur,ext{i})));
+end
+
+folderResultCur       =  fullfile(folderResult, [setTestCur,'_',int2str(noiseSigma)]);
+if ~isdir(folderResultCur)
+    mkdir(folderResultCur)
+end
+
+
+% PSNR and SSIM
+PSNRs = zeros(1,length(filePaths));
+SSIMs = zeros(1,length(filePaths));
+
+
+for i = 1 : length(filePaths)
+    
+    % read image
+    label = imread(fullfile(folderTest,setTestCur,filePaths(i).name));
+    [~,nameCur,extCur] = fileparts(filePaths(i).name);
+    [w,h,c]=size(label);
+    if c==3
+        label = rgb2gray(label);
+    end
+    
+    % add additive Gaussian noise
+    randn('seed',0);
+    noise = noiseSigma/255.*randn(size(label));
+    input = im2single(label) + single(noise);
+    
+    if gpu
+        input = gpuArray(input);
+    end
+    net.eval({'input', input}) ;
+    % output (single)
+    output = gather(squeeze(gather(net.vars(out1).value)));
+    
+    
+    % calculate PSNR and SSIM
+    [PSNRCur, SSIMCur] = Cal_PSNRSSIM(label,im2uint8(output),0,0);
+    if showresult
+        imshow(cat(2,im2uint8(input),im2uint8(label),im2uint8(output)));
+        title([filePaths(i).name,'    ',num2str(PSNRCur,'%2.2f'),'dB','    ',num2str(SSIMCur,'%2.4f')])
+        imwrite(im2uint8(output), fullfile(folderResultCur, [nameCur, '_' int2str(noiseSigma),'_PSNR_',num2str(PSNRCur*100,'%4.0f'), extCur] ));
+        drawnow;
+       % pause()
+    end
+    PSNRs(i) = PSNRCur;
+    SSIMs(i) = SSIMCur;
+end
+
+
+disp([mean(PSNRs),mean(SSIMs)]);
+
+
+
+
diff --git a/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/Demo_Train_DnCNN_DAG.m b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/Demo_Train_DnCNN_DAG.m
new file mode 100644
index 00000000..55f911a0
--- /dev/null
+++ b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/Demo_Train_DnCNN_DAG.m
@@ -0,0 +1,61 @@
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% @article{zhang2017beyond,
+%   title={Beyond a {Gaussian} denoiser: Residual learning of deep {CNN} for image denoising},
+%   author={Zhang, Kai and Zuo, Wangmeng and Chen, Yunjin and Meng, Deyu and Zhang, Lei},
+%   journal={IEEE Transactions on Image Processing},
+%   year={2017},
+%   volume={26}, 
+%   number={7}, 
+%   pages={3142-3155}, 
+% }
+
+% by Kai Zhang (1/2018)
+% cskaizhang@gmail.com
+% https://github.com/cszn
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% The training data is generated by '[imdb] = generatepatches;' in line 126 of 'DnCNN_train_dag.m'.
+
+rng('default')
+
+addpath('utilities');
+%-------------------------------------------------------------------------
+% Configuration
+%-------------------------------------------------------------------------
+opts.modelName        = 'DnCNN'; % model name
+opts.learningRate     = [logspace(-3,-3,22) logspace(-4,-4,105)];% you can change the learning rate
+opts.batchSize        = 128; % 
+opts.gpus             = [1]; 
+opts.numSubBatches    = 2;
+
+% solver
+opts.solver           = 'Adam'; % global
+opts.derOutputs       = {'objective',1} ;
+
+opts.backPropDepth    = Inf;
+%-------------------------------------------------------------------------
+%   Initialize model
+%-------------------------------------------------------------------------
+
+net  = feval([opts.modelName,'_Init']);
+
+%-------------------------------------------------------------------------
+%   Train
+%-------------------------------------------------------------------------
+
+[net, info] = DnCNN_train_dag(net,  ...
+    'learningRate',opts.learningRate, ...
+    'derOutputs',opts.derOutputs, ...
+    'numSubBatches',opts.numSubBatches, ...
+    'backPropDepth',opts.backPropDepth, ...
+    'solver',opts.solver, ...
+    'batchSize', opts.batchSize, ...
+    'modelname', opts.modelName, ...
+    'gpus',opts.gpus) ;
+
+
+
+
+
+
diff --git a/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/DnCNN_Init.m b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/DnCNN_Init.m
new file mode 100644
index 00000000..62eff08a
--- /dev/null
+++ b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/DnCNN_Init.m
@@ -0,0 +1,213 @@
+function net = DnCNN_Init()
+
+% by Kai Zhang (1/2018)
+% cskaizhang@gmail.com
+% https://github.com/cszn
+
+% Create DAGNN object
+net = dagnn.DagNN();
+
+% conv + relu
+blockNum = 1;
+inVar = 'input';
+channel= 1; % grayscale image
+dims   = [3,3,channel,64];
+pad    = [1,1];
+stride = [1,1];
+lr     = [1,1];
+[net, inVar, blockNum] = addConv(net, blockNum, inVar, dims, pad, stride, lr);
+[net, inVar, blockNum] = addReLU(net, blockNum, inVar);
+
+for i = 1:15
+    % conv + bn + relu
+    dims   = [3,3,64,64];
+    pad    = [1,1];
+    stride = [1,1];
+    lr     = [1,0];
+    [net, inVar, blockNum] = addConv(net, blockNum, inVar, dims, pad, stride, lr);
+    n_ch   = dims(4);
+    [net, inVar, blockNum] = addBnorm(net, blockNum, inVar, n_ch);
+    [net, inVar, blockNum] = addReLU(net, blockNum, inVar);
+end
+
+% conv
+dims   = [3,3,64,channel];
+pad    = [1,1];
+stride = [1,1];
+lr     = [1,0]; % or [1,1], it does not influence the results
+[net, inVar, blockNum] = addConv(net, blockNum, inVar, dims, pad, stride, lr);
+
+% sum
+inVar = {inVar,'input'};
+[net, inVar, blockNum] = addSum(net, blockNum, inVar);
+
+outputName = 'prediction';
+net.renameVar(inVar,outputName)
+
+% loss
+net.addLayer('loss', dagnn.Loss('loss','L2'), {'prediction','label'}, {'objective'},{});
+net.vars(net.getVarIndex('prediction')).precious = 1;
+
+
+end
+
+
+
+
+% Add a Concat layer
+function [net, inVar, blockNum] = addConcat(net, blockNum, inVar)
+
+outVar   = sprintf('concat%d', blockNum);
+layerCur = sprintf('concat%d', blockNum);
+
+block = dagnn.Concat('dim',3);
+net.addLayer(layerCur, block, inVar, {outVar},{});
+
+inVar = outVar;
+blockNum = blockNum + 1;
+end
+
+
+% Add a loss layer
+function [net, inVar, blockNum] = addLoss(net, blockNum, inVar)
+
+outVar   = 'objective';
+layerCur = sprintf('loss%d', blockNum);
+
+block    = dagnn.Loss('loss','L2');
+net.addLayer(layerCur, block, inVar, {outVar},{});
+
+inVar    = outVar;
+blockNum = blockNum + 1;
+end
+
+
+% Add a sum layer
+function [net, inVar, blockNum] = addSum(net, blockNum, inVar)
+
+outVar   = sprintf('sum%d', blockNum);
+layerCur = sprintf('sum%d', blockNum);
+
+block    = dagnn.Sum();
+net.addLayer(layerCur, block, inVar, {outVar},{});
+
+inVar    = outVar;
+blockNum = blockNum + 1;
+end
+
+
+% Add a relu layer
+function [net, inVar, blockNum] = addReLU(net, blockNum, inVar)
+
+outVar   = sprintf('relu%d', blockNum);
+layerCur = sprintf('relu%d', blockNum);
+
+block    = dagnn.ReLU('leak',0);
+net.addLayer(layerCur, block, {inVar}, {outVar},{});
+
+inVar    = outVar;
+blockNum = blockNum + 1;
+end
+
+
+% Add a bnorm layer
+function [net, inVar, blockNum] = addBnorm(net, blockNum, inVar, n_ch)
+
+trainMethod = 'adam';
+outVar   = sprintf('bnorm%d', blockNum);
+layerCur = sprintf('bnorm%d', blockNum);
+
+params={[layerCur '_g'], [layerCur '_b'], [layerCur '_m']};
+net.addLayer(layerCur, dagnn.BatchNorm('numChannels', n_ch), {inVar}, {outVar},params) ;
+
+pidx = net.getParamIndex({[layerCur '_g'], [layerCur '_b'], [layerCur '_m']});
+b_min                           = 0.025;
+net.params(pidx(1)).value       = clipping(sqrt(2/(9*n_ch))*randn(n_ch,1,'single'),b_min);
+net.params(pidx(1)).learningRate= 1;
+net.params(pidx(1)).weightDecay = 0;
+net.params(pidx(1)).trainMethod = trainMethod;
+
+net.params(pidx(2)).value       = zeros(n_ch, 1, 'single');
+net.params(pidx(2)).learningRate= 1;
+net.params(pidx(2)).weightDecay = 0;
+net.params(pidx(2)).trainMethod = trainMethod;
+
+net.params(pidx(3)).value       = [zeros(n_ch,1,'single'), 0.01*ones(n_ch,1,'single')];
+net.params(pidx(3)).learningRate= 1;
+net.params(pidx(3)).weightDecay = 0;
+net.params(pidx(3)).trainMethod = 'average';
+
+inVar    = outVar;
+blockNum = blockNum + 1;
+end
+
+
+% add a ConvTranspose layer
+function [net, inVar, blockNum] = addConvt(net, blockNum, inVar, dims, crop, upsample, lr)
+opts.cudnnWorkspaceLimit = 1024*1024*1024*2; % 2GB
+convOpts = {'CudnnWorkspaceLimit', opts.cudnnWorkspaceLimit} ;
+trainMethod = 'adam';
+
+outVar      = sprintf('convt%d', blockNum);
+
+layerCur    = sprintf('convt%d', blockNum);
+
+convBlock = dagnn.ConvTranspose('size', dims, 'crop', crop,'upsample', upsample, ...
+    'hasBias', true, 'opts', convOpts);
+
+net.addLayer(layerCur, convBlock, {inVar}, {outVar},{[layerCur '_f'], [layerCur '_b']});
+
+f  = net.getParamIndex([layerCur '_f']) ;
+sc = sqrt(2/(dims(1)*dims(2)*dims(4))) ; %improved Xavier
+net.params(f).value        = sc*randn(dims, 'single');
+net.params(f).learningRate = lr(1);
+net.params(f).weightDecay  = 1;
+net.params(f).trainMethod  = trainMethod;
+
+f = net.getParamIndex([layerCur '_b']) ;
+net.params(f).value        = zeros(dims(3), 1, 'single');
+net.params(f).learningRate = lr(2);
+net.params(f).weightDecay  = 1;
+net.params(f).trainMethod  = trainMethod;
+
+inVar    = outVar;
+blockNum = blockNum + 1;
+end
+
+
+% add a Conv layer
+function [net, inVar, blockNum] = addConv(net, blockNum, inVar, dims, pad, stride, lr)
+opts.cudnnWorkspaceLimit = 1024*1024*1024*2; % 2GB
+convOpts    = {'CudnnWorkspaceLimit', opts.cudnnWorkspaceLimit} ;
+trainMethod = 'adam';
+
+outVar      = sprintf('conv%d', blockNum);
+layerCur    = sprintf('conv%d', blockNum);
+
+convBlock   = dagnn.Conv('size', dims, 'pad', pad,'stride', stride, ...
+    'hasBias', true, 'opts', convOpts);
+
+net.addLayer(layerCur, convBlock, {inVar}, {outVar},{[layerCur '_f'], [layerCur '_b']});
+
+f = net.getParamIndex([layerCur '_f']) ;
+sc = sqrt(2/(dims(1)*dims(2)*max(dims(3), dims(4)))) ; %improved Xavier
+net.params(f).value        = sc*randn(dims, 'single') ;
+net.params(f).learningRate = lr(1);
+net.params(f).weightDecay  = 1;
+net.params(f).trainMethod  = trainMethod;
+
+f = net.getParamIndex([layerCur '_b']) ;
+net.params(f).value        = zeros(dims(4), 1, 'single');
+net.params(f).learningRate = lr(2);
+net.params(f).weightDecay  = 1;
+net.params(f).trainMethod  = trainMethod;
+
+inVar    = outVar;
+blockNum = blockNum + 1;
+end
+
+
+function A = clipping(A,b)
+A(A>=0&A<b) = b;
+A(A<0&A>-b) = -b;
+end
diff --git a/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/DnCNN_train_dag.m b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/DnCNN_train_dag.m
new file mode 100644
index 00000000..a9a86542
--- /dev/null
+++ b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/DnCNN_train_dag.m
@@ -0,0 +1,536 @@
+function [net,stats] = DnCNN_train_dag(net, varargin)
+
+%CNN_TRAIN_DAG Demonstrates training a CNN using the DagNN wrapper
+%    CNN_TRAIN_DAG() is similar to CNN_TRAIN(), but works with
+%    the DagNN wrapper instead of the SimpleNN wrapper.
+
+% Copyright (C) 2014-16 Andrea Vedaldi.
+% All rights reserved.
+%
+% This file is part of the VLFeat library and is made available under
+% the terms of the BSD license (see the COPYING file).
+
+
+%%%-------------------------------------------------------------------------
+%%% solvers: SGD(default) and Adam with(default)/without gradientClipping
+%%%-------------------------------------------------------------------------
+
+%%% solver: Adam
+%%% opts.solver = 'Adam';
+opts.beta1   = 0.9;
+opts.beta2   = 0.999;
+opts.alpha   = 0.01;
+opts.epsilon = 1e-8;
+
+%%% solver: SGD
+opts.solver = 'SGD';
+opts.learningRate = 0.01;
+opts.weightDecay  = 0.0005;
+opts.momentum     = 0.9 ;
+
+%%% GradientClipping
+opts.gradientClipping = false;
+opts.theta            = 0.005;
+
+%%%-------------------------------------------------------------------------
+%%%  setting for dag
+%%%-------------------------------------------------------------------------
+
+opts.conserveMemory = true;
+opts.mode = 'normal';
+opts.cudnn = true ;
+opts.backPropDepth = +inf ;
+opts.skipForward = false;
+opts.numSubBatches = 1;
+%%%-------------------------------------------------------------------------
+%%%  setting for model
+%%%-------------------------------------------------------------------------
+
+opts.batchSize = 128 ;
+opts.gpus = [];
+opts.numEpochs = 300 ;
+opts.modelName   = 'model';
+opts.expDir = fullfile('data',opts.modelName) ;
+opts.numberImdb   = 1;
+opts.imdbDir      = opts.expDir;
+opts.derOutputs = {'objective', 1} ;
+
+%%%-------------------------------------------------------------------------
+%%%  update settings
+%%%-------------------------------------------------------------------------
+
+opts = vl_argparse(opts, varargin);
+opts.numEpochs = numel(opts.learningRate);
+
+if ~exist(opts.expDir, 'dir'), mkdir(opts.expDir) ; end
+
+%%% load training data
+% opts.imdbPath    = fullfile(opts.imdbDir, 'imdb.mat');
+% imdb = load(opts.imdbPath) ;
+
+% if  mod(epoch,5)~=1 && isfield(imdb,'set') ~= 0
+%
+% else
+%     clear imdb;
+%     [imdb] = generatepatches;
+% end
+%
+% opts.train = find(imdb.set==1);
+
+opts.continue = true;
+opts.prefetch = true;
+opts.saveMomentum = false;
+opts.nesterovUpdate = false ;
+
+opts.profile = false ;
+opts.parameterServer.method = 'mmap' ;
+opts.parameterServer.prefix = 'mcn' ;
+
+
+opts.extractStatsFn = @extractStats ;
+opts = vl_argparse(opts, varargin) ;
+
+% -------------------------------------------------------------------------
+%                                                            Initialization
+% -------------------------------------------------------------------------
+opts.train = true;
+evaluateMode = isempty(opts.train) ;
+
+% -------------------------------------------------------------------------
+%                                                        Train
+% -------------------------------------------------------------------------
+
+
+modelPath = @(ep) fullfile(opts.expDir, sprintf([opts.modelName,'-epoch-%d.mat'], ep));
+start = findLastCheckpoint(opts.expDir,opts.modelName) ;
+
+if start>=1
+    fprintf('%s: resuming by loading epoch %d\n', mfilename, start) ;
+    [net] = loadState(modelPath(start)) ;
+end
+state = [] ;
+
+
+% for iobj = numel(opts.derOutputs)
+net.vars(net.getVarIndex(opts.derOutputs{1})).precious = 1;
+% end
+
+imdb = [];
+
+for epoch=start+1:opts.numEpochs
+    
+    if  mod(epoch,10)~=1 && isfield(imdb,'set') ~= 0
+        
+    else
+        clear imdb;
+        [imdb] = generatepatches;
+    end
+    
+    opts.train = find(imdb.set==1);
+    
+    prepareGPUs(opts, epoch == start+1) ;
+    
+    % Train for one epoch.
+    params = opts;
+    params.epoch = epoch ;
+    params.learningRate = opts.learningRate(min(epoch, numel(opts.learningRate))) ;
+    params.thetaCurrent = opts.theta(min(epoch, numel(opts.theta)));
+    params.train = opts.train(randperm(numel(opts.train))) ; % shuffle
+    params.getBatch = getBatch ;
+    
+    if numel(opts.gpus) <= 1
+        [net,~] = processEpoch(net, state, params, 'train',imdb) ;
+        if ~evaluateMode
+            saveState(modelPath(epoch), net) ;
+        end
+        %  lastStats = state.stats ;
+    else
+        spmd
+            [net, ~] = processEpoch(net, state, params, 'train',imdb) ;
+            if labindex == 1 && ~evaluateMode
+                saveState(modelPath(epoch), net) ;
+            end
+            %  lastStats = state.stats ;
+        end
+        %lastStats = accumulateStats(lastStats) ;
+    end
+    
+    % stats.train(epoch) = lastStats.train ;
+    % stats.val(epoch) = lastStats.val ;
+    % clear lastStats ;
+    % saveStats(modelPath(epoch), stats) ;
+    
+end
+
+% With multiple GPUs, return one copy
+if isa(net, 'Composite'), net = net{1} ; end
+
+% -------------------------------------------------------------------------
+function [net, state] = processEpoch(net, state, params, mode, imdb)
+% -------------------------------------------------------------------------
+% Note that net is not strictly needed as an output argument as net
+% is a handle class. However, this fixes some aliasing issue in the
+% spmd caller.
+% initialize with momentum 0
+if isempty(state) || isempty(state.momentum)
+    state.momentum = num2cell(zeros(1, numel(net.params))) ;
+    state.m = num2cell(zeros(1, numel(net.params))) ;
+    state.v = num2cell(zeros(1, numel(net.params))) ;
+    state.t = num2cell(zeros(1, numel(net.params))) ;
+end
+
+% move CNN  to GPU as needed
+numGpus = numel(params.gpus) ;
+if numGpus >= 1
+    net.move('gpu') ;
+    state.momentum = cellfun(@gpuArray, state.momentum, 'uniformoutput', false) ;
+    state.m = cellfun(@gpuArray,state.m,'UniformOutput',false) ;
+    state.v = cellfun(@gpuArray,state.v,'UniformOutput',false) ;
+    state.t = cellfun(@gpuArray,state.t,'UniformOutput',false) ;
+end
+
+
+if numGpus > 1
+    parserv = ParameterServer(params.parameterServer) ;
+    net.setParameterServer(parserv) ;
+else
+    parserv = [] ;
+end
+
+% profile
+if params.profile
+    if numGpus <= 1
+        profile clear ;
+        profile on ;
+    else
+        mpiprofile reset ;
+        mpiprofile on ;
+    end
+end
+
+num = 0 ;
+epoch = params.epoch ;
+subset = params.(mode) ;
+%adjustTime = 0 ;
+
+stats.num = 0 ; % return something even if subset = []
+stats.time = 0 ;
+count = 0;
+%start = tic ;
+for t=1:params.batchSize:numel(subset)
+    %     fprintf('%s: epoch %02d: %3d/%3d:', mode, epoch, ...
+    %         fix((t-1)/params.batchSize)+1, ceil(numel(subset)/params.batchSize)) ;
+    
+    batchSize = min(params.batchSize, numel(subset) - t + 1) ;
+    count = count + 1;
+    for s=1:params.numSubBatches
+        % get this image batch and prefetch the next
+        batchStart = t + (labindex-1) + (s-1) * numlabs ;
+        batchEnd = min(t+params.batchSize-1, numel(subset)) ;
+        batch = subset(batchStart : params.numSubBatches * numlabs : batchEnd) ;
+        num = num + numel(batch) ;
+        if numel(batch) == 0, continue ; end
+        
+        inputs = params.getBatch(imdb, batch) ;
+        
+        if params.prefetch
+            if s == params.numSubBatches
+                batchStart = t + (labindex-1) + params.batchSize ;
+                batchEnd = min(t+2*params.batchSize-1, numel(subset)) ;
+            else
+                batchStart = batchStart + numlabs ;
+            end
+            nextBatch = subset(batchStart : params.numSubBatches * numlabs : batchEnd) ;
+            params.getBatch(imdb, nextBatch) ;
+        end
+        
+        if strcmp(mode, 'train')
+            net.mode = 'normal' ;
+            net.accumulateParamDers = (s ~= 1) ;
+            net.eval(inputs, params.derOutputs, 'holdOn', s < params.numSubBatches) ;
+        else
+            net.mode = 'test' ;
+            net.eval(inputs) ;
+        end
+    end
+    
+    % Accumulate gradient.
+    if strcmp(mode, 'train')
+        if ~isempty(parserv), parserv.sync() ; end
+        state = accumulateGradients(net, state, params, batchSize, parserv) ;
+    end
+    
+    
+    
+    %%%--------add your code here------------------------
+    
+    
+    %%%--------------------------------------------------
+    loss2 = squeeze(gather(net.vars(net.getVarIndex(params.derOutputs{1})).value));
+   
+    fprintf('%s: epoch %02d : %3d/%3d:', mode, epoch,  ...
+        fix((t-1)/params.batchSize)+1, ceil(numel(subset)/params.batchSize)) ;
+    fprintf('error: %f \n', loss2) ;
+    
+end
+
+% Save back to state.
+state.stats.(mode) = stats ;
+if params.profile
+    if numGpus <= 1
+        state.prof.(mode) = profile('info') ;
+        profile off ;
+    else
+        state.prof.(mode) = mpiprofile('info');
+        mpiprofile off ;
+    end
+end
+if ~params.saveMomentum
+    state.momentum = [] ;
+    state.m = [] ;
+    state.v = [] ;
+    state.t = [] ;
+else
+    state.momentum = cellfun(@gather, state.momentum, 'uniformoutput', false) ;
+    state.m = cellfun(@gather, state.m, 'uniformoutput', false) ;
+    state.v = cellfun(@gather, state.v, 'uniformoutput', false) ;
+    state.t = cellfun(@gather, state.t, 'uniformoutput', false) ;
+end
+
+net.reset() ;
+net.move('cpu') ;
+
+% -------------------------------------------------------------------------
+function state = accumulateGradients(net, state, params, batchSize, parserv)
+% -------------------------------------------------------------------------
+% numGpus = numel(params.gpus) ;
+% otherGpus = setdiff(1:numGpus, labindex) ;
+for p=1:numel(net.params)
+    
+    if ~isempty(parserv)
+        parDer = parserv.pullWithIndex(p) ;
+    else
+        parDer = net.params(p).der ;
+    end
+    
+    switch params.solver
+        
+        case 'SGD' %%% solver: SGD
+            
+            switch net.params(p).trainMethod
+                
+                case 'average' % mainly for batch normalization
+                    thisLR = net.params(p).learningRate;
+                    net.params(p).value = vl_taccum(...
+                        1 - thisLR, net.params(p).value, ...
+                        (thisLR/batchSize/net.params(p).fanout),  parDer) ;
+                    
+                otherwise
+                    thisDecay = params.weightDecay * net.params(p).weightDecay ;
+                    thisLR = params.learningRate * net.params(p).learningRate ;
+                    
+                    % Normalize gradient and incorporate weight decay.
+                    parDer = vl_taccum(1/batchSize, parDer, ...
+                        thisDecay, net.params(p).value) ;
+                    
+                    theta  = params.thetaCurrent/lr;
+                    parDer = gradientClipping(parDer,theta,params.gradientClipping);
+                    
+                    % Update momentum.
+                    state.momentum{p} = vl_taccum(...
+                        params.momentum, state.momentum{p}, ...
+                        -1, parDer) ;
+                    
+                    % Nesterov update (aka one step ahead).
+                    if params.nesterovUpdate
+                        delta = vl_taccum(...
+                            params.momentum, state.momentum{p}, ...
+                            -1, parDer) ;
+                    else
+                        delta = state.momentum{p} ;
+                    end
+                    
+                    % Update parameters.
+                    net.params(p).value = vl_taccum(...
+                        1,  net.params(p).value, thisLR, delta) ;
+            end
+            
+        case  'Adam'
+            
+            switch net.params(p).trainMethod
+                
+                case 'average' % mainly for batch normalization
+                    thisLR = net.params(p).learningRate;
+                    net.params(p).value = vl_taccum(...
+                        1 - thisLR, net.params(p).value, ...
+                        (thisLR/batchSize/net.params(p).fanout),  parDer) ;
+
+                otherwise
+                    
+                    thisLR = params.learningRate * net.params(p).learningRate ;
+                    state.t{p} = state.t{p} + 1;
+                    t = state.t{p};
+                    alpha = thisLR; % opts.alpha; 
+                    lr = alpha * sqrt(1 - params.beta2^t) / (1 - params.beta1^t);
+                    
+                    state.m{p} = state.m{p} + (1 - params.beta1) .* (net.params(p).der - state.m{p});
+                    state.v{p} = state.v{p} + (1 - params.beta2) .* (net.params(p).der .* net.params(p).der - state.v{p});
+                    net.params(p).value = net.params(p).value - lr * state.m{p} ./ (sqrt(state.v{p}) + params.epsilon);%   - thisLR * 0.0005 * net.params(p).value;
+            end
+    end
+end
+
+
+%%%-------------------------------------------------------------------------
+function A = smallClipping(A, theta)
+%%%-------------------------------------------------------------------------
+A(A>theta)  = A(A>theta) -0.0001;
+A(A<-theta) = A(A<-theta)+0.0001;
+%%%-------------------------------------------------------------------------
+function A = smallClipping2(A, theta1,theta2)
+%%%-------------------------------------------------------------------------
+A(A>theta1)  = A(A>theta1)-0.02;
+A(A<theta2)  = A(A<theta2)+0.02;
+
+function A = smallClipping3(A, theta1,theta2)
+%%%-------------------------------------------------------------------------
+A(A>theta1)  = A(A>theta1) -0.1;
+A(A<theta2)  = A(A<theta2) +0.1;
+% % -------------------------------------------------------------------------
+% function stats = accumulateStats(stats_)
+% % -------------------------------------------------------------------------
+%
+% for s = {'train', 'val'}
+%     s = char(s) ;
+%     total = 0 ;
+%
+%     % initialize stats stucture with same fields and same order as
+%     % stats_{1}
+%     stats__ = stats_{1} ;
+%     names = fieldnames(stats__.(s))' ;
+%     values = zeros(1, numel(names)) ;
+%     fields = cat(1, names, num2cell(values)) ;
+%     stats.(s) = struct(fields{:}) ;
+%
+%     for g = 1:numel(stats_)
+%         stats__ = stats_{g} ;
+%         num__ = stats__.(s).num ;
+%         total = total + num__ ;
+%
+%         for f = setdiff(fieldnames(stats__.(s))', 'num')
+%             f = char(f) ;
+%             stats.(s).(f) = stats.(s).(f) + stats__.(s).(f) * num__ ;
+%
+%             if g == numel(stats_)
+%                 stats.(s).(f) = stats.(s).(f) / total ;
+%             end
+%         end
+%     end
+%     stats.(s).num = total ;
+% end
+
+% -------------------------------------------------------------------------
+function stats = extractStats(stats, net)
+% -------------------------------------------------------------------------
+sel = find(cellfun(@(x) isa(x,'dagnn.Loss'), {net.layers.block})) ;
+for i = 1:numel(sel)
+    stats.(net.layers(sel(i)).outputs{1}) = net.layers(sel(i)).block.average ;
+end
+
+% -------------------------------------------------------------------------
+function saveState(fileName, net_)
+% -------------------------------------------------------------------------
+net = net_.saveobj() ;
+save(fileName, 'net') ;
+
+% -------------------------------------------------------------------------
+function saveStats(fileName, stats)
+% -------------------------------------------------------------------------
+if exist(fileName)
+    save(fileName, 'stats', '-append') ;
+else
+    save(fileName, 'stats') ;
+end
+
+% -------------------------------------------------------------------------
+function [net] = loadState(fileName)
+% -------------------------------------------------------------------------
+load(fileName, 'net') ;
+net = dagnn.DagNN.loadobj(net) ;
+% if isempty(whos('stats'))
+%     error('Epoch ''%s'' was only partially saved. Delete this file and try again.', ...
+%         fileName) ;
+% end
+
+%%%-------------------------------------------------------------------------
+function epoch = findLastCheckpoint(modelDir,modelName)
+%%%-------------------------------------------------------------------------
+list = dir(fullfile(modelDir, [modelName,'-epoch-*.mat'])) ;
+tokens = regexp({list.name}, [modelName,'-epoch-([\d]+).mat'], 'tokens') ;
+epoch = cellfun(@(x) sscanf(x{1}{1}, '%d'), tokens) ;
+epoch = max([epoch 0]) ;
+
+
+% -------------------------------------------------------------------------
+function clearMex()
+% -------------------------------------------------------------------------
+clear vl_tmove vl_imreadjpeg ;
+
+% -------------------------------------------------------------------------
+function prepareGPUs(opts, cold)
+% -------------------------------------------------------------------------
+numGpus = numel(opts.gpus) ;
+if numGpus > 1
+    % check parallel pool integrity as it could have timed out
+    pool = gcp('nocreate') ;
+    if ~isempty(pool) && pool.NumWorkers ~= numGpus
+        delete(pool) ;
+    end
+    pool = gcp('nocreate') ;
+    if isempty(pool)
+        parpool('local', numGpus) ;
+        cold = true ;
+    end
+    
+end
+if numGpus >= 1 && cold
+    fprintf('%s: resetting GPU\n', mfilename)
+    clearMex() ;
+    if numGpus == 1
+        gpuDevice(opts.gpus)
+    else
+        spmd
+            clearMex() ;
+            gpuDevice(opts.gpus(labindex))
+        end
+    end
+end
+
+
+%%%-------------------------------------------------------------------------
+function A = gradientClipping(A, theta,gradientClip)
+%%%-------------------------------------------------------------------------
+if gradientClip
+    A(A>theta)  = theta;
+    A(A<-theta) = -theta;
+else
+    return;
+end
+
+% -------------------------------------------------------------------------
+function fn = getBatch()
+% -------------------------------------------------------------------------
+fn = @(x,y) getDagNNBatch(x,y) ;
+
+% -------------------------------------------------------------------------
+function [inputs2] = getDagNNBatch(imdb, batch)
+% -------------------------------------------------------------------------
+noiselevel = 25;
+label      = imdb.labels(:,:,:,batch);
+label      = data_augmentation(label,randi(8));
+input      = label + noiselevel/255*randn(size(label),'single');  % add AWGN with noise level noiselevel
+input      = gpuArray(input);
+label      = gpuArray(label);
+inputs2    = {'input', input, 'label', label} ;
+
+
diff --git a/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/generatepatches.m b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/generatepatches.m
new file mode 100644
index 00000000..6d3e58b4
--- /dev/null
+++ b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/generatepatches.m
@@ -0,0 +1,100 @@
+function [imdb] = generatepatches
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% @article{zhang2017beyond,
+%   title={Beyond a {Gaussian} denoiser: Residual learning of deep {CNN} for image denoising},
+%   author={Zhang, Kai and Zuo, Wangmeng and Chen, Yunjin and Meng, Deyu and Zhang, Lei},
+%   journal={IEEE Transactions on Image Processing},
+%   year={2017},
+%   volume={26}, 
+%   number={7}, 
+%   pages={3142-3155}, 
+% }
+
+% by Kai Zhang (1/2018)
+% cskaizhang@gmail.com
+% https://github.com/cszn
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+addpath('utilities');
+batchSize     = 128;        % batch size
+folder        = 'Train400';  %
+nchannel      = 1;           % number of channels
+patchsize     = 40;
+
+stride        = 9;
+
+step1         = randi(stride)-1;
+step2         = randi(stride)-1;
+count         = 0;
+ext           =  {'*.jpg','*.png','*.bmp','*.jpeg'};
+filepaths     =  [];
+
+for i = 1 : length(ext)
+    filepaths = cat(1,filepaths, dir(fullfile(folder, ext{i})));
+end
+
+% count the number of extracted patches
+scales  = [1 0.9 0.8 0.7]; % scale the image to augment the training data
+
+for i = 1 : length(filepaths)
+    
+    image = imread(fullfile(folder,filepaths(i).name)); % uint8
+    if size(image,3)==3
+        image = rgb2gray(image);
+    end
+    %[~, name, exte] = fileparts(filepaths(i).name);
+    if mod(i,100)==0
+        disp([i,length(filepaths)]);
+    end
+    for s = 1:4
+        image = imresize(image,scales(s),'bicubic');
+        [hei,wid,~] = size(image);
+        for x = 1+step1 : stride : (hei-patchsize+1)
+            for y = 1+step2 :stride : (wid-patchsize+1)
+                count = count+1;
+            end
+        end
+    end
+end
+
+numPatches  = ceil(count/batchSize)*batchSize;
+diffPatches = numPatches - count;
+disp([int2str(numPatches),' = ',int2str(numPatches/batchSize),' X ', int2str(batchSize)]);
+
+
+count = 0;
+imdb.labels  = zeros(patchsize, patchsize, nchannel, numPatches,'single');
+
+for i = 1 : length(filepaths)
+    
+    image = imread(fullfile(folder,filepaths(i).name)); % uint8
+    %[~, name, exte] = fileparts(filepaths(i).name);
+    if size(image,3)==3
+        image = rgb2gray(image);
+    end
+    if mod(i,100)==0
+        disp([i,length(filepaths)]);
+    end
+    for s = 1:4
+        image = imresize(image,scales(s),'bicubic');
+        for j = 1:1
+            image_aug   = data_augmentation(image, j);  % augment data
+            im_label    = im2single(image_aug);         % single
+            [hei,wid,~] = size(im_label);
+            
+            for x = 1+step1 : stride : (hei-patchsize+1)
+                for y = 1+step2 :stride : (wid-patchsize+1)
+                    count       = count+1;
+                    imdb.labels(:, :, :, count)   = im_label(x : x+patchsize-1, y : y+patchsize-1,:);
+                    if count<=diffPatches
+                        imdb.labels(:, :, :, end-count+1)   = im_label(x : x+patchsize-1, y : y+patchsize-1,:);
+                    end
+                end
+            end
+        end
+    end
+end
+
+imdb.set    = uint8(ones(1,size(imdb.labels,4)));
+
diff --git a/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/readme.png b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/readme.png
new file mode 100644
index 00000000..ad0fb8b3
Binary files /dev/null and b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/readme.png differ
diff --git a/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/01.png b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/01.png
new file mode 100644
index 00000000..795f9872
Binary files /dev/null and b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/01.png differ
diff --git a/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/02.png b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/02.png
new file mode 100644
index 00000000..ecc66936
Binary files /dev/null and b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/02.png differ
diff --git a/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/03.png b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/03.png
new file mode 100644
index 00000000..bc20886a
Binary files /dev/null and b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/03.png differ
diff --git a/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/04.png b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/04.png
new file mode 100644
index 00000000..2a80fdb3
Binary files /dev/null and b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/04.png differ
diff --git a/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/05.png b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/05.png
new file mode 100644
index 00000000..ece8444d
Binary files /dev/null and b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/05.png differ
diff --git a/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/06.png b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/06.png
new file mode 100644
index 00000000..72a8a6be
Binary files /dev/null and b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/06.png differ
diff --git a/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/07.png b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/07.png
new file mode 100644
index 00000000..ea9f17d9
Binary files /dev/null and b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/07.png differ
diff --git a/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/08.png b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/08.png
new file mode 100644
index 00000000..df777bd8
Binary files /dev/null and b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/08.png differ
diff --git a/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/09.png b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/09.png
new file mode 100644
index 00000000..3b6a71bc
Binary files /dev/null and b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/09.png differ
diff --git a/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/10.png b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/10.png
new file mode 100644
index 00000000..931cfa1f
Binary files /dev/null and b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/10.png differ
diff --git a/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/11.png b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/11.png
new file mode 100644
index 00000000..d65b82ce
Binary files /dev/null and b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/11.png differ
diff --git a/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/12.png b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/12.png
new file mode 100644
index 00000000..b9019ad3
Binary files /dev/null and b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/testsets/Set12/12.png differ
diff --git a/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/utilities/Cal_PSNRSSIM.m b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/utilities/Cal_PSNRSSIM.m
new file mode 100644
index 00000000..7307e739
--- /dev/null
+++ b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/utilities/Cal_PSNRSSIM.m
@@ -0,0 +1,221 @@
+function [psnr_cur, ssim_cur] = Cal_PSNRSSIM(A,B,row,col)
+
+
+[n,m,ch]=size(B);
+A = A(row+1:n-row,col+1:m-col,:);
+B = B(row+1:n-row,col+1:m-col,:);
+A=double(A); % Ground-truth
+B=double(B); %
+
+e=A(:)-B(:);
+mse=mean(e.^2);
+psnr_cur=10*log10(255^2/mse);
+
+if ch==1
+    [ssim_cur, ~] = ssim_index(A, B);
+else
+    ssim_cur = (ssim_index(A(:,:,1), B(:,:,1)) + ssim_index(A(:,:,2), B(:,:,2)) + ssim_index(A(:,:,3), B(:,:,3)))/3;
+end
+
+
+function [mssim, ssim_map] = ssim_index(img1, img2, K, window, L)
+
+%========================================================================
+%SSIM Index, Version 1.0
+%Copyright(c) 2003 Zhou Wang
+%All Rights Reserved.
+%
+%The author is with Howard Hughes Medical Institute, and Laboratory
+%for Computational Vision at Center for Neural Science and Courant
+%Institute of Mathematical Sciences, New York University.
+%
+%----------------------------------------------------------------------
+%Permission to use, copy, or modify this software and its documentation
+%for educational and research purposes only and without fee is hereby
+%granted, provided that this copyright notice and the original authors'
+%names appear on all copies and supporting documentation. This program
+%shall not be used, rewritten, or adapted as the basis of a commercial
+%software or hardware product without first obtaining permission of the
+%authors. The authors make no representations about the suitability of
+%this software for any purpose. It is provided "as is" without express
+%or implied warranty.
+%----------------------------------------------------------------------
+%
+%This is an implementation of the algorithm for calculating the
+%Structural SIMilarity (SSIM) index between two images. Please refer
+%to the following paper:
+%
+%Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, "Image
+%quality assessment: From error measurement to structural similarity"
+%IEEE Transactios on Image Processing, vol. 13, no. 1, Jan. 2004.
+%
+%Kindly report any suggestions or corrections to zhouwang@ieee.org
+%
+%----------------------------------------------------------------------
+%
+%Input : (1) img1: the first image being compared
+%        (2) img2: the second image being compared
+%        (3) K: constants in the SSIM index formula (see the above
+%            reference). defualt value: K = [0.01 0.03]
+%        (4) window: local window for statistics (see the above
+%            reference). default widnow is Gaussian given by
+%            window = fspecial('gaussian', 11, 1.5);
+%        (5) L: dynamic range of the images. default: L = 255
+%
+%Output: (1) mssim: the mean SSIM index value between 2 images.
+%            If one of the images being compared is regarded as
+%            perfect quality, then mssim can be considered as the
+%            quality measure of the other image.
+%            If img1 = img2, then mssim = 1.
+%        (2) ssim_map: the SSIM index map of the test image. The map
+%            has a smaller size than the input images. The actual size:
+%            size(img1) - size(window) + 1.
+%
+%Default Usage:
+%   Given 2 test images img1 and img2, whose dynamic range is 0-255
+%
+%   [mssim ssim_map] = ssim_index(img1, img2);
+%
+%Advanced Usage:
+%   User defined parameters. For example
+%
+%   K = [0.05 0.05];
+%   window = ones(8);
+%   L = 100;
+%   [mssim ssim_map] = ssim_index(img1, img2, K, window, L);
+%
+%See the results:
+%
+%   mssim                        %Gives the mssim value
+%   imshow(max(0, ssim_map).^4)  %Shows the SSIM index map
+%
+%========================================================================
+
+
+if (nargin < 2 || nargin > 5)
+    ssim_index = -Inf;
+    ssim_map = -Inf;
+    return;
+end
+
+if (size(img1) ~= size(img2))
+    ssim_index = -Inf;
+    ssim_map = -Inf;
+    return;
+end
+
+[M N] = size(img1);
+
+if (nargin == 2)
+    if ((M < 11) || (N < 11))
+        ssim_index = -Inf;
+        ssim_map = -Inf;
+        return
+    end
+    window = fspecial('gaussian', 11, 1.5);	%
+    K(1) = 0.01;								      % default settings
+    K(2) = 0.03;								      %
+    L = 255;                                  %
+end
+
+if (nargin == 3)
+    if ((M < 11) || (N < 11))
+        ssim_index = -Inf;
+        ssim_map = -Inf;
+        return
+    end
+    window = fspecial('gaussian', 11, 1.5);
+    L = 255;
+    if (length(K) == 2)
+        if (K(1) < 0 || K(2) < 0)
+            ssim_index = -Inf;
+            ssim_map = -Inf;
+            return;
+        end
+    else
+        ssim_index = -Inf;
+        ssim_map = -Inf;
+        return;
+    end
+end
+
+if (nargin == 4)
+    [H W] = size(window);
+    if ((H*W) < 4 || (H > M) || (W > N))
+        ssim_index = -Inf;
+        ssim_map = -Inf;
+        return
+    end
+    L = 255;
+    if (length(K) == 2)
+        if (K(1) < 0 || K(2) < 0)
+            ssim_index = -Inf;
+            ssim_map = -Inf;
+            return;
+        end
+    else
+        ssim_index = -Inf;
+        ssim_map = -Inf;
+        return;
+    end
+end
+
+if (nargin == 5)
+    [H W] = size(window);
+    if ((H*W) < 4 || (H > M) || (W > N))
+        ssim_index = -Inf;
+        ssim_map = -Inf;
+        return
+    end
+    if (length(K) == 2)
+        if (K(1) < 0 || K(2) < 0)
+            ssim_index = -Inf;
+            ssim_map = -Inf;
+            return;
+        end
+    else
+        ssim_index = -Inf;
+        ssim_map = -Inf;
+        return;
+    end
+end
+
+C1 = (K(1)*L)^2;
+C2 = (K(2)*L)^2;
+window = window/sum(sum(window));
+img1 = double(img1);
+img2 = double(img2);
+
+mu1   = filter2(window, img1, 'valid');
+mu2   = filter2(window, img2, 'valid');
+mu1_sq = mu1.*mu1;
+mu2_sq = mu2.*mu2;
+mu1_mu2 = mu1.*mu2;
+sigma1_sq = filter2(window, img1.*img1, 'valid') - mu1_sq;
+sigma2_sq = filter2(window, img2.*img2, 'valid') - mu2_sq;
+sigma12 = filter2(window, img1.*img2, 'valid') - mu1_mu2;
+
+if (C1 > 0 & C2 > 0)
+    ssim_map = ((2*mu1_mu2 + C1).*(2*sigma12 + C2))./((mu1_sq + mu2_sq + C1).*(sigma1_sq + sigma2_sq + C2));
+else
+    numerator1 = 2*mu1_mu2 + C1;
+    numerator2 = 2*sigma12 + C2;
+    denominator1 = mu1_sq + mu2_sq + C1;
+    denominator2 = sigma1_sq + sigma2_sq + C2;
+    ssim_map = ones(size(mu1));
+    index = (denominator1.*denominator2 > 0);
+    ssim_map(index) = (numerator1(index).*numerator2(index))./(denominator1(index).*denominator2(index));
+    index = (denominator1 ~= 0) & (denominator2 == 0);
+    ssim_map(index) = numerator1(index)./denominator1(index);
+end
+
+mssim = mean2(ssim_map);
+
+return
+
+
+
+
+
+
+
diff --git a/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/utilities/Demo_DagNN_Merge_Bnorm.m b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/utilities/Demo_DagNN_Merge_Bnorm.m
new file mode 100644
index 00000000..c4c33046
--- /dev/null
+++ b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/utilities/Demo_DagNN_Merge_Bnorm.m
@@ -0,0 +1,200 @@
+function [] = Demo_DagNN_Merge_Bnorm()
+
+% merge bnorm:  'DnCNN-epoch-50.mat' ------> 'DnCNN-epoch-0.mat'
+
+inputfileName  = 'DnCNN-epoch-50.mat';
+targetfileName = 'DnCNN-epoch-0.mat';
+
+% Merge Bnorm to (1) accelerate the testing inference; and (2) fine-tune the model with small learning rate for better PSNR.
+
+load(inputfileName);
+net = dagnn.DagNN.loadobj(net) ;
+
+%CNN_IMAGENET_DEPLOY  Deploy a CNN
+
+isDag = isa(net, 'dagnn.DagNN');
+
+% if isDag
+%   dagRemoveLayersOfType(net, 'dagnn.Loss') ;
+%   dagRemoveLayersOfType(net, 'dagnn.DropOut') ;
+% else
+%   net = simpleRemoveLayersOfType(net, 'softmaxloss') ;
+%   net = simpleRemoveLayersOfType(net, 'dropout') ;
+% end
+
+if isDag
+  dagMergeBatchNorm(net) ;
+  dagRemoveLayersOfType(net, 'dagnn.BatchNorm') ;
+else
+  net = simpleMergeBatchNorm(net) ;
+  net = simpleRemoveLayersOfType(net, 'bnorm') ;
+end
+
+net = net.saveobj() ;
+save(targetfileName, 'net') ;
+
+
+
+
+% Switch to use MatConvNet default memory limit for CuDNN (512 MB)
+% if ~isDag
+%   for l = simpleFindLayersOfType(net, 'conv')
+%     net.layers{l}.opts = removeCuDNNMemoryLimit(net.layers{l}.opts) ;
+%   end
+% else
+%   for name = dagFindLayersOfType(net, 'dagnn.Conv')
+%     l = net.getLayerIndex(char(name)) ;
+%     net.layers(l).block.opts = removeCuDNNMemoryLimit(net.layers(l).block.opts) ;
+%   end
+% end
+
+% -------------------------------------------------------------------------
+function opts = removeCuDNNMemoryLimit(opts)
+% -------------------------------------------------------------------------
+remove = false(1, numel(opts)) ;
+for i = 1:numel(opts)
+  if isstr(opts{i}) && strcmp(lower(opts{i}), 'CudnnWorkspaceLimit')
+    remove([i i+1]) = true ;
+  end
+end
+opts = opts(~remove) ;
+
+
+
+% -------------------------------------------------------------------------
+function net = simpleRemoveMomentum(net)
+% -------------------------------------------------------------------------
+for l = 1:numel(net.layers)
+  if isfield(net.layers{l}, 'momentum')
+    net.layers{l} = rmfield(net.layers{l}, 'momentum') ;
+  end
+end
+
+% -------------------------------------------------------------------------
+function layers = simpleFindLayersOfType(net, type)
+% -------------------------------------------------------------------------
+layers = find(cellfun(@(x)strcmp(x.type, type), net.layers)) ;
+
+% -------------------------------------------------------------------------
+function net = simpleRemoveLayersOfType(net, type)
+% -------------------------------------------------------------------------
+layers = simpleFindLayersOfType(net, type) ;
+net.layers(layers) = [] ;
+
+% -------------------------------------------------------------------------
+function layers = dagFindLayersWithOutput(net, outVarName)
+% -------------------------------------------------------------------------
+layers = {} ;
+for l = 1:numel(net.layers)
+  if any(strcmp(net.layers(l).outputs, outVarName))
+    layers{1,end+1} = net.layers(l).name ;
+  end
+end
+
+% -------------------------------------------------------------------------
+function layers = dagFindLayersOfType(net, type)
+% -------------------------------------------------------------------------
+layers = [] ;
+for l = 1:numel(net.layers)
+  if isa(net.layers(l).block, type)
+    layers{1,end+1} = net.layers(l).name ;
+  end
+end
+
+% -------------------------------------------------------------------------
+function dagRemoveLayersOfType(net, type)
+% -------------------------------------------------------------------------
+names = dagFindLayersOfType(net, type) ;
+for i = 1:numel(names)
+  layer = net.layers(net.getLayerIndex(names{i})) ;
+  net.removeLayer(names{i}) ;
+  net.renameVar(layer.outputs{1}, layer.inputs{1}, 'quiet', true) ;
+end
+
+% -------------------------------------------------------------------------
+function dagMergeBatchNorm(net)
+% -------------------------------------------------------------------------
+names = dagFindLayersOfType(net, 'dagnn.BatchNorm') ;
+for name = names
+  name = char(name) ;
+  layer = net.layers(net.getLayerIndex(name)) ;
+
+  % merge into previous conv layer
+  playerName = dagFindLayersWithOutput(net, layer.inputs{1}) ;
+  playerName = playerName{1} ;
+  playerIndex = net.getLayerIndex(playerName) ;
+  player = net.layers(playerIndex) ;
+  if ~isa(player.block, 'dagnn.Conv')
+    error('Batch normalization cannot be merged as it is not preceded by a conv layer.') ;
+  end
+
+  % if the convolution layer does not have a bias,
+  % recreate it to have one
+  if ~player.block.hasBias
+    block = player.block ;
+    block.hasBias = true ;
+    net.renameLayer(playerName, 'tmp') ;
+    net.addLayer(playerName, ...
+                 block, ...
+                 player.inputs, ...
+                 player.outputs, ...
+                 {player.params{1}, sprintf('%s_b',playerName)}) ;
+    net.removeLayer('tmp') ;
+    playerIndex = net.getLayerIndex(playerName) ;
+    player = net.layers(playerIndex) ;
+    biases = net.getParamIndex(player.params{2}) ;
+    net.params(biases).value = zeros(block.size(4), 1, 'single') ;
+  end
+
+  filters = net.getParamIndex(player.params{1}) ;
+  biases = net.getParamIndex(player.params{2}) ;
+  multipliers = net.getParamIndex(layer.params{1}) ;
+  offsets = net.getParamIndex(layer.params{2}) ;
+  moments = net.getParamIndex(layer.params{3}) ;
+
+  [filtersValue, biasesValue] = mergeBatchNorm(...
+    net.params(filters).value, ...
+    net.params(biases).value, ...
+    net.params(multipliers).value, ...
+    net.params(offsets).value, ...
+    net.params(moments).value) ;
+
+  net.params(filters).value = filtersValue ;
+  net.params(biases).value = biasesValue ;
+  net.params(biases).learningRate = 1;
+end
+
+% -------------------------------------------------------------------------
+function net = simpleMergeBatchNorm(net)
+% -------------------------------------------------------------------------
+
+for l = 1:numel(net.layers)
+  if strcmp(net.layers{l}.type, 'bnorm')
+    if ~strcmp(net.layers{l-1}.type, 'conv')
+      error('Batch normalization cannot be merged as it is not preceded by a conv layer.') ;
+    end
+    [filters, biases] = mergeBatchNorm(...
+      net.layers{l-1}.weights{1}, ...
+      net.layers{l-1}.weights{2}, ...
+      net.layers{l}.weights{1}, ...
+      net.layers{l}.weights{2}, ...
+      net.layers{l}.weights{3}) ;
+    net.layers{l-1}.weights = {filters, biases} ;
+  end
+  f = net.getParamIndex(net.layers{l-1}.params) ;
+  net.params(f(2)).learningRate = 1;
+  
+  
+end
+
+% -------------------------------------------------------------------------
+function [filters, biases] = mergeBatchNorm(filters, biases, multipliers, offsets, moments)
+% -------------------------------------------------------------------------
+% wk / sqrt(sigmak^2 + eps)
+% bk - wk muk / sqrt(sigmak^2 + eps)
+a = multipliers(:) ./ moments(:,2) ;
+b = offsets(:) - moments(:,1) .* a ;
+biases(:) = biases(:) + b(:) ;
+sz = size(filters) ;
+numFilters = sz(4) ;
+filters = reshape(bsxfun(@times, reshape(filters, [], numFilters), a'), sz) ;
diff --git a/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/utilities/DnCNN-epoch-50.mat b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/utilities/DnCNN-epoch-50.mat
new file mode 100644
index 00000000..261c9cc7
Binary files /dev/null and b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/utilities/DnCNN-epoch-50.mat differ
diff --git a/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/utilities/data_augmentation.m b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/utilities/data_augmentation.m
new file mode 100644
index 00000000..349c61a8
--- /dev/null
+++ b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/utilities/data_augmentation.m
@@ -0,0 +1,60 @@
+function image = data_augmentation(image, mode)
+
+if mode == 1
+    return;
+end
+
+if mode == 2 % flipped
+    image = flipud(image);
+    return;
+end
+
+if mode == 3 % rotation 90
+    image = rot90(image,1);
+    return;
+end
+
+if mode == 4 % rotation 90 & flipped
+    image = rot90(image,1);
+    image = flipud(image);
+    return;
+end
+
+if mode == 5 % rotation 180
+    image = rot90(image,2);
+    return;
+end
+
+if mode == 6 % rotation 180 & flipped
+    image = rot90(image,2);
+    image = flipud(image);
+    return;
+end
+
+if mode == 7 % rotation 270
+    image = rot90(image,3);
+    return;
+end
+
+if mode == 8 % rotation 270 & flipped
+    image = rot90(image,3);
+    image = flipud(image);
+    return;
+end
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/utilities/modcrop.m b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/utilities/modcrop.m
new file mode 100644
index 00000000..bfc6a3fe
--- /dev/null
+++ b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/utilities/modcrop.m
@@ -0,0 +1,12 @@
+function imgs = modcrop(imgs, modulo)
+if size(imgs,3)==1
+    sz = size(imgs);
+    sz = sz - mod(sz, modulo);
+    imgs = imgs(1:sz(1), 1:sz(2));
+else
+    tmpsz = size(imgs);
+    sz = tmpsz(1:2);
+    sz = sz - mod(sz, modulo);
+    imgs = imgs(1:sz(1), 1:sz(2),:);
+end
+
diff --git a/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/utilities/shave.m b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/utilities/shave.m
new file mode 100644
index 00000000..9790fe39
--- /dev/null
+++ b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/utilities/shave.m
@@ -0,0 +1,3 @@
+function I = shave(I, border)
+I = I(1+border(1):end-border(1), ...
+      1+border(2):end-border(2), :, :);
diff --git a/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/utilities/simplenn_matlab.m b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/utilities/simplenn_matlab.m
new file mode 100644
index 00000000..7dafab8a
--- /dev/null
+++ b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/utilities/simplenn_matlab.m
@@ -0,0 +1,27 @@
+function res = simplenn_matlab(net, input)
+
+%% If you did not install the matconvnet package, you can use this for testing.
+
+n = numel(net.layers);
+
+res = struct('x', cell(1,n+1));
+res(1).x = input;
+
+for ilayer = 1 : n
+    l = net.layers{ilayer};
+    switch l.type
+        case 'conv'
+            for noutmaps = 1 : size(l.weights{1},4)
+                z = zeros(size(res(ilayer).x,1),size(res(ilayer).x,2),'single');
+                for ninmaps = 1 : size(res(ilayer).x,3)
+                    z = z + convn(res(ilayer).x(:,:,ninmaps), rot90(l.weights{1}(:,:,ninmaps,noutmaps),2),'same');
+                end
+                res(ilayer+1).x(:,:,noutmaps) = z + l.weights{2}(noutmaps);
+            end
+        case 'relu'
+            res(ilayer+1).x = max(res(ilayer).x,0);
+    end
+    res(ilayer).x = [];
+end
+
+end
diff --git a/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/vl_nnloss.m b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/vl_nnloss.m
new file mode 100644
index 00000000..3c42e954
--- /dev/null
+++ b/TrainingCodes/DnCNN_TrainingCodes_DagNN_v1.1/vl_nnloss.m
@@ -0,0 +1,12 @@
+function Y = vl_nnloss(X,c,dzdy,varargin)
+
+% --------------------------------------------------------------------
+% pixel-level L2 loss
+% --------------------------------------------------------------------
+if nargin <= 2 || isempty(dzdy)
+    t = ((X-c).^2)/2;
+    Y = sum(t(:))/size(X,4); % reconstruction error per sample;
+else
+    Y = bsxfun(@minus,X,c).*dzdy;
+end
+