!Tweak-Begin.F95
!2018.10.09.1015cdt JMS- Introducing:  "Tweak-Engine",
!                                   *******************
!                                  &: Global Variables.

!  Jeffrey M. Setterholm, 8095 230th St. E., Lakeville, Minnesota 55044, USA
!   I have authored the four Fortran *.f95 source code files listed below. 
!                       I hereby place these four files:
!       Tweak-Begin.f95, Tweak-Engine.F95, Tweak-User.F95, & Tweak-Vis.f95
!               and the algorithms which are demonstrated therein,
!                     in the public domain (a.k.a.: "free").
!Disclaimer:
!       ***************************************************************
!       **********   Individual cognition is always flawed,  **********
!       **********           including yours & mine.         **********
!       **********                  - So: -                  **********
!       **********       Use this code at your own risk.     **********
!       ***************************************************************

!This code realizes my present understanding of:

!                  Best-Fitting the Continuous Parameters of
!   Your   :       Linear/Non-Linear Algorithmic Models of
!   Your   :             Arbitrarily Large Datasets
!                          with Extreme Accuracy.
 
!When YouTweak.f95 is tailored for your purposes:

!               You can determine your math model's adequacy,
!            or need for improvement,in accounting for your data.

!If you enjoy searching,    -and-
! when these algorithms are subsequently perfected,
!   you'll find that:
 
!             Mother nature's error evaluations of math models
!      are delightfully impartial, non-intrusive, and deeply insightful...
!     much like the smiles & frowns of the great teachers in the Himalayas.

!******************************************************************************

!Tweak-Begin.f95 organizes the key variables within Tweak.exe
!An overview of the Engine (timestamps & sizes are approximate):

!                        Lines    Bytes
!10/08/2018  05:30 AM             3,560 Tweak.gui       <- Absoft's GUI
!10/08/2018  11:54 AM     ~470   23,907 Tweak-Begin.f95  <- This file
!                                         Module TweakRec 
!                                         Module KPrec
!                                         Module LZrec 
!                                         Module MDrec
!                                         Module MDrec
!                                         Subroutine TweakAllocateAll 
!10/08/2018  08:11 AM    ~1260   60,540 Tweak-Engine.F95    Includes:
!                                         Program Tweak
!                                         Subroutine EvalFit(RSSL,iP)
!                                         Subroutine PrintIter(iP)
!                                         Subroutine RSSPartials(jpUL)
!                                         Subroutine DatapointPartials
!                                         Subroutine Invert(N,A,ValMin,iRank...
!                                         Subroutine PrintA(N,A,Noise,iRank,...
!                                         Subroutine SelectStepMult(iP)
!                                         Function om(Value1,Value0,iP)
!                                         Subroutine FloatWrite(R16In,a40out)
!                                         Subroutine FDate23(DaTime)
!                                         Subroutine Beamer(n,nTot)
!10/08/2018  08:14 AM     ~400   19,167 Tweak-User.F95  <- Your Model goes here
!                                         Subroutine YouTweak(Mode)
!                                         Subroutine YouTweak3D

!10/08/2018  07:22 AM     ~850   50,436 Tweak-Vis.f95    <- 3D Visualization
!                                         Subroutine Draw3dJTL(iColor,...
!                                         Subroutine Morph3dJTL(TLin,iP)
!                                         Function Log10LF(Value1,iP)
!                                         Module CharDef
!                                          ~an English character vector font
!                                         Subroutine AlphaJS(cLabelL,PosLLCq...
!                                         Subroutine Model7DoF(Xyzh,Rpyh,S,...

!10/08/2018  12:27 PM         1,037,824 Tweak.exe
!10/08/2018  11:58 AM             1,621 Tweak-User.nml   <- For mods @ runtime

!10/08/2018  12:03 PM           141,391 Tweak.txt        <- Reports the process
!10/08/2018  06:18 AM             1,301 Tweak-Params.txt <- Reports parameters

!10/08/2018  12:03 PM         1,412,272 Tweak-3dDraw.3dv <- 3D output.

!             ~9 File(s)     ~3,089,439 bytes

!Subdirectory /Images has 23 selfies of five solution trajectories.

!******************************************************************************

!Five "modules" follow herein. The summarizing 'use' statements are:
!---                                                                  globals
!use Tweakrec,only: jPhase,jMode  ,cVersion,cDateTime & !Tweak's FYI
!                  ,jBOn,jUnClamp,jStepMult,jPrev     & !Solution strategy
!                  ,jpU10,jpU,jpD           ,cFloat40 & !Printout-Alphanumeric
!                  ,jpU3d,TLrec,TL,TLprev,TL2,TLsave  & !Printout-3D
!                  ,jTLmorph,TLiter,omjSave           & !         "
!                  ,jItertot,jIter,jDone              & !Iteration control 
!                  ,RSS,RSSbase,Weight,Delta,offon    & !Tweaking & errors
!                  ,AbsDet,NoiseFloor,iRank,kPChanged & !Inverter outputs
!                  ,jMintot,jMin,StepMult,omj         & !Minimization passes
!                  ,B,BtB,BtZ                         & !Allocated matrices
!                  ,What,Why,How,Who,When,Where1      & !Project context
!                  ,jUserPhase,jUserConfig,cjUserFile & !Use in YouTweak()
!                  ,TweakNml                            !Runtime reconfig.
!Use KPrec,   only: Kptot,Kp,Kp2,  Pr   ,Pu,PstepMag    !Parameters -to fit to-
!Use LZrec,   only: Lztot,Lz,      Zr   ,Zu,Zu2,Z0      !Outputs    - of your -
!use MDrec,   only: Mdtot,Md,MdMax,Datae                !Dataset
!use NXrec,   only: Nxtot,Nx,      Xr                   !Independent variables
!                                  ^^:These are allocated records.
!                    ^^:At the outset you size these four values.
!---

!Table of Contents:
!Module TweakRec !Program control and various summary values:
!Module KPrec    !Parameters:
!Module LZrec    !Outputs:
!Module MDrec    !dataset storage & access:
!Module MDrec    !dataset storage & access:
!      &
!Subroutine TweakAllocateAll !Allocates Tweak's arrays.
!---------------------------------------------------------------------------7 9
!//////////////////////////////////////////////////////////////////////////////

Module TweakRec !Program control and various summary values:      !use TweakRec
!              = j- integers & various reals

!Tweak's FYI
  character::jPhase*80 !Marks key phases of the fitting process.    Purely: FYI
!                                                          For Your Information
!              ="0: Allocates records based on Kptot,Lztot,Mdtot,& Nxtot"c
!                                                in YouTweak(0)
!              ="1: You quantify your system model in the allocated arrays"c
!                                                in YouTweak(1)
!                            YouTweak(2,3,& 4) - exercised by the Tweak engine.
!              ="2: Top-of-each iteration: establish baseline residual errors"c
!              ="3: Populate BtB(Kptot,Kptot) and BtZ(Kptot)"c
!              ="4: Normalize and invert BtB(Kptot,Kptot)"c 
!              ="5: Inverters stepsize: Pr.Pistep = BtBinv(denormalized)*BtZ"c

!              ="3,4,& 5: RSSPartials uses RSS error gradient instead."c

!              ="6: Find an error-reducing actual stepsize"c
!              ="7: Iterative pass done. Your intervention opportunity:"c
!                                                in YouTweak(5)
!              ="8: Tweak is done interating. Export the results & close"c
!                                                in YouTweak(6)

  integer*4::jMode !Interface of Tweak() to YouTweak().             Purely: FYI
!              = 0: Define Nxtot, Lztot, Mdtot, & Kptot:
!              = 1: Initialize Xr(), Zr(), Dr(), Pr(), & Datae():
!              = 2: Md=0  ~ "Rewind" your dataset:
!              = 3: Md=Md+1 - Datae(*,Md) -> Xr.Xdata and Zr.Zdata:
!              = 4: Exercise your system. Compute Zr.Z using Pr.P:
!              = 5: End of iterative pass -intervention opportunity:
!              = 6: Iterating done. Final printout(s) opportunity:

  character::cVersion*60  & !Present software version
               ="Tweak.Exe, version:1.0, 2018.10.09, JMS, Public Domain"c

  character::cDateTime*23!Current date & time. 

!Solution strategy:
  integer*4::jBOn     = 3 !Iteration# to start calling DatapointPartials()
!                          For jIter<jBOn the call is to RSSPartials() instead
!                     = 1 :Will probably work well for linear models. 
!                     > 1 :A few interations using RSSPartials()
!                          hastens the convergence of non-linear models.

  integer*4::jUnClamp = 3  !Number of iterations to restrict Pr.Pistep to <= 1.

  integer*4::jStepMult= 2!Minimization algorithm - computes StepMult.
!                     = 1!   - StepMult=1.0    .e.g.: using inverter's result
!                     = 2!   - calls SelectStepMult()

  character*9::omj  &    !Expresses the value1/Value0 ratio of positive numbers
             ='undefined'! in exponential form, & places the signed exponent
                         ! on the left side of the result.
                         !Exponents > 0 represent ratios >= 1.
                         !Exponents <=0 represent ratios <  1.
                         !1. > mantissa >= .1, except when Value1 = 0.
                         !"__sme.nnn"
                         !"_smme.nnn"

  integer*4::jPrev    =-1!Use previous parameter values in Tweak-Params.txt?
!                     =-1: no
!                     = 0: ask
!                     = 1: yes
!                         Ordinarily you'll solve each problem from scratch.
!                         For problems whose iteration time is more than a few
!                         minutes, easy resumption is convenient.

!Printout-Alphanumeric:
  integer*4::jpU10   =10 !Unit number- for Tweak's parameter reads & writes
!                                             to/from "Tweak-Params.txt"
  integer*4::jpU     =11 !Unit number- Output of general information
!                    = 6 : outputs to DOS screen                  
!                   >=11 : writes/appends "Tweak.txt"
  integer*4::jpD     = 4 !Print Detail level, use: "if(jpD >= ___)"
!                    = 0 : Disabled, e.g: when jpU<6
!                    = 1 : User selected important information/program controls
!                    = 2 : Enables screen/keyboard/program interactions
!                    = 3 : Report/summarize- final iteration
!                    = 4 :                 - each iteration
!                    = 5 :                 - some main program info.
!                    = 6 :                 - completed matrices & minimization
!                    = 7 :                 - some subroutine details
!                    = 8 :                 - in-depth details
  character::cFloat40*40 !used for: call FloatWrite(R16In,cFloat40) printing

!Printout-3D: -----
!YouTweak can interface with 3D visual data collector Draw3dJTL().
!  A minimal interface follows. 
!If no calls to Draw3dJTL() are made, then TweakVis.f95 isn't needed
!  to compile & link Tweak.exe.
  integer*4::jpU3d  = 13 !Unit number- for exporting 3D viewing coordinates.

  integer*4::jTLmorph = 1!Morph3dJTL() mode control.
!                           Morph3dJTL() centralizes the arbitrary rescale
!                             (~morphing) of TL.XYZ plot points.
!                         Example uses:
!                     = 0: Z-axis- no change
!                                     (<,-1.) [-1.,1.]   ((1.,>)   :Z
!                     = 1:       - -log10(-Z)| linear |+log10(+Z)
!                     = 2:       -   2.-1./Z | linear | -2.+1./Z
!                     = 3:       - User defined morph of X,Y, & Z.

 type:: TLrec            !"Thick Line" record
  integer*4:: iC           ! DOS color of the line [0,15] to the next dot
!                          =  0: move-to without drawing
!                          = -1: closes the file
  real*4   :: Th           ! Line thickness [1.0 - 20.0] & width      in pixels
  real*16  :: XYZ(3)       ! (X,Y,Z) "next dot"
 end type TLrec; 
  type(TLrec):: TL       !Draw3dJTL- Your "next dot" input
  type(TLrec):: TLprev   !         - the most recent dot added
  type(TLrec):: TL2      !         -      "next dot" used by AlphaJS
  type(TLrec):: TLsave   !A place to save TLprev during dot digressions.
  type(TLrec):: TLiter(0:100)
  character*9:: omjSave(0:100) !omj values of the first 100 iterations
!------------------ 

!Iteration control:
  integer*4::jItertot=30!Fit iterations- total
  integer*4::jIter      !              - current #
  integer*4::jDone      != 0: will do another iteration
                        !> 0: reduced error not found.

!Evaluation:
  real*16  ::RSS        !Root-Sum-Squared- present
  real*16  ::RSSbase    !                - previous iteration
!Numerical Partial Derivative default values:
  real*16  ::Weight     =  1._16     !Nominal output weighting factor
!                       =  0._16     !                       - set individually
  real*16  ::Delta      =  1.e-6_16  !Numerical partial tweak-
!                                    -adjusted in YouTweak(5)
!                       =  0._16     !                       - set individually
  real*16  ::offon      =  1._16     !Defaults to all- on
!                       =  0._16     !               - off- set individually

!Inversion output values:
  real*16  ::AbsDet                  !Determinant- absolute value
  real*16  ::NoiseFloor              !Inversion noise floor
  integer*4::iRank                   !BtB's rank, < Kptot if linearly dependent
  integer*4::kPChanged               !# of parameters whose values changed.

!Minimization passes to determine the step multiplier:
  integer*4::jMintot =20
  integer*4::jMin
  real*16  ::StepMult    !Pr.Pnew = Pr.Pbase + SetpMult * Pr.Pistep

  real*16   ,allocatable::B(  :,:) !(Lztot,Kptot) Numerical partials array
  real*16   ,allocatable::BtB(:,:) !(Kptot,Kptot) Outer product array
  real*16   ,allocatable::BtZ(:  ) !(Kptot      )

!Project Context: Consider having your dataset(s) define these:
  character:: What  *60 ! & !Title    = "Title this particular parameter fit."c
  character:: Why   *60 ! & !Objective= "Briefly describe your objective(s)."c
  character:: How   *60 ! & !Approach = "Briefly describe your approach."c
  character:: Who   *60 ! & !Author(s)= "The model & dataset authors."c
  character:: When  *60 ! & !date/time= "The model/dataset date?"c
  character:: Where1*60 ! & !Location = "Describe the location."c
!                    ^^ :60 characters max

!User-assignable variables, for use in YouTweak()
  integer*4:: jUserPhase     
  integer*4:: jUserConfig
  character:: cjUserFile*60

!Namelist interface to TweakRec - for runtime reconfiguration of the .exe:
NAMELIST / TweakNml /                 & 
   jBOn,jUnClamp,jStepMult            &  !Solution strategy
  ,jpU,jpD                            &  !Printout-Alphanumeric
  ,jTLmorph                           &  !Printout-3D
  ,jItertot                           &  !Iteration control
  ,Delta                              &  !Numerical Partials
  ,jMintot                            &  !Minimization passes
  ,What,Why,How,Who,When,Where1       &  !Project Context
  ,jUserPhase,jUserConfig,cjUserFile     !User-assigns, for use in YouTweak()

!Contains
End Module TweakRec
!//////////////////////////////////////////////////////////////////////////////

Module KPrec    !Parameters:                                         !Use KPrec
!                ... the "parameters" of "parameter identification".
!              = P's with K indices
!Size & indexes:
  integer*4::Kptot    !Parameter- Total number
  integer*4::Kp       !         - # in use
  integer*4::Kp2      !         - # in use

!Parameter properties:
 type:: Prec          !         - type:
  real*16  ::P        !           - value- in use
  real*16  ::Pbase    !           -      - baseline at iteration start
  real*16  ::offon  & !           = 0._16: disables fitting
                =1._16!           = 1._16:  enables fitting    
  real*16  ::Delta    !           - numerical differentiation increment
! real*16  ::Pistep   !           - Initial step prediction

  real*16  ::Pstep(4,2)!          - Step size- (1,1): from RSS 
!                                              (1,2):         scaled/bounder
!                                              (2,1): from BtBinv*BtZ   "Inv"
!                                              (2,2):         scaled/bounder
!                                              (3,1):                <<unused>>
!                                              (3,2):                <<unused>>
!                                              (4,1): input to minimizer
!                                              (4,2): Minimizer's net tweak
  real*16  ::Pnew     !           - New coefficient value
  real*16  ::Vstep    !           - V-minimizer output
  integer*4::Inverted !           = 1: linearly independent
!                     !           = 0: linearly   dependent, Pr().Pistep=0.
  integer*4::Kp       !           - # of this Pr(Kp)
  integer*4::ispare(2)!
  character::Pname*16 !		    - Name
  real*16  ::BtBnorm  !y Diagonal of BtB()- used as a normalizing factor
 end type Prec
 type(Prec)            ::Pu,P0
 type(Prec),allocatable::Pr(:)    !(KpTot)

  real*16  ::PstepMag !Magnitude of the Pr.Pistep(s) vector
!Contains
End Module KPrec

!Adjustment of Delta is shown in YouTweak(5)
!//////////////////////////////////////////////////////////////////////////////

Module LZrec    !Outputs:                                            !Use LZrec
!              = Z's with L indices:
!Size & indexes:
  integer*4::Lztot    !Output- Total number
  integer*4::Lz       !      - # in use
!Output properties:
 type:: Zrec          !      - type:
  real*16  ::Z        !y        - value- in use
  real*16  ::Zbase    !                - baseline at iteration start
  real*16  ::Wt       !y        - Weight of output in the fit, Nominally 1.0
  real*16  ::WZ       !y        - Weight*value
  real*16  ::Zdata    !         - value- of/from the current Md datapoint
!                                        Managed & used by YouTweak()
  real*16  ::RMS      !y        - error- RMS- present
  real*16  ::RMSbase  !                     - at iteration start    
  real*16  ::Zmax     !y               - max- value- computed by EvalFit()
  integer*4::Mdmax    !y                    - @ Md#
  integer*4::Lz        !         - # of this Zr(Lz)
  integer*4::ispare(2)!
  character::Zname*16 !         - Name
 !integer*4::L        !         - #
 end type Zrec
 type(Zrec)            ::Zu,Zu2,Z0
 type(Zrec),allocatable::Zr(:) !(Lztot)
!Contains
End Module LZrec
!//////////////////////////////////////////////////////////////////////////////

Module MDrec    !Dataset storage & access:                           !use MDrec
!              = Datae with M indices:
!Size & indexes:
  integer*4::Mdtot   !Datapoints- total number in the Dataset
  integer*4::Md      !          - # in use
  integer*4::Mdmax
  real*16,allocatable::Datae(:,:) !( 0 :Nxtot+Lztot,Mdtot)
!                                   ^0 stores the dataset index as a real*16
!Contains
End Module MDrec
!//////////////////////////////////////////////////////////////////////////////

Module NXrec    !Independent variables of your dataset:              !use NXrec
!              = X's with N indices:

!Size & indexes:
  integer*4::Nxtot    !Variables- Total number
  integer*4::Nx       !         - # in use
!Variables properties:
 type:: Xrec          !         - type:
  real*16  ::Xdata    !           - value- of/from the current Md datapoint
!                                          Managed & used by YouTweak()
  integer*4::Kcoes    !           - number of (e.g. polynomial) coefficients.
  integer*4::Nx       !           -  # of this Xr(Nx)
  integer*4::ispare(2)
  character::Xname*16 !		      
 end type Xrec
 type(Xrec)            ::Xu,X0
 type(Xrec),allocatable::Xr(:) !(Nxtot)
!Contains
End Module NXrec
!//////////////////////////////////////////////////////////////////////////////
!---------------------------------------------------------------------------7 9

Subroutine TweakAllocateAll  !Allocates Tweak's arrays.
!2018.09.25.2010cdt JMS- Traveler2/Athlon64/WinXPPro-32/APF9.0-32
 use TweakRec, only:B,BtB,BtZ                                        !Globals
 Use KPrec
 Use LZrec
 use MDrec
 use NXrec
 !---
 implicit none
  !---                                                               !internals
   integer*4::iAlloc
   !-----                                                            !end defs
   if(Kptot<=0) then; write(6,"(1x,'Tweak halt!  Kptot=',i9)") Kptot ;stop;endif
   if(Lztot<=0) then; write(6,"(1x,'Tweak halt!  Lztot=',i9)") Lztot ;stop;endif
   if(Mdtot<=0) then; write(6,"(1x,'Tweak halt!  Mdtot=',i9)") Mdtot ;stop;endif
   if(Nxtot<=0) then; write(6,"(1x,'Tweak halt!  Nxtot=',i9)") Nxtot ;stop;endif
!  --------  
   allocate(Pr(Kptot),stat=iAlloc)          !Parameters array
     if(iAlloc /= 0) then
       pause 'TweakAllocateAll: Pr(Kptot) allocation error. `Enter` to exit.'
       stop
     endif!(iAlloc/=0)
   Pr(1:Kptot)=P0; Pr.Pname=char(32)
!  --------
   allocate(Zr(Lztot),stat=iAlloc)          !Outputs array
     if(iAlloc /= 0) then
       pause 'TweakAllocateAll: Zr(Lztot) allocation error. `Enter` to exit.'
       stop
     endif!(iAlloc/=0)
   Zr(1:Lztot)=Z0; Zr.Zname=char(32)
!  --------
   allocate(Datae(0:Nxtot+Lztot,Mdtot),stat=iAlloc) !Data array
     if(iAlloc /= 0) then
       pause 'TweakAllocateAll: Datae() allocation error. `Enter` to exit.'
       stop
     endif!(iAlloc/=0)
   Datae=0._16
!  --------
   allocate(Xr(0:Nxtot),stat=iAlloc)        !Independent variables array
     if(iAlloc /= 0) then
       pause 'TweakAllocateAll: Xr(Nxtot) allocation error. `Enter` to exit.'
       stop
     endif!(iAlloc/=0)
   Xr(0:NxTot)=X0; Xr.Xname=char(32)
!  --------
   allocate(B(Lztot,Kptot),stat=iAlloc)     !B: Numerical partials array
     if(iAlloc /= 0) then
       pause 'TweakAllocateAll: B(Lztot,Kptot) alloc. error.`Enter` to exit.'
       stop
     endif!(iAlloc/=0)
   B=0._16
!  --------
   allocate(BtB(Kptot,Kptot),stat=iAlloc)   !Bt*B: Outer product array
     if(iAlloc /= 0) then
       pause 'TweakAllocateAll: B(Lztot,Kptot) alloc. error.`Enter`to exit.'
       stop
     endif!(iAlloc/=0)
   BtB=0._16
!  --------
   allocate(BtZ(Kptot),stat=iAlloc)         !Bt*Z array
     if(iAlloc /= 0) then
       pause 'TweakAllocateAll: BtZ(Kptot) allocation error. `Enter` to exit.'
       stop
     endif!(iAlloc/=0)
   BtZ=0._16
!  -------- Array allocations completed.
                                                                         return
End Subroutine TweakAllocateAll

