!************************************************************************* !* * !* A-PWLSVR - Adaptive Piecewise Linear Support Vector Regression * !* * !* by Napsu Karmitsa 2019, Adil Bagirov and Kaisa Joki * !* (last modified 01.02.2022). * !* * !* The software is free for academic teaching and research * !* purposes but we ask you to refer the appropriate references * !* given below, if you use it. * !* * !************************************************************************* !* !* !* Codes included: !* !* !* spr.f03 - Mainprogram (this file). !* constants.f03 - Parameters and constants. !* initspr.f03 - Initialization of parameters for A-PWLSVR and DBDC. !* functions.f03 - Computation of DC components f_1 and f_2 and their subgradients for the PWLSVR problem. !* bundle1.f03 - Bundle of DC component f_1. !* bundle2.f03 - Bundle of DC component f_2. !* dbdc.f03 - DBDC method. !* plqdf1.f - Quadratic solver by L. Luksan. !* !* Makefile - makefile. !* !* !* Calling sequence: !* !* ./spr infile nrecord ntrain nft maxmin [noutcom maxlin eps] !* !* with the following arguments: !* !* infile - the data set; !* nrecord - number of records in the data set; !* ntrain - size of the training set, ntrain <= nrecord; !* nft - number of features in the data set; !* maxmin - maximum number of min functions under maximum in the model; !* noutcom - number of output column, optional, nft by default; !* maxlin - number of linear functions under each minimum in the model, optional, !* Default value = 3; !* eps - stopping tolerance, optional, 10^(-2) by default. !* !* If you want to tune other parameters, modify initspr.f03. !* As default the results are written in !* !* predictions.txt - Evaluation criteria (e.g. RMSE, MAE,...) with training and predictions; !* reg_results.txt - Regression results with hyperplanes. !* pwl_points.txt - Points of original data and pwl-model obtained !* for drawing purposis. Works only with nft = 2 and maxlin = 3 (otherwise empty). !* !* !* References: !* A. Bagirov, S. Taheri, N. Karmitsa, K. Joki, M.M. Mäkelä, "Adaptive piecewise linear support vector regression", submitted 2020. !* !* !* Acknowledgements: !* !* The work was financially supported by the Academy of Finland (Project No. 289500 and 319274). !* !* !***************************************************************************************** !* !* * PROGRAM spr * !* !* Main program for A-PWLSVR piecewise linear support vector regression. !* !***************************************************************************************** PROGRAM spr USE constants, ONLY : dp ! Precision for reals. USE initspr, ONLY : & ! Initialization of parameters. infile, & ! Name of the dataset file. outfile1, & ! Regression results with hyperplanes. outfile2, & ! Regression results with predictions. outfile3, & ! Data points and pwl-model (only usable for two dim data). maxv, & ! Maximum value of data points. minv, & ! Minimum value of data points. nrecord, & ! Number of instances in data, from user. ntrain, & ! Size of the training set, ntrain <= nrecord, from user. nft, & ! Number of features in data, from user. noutcom, & ! Output column. iscale, & ! Scaling is ! 0 - not used, ! 1 - used. inc, & ! Incremental algorithm is ! 0 - not used, ! 1 - used. ipredi, & ! Prediction is computed ! 0 - only at the termination with complite model, ! 1 - at every step. maxmin, & ! Maximum number of min functions under maximum in the model. maxmin2, & ! maxmin2 = 1,...,maxmin. maxlin, & ! Number of linear functions under each minimum in the model. mk, & ! Vector of the number of linear functions under each ! minimum functions (for the future, now always the same number). maxdim, & ! Maximum dimension of the problem. tol, & ! Stopping tolerance. eps, & ! Epsilon used in SVM. a, & ! Data matrix a(nrecord,nft), from input file. a4, & ! Auxiliary data matrix a4(ntrain,nft). x_0=>x0, & ! The starting point. x0reg, & ! The SVR solution. x_solution, & ! The solution. x_solution2, & ! The solution. x_unscaled, & ! The unscaled solution. f_solution, & ! The objective function value at the solution 'x_solution' f_solution2, & ! The objective function value at the solution 'x_solution2' pred, & ! the values of predicted outputs dlin, & clin, & init_par, & scaling, & rescaling USE initdbdc, ONLY : & ! Initialization of parameters. problem1, & problem2, & agg_used, & stepsize_used, & termination, & counter, & iprint, & mit, & mrounds, & mrounds_clarke USE bundle1 ! The BUNDLE of the DC component f_1 USE bundle2 ! The BUNDLE of the DC component f_2 USE dbdc_non ! DBDC method IMPLICIT NONE !local REAL(KIND=dp) :: & h,h2,tmp,tmp2, & atmp, & rmse1, & mae1, & mse1, & mape1, & mape2, & ce1, & r1, & e0_test,e0_tr,e1_test,e1_tr, & f_1,f_2, & sx,sy,rcor,ytest, & large,piece,piece_best,piece_min, & start_time, finish_time ! start and finish CPU time INTEGER :: & count, & ! Number of commandline arguments ntest, & ! Size of the test set. neval, & ! Total number of linear function evaluations. i,i2,i3,j,j1,j2,j3, & ipiecemin,ipiecebest, & pind, & mtmp CHARACTER(len=50) :: arg ! Intrinsic Functions INTRINSIC ABS,MAX count = command_argument_count() IF (count < 1) THEN PRINT* PRINT*,'A-PWLSVR - Adaptive Piecewise Linear Support Vector Regression Method' PRINT* PRINT*,'Calling sequence:' PRINT* PRINT*,' ./spr infile nrecord ntrain nft maxmin [noutcom maxlin tol]' PRINT* PRINT*,'with the following arguments:' PRINT* PRINT*,' infile - the data set;' PRINT*,' nrecord - number of records in the data set;' PRINT*,' ntrain - number of records in the training set;' PRINT*,' nft - number of features in the data set;' PRINT*,' maxmin - maximum number of min functions under maximum in the model;' PRINT*,' noutcom - optional, number of output column, nft by default;' PRINT*,' maxlin - optional, number of linear functions under each minimum in the model,' PRINT*,' 2 by default;' PRINT*,' tol - optional, stopping tolerance, 10^(-2) by default.' PRINT* PRINT*,'If you want to tune other parameters, modify initspr.f03.' PRINT*,'As default the results are written in' PRINT* PRINT*,' predictions.txt - Evaluation criteria (e.g. RMSE, MAE,...) with training and predictions;' PRINT*,' reg_results.txt - Regression results with hyperplanes.' PRINT*,' pwl_points.txt - Points of original data and pwl-model obtained for drawing purposis.' PRINT*,' Works only with nft = 2 and maxlin = 3 (otherwise empty).' PRINT* STOP END IF IF (count < 5) STOP ' Too few commandline arguments. Type ./spr for help.' CALL get_command_argument(2, arg) READ (arg,*) nrecord CALL get_command_argument(3, arg) READ (arg,*) ntrain CALL get_command_argument(4, arg) READ (arg,*) nft CALL get_command_argument(5, arg) READ (arg,*) maxmin IF (count >= 6) THEN CALL get_command_argument(6, arg) READ (arg,*) noutcom ELSE noutcom=nft END IF IF (count >= 7) THEN CALL get_command_argument(7, arg) READ (arg,*) maxlin END IF IF (count >= 8) THEN CALL get_command_argument(8, arg) READ (arg,*) tol END IF CALL get_command_argument(1, infile) PRINT*,' Input file: ',infile PRINT*,' Number of records = ',nrecord PRINT*,' Size of the training set = ',ntrain PRINT*,' Number of features = ',nft PRINT*,' Maximum number of min functions = ',maxmin PRINT*,' Maximum number of linear functions = ',maxlin CALL init_par() ! sets dimension for optimization problem etc. allocate(a(nrecord,nft),a4(ntrain,nft)) allocate(x_0(maxdim),x0reg(maxdim),x_solution(maxdim),x_solution2(maxdim),x_unscaled(maxdim),clin(nft),dlin(nft)) ! maxdim= nft*q, where q=sum_k=1,...K (M_k) allocate(pred(MAX(nrecord-ntrain,ntrain))) allocate(mk(maxmin)) neval=0 large=1.0E+36_dp mk=maxlin ! each min function has maxlin linear functions OPEN(43,file=outfile3) OPEN(42,file=outfile1) OPEN(41,file=outfile2) OPEN(78,file=infile,status='old',form='formatted') DO i=1,nrecord READ(78,*) (a(i,j),j=1,nft) END DO SELECT CASE (inc) CASE(0) WRITE(42,*) 'Piecewise linear regression without incremental approach.' WRITE(42,*) 'Data set: ',infile WRITE(42,100) maxmin,maxlin,eps IF(iscale == 1) WRITE(42,*) 'Scaling is used. ' IF(iscale == 0) WRITE(42,*) 'No scaling. ' WRITE(41,*) 'Piecewise linear regression without incremental approach.' WRITE(41,*) 'Data set: ',infile WRITE(41,100) maxmin,maxlin,eps IF(iscale == 1) WRITE(41,*) 'Scaling is used. ' IF(iscale == 0) WRITE(41,*) 'No scaling. ' CASE(1) WRITE(42,*) 'A-PWLSVR - Adaptive Piecewise linear support vector regression.' WRITE(42,*) 'Data set: ',infile WRITE(42,100) maxmin,maxlin,eps IF(iscale == 1) WRITE(42,*) 'Scaling is used. ' IF(iscale == 0) WRITE(42,*) 'No scaling. ' WRITE(41, *) 'A-PWLSVR - Adaptive Piecewise linear support vector regression.' WRITE(41, *) 'Data set: ',infile WRITE(41,100) maxmin,maxlin,eps IF(iscale == 1) WRITE(41,*) 'Scaling is used. ' IF(iscale == 0) WRITE(41,*) 'No scaling. ' CASE(2) WRITE(42,*) 'A-PWLSVR - Adaptive Piecewise linear support vector regression with double starting points.' WRITE(42,*) 'Data set: ',infile WRITE(42,100) maxmin,maxlin,eps IF(iscale == 1) WRITE(42,*) 'Scaling is used. ' IF(iscale == 0) WRITE(42,*) 'No scaling. ' WRITE(41, *) 'A-PWLSVR - Adaptive Piecewise linear support vector regression with double starting points.' WRITE(41, *) 'Data set: ',infile WRITE(41,100) maxmin,maxlin,eps IF(iscale == 1) WRITE(41,*) 'Scaling is used. ' IF(iscale == 0) WRITE(41,*) 'No scaling. ' CASE DEFAULT END SELECT 100 FORMAT(' Model with maxmin = ',i3,', maxlin = ', i3,', and epsilon = ',f5.3) IF (noutcom < nft) THEN ! swap features DO i = 1,nrecord atmp = a(i,noutcom) a(i,noutcom) = a(i,nft) a(i,nft) = atmp END DO END IF a4=a CALL cpu_time(start_time) IF (iscale == 1) CALL scaling WRITE(42,*) WRITE(42,*) 'Number of records: ',nrecord WRITE(42,*) 'Size of the training set: ',ntrain WRITE(42,*) 'Number of features: ',nft WRITE(41,*) WRITE(41,*) 'Number of records: ',nrecord WRITE(41,*) 'Size of the training set: ',ntrain WRITE(41,*) 'Number of features: ',nft x_0 = 2_dp ! the very first x_0 q=1 ! for linear regression maxdim=nft ! for linear regression maxmin2 = 1 ! for linear regression mk=1 ! for linear regression counter=0 ! First compute just linear regression ! That is, we have only one linear function under min and only one min under max. ! The number of variables in the optimization problem is nft. CALL DBDC_algorithm( x_0(1:maxdim), x_solution(1:maxdim), f_solution, mit, mrounds, & & mrounds_clarke, termination, counter, agg_used, & & stepsize_used, iprint ) neval=neval+counter(3)+counter(7) WRITE(42,*) WRITE(42,*) 'Fit function value for linear regression: ',f_solution WRITE(42,*) 'Number of linear function evaluations: ',neval IF (iscale == 1) THEN CALL rescaling !(x_solution,x_unscaled) ELSE x_unscaled=x_solution END IF WRITE(42,*) ! Print the solution WRITE(42,*) 'x: ',x_unscaled(1:maxdim) ! x0reg=x_solution ! If scaling is used this is the scaled starting point. f_1=f_solution ! Indices for training data for linear regression IF (nrecord > ntrain .AND. ipredi == 1) THEN ! there are training and test sets e0_tr=0_dp ytest=0_dp DO i = 1, ntrain e0_tr=e0_tr+a(i,nft) pred(i)=x_unscaled(nft) DO j=1, nft-1 pred(i) = pred(i)+x_unscaled(j)*a(i,j) END DO ytest=ytest+pred(i) END DO ytest=ytest/REAL(ntrain,dp) e0_tr=e0_tr/REAL(ntrain,dp) rmse1=0_dp mae1=0_dp e1_tr=0_dp mape1=0_dp DO i=1,ntrain rmse1=rmse1+(pred(i)-a(i,nft))**2 mae1=mae1+ABS(pred(i)-a(i,nft)) if(abs(a(i,nft)).le.0.01_dp) then mape1=mape1+abs(pred(i)-a(i,nft)) else mape1=mape1+abs((pred(i)-a(i,nft))/a(i,nft)) end if e1_tr=e1_tr+(a(i,nft)-e0_tr)**2 !e0? END DO mse1=rmse1/REAL(ntrain-nft,dp) ce1=1_dp-(rmse1/e1_tr) rmse1=SQRT(rmse1/REAL(ntrain,dp)) mae1=mae1/REAL(ntrain,dp) mape1=mape1/REAL(ntrain,dp) IF(ntrain > 1) THEN sx=0_dp sy=0_dp rcor=0_dp DO i=1,ntrain sx=sx+(pred(i)-ytest)**2 sy=sy+(a(i,nft)-e0_tr)**2 rcor=rcor+(pred(i)-ytest)*(a(i,nft)-e0_tr) END DO sx=SQRT(sx/REAL(ntrain-1)) sy=SQRT(sy/REAL(ntrain-1)) r1=rcor/(REAL(ntrain-1,dp)*sx*sy) END IF WRITE(41,*) WRITE(41,*) 'Evaluation criteria for training set in linear regression (rmse, mae, mse, mape,ce, r): ' WRITE(41,*) rmse1!,mae1,mse1,ce1,r1 WRITE(41,*) mae1 WRITE(41,*) mse1 WRITE(41,*) mape1 WRITE(41,*) ce1 WRITE(41,*) r1 END IF ! Indices for test data for linear regression IF (nrecord > ntrain .AND. ipredi == 1) THEN ! prediction ntest=nrecord-ntrain e0_test=0_dp ytest=0_dp DO i = 1, ntest e0_test=e0_test+a(ntrain+i,nft) pred(i)=x_unscaled(nft) DO j=1, nft-1 pred(i) = pred(i)+x_unscaled(j)*a(ntrain+i,j) END DO ytest=ytest+pred(i) END DO ytest=ytest/REAL(ntest,dp) e0_test=e0_test/REAL(ntest,dp) rmse1=0_dp mae1=0_dp mape2=0_dp e1_test=0_dp DO i=1,ntest rmse1=rmse1+(pred(i)-a(ntrain+i,nft))**2 mae1=mae1+ABS(pred(i)-a(ntrain+i,nft)) if(abs(a(ntrain+i,nft)).le.0.01_dp) then mape2=mape2+abs(pred(i)-a(ntrain+i,nft)) else mape2=mape2+abs((pred(i)-a(ntrain+i,nft))/a(ntrain+i,nft)) end if e1_test=e1_test+(a(ntrain+i,nft)-e0_test)**2 END DO mse1=rmse1/REAL(ntest-nft,dp) ce1=1_dp-(rmse1/e1_test) rmse1=SQRT(rmse1/REAL(ntest,dp)) mae1=mae1/REAL(ntest,dp) mape2=mape2/REAL(ntest,dp) IF(ntest > 1) THEN sx=0_dp sy=0_dp rcor=0_dp DO i=1,ntest sx=sx+(pred(i)-ytest)**2 sy=sy+(a(ntrain+i,nft)-e0_test)**2 rcor=rcor+(pred(i)-ytest)*(a(ntrain+i,nft)-e0_test) END DO sx=SQRT(sx/REAL(ntest-1)) sy=SQRT(sy/REAL(ntest-1)) r1=rcor/(REAL(ntest-1,dp)*sx*sy) END IF WRITE(41,*) WRITE(41,*) 'Evaluation criteria for linear regression (rmse, mae, mse, mape, ce, r): ' WRITE(41,*) rmse1!,mae1,mse1,ce1,r1 WRITE(41,*) mae1 WRITE(41,*) mse1 WRITE(41,*) mape2 WRITE(41,*) ce1 WRITE(41,*) r1 END IF f_2=f_1 IF (inc >= 1) THEN ! Compute starting points and structure incremental. ! Computation of min minloop: DO j=2,maxlin PRINT*,j,'linear functions under min.' ! Set a new starting point x_0(1:maxdim)=x_solution(1:maxdim) ! maxdim first elements, note that maxdim is not yet updated. x_0(maxdim+1:j*nft)=x_solution(maxdim-nft+1:maxdim) ! the last nft elements mk=mk+1 maxdim=maxdim+nft q=q+1 CALL DBDC_algorithm( x_0(1:maxdim), x_solution(1:maxdim), f_solution, mit, mrounds, & & mrounds_clarke, termination, counter, agg_used, & & stepsize_used, iprint ) neval=neval + j*(counter(3)+counter(7)) WRITE(42,*) WRITE(42,101) 1,j WRITE(42,*) 'Fit function value for PWLR is ',f_solution WRITE(42,*) 'Number of linear function evaluations: ',neval WRITE(42,*) 'The difference between function values is ',(f_2-f_solution)/(f_1+1_dp) WRITE(42,*) 'The 2. diff. between function values is ',(f_2-f_solution)/(f_2+1_dp) WRITE(42,*) 101 FORMAT(' Model with maxmin = ',i3,' and maxlin = ', i3) IF (iscale == 1) THEN CALL rescaling !(x_solution,x_unscaled) ELSE x_unscaled=x_solution END IF ! Print the solution WRITE(42,*) 'x: ',x_unscaled(1:maxdim) ! Stop here if f is small enough IF (f_solution < tol) EXIT minloop ! Commented, causes problems with next loop arrays. ! If you add this you should adjust maxlin etc corresbondingly. !IF ((f_2-f_solution)/(f_1+1_dp) < tol) THEN !PRINT*,' Early stopping from the first loop. ' !EXIT minloop !END IF f_2=f_solution ! Indices for training data for PWL-regression IF (nrecord > ntrain .AND. ipredi == 1) THEN ! training and test sets are used WRITE(41,*) WRITE(41,101) 1,j mtmp=nft*j ytest=0_dp DO i = 1, ntrain ! computation of correct piece pind=1 ipiecebest=1 piece_best=-large ipiecemin=pind piece_min=large DO j2=1,j piece=x_unscaled(pind*nft) ! DO j3=1, nft-1 piece = piece+x_unscaled((pind-1)*nft+j3)*a(i,j3) END DO IF (piece piece_best) THEN ipiecebest=ipiecemin piece_best=piece_min END IF pred(i)=x_unscaled(ipiecebest*nft) DO j1=1, nft-1 pred(i) = pred(i)+x_unscaled((ipiecebest-1)*nft+j1)*a(i,j1) END DO ytest=ytest+pred(i) END DO ytest=ytest/REAL(ntrain,dp) rmse1=0_dp mae1=0_dp mape1=0_dp DO i=1,ntrain rmse1=rmse1+(pred(i)-a(i,nft))**2 mae1=mae1+ABS(pred(i)-a(i,nft)) if(abs(a(i,nft)).le.0.01_dp) then mape1=mape1+abs(pred(i)-a(i,nft)) else mape1=mape1+abs((pred(i)-a(i,nft))/a(i,nft)) end if END DO mse1=rmse1/REAL(ntrain-mtmp,dp) ce1=1_dp-(rmse1/e1_tr) rmse1=SQRT(rmse1/REAL(ntrain,dp)) mae1=mae1/REAL(ntrain,dp) mape1=mape1/REAL(ntrain,dp) IF(ntrain > 1) THEN sx=0_dp sy=0_dp rcor=0_dp DO i=1,ntrain sx=sx+(pred(i)-ytest)**2 sy=sy+(a(i,nft)-e0_tr)**2 rcor=rcor+(pred(i)-ytest)*(a(i,nft)-e0_tr) END DO sx=SQRT(sx/REAL(ntrain-1)) sy=SQRT(sy/REAL(ntrain-1)) r1=rcor/(REAL(ntrain-1,dp)*sx*sy) END IF WRITE(41,*) WRITE(41,*) 'Evaluation criteria for training set in SVM piecewise linear regression (rmse, mae, mse, mape, ce, r): ' WRITE(41,*) rmse1!,mae1,mse1,mape1,ce1,r1 WRITE(41,*) mae1 WRITE(41,*) mse1 WRITE(41,*) mape1 WRITE(41,*) ce1 WRITE(41,*) r1 END IF ! Indices for test data for PWL-regression !!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF (nrecord > ntrain .AND. ipredi == 1) THEN ! prediction WRITE(41,*) WRITE(41,101) 1,j mtmp=nft*j ytest=0_dp DO i = 1, ntest ! computation of correct piece pind=1 ipiecebest=1 piece_best=-large ipiecemin=pind piece_min=large DO j2=1,j piece=x_unscaled(pind*nft) ! DO j3=1, nft-1 piece = piece+x_unscaled((pind-1)*nft+j3)*a(ntrain+i,j3) END DO IF (piece piece_best) THEN ipiecebest=ipiecemin piece_best=piece_min END IF pred(i)=x_unscaled(ipiecebest*nft) DO j1=1, nft-1 pred(i) = pred(i)+x_unscaled((ipiecebest-1)*nft+j1)*a(ntrain+i,j1) END DO ytest=ytest+pred(i) END DO ytest=ytest/REAL(ntest,dp) rmse1=0_dp mae1=0_dp mape2=0_dp DO i=1,ntest rmse1=rmse1+(pred(i)-a(ntrain+i,nft))**2 mae1=mae1+ABS(pred(i)-a(ntrain+i,nft)) if(abs(a(ntrain+i,nft)).le.0.01_dp) then mape2=mape2+abs(pred(i)-a(ntrain+i,nft)) else mape2=mape2+abs((pred(i)-a(ntrain+i,nft))/a(ntrain+i,nft)) end if END DO mse1=rmse1/REAL(ntest-mtmp,dp) ce1=1_dp-(rmse1/e1_test) rmse1=SQRT(rmse1/REAL(ntest,dp)) mae1=mae1/REAL(ntest,dp) mape2=mape2/REAL(ntest,dp) IF(ntest > 1) THEN sx=0_dp sy=0_dp rcor=0_dp DO i=1,ntest sx=sx+(pred(i)-ytest)**2 sy=sy+(a(ntrain+i,nft)-e0_test)**2 rcor=rcor+(pred(i)-ytest)*(a(ntrain+i,nft)-e0_test) END DO sx=SQRT(sx/REAL(ntest-1)) sy=SQRT(sy/REAL(ntest-1)) r1=rcor/(REAL(ntest-1,dp)*sx*sy) END IF WRITE(41,*) WRITE(41,*) 'Evaluation criteria for SVM piecewise linear regression (rmse, mae, mse, mape, ce, r): ' WRITE(41,*) rmse1!,mae1,mse1,ce1,r1 WRITE(41,*) mae1 WRITE(41,*) mse1 WRITE(41,*) mape2 WRITE(41,*) ce1 WRITE(41,*) r1 END IF END DO minloop IF (f_solution > tol) THEN maxloop: DO j=2, maxmin PRINT*,j,'min under max' f_2=f_solution ! Set a new starting point x_0(1:maxdim)=x_solution(1:maxdim) ! maxdim first elements, note that maxdim is not yet updated. x_0(maxdim+1:j*maxlin*nft)=x_solution(maxdim-nft*maxlin+1:maxdim) ! the last nft elements q=maxlin*j maxdim=nft*q maxmin2=j CALL DBDC_algorithm( x_0(1:maxdim), x_solution(1:maxdim), f_solution, mit, mrounds, & & mrounds_clarke, termination, counter, agg_used, & & stepsize_used, iprint ) neval = neval+q*(counter(3)+counter(7)) IF (inc==2) THEN ! Here we need to compute more if inc==2 ! Set a new starting point DO i=1,j*maxlin x_0((i-1)*nft+1:i*nft)=x0reg(1:nft) END DO CALL DBDC_algorithm( x_0(1:maxdim), x_solution2(1:maxdim), f_solution2, mit, mrounds, & & mrounds_clarke, termination, counter, agg_used, & & stepsize_used, iprint ) neval = neval+q*(counter(3)+counter(7)) IF (f_solution2 < f_solution) THEN x_solution = x_solution2 f_solution = f_solution2 END IF END IF WRITE(42,*) WRITE(42,101) j,maxlin WRITE(42,*) 'Fit function value1 for PWLR is ',f_solution !WRITE(42,*) 'Fit function value2 for PWLR is ',f_solution2 WRITE(42,*) 'Number of linear function evaluations: ',neval,q WRITE(42,*) 'The difference between function values is ',(f_2-f_solution)/(f_1+1_dp) !WRITE(42,*) 'The difference between function values is ',(f_2-f_solution2)/(f_1+1_dp) WRITE(42,*) 'The 2. diff. between function values is ',(f_2-f_solution)/(f_2+1_dp) WRITE(42,*) !IF (inc==2) THEN ! This could be given above but then we cannot print. Fix to final version. ! IF (f_solution2 < f_solution) THEN ! x_solution = x_solution2 ! f_solution = f_solution2 ! END IF !END IF IF (iscale == 1) THEN CALL rescaling !(x_solution,x_unscaled) ELSE x_unscaled=x_solution END IF !result can be printed here WRITE(42,*) 'x: ',x_unscaled(1:maxdim) IF (nft == 2 .AND. maxlin == 3) THEN ! Note works only for two dimensional data and 3 linear functions under maximum WRITE(43,101) j,maxlin h=minv h2=(maxv-minv)/1000_dp DO i=1,1000 tmp2 = - large h=h+h2 i3=1 DO i2=1,j tmp = MIN(x_unscaled(i3)*h+x_unscaled(i3+1),x_unscaled(i3+2)*h+x_unscaled(i3+3)) tmp = MIN(tmp,x_unscaled(i3+4)*h+x_unscaled(i3+5)) i3 = i3 + 6 IF (tmp > tmp2) tmp2 = tmp END DO WRITE(43,*) h, tmp2 END DO END IF mtmp=j ! Stop here if no progress !IF (f_solution < tol) EXIT maxloop !IF ((f_2-f_solution)/(f_1+1_dp) < tol) EXIT maxloop ! Indices for training IF (nrecord > ntrain .AND. ipredi == 1) THEN ! prediction WRITE(41,*) WRITE(41,101) mtmp,maxlin mtmp=mtmp*nft*maxlin ytest=0_dp DO i = 1, ntrain ! computation of correct piece pind=1 ipiecebest=1 piece_best=-large DO j1=1,j ipiecemin=pind piece_min=large DO j2=1,maxlin piece=x_unscaled(pind*nft) ! DO j3=1, nft-1 piece = piece+x_unscaled((pind-1)*nft+j3)*a(i,j3) END DO IF (piece piece_best) THEN ipiecebest=ipiecemin piece_best=piece_min END IF END DO pred(i)=x_unscaled(ipiecebest*nft) DO j1=1, nft-1 pred(i) = pred(i)+x_unscaled((ipiecebest-1)*nft+j1)*a(i,j1) END DO ytest=ytest+pred(i) END DO ytest=ytest/REAL(ntrain,dp) rmse1=0_dp mae1=0_dp mape1=0_dp DO i=1,ntrain rmse1=rmse1+(pred(i)-a(i,nft))**2 mae1=mae1+ABS(pred(i)-a(i,nft)) if(abs(a(i,nft)).le.0.01_dp) then mape1=mape1+abs(pred(i)-a(i,nft)) else mape1=mape1+abs((pred(i)-a(i,nft))/a(i,nft)) end if END DO mse1=rmse1/REAL(ntrain-mtmp,dp) ce1=1_dp-(rmse1/e1_tr) ! rmse1=SQRT(rmse1/REAL(ntrain,dp)) mae1=mae1/REAL(ntrain,dp) mape1=mape1/REAL(ntrain,dp) IF(ntrain > 1) THEN sx=0_dp sy=0_dp rcor=0_dp DO i=1,ntrain sx=sx+(pred(i)-ytest)**2 sy=sy+(a(i,nft)-e0_tr)**2 rcor=rcor+(pred(i)-ytest)*(a(i,nft)-e0_tr) END DO sx=SQRT(sx/REAL(ntrain-1)) sy=SQRT(sy/REAL(ntrain-1)) r1=rcor/(REAL(ntrain-1,dp)*sx*sy) END IF WRITE(41,*) WRITE(41,*) 'Evaluation criteria for training in SVM piecewise linear regression (rmse, mae, mse, mape, ce, r): ' WRITE(41,*) rmse1!,mae1,mse1,ce1,r1 WRITE(41,*) mae1 WRITE(41,*) mse1 WRITE(41,*) mape1 WRITE(41,*) ce1 WRITE(41,*) r1 END IF ! Indices for test set. mtmp=j IF (nrecord > ntrain .AND. ipredi == 1) THEN ! prediction WRITE(41,*) WRITE(41,101) mtmp,maxlin mtmp=mtmp*nft*maxlin ytest=0_dp DO i = 1, ntest ! computation of correct piece pind=1 ipiecebest=1 piece_best=-large DO j1=1,j ipiecemin=pind piece_min=large DO j2=1,maxlin piece=x_unscaled(pind*nft) ! DO j3=1, nft-1 piece = piece+x_unscaled((pind-1)*nft+j3)*a(ntrain+i,j3) END DO IF (piece piece_best) THEN ipiecebest=ipiecemin piece_best=piece_min END IF END DO pred(i)=x_unscaled(ipiecebest*nft) DO j1=1, nft-1 pred(i) = pred(i)+x_unscaled((ipiecebest-1)*nft+j1)*a(ntrain+i,j1) END DO ytest=ytest+pred(i) END DO ytest=ytest/REAL(ntest,dp) rmse1=0_dp mae1=0_dp mape2=0_dp DO i=1,ntest rmse1=rmse1+(pred(i)-a(ntrain+i,nft))**2 mae1=mae1+ABS(pred(i)-a(ntrain+i,nft)) if(abs(a(ntrain+i,nft)).le.0.01_dp) then mape2=mape2+abs(pred(i)-a(ntrain+i,nft)) else mape2=mape2+abs((pred(i)-a(ntrain+i,nft))/a(ntrain+i,nft)) end if END DO mse1=rmse1/REAL(ntest-mtmp,dp) ce1=1_dp-(rmse1/e1_test) rmse1=SQRT(rmse1/REAL(ntest,dp)) mae1=mae1/REAL(ntest,dp) mape2=mape2/REAL(ntest,dp) IF(ntest > 1) THEN sx=0_dp sy=0_dp rcor=0_dp DO i=1,ntest sx=sx+(pred(i)-ytest)**2 sy=sy+(a(ntrain+i,nft)-e0_test)**2 rcor=rcor+(pred(i)-ytest)*(a(ntrain+i,nft)-e0_test) END DO sx=SQRT(sx/REAL(ntest-1)) sy=SQRT(sy/REAL(ntest-1)) r1=rcor/(REAL(ntest-1,dp)*sx*sy) END IF WRITE(41,*) WRITE(41,*) 'Evaluation criteria for SVM piecewise linear regression (rmse, mae, mse, mape, ce, r): ' WRITE(41,*) rmse1!,mae1,mse1,ce1,r1 WRITE(41,*) mae1 WRITE(41,*) mse1 WRITE(41,*) mape2 WRITE(41,*) ce1 WRITE(41,*) r1 END IF mtmp=j IF (f_solution < tol) EXIT maxloop IF ((f_2-f_solution)/(f_1+1_dp) < tol) EXIT maxloop END DO maxloop IF (ipredi == 0) THEN ! Print only the final result IF (nrecord > ntrain) THEN ! prediction WRITE(41,*) WRITE(41,*) 'Final model with maxmin = ',mtmp,' and maxlin = ', maxlin mtmp=mtmp*nft*maxlin ytest=0_dp DO i = 1, ntest ! computation of correct piece pind=1 ipiecebest=1 piece_best=-large DO j=1,maxmin ipiecemin=pind piece_min=large DO j2=1,maxlin piece=x_unscaled(pind*nft) ! DO j3=1, nft-1 piece = piece+x_unscaled((pind-1)*nft+j3)*a(ntrain+i,j3) END DO IF (piece piece_best) THEN ipiecebest=ipiecemin piece_best=piece_min END IF END DO pred(i)=x_unscaled(ipiecebest*nft) DO j=1, nft-1 pred(i) = pred(i)+x_unscaled((ipiecebest-1)*nft+j)*a(ntrain+i,j) END DO ytest=ytest+pred(i) END DO ytest=ytest/REAL(ntest,dp) rmse1=0_dp mae1=0_dp mape2=0_dp DO i=1,ntest rmse1=rmse1+(pred(i)-a(ntrain+i,nft))**2 mae1=mae1+ABS(pred(i)-a(ntrain+i,nft)) if(abs(a(ntrain+i,nft)).le.0.01_dp) then mape2=mape2+abs(pred(i)-a(ntrain+i,nft)) else mape2=mape2+abs((pred(i)-a(ntrain+i,nft))/a(ntrain+i,nft)) end if END DO mse1=rmse1/REAL(ntest-mtmp,dp) ce1=1_dp-(rmse1/e1_test) rmse1=SQRT(rmse1/REAL(ntest,dp)) mae1=mae1/REAL(ntest,dp) mape2=mape2/REAL(ntest,dp) IF(ntest > 1) THEN sx=0_dp sy=0_dp rcor=0_dp DO i=1,ntest sx=sx+(pred(i)-ytest)**2 sy=sy+(a(ntrain+i,nft)-e0_test)**2 rcor=rcor+(pred(i)-ytest)*(a(ntrain+i,nft)-e0_test) END DO sx=SQRT(sx/REAL(ntest-1)) sy=SQRT(sy/REAL(ntest-1)) r1=rcor/(REAL(ntest-1,dp)*sx*sy) END IF WRITE(41,*) WRITE(41,*) 'Evaluation criteria for SVM piecewise linear regression (rmse, mae, mse, mape, ce, r): ' WRITE(41,*) rmse1!,mae1,mse1,ce1,r1 WRITE(41,*) mae1 WRITE(41,*) mse1 WRITE(41,*) mape2 WRITE(41,*) ce1 WRITE(41,*) r1 END IF END IF ELSE ! f_solution <= tol with only on min function under the maximum WRITE(42,*) WRITE(42,*) 'Fit function value for piecewise linear regression with ', j WRITE(42,*) 'linear functions under the minimum and ',1,' min under the maximum' WRITE(42,*) 'is ',f_solution IF (iscale == 1) THEN CALL rescaling !(x_solution,x_unscaled) ELSE x_unscaled=x_solution END IF ! Print the solution WRITE(42,*) 'x: ',x_unscaled WRITE(41,*) WRITE(41,*) 'Final model with maxmin = ',1,' and maxlin = ', j WRITE(41,*) 'Write prediction to "else" part.' WRITE(41,*) 'This never happened in our numerical experiments.' WRITE(41,*) 'If you get this result, please, send message to napsu@karmitsa.fi' END IF ELSE ! Compute the known model with maxmin and maxlin functions ! Set a new starting point DO i=1,maxmin*maxlin x_0((i-1)*maxdim+1:i*maxdim)=x_solution(1:maxdim) END DO CALL init_par() ! sets dimension for optimization problem and q. mk=maxlin ! each min function has maxlin linear functions maxmin2=maxmin CALL DBDC_algorithm( x_0, x_solution, f_solution, mit, mrounds, & & mrounds_clarke, termination, counter, agg_used, & & stepsize_used, iprint ) WRITE(42,*) WRITE(42,*) 'Fit function value for PWLR: ',f_solution IF (iscale == 1) THEN CALL rescaling !(x_solution,x_unscaled) ELSE x_unscaled=x_solution END IF ! print the results WRITE(42,*) 'x: ',x_unscaled IF (nrecord > ntrain) THEN ! prediction WRITE(41,*) WRITE(41,*) 'Final model with maxmin = ',maxmin,' and maxlin = ', maxlin ytest=0_dp DO i = 1, ntest ! computation of correct piece pind=1 ipiecebest=1 piece_best=-large DO j=1,maxmin ipiecemin=pind piece_min=large DO j2=1,maxlin piece=x_unscaled(pind*nft) ! DO j3=1, nft-1 piece = piece+x_unscaled((pind-1)*nft+j3)*a(ntrain+i,j3) END DO IF (piece piece_best) THEN ipiecebest=ipiecemin piece_best=piece_min END IF END DO pred(i)=x_unscaled(ipiecebest*nft) DO j=1, nft-1 pred(i) = pred(i)+x_unscaled((ipiecebest-1)*nft+j)*a(ntrain+i,j) END DO ytest=ytest+pred(i) END DO ytest=ytest/REAL(ntest,dp) rmse1=0_dp mae1=0_dp mape2=0_dp DO i=1,ntest rmse1=rmse1+(pred(i)-a(ntrain+i,nft))**2 mae1=mae1+ABS(pred(i)-a(ntrain+i,nft)) if(abs(a(ntrain+i,nft)).le.0.01_dp) then mape2=mape2+abs(pred(i)-a(ntrain+i,nft)) else mape2=mape2+abs((pred(i)-a(ntrain+i,nft))/a(ntrain+i,nft)) end if END DO mse1=rmse1/REAL(ntest-nft*maxmin*maxlin,dp) ce1=1_dp-(rmse1/e1_test) rmse1=SQRT(rmse1/REAL(ntest,dp)) mae1=mae1/REAL(ntest,dp) mape2=mape2/REAL(ntest,dp) IF(ntest > 1) THEN sx=0_dp sy=0_dp rcor=0_dp DO i=1,ntest sx=sx+(pred(i)-ytest)**2 sy=sy+(a(ntrain+i,nft)-e0_test)**2 rcor=rcor+(pred(i)-ytest)*(a(ntrain+i,nft)-e0_test) END DO sx=SQRT(sx/REAL(ntest-1)) sy=SQRT(sy/REAL(ntest-1)) r1=rcor/(REAL(ntest-1,dp)*sx*sy) END IF WRITE(41,*) WRITE(41,*) 'Evaluation criteria for SVM piecewise linear regression (rmse, mae, mse, mape, ce, r): ' WRITE(41,*) rmse1!,mae1,mse1,ce1,r1 WRITE(41,*) mae1 WRITE(41,*) mse1 WRITE(41,*) ce1 WRITE(41,*) r1 END IF END IF CALL cpu_time(finish_time) WRITE(41,*) WRITE(41,*) 'CPU time used: ',finish_time-start_time WRITE(42,*) WRITE(42,*) 'CPU time used: ',finish_time-start_time ! Changing the order back ! Maybe no need to do this. Instead do prediction. IF (noutcom < nft) THEN DO i=1,nrecord atmp = a(i,noutcom) a(i,noutcom)=a(i,nft) a(i,nft)=atmp END DO END IF CLOSE(42) CLOSE(41) CLOSE(43) CLOSE(78) deallocate(x_0,x0reg,x_solution,x_unscaled) deallocate(a,a4) deallocate(pred,dlin,clin) deallocate(mk) STOP END PROGRAM spr