Subroutine optsensitivity

633 : 
634 : !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
635 : !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
636 : !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
637 : ! ADDED BY STAN TO ASSESS SENSITIVITY OF OBJECTIVE FUNCTION TO EACH OPTIMISED
638 : ! PARAMETER:
639 : !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
640 : 
641 :       subroutine optsensitivity ()
642 :       use vom_sce_mod
643 :       implicit none
644 : 
645 :       INTEGER              :: ii, jj, kk, failed10, pos
646 :       REAL*8               :: distmin, distmax
647 :       REAL*8               :: oldpar, newpar
648 :       REAL*8               :: ofvec2, ofchange, parchange
649 :       CHARACTER(300)       :: writeformat
650 :       CHARACTER(len=135)   :: msg
651 :       REAL*8, ALLOCATABLE  :: tmp_8(:)
652 : 
653 :         shufflevar2(:) = shufflevar(:,1)
654 :         ofvec2 = ofvec(1)
655 :         failed10 = 0
656 : 
657 :         write(kfile_progress,*) 'SENSITIVITY ANALYSIS'
658 :         write(kfile_progress,*) 'changes in % of feasible range)'
659 :         dataarray(1,1:nopt) = shufflevar2(optid(:))
660 :         dataarray(1,nopt+1) = ofvec2
661 :         evolution = 'test'
662 :         pos = 0
663 : 
664 :       do ii = 1, nopt
665 : 
666 :           oldpar = shufflevar2(optid(ii))
667 : 
668 : !         * TO MAKE SURE THAT PERTURBATIONS DO NOT EXCEED THE FEASIBLE RANGE
669 : !         * AND THE PARAMETER VALUES THEMSELVES
670 : 
671 :           distmin = oldpar - parmin(optid(ii))
672 :           distmax = parmax(optid(ii)) - oldpar
673 :           write(kfile_progress,*) " "
674 :           writeformat = '("change of var: ",A9, "(",E9.3,")")'
675 :           write(msg,writeformat) parname(optid(ii)), oldpar
676 :           write(kfile_progress,*) TRIM(msg)
677 :           writeformat = '(f7.3,"% (",e14.7,")",": change of OF by ",'
678 :           writeformat(44:68) = 'e10.3,"%"," (",e10.3,")")'
679 : 
680 :           jj = 0                        ! counts from 0 to 3 and than from 3 to 0
681 :           kk = 1                        ! counter: first +1 after 3 steps it becomes -1
682 : 
683 :         do while (jj .ge. 0)
684 : 
685 :             if (kk .eq. 1) then
686 :               newpar = oldpar - distmin * 10.d0 ** (-jj)
687 :             else
688 :               newpar = oldpar + distmax * 10.d0 ** (-jj)
689 :             endif
690 :             nrun = nrun + 1
691 :             shufflevar2(optid(ii)) = newpar
692 : 
693 :           call transpmodel(shufflevar2(:), SIZE(shufflevar2(:)), ofvec2, 1)
694 : 
695 :             if (ofvec2 .gt. bestobj) then
696 :               bestobj = ofvec2
697 :               call write_lastbest(shufflevar2(:), vom_npar, bestobj, 0)
698 :             endif
699 : 
700 : !             * use temporary variable to prevent warning in ifort
701 :               allocate(tmp_8(nopt))
702 :               tmp_8(:) = shufflevar2(optid(:))
703 :               write(kfile_sceout,outformat) tmp_8(:), ofvec2
704 :               deallocate(tmp_8)
705 : 
706 :             if (kk .eq. 1) then
707 :               dataarray((ii-1)*8+jj+2,1:nopt) = shufflevar2(optid(:))
708 :               dataarray((ii-1)*8+jj+2,nopt+1) = ofvec2
709 :             else
710 :               dataarray((ii-1)*8+jj+6,1:nopt) = shufflevar2(optid(:))
711 :               dataarray((ii-1)*8+jj+6,nopt+1) = ofvec2
712 :             endif
713 :             ofchange = (ofvec2 - dataarray(1,nopt + 1))                &
714 :      &               / ABS(dataarray(1,nopt + 1)) * 100.d0
715 :             parchange = (newpar - oldpar) / (parmax(optid(ii))         &
716 :      &                - parmin(optid(ii))) * 100.d0  ! the change of the parameter in % of feasible range
717 :             write(msg,writeformat) parchange, newpar, ofchange, ofvec2
718 :             write(kfile_progress,*) TRIM(msg)
719 :             flush(kfile_progress)
720 : 
721 :             if (ofchange .gt. 1.0d-10) then
722 :               shufflevar(:,sopt-pos) = shufflevar2(:)
723 :               ofvec(sopt-pos) = ofvec2
724 :               pos = pos + 1
725 :               if (ABS(parchange) .gt. i_resolution) then
726 :                 failed10 = failed10 + 1
727 :               endif
728 :             endif
729 : 
730 :             jj = jj + kk
731 :             if (jj .eq. 4) then
732 :               jj = 3
733 :               kk = -1
734 :             endif
735 :         enddo
736 : 
737 :           shufflevar2(optid(ii)) = oldpar
738 : 
739 :       enddo
740 : 
741 : !     * ASSESS SECOND CONVERGENCE CRITERIUM: NO PARAMETER CHANGE BY MORE THAN
742 : !     * 10% LEADS TO AN INCREASE IN OBJECTIVE FUNCTION
743 : 
744 :         if (failed10 .gt. 0) then
745 :           writeformat = '(I2," parameter(s) more than ",F6.3,"% out of optimum.")'
746 :           write(msg,writeformat) failed10, i_resolution
747 :           write(kfile_progress,*) TRIM(msg)
748 :           write(kfile_progress,*) "Optimisation continued..."
749 :         else
750 :           write(kfile_progress,*) " "
751 :           write(kfile_progress,*) " "
752 :           write(kfile_progress,*) "Second convergence criterion satisfied:"
753 :           writeformat = '("no parameter shift by more than ",F6.3,"% of max distance leads to")'
754 :           write(msg,writeformat) i_resolution
755 :           write(kfile_progress,*) TRIM(msg)
756 :           write(kfile_progress,*) "an increase in the objective function."
757 :           success = 1
758 :         endif
759 : 
760 : !     * IF CRITERIUM NOT SATISFIED, APPEND NEWLY CREATED DATA TO THE END OF
761 : !     * shufflevar AND ofvec ARRAYS AND RETURN TO MAIN LOOP TO RE-RUN SCE-LOOP
762 : !     * --> deleted
763 : 
764 :       return
765 :       end subroutine optsensitivity