module fit_mod implicit none contains ! > Routine to controll the genetic fitting algorithm subroutine fit(q,x1,x2,y,frms,difper,wt,par,p_spread,mut, > npar,p_act,seed,gtype,nset,nsel,chkpnt,old,iter,maxit, $ micit,ny,filename) use idxsrt_mod, only: idxsrt use dim_parameter,only: qn,numdatpt,ntot use init_mod,only: actinit use write_mod,only: write_output #ifdef mpi_version use mpi #endif implicit none ! MPI Variables #ifdef mpi_version integer ierror,my_rank,workernum,mpi_control_data(4) #endif ! Input variables (not changed within this subroutine). double precision q(qn,numdatpt),x1(qn,numdatpt),x2(qn,numdatpt)!< coordinates/input values double precision y(ntot,numdatpt),ny(ntot,numdatpt) !< Output/(energie and ci) values double precision wt(ntot,numdatpt) !< weights integer npar !< number of parameters integer nset !< number of parameter sets integer maxit !< maximum number of macroiterations integer micit !< maximum number of microiterations (i.e. LM iterations) !Fabian 15.03.2022: Used for babies or parent generation integer gtype!< type of random number generator --> is this ever really used?? integer nsel !< number of parents selected for having babies integer seed !< random seed for babies generation double precision p_spread(npar) double precision difper, mut !Fabian 15.03.2022: Used for checkfile character(len=80) :: chkpnt character(len=10) :: writer integer iter double precision old !< old rms !Fabian 15.03.2022: Used for wrout character(len=80) filename !Fabian 15.03.2022: Used in parameter initialization integer p_act(npar) !< array that contains info if a parameter is active or not == equivalent to ma in mrqmin ! Input/output variables (changed/updated within this subroutine) double precision par(npar,nset) !< parameters ! Output variables double precision frms !< best rms after macro iteration ! Internal variables integer i ! logical conver, ldum !< logicals for checking if calculation is converged logical ldum !< logicals for checking if calculation is converged integer start logical enough_parents integer mfit !< number of active parameters integer flag !< flag for write routine for fitting status(converged,maxiterationsreach,no convergence) ! Fabian 12.04. These are automatic arrays, maybe make them allocated or static integer idx(nset) !< array for sorting the parameter sets after their rms double precision rms(nset) !< array that contains rms of all parameter sets integer lauf !< counter for macroiteration double precision newpar(npar,nset) !< temporary storage array before parents&babies integer iact(npar) !< array pointing to the position of active parameters cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc #ifdef mpi_version call MPI_Comm_rank(MPI_COMM_WORLD, my_rank, ierror) #endif !> Initialize mfit,iact call actinit(p_act,iact,mfit,npar) #ifdef mpi_version call bcastact(p_act,npar) call MPI_Bcast(mfit, 1, MPI_INT, 0, MPI_COMM_WORLD,ierror) #endif !> Initialize rms vector rms=0.d0 rms(1:nset)=1d10 !> Write number of the present iteration and increase start to iter, if it is a restarted fit if (iter.ne.1) then write(6,*) 'Genetic restart, proceed with iteration', iter endif start=iter !> Start the genetic algorithm that consists of maxit macroiterations do lauf=start,maxit write(6,*) '' write(6,'(150("#"))') write(6,*) '' write(6,'(''Iteration:'',i5)') lauf !ATTENTION: THIS SUBROUTINE IS THE PARALLIZED SECTION !!! !Perform optimization for the parameter sets of generation lauf call fit_sets(lauf,nset,npar,par,rms, $ p_act,mfit,micit) !------------------------------------------------------------------------- !Sort the rms vector and assign the set number to each rms call idxsrt(rms,idx,nset) !write out sorted errors and indicate with which set each error was obtained do i=1,nset write(6,'(A8,I3,A8,F12.8,A8,I3)') 'Rank:', i,'RMS:', rms(i), $ 'Set',idx(i) enddo !write best rms onto the output variable frms frms=rms(1) !------------------------------------------------------------------------- !Resort the parameter array sucht that the parameter sets with the lowest rms are listed first newpar(1:npar,1:nset)=par(1:npar,idx(1:nset)) par(1:npar,1:nset)=newpar(1:npar,1:nset) !------------------------------------------------------------------------- !Return if maximum number of macro iterations is reached if (lauf.ge.maxit) return !------------------------------------------------------------------------- !Prepare next iteration of the genetic algorithm !Select the best parameter sets and sufficiently distinct sets as parents for the next iteration !Note: After parents, the first nsel entries of par and rms contain the parents !Note: However, rms is not strictly sorted after this (especially if the best parameter set were too similar) call parents(par,rms,difper,npar,idx,nset,nsel,p_act,mfit $ ,enough_parents) !Check for convergence of genetic algorithm, i.e. whether the generation of new parents leads to !a decrease of the rms as well as sufficiently distinct parameter set; return if convergence is reached ldum=conver(old,rms,idx,nsel) ! initialize flag for write routine flag=1 ! set converged flag for write routine if (ldum) flag=2 ! write intermediate output call write_output(q,x1,x2,y,wt,par,p_act,p_spread, > nset,npar,flag,lauf) if (ldum) return ! call flush ! flush(6) !Check if there are enough parents for next macro iteration if (enough_parents .eqv. .false.) then write(6,*) "Warning: Found too few different parents $ for next macroiteration, exit genetic algorithm" exit endif !Generate new parameter sets and proceed to the next iteration call babies(par,p_spread,mut,npar,mfit,nset,nsel,iact, $ seed,gtype) iter=iter+1 !------------------------------------------------------------------------- !write checkpoint: ! writer='write' ! call chkfile(chkpnt,par,npar,p_act,seed,gtype,nset,iter, ! & old,writer) !------------------------------------------------------------------------- enddo write(6,*) "Finished fit, return to main program" end subroutine C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% subroutine fit_sets(lauf,nset,npar,par,rms, $ p_act,mfit,micit) use dim_parameter,only: lbfgs use marq_mod,only: mrqmin use lbfgsb_mod,only: lbfgs_driver #ifndef mpi_version use omp_lib #else use mpi integer ierror,my_rank integer workernum #endif ! Input variables integer lauf !number of the current macroiteration integer nset !number of parameter sets integer npar !number of parameters !Input / output variables double precision par(npar,nset) !< parameters double precision rms(nset) !< array that contains rms of all parameter sets ! Input variables (necessary solely for mrqmin) integer p_act(npar) !< array that contains info if a parameter is active or not == equivalent to ma in mrqmin integer mfit !< number of active parameters integer micit ! number of microiterations ! Internal variables in parallel section double precision lrms !< rms for one parameter set double precision lpar(npar) !array for one parameter set !Fabian 31.03.2022: New test to reduce sice of parameters integer i,j ! Internal variables for OpenMP double precision startzeit,endzeit,start_totzeit,end_totzeit integer thread_num cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !> ATTENTION: THIS IS THE PARALLIZED SECTION !!! !> Perform non-linear least squares fit for each parameter set: #ifdef mpi_version ! write(*,*) 'mpi_version' start_totzeit=MPI_Wtime() call MPI_Comm_size(MPI_COMM_WORLD, workernum, ierror) call workshare(workernum, par, rms, npar, nset) end_totzeit=MPI_Wtime() #else start_totzeit=omp_get_wtime() !$omp parallel do schedule(dynamic) !$omp& default(shared) !$omp& private(i,j,lpar,lrms,thread_num,startzeit,endzeit) do i=lauf,nset ! > Fabian 15.03.2022: Variable for timing the duration of optimizing one parameter set startzeit=omp_get_wtime() !Fabian !> Write the parameters and the initial rms for this set onto private variables lpar(1:npar)=par(1:npar,i) lrms=rms(i) !Fabian 05.04.2022: Here I could separate the active and inactive parameters and perform the LM optimization purely with the active params !Fabian 05.04.2022: However, this would require to store the inactive parameter and the vector that decides if a variable is active onto a module since I need it in funcs then! !> Levenberg-Marquardt-Optimization of the current parameter set !Fabian 16.03.2022: This version might be MPI compatible since it contains purely of private variables !Fabian 16.03.2022: Use this instead of the above, if the data is declared global via a module and pst is only then used when necessary! if(lbfgs) then call lbfgs_driver(lpar,npar,p_act,mfit, & lrms,micit,i) else call mrqmin(lpar,npar,p_act,mfit, & lrms,micit,i) endif !> Write the optimized parameters and the optimized rms back onto the arrays that collect all parameters and rms par(1:npar,i)=lpar(1:npar) rms(i)=lrms !> Fabian 15.03.2022: Some output for timing the duration of optimizing one parameter set thread_num = omp_get_thread_num() endzeit=omp_get_wtime() write(6,*) 'Thread', thread_num ,'Time:', endzeit-startzeit !> Write output for the spezific set of parameters write(6,99) i, rms(i), rms(i)*219474.69d0 99 format('Set:',i5,5x,'RMS:',g16.6,3x,'(',f16.6,' cm-1)') enddo !$omp end parallel do end_totzeit=omp_get_wtime() #endif write(6,*) 'Total time for Macroiteration: ' > ,end_totzeit-start_totzeit write(6,*) 'Finished parallel fit for Iteration', lauf end subroutine C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C% SUBROUTINE PARENTS(...) C% C% subroutine to select the parent parameter sets according to their C% RMS error C% C % variables: C % par: parameter vector (double[npar,nset]) C % rms: error for each set (double[nset]) C % difper: C % npar: number of parameters (int) C % idx: sorted indeces according to rms(1..nset) (int[nset]) C % nset: number of sets C % nsel: number of selected parents C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% subroutine parents(par,rms,difper,npar,idx,nset,nsel,p_act,mfit $ ,enough_parents) implicit none integer i, j, k, nset, idx(nset), npar, nsel, p_act(npar), mfit double precision par(npar,nset), dum(npar,nset), rms(nset), last double precision thr double precision difper, drms(nset) integer dum_idx(nset), rank_parent(nsel) ! logical difchk logical enough_parents thr=1.d-8 dum=0.d0 dum_idx = 0 rank_parent = 0 drms=0.d0 c write the best parameter set on the dummy dum(1:npar,1)=par(1:npar,1) dum_idx(1)=idx(1) rank_parent(1) = 1 !Choose exactly (beside the best parameter set) nsel-1 parameter sets as new parents and write them on dum !These parents are selected according to the lowest possible rms AND sufficient dissimilarity !to the overall best parameter sets last=1.d14 k=1 do i=1,nset if (difchk(dum,par(1:npar,i),difper,k,npar,p_act,mfit,nset)) > then k=k+1 dum(1:npar,k)=par(1:npar,i) drms(k)=rms(i) dum_idx(k) = idx(i) rank_parent(k) = i endif if (k.eq.nsel) exit enddo !Terminate programm if too few parents are found enough_parents=.true. if(k.lt.nsel) then enough_parents=.false. endif !Copy the selected parent parameter sets back to the array par do i=2,nsel par(1:npar,i)=dum(1:npar,i) rms(i)=drms(i) enddo !Write out some information on the chosen parent parameter sets write(6,*) 'nsel:', nsel write(6,*) write(6,*) 'Selected parents:' do j=1,nsel write(6,201) rank_parent(j), rms(j), dum_idx(j) write(6,200) (par(k,j), k=1,npar) enddo 200 format('Par:',6g16.7) 201 format('>>> Rank:',i5,' RMS:' ,g14.4,' set:',i5,' <<<' ) ! call flush ! flush(6) end subroutine !---------------------------------------------------------------------- ! function to check whether new parameter set is sufficiently different ! from already selected sets: logical function difchk(dum,par,difper,k,npar,p_act,mfit,nset) implicit none integer i, j, k, npar, p_act(npar), mfit,nset double precision dum(npar,nset), par(npar), per, thr, difper double precision epsilon parameter(epsilon=1d-8) !.. this threshold specifies that parameter set must have an average ! difference of at least 1% with respect to any other selected set. thr=1.d0-difper if (thr.gt.0.99d0) thr=0.99d0 !avoids no difference difchk=.true. do i=1,k per=0.d0 !Calculate relative difference between between current set (par) and the already selected sets (dum) do j=1,npar if (p_act(j).ge.1) then !Added flexible value for p_act; Nicole 15.12.2022; only active parameters are counted per=per+(min(dum(j,i),par(j))+epsilon) $ /(max(dum(j,i),par(j))+epsilon) endif enddo per=per/mfit !Modified Version that only active parameters are counted; Fabian 14.12.2021 !Discard the current set if it is too similar to one already selected if (per.gt.thr) then difchk=.false. return endif enddo end function !-------------------------------------------------------------------- ! subroutine to create the baby sets of parameters from the selected ! parent sets subroutine babies(par,p_spread,mut,npar,mfit,nset,nsel,iact, $ seed,gtype) implicit none c functions double precision rn !gets one random number integer i, j, k, npar, nset, nsel, mfit, iact(npar) double precision par(npar,nset), p_spread(npar), mut, dum integer seed,gtype !loop over all dieing sets (only the nsel parent sets survive) do i=nsel+1,nset !loop over all active parameters do j=1,mfit !picking a random parameter set of the first nsel parent sets !(Fabian 16.03.2022: Add feature, to ensure that at least one baby is generated from each parent?) k=int(rn(seed,gtype,0)*nsel)+1 !Fabian 08.04.2022: Even though seed isnt passed here, the rn call is dependent on the earlier initialized seed !writing the j'th parameter of the selected parent set onto the j'th parameter of the i'th of the remaining sets (only the active parameters are copied) !(Fabian 16.03.2022: This way, I recombinate a number of parents to new babies. However, recombination might not be good, if these parent sets are relatively distinct; maybe use only two parent sets for recombination?) par(iact(j),i)=par(iact(j),k) !select whether the j'th parameter of this new set is mutated !(Fabian 16.03.2022: Add feature, to ensure that at least one parameter is mutated?) if (rn(seed,gtype,0).lt.mut) then dum=rn(seed,gtype,0) - 0.5d0 dum=dum*p_spread(iact(j)) par(iact(j),i)=par(iact(j),i)*(1.d0+dum) endif enddo enddo end subroutine !----------------------------------------------------------------------- ! check convergence of genetic algorithm function conver(old,rms,idx,nsel) implicit none integer i, j, nsel, idx(*), baby double precision rms(*), new, old, thresh, percent, thrper logical conver !Thresholds and initializiation conver=.false. thresh=old*1.d-3 thrper=0.2d0 ! Lets use all values in the selected subset: j=nsel baby=0 ! Calculate average error for the nsel best parameter sets new=0.d0 do i=1,j new=new+rms(i) enddo new=new/dble(j) ! calculate the number of selected parent sets that were originally babies in the previous iteration do i=1,nsel if (idx(i).gt.nsel) baby=baby+1 enddo ! calculate the percentage percent=dble(baby)/dble(nsel) ! some output write(6,100) baby write(6,101) new, j write(6,*) 100 format('Number of babies in chosen subsets:', i3) 101 format('Average RMS error of chosen subsets:', g12.4, $ ' / averaged values:', i4) write(6,110) percent*100.d0 write(6,111) old, new, old-new 110 format('Percent babies:',f6.1) 111 format('Old RMS:',d12.4,' New RMS:',d12.4,' Diff:',d12.4) !Set convergence to true if !1. too few previous babies are among the new parents !2. or the average rms of the selected parents between the current & previous macro iteration is sufficiently small conver=(percent.le.thrper).and.(abs(new-old).lt.thresh) write(6,*) 'Convergence:', conver !Set average rms of this iteration to the comparison variable old for the next iteration old=new end function end module fit_mod