Upgrade of a server host from Opensuse Leap 15.1 to 15.2 – II – smartd configuration, 3ware controllers, mdadm RAID

During the upgrade of one of my server systems from Opensuse Leap 15.1 to 15.2 I noticed that the smartd-daemon did not start successfully. Logs revealed that the upgrade itself was NOT the cause of this failure; it is nevertheless a good occasion to briefly point out how one can configure the smartd-daemon to keep an overview over the evolution of disk properties on a Leap 15.X system – even when the disks are an integral part of one or more Raid-arrays. “Overview” is also a reason why I use “rsyslog” in addition to a pure systemd journal log to track the SMART data.

A problem with a 3ware based Raid array

The main cause of the smartd-failure on the upgraded server was an active line with the command “DEVICESCAN” in Opensuse’s version of the configuration file “/etc/smartd.conf”. As the command’s name indicates, the system is scanned for disks to monitor. Unfortunately, smartd did not find any SMART enabled devices on my server.

Feb 05 17:01:55 MySRV smartd[14056]: Drive: DEVICESCAN, implied '-a' Directive on line 32 of file /etc/smartd.conf
Feb 05 17:01:55 MySRV smartd[14056]: Configuration file /etc/smartd.conf was parsed, found DEVICESCAN, scanning devices
Feb 05 17:01:55 MySRV smartd[14056]: Unable to monitor any SMART enabled devices. Try debug (-d) option. Exiting...
...

Well, this finding was correct. All internal disks of the server were members of Raid arrays on elderly 3ware controllers. 3ware controllers are supported by the Linux kernel – in my case via a module “3w_9xxx”. During the installation an Opensuse Leap system the driver module was automatically added to the configurations of GRUB, the initramfs and the list of modules loaded by the eventually booted kernel. Upgrades have respected the configuration so far.

But with hardware RAID controllers the attached individual disks are NOT directly accessible to standard “smartd” directives or a standard plain “smartctl” command. You need to a special options. The HW-RAID-array is seen by the Linux system and smartd as a block device – without an inner structure. And the block device corresponding to the HW-controlled RAID-array has no SMART capabilities – which is, of course, very reasonable. Thus DEVICESCAN fails.

I wondered a bit why I had not seen the problem on the same system with Leap 15.1. Actually, this was my fault: Looking into a backup I saw in the last logs of the host during its time with Leap 15.1 that the problem had existed already before the upgrade – and was caused by an old, original “smartd.conf”-file from SuSE with which I had accidentally overwritten my own version of “smartd.conf”. Stupid me! But the really concerning thing was that I should have noticed the consequences in some logs. But, since I became a almost full-time employee I have been too lazy to work regularly with specific filter-options of systemd’s command “journalctl” – so I missed the relevant messages within the flood of messages in the journal. Time to reactivate rsyslog to get specific information into specific files – see below.

Anyway … It is reasonable that the failure of a standard DEVICESCAN should lead to a failure of the smartd service. The attentive admin thus gets a clear hint that his intention to watch disk health variables is not working as expected – even if the “smartd.conf” should contain other valid lines. The comments in “/etc/smartd.conf” clearly say:

# The word DEVICESCAN will cause any remaining lines in this configuration file to be ignored: 
# it tells smartd to scan for all ATA and SCSI devices. DEVICESCAN may be followed by any of the
# Directives listed below, which will be applied to all devices that are found.
#  
# Most users should comment out DEVICESCAN and explicitly list the devices that they wish to monitor.

I changed the line breaks a bit to stress the last statement. The obvious
and trivial solution in my case was to comment out the DEVICESCAN-line and activate some explicit statements for the individual disks in the 3ware array. To make things a bit more interesting, let us look at a different system MySYS which uses a combination of 3ware controlled hard-disk arrays and mdadm controlled SSD-arrays.

Off-topic:
The server host discussed above does not require a particular high disk throughput. Storage capacity in TB is more important. Actually, I would never buy a 3ware or LSI controller again for a semi-professional Linux machine. Too much money for a too low performance. With Linux you are certainly better off with mdadm and SW-RAID on a modern multi-core-processor. The general impact on an i7 or i9 on the overall performance is negligible in my experience. Especially, when you have a lot of fast RAM. In addition mdadm gives you much (!) more flexibility. You can, for example, build a set of multiple RAID-arrays working in parallel with one and the same stack of disks/SSDs by assigning different partitions to different distinct arrays – which may even be of a different type. But this is another story …

Entries for disks attached to 3ware controllers and SSDs in mdadm-Raid-arrays in the “/etc/smartd.conf” file

In the system MySRV I use a configuration

  • of multiple mdadm RAID-arrays consisting of assigned SSD-partitions of a stack of 4 SSDs
  • and 4 conventional hard-disks attached to a 3ware controller, forming a RAID 10-array.

In addition the host contains a further stand alone SSD.

For disks attached to a 3ware-controller we need a special form of the directives in the “smartd.conf”-file. Actually, some examples for such directives are given in a 3ware-related commentary section of “/etc/smartd.conf”. You may find something directly suitable for your purposes there.

What about devices which are members of mdadm-controlled SW-Raid arrays? Well, the nice thing is that in contrast to HW-controlled Raid-disks the SMART-daemon can access SSDs (or hard disks) on a mdadm-controlled SW-Raid-array directly! Without any special options! Actually, folks familiar with “mdadm” would have expected this for very basic reasons … But this is yet another story … We just accept the fact as one example ofthe power of Linux.

Thus, for my special system MySYS I just disabled the line with “Devicescan” and inserted the following lines into “/etc/smartd.conf”:

...
DEFAULT -d removable -s (S/../../7/03)
...
#DEVICESCAN
...
# Monitor SSDs contributing to different mdadm based SSD-Raid-arrays via (vertical) stacks of partitions 
/dev/sda -d sat -a -o off
/dev/sdb -d sat -a -o off
/dev/sdc -d sat -a -o off
/dev/sdd -d sat -a -o off

# Monitor additional stand alone SSD 
/dev/sde -d sat -a -o off

# Monitor HDs attached to 3ware controller with a Raid 10 configuration 
/dev/twa0 -d 3ware,0 -F samsung -a -o off
/dev/twa0 -d 3ware,1 -F samsung -a -o off
/dev/twa0 -d 3ware,2 -F samsung -a -o off
/dev/twa0 -d 3ware,3 -F samsung -a -o off
...

Disclaimer: Never copy the statements and later discussed smartctl-commands above without a thorough study of the literature and documentation on smart, smartctl and smartd. Check the preconditions of applying the directives and commands carefully with respect to your actual HW-/SW-configuration and the type of SSDs and hard disks used there! You have to find out about the correct settings for your system-configuration on your own.
Warning: The effects of the “-a” option and especially of the planned automatic disk tests initiated via the “DEFAULT”-statement may be different for SCSI or SATA-disks. Some options may even lead to data loss on older Samsung disks. I will not take any responsibility
for any problems caused by applying the settings displayed above to YOUR systems. Get advice of a professional regarding RAID-configurations and smartctl ahead of any experiments.

Having this said, the expert sees that I plan short self tests of the disks on every Sunday at 3 AM. Such tests will temporarily reduce the performance of the affected disks and as a consequence also of the various RAID-configurations the disk or any of its partitions might be a part of. At least on my system such a test will, however, not render the defined RAID-arrays unusable.

I could test this for a SATA-device and mdadm e.g. by using a suitable smartctl-command as

smartctl -t short /dev/sdc

and watching the continuing operation of a related RAID-array by

mdadm -D /dev/mdxxx

The corresponding smartcl-command for a 3ware RAID-disk would e.g. be

smartctl -t short /dev/twa0 -d 3ware,3 -F samsung.

As said above: Be careful with such tests on productive systems.

By the way: Only for the 3ware disks you may get an information about the starting of the test in a smartd related log-file (see below). If no errors are found, the logs will NOT show any special message.

At least on an Opensuse system the other settings and options for the various disks are discussed and explained in the commentary sections of the file:

#   -d TYPE Set the device type: ata, scsi, marvell, removable, 3ware,N, hpt,L/M/N
#   -o VAL  Enable/disable automatic offline tests (on/off)
#   -H      Monitor SMART Health Status, report if failed
#   -l TYPE Monitor SMART log.  Type is one of: error, selftest
#   -f      Monitor for failure of any 'Usage' Attributes
#   -p      Report changes in 'Prefailure' Normalized Attributes
#   -u      Report changes in 'Usage' Normalized Attributes
#   -t      Equivalent to -p and -u Directives
#   -a      Default: equivalent to -H -f -t -l error -l selftest -C 197 -U 198
#   -F TYPE Use firmware bug workaround. Type is one of: none, samsung

and

# Monitor 4 ATA disks connected to a 3ware 6/7/8000 controller which uses
# the 3w-xxxx driver. Start long self-tests Sundays between 1-2, 2-3, 3-4, 
# and 4-5 am.
# NOTE: starting with the Linux 2.6 kernel series, the /dev/sdX interface
# is DEPRECATED.  Use the /dev/tweN character device interface instead.
# For example /dev/twe0, /dev/twe1, and so on.
#/dev/sdc -d 3ware,0 -a -s L/../../7/01
#/dev/sdc -d 3ware,1 -a -s L/../../7/02
#/dev/sdc -d 3ware,2 -a -s L/../../7/03
#/dev/sdc -d 3ware,3 -a -s L/../../7/04

On my system (with one a 3ware controller and mdadm Raid-arrays) the smartd daemon started without any problems afterwards.

Configure rsyslog to write smartd-information into a distinct separate log-file

Old fashioned as I am, I use rsyslog aside systemd’s journal-logging. rsyslog can easily be configured to write log entries for selected daemons/processes into predefined separate log-files. Such files can directly be opened as they are plain text-files. Let us configure rsyslog for the smartd daemon:

  • Step 1: We create an empty file “/var/log/smartd.log” (as root) on the host in question by the command

    “touch /var/log(smartd.conf”

  • Step 2: Change rights by “chmod 640 /var/log/smartd.log”
  • Step 3: Add the following line to the file “/etc/rsyslog.conf” (close to the end)

    :programname, isequal, “smartd” /var/log/smartd.log

  • Step 4: Restart the rsyslog-daemon by “systemctl restart rsyslog”

You
surely have noticed that this is for local logging. But it is easy to adapt the settings to logging on remote systems.

What do smartd-logs look like?

Depending on your hardware you may get specific messages for your devices at the start of the smartd-daemon. Depending on the execution of scheduled self-tests you may get some messages about starting such test. After some time you should (hopefully) see successful messages regarding short or long self-tests – depending on your settings. In my case:

2021-02-11T14:29:43.511950+01:00 MySYS smartd[1723]: Device: /dev/sda [SAT], previous self-test completed without error
...
2021-02-11T14:29:44.026256+01:00 MySYS smartd[1723]: Device: /dev/twa0 [3ware_disk_00], previous self-test completed without error
...

In addition that you will get continuous information about changes of variables describing the health-status of your disks. In my case only the variation of the temperature due to changing airflow is shown:

...
2021-02-11T14:59:43.340341+01:00 MySYS smartd[1723]: Device: /dev/sda [SAT], SMART Usage Attribute: 190 Airflow_Temperature_Cel changed from 76 to 78
2021-02-11T14:59:43.350339+01:00 MySYS smartd[1723]: Device: /dev/sdb [SAT], SMART Usage Attribute: 190 Airflow_Temperature_Cel changed from 77 to 78
2021-02-11T14:59:43.353660+01:00 MySYS smartd[1723]: Device: /dev/sdc [SAT], SMART Usage Attribute: 190 Airflow_Temperature_Cel changed from 77 to 79
2021-02-11T14:59:43.357007+01:00 MySYS smartd[1723]: Device: /dev/sdd [SAT], SMART Usage Attribute: 190 Airflow_Temperature_Cel changed from 77 to 79
2021-02-11T14:59:43.360136+01:00 MySYS smartd[1723]: Device: /dev/sde [SAT], SMART Usage Attribute: 190 Airflow_Temperature_Cel changed from 77 to 78
...

Conclusion

If for some reason the “smartd”-daemon and the related service does not start on a Linux host with Opensuse Leap 15.2 it is easy to get it running again. Even, when your disk devices are members of RAID arrays. The configuration file “/etc/smartd.log” contains a lot of hints how to set up reasonable directives regarding various HW-RAID-controllers supported by the kernel. A really good thing is that we do not need any special commands or options regarding disk devices whose partitions are members of mdamd-controlled SW-RAID arrays.

Links

https://linux.die.net/man/5/smartd.conf
https://linux.die.net/man/8/smartd

https://www.linux.com/topic/desktop/monitor-sata-and-ssd-health-smart/
https://linuxconfig.org/introduction-to-the-systemd-journal
https://www.debugpoint.com/2020/12/systemd-journalctl/

Upgrade of a server host from Opensuse Leap 15.1 to 15.2 – I – a problem with Apache2/MPM and PHP7

The official support of Opensuse Leap 15.1 is about to end. Time to upgrade some servers to Opensuse Leap 15.2 in my home-office. While the basic upgrade of the OS was a smooth process I stumbled across some problems and challenges with a few services. I shall cover these problems step by step in this series – and hope that my remarks and hints will help some of my readers.

A central server system with network related services and with virtualized guests

One of the two central server systems in my LAN provides DHCP-, DNS-, OpenLDAP-services plus an internal HTTP-service (based on Apache2) for development work. Its list of services in addition comprises NFS, Samba and some SVN-services for multiple purposes, among others as backup space for other systems. The host also works as a CUPS server – offering some pre-defined queues for other Linux-systems in the LAN.

The system also serves as a KVM- and LXC-host for some virtualization guests: A KVM based LAMP/GIT-server for production data, a KVM based mail-server, a KVM based Kali-system plus some LXC-containers for databases. Some LVM volumes on this server are encrypted and all disk-space is based on Raid configurations – both with mdadm for different SSD-raids and with a 3ware controller for a standard hard disk based raid10. This server is also critical as it separates its internal virtual network structure from real LAN segments by a packet-filter. It must cover routing tasks.

I run this Linux host with Opensuse Leap – not SLES. So, it is a good test object for the cleanliness of Opensuse’s upgrade practices. Actually, I am always a bit nervous during upgrade. Fortunately, the HW is from the time of my freelance business – and thus safe and proven to use. It contains a passively cooled Nvidia card, for which the community repositories still provides kernel modules. The 3ware drivers also have worked reasonably well over the last 4 years. So, drivers at least should not pose any problems.

To make a long story short: The basic upgrade did not cost me much headache. I outline the upgrade procedure below – and then turn to the problems I have seen in its wake.

Backup and basic upgrade procedure

It should be clear that a server as the described one requires some regular backup routine for important data. I shall not elaborate on this topic here. Regarding the backup of the “/”-partition I use an approach described in the article
Upgrade workstation from Opensuse Leap 15.1 to Leap 15.2:

I.e. I copied the whole partition to a different, external hard-disk and afterwards in addition to a free partition on the server system itself – all with the help of the “dd”-command. You may do the same with other partitions which you want to re-create quickly in case of an emergency during upgrade.

Regarding the upgrade itself I cling to the procedure described by the Linux Kamarada; see the article mentioned above. (The required work with the Opensuse repositories is much simpler for my servers than for my workstations and laptops as I use much less repositories; the only repositories aside the Update repository for the present Leap version are SUSE’s “Security” and “Network” repositories plus the Nvidia community repository. I have made good experiences with this habit over the last 10 years.)

The upgrade worked almost perfectly. Some quick checks verified that my virtualized (and not yet upgraded) virtual guest systems could be started and worked as expected. The same was true for the slapd-, smb-, the wsdd- and the dns- (=named-) service.

But, there were also exceptions which required additional work:

  • Graphics Failure: The X-system did not start and the “graphical target” of systemd could not be reached directly after the system
    upgrade (expected).
  • Service failure: The Apache2 service did not start (unexpected).
  • Service failure: The smartd-service did not start (unexpected, but probably due to an overwritten configuration file).
  • Required Adjustment: The CUPS libraries had to be changed to Opensuse’s “Printing” repository to become consistent with other Linux hosts which use printer queues on the server (expected).
  • Required Adjustment: The rkhunter-configuration had to be adjusted (expected).
  • Required Adjustment: A specific KDE-settings had to be adjusted to allow for screen-saving functions – even for SDDM screens without started KDE-sessions (unexpected).
  • Required Adjustment: “virt-viewer” access via “qemu+ssh” to some clients posed grave problems (unexpected). These problems may have existed already before, but are worth a description.

I cover the first 2 points in this article. Solving the graphics problem is a short story. The explanation and solution of the problem with the Apache service, however, requires some space. But I hope the information provided in the related text is interesting for some readers who are no experts in Opensuse’s handling of the Apache service.

Nvidia – X-Windows

You may ask: Why X11 on a server at all? Good question; well, as a semi-professional in server-administration, I think that one can do some things faster than without. Normally, I access the server remotely via ssh or “ssh -X”. But during my years with SLES and Opensuse I experienced some situations during which I was happy to be able to access the server directly on a graphical console; certainly not a popular view among profs. Anyway – the Nvidia/X-Window topic is a standard problem and easy to solve:

I just had to reinstall the required drivers from the Nvidia repository (via the CLI-based YaST), perform a “mkinitrd” afterwards and restart the system.

However, I nevertheless recommend not to boot the server system into the graphical systemd target, if not required. The relevant command to set the default target is

systemctl set-default multi-user.target

See the related link at the end of this post.

Off topic: If you prefer KDE as a graphical desktop environment do not even try to switch to Wayland on Leap 15.2. In my experience this is only a frustrating endeavor – as we still do not have Plasma 5.20 available but Plasma 5.18 (see a related link at the end).

Apache does not start due to a problem with “mod_php7” and MPM

Let us turn to something more interesting: An Apache problem with the current PHP-module and MPM-settings.

I admit, I am a bit old-fashioned. Most developers and admins, today, probably use Nginx and in addition Fast-CGI in combination with PHP. I, instead, cling to an old Apache installation for LAMP-services on the server named above. One reason is that I need some compatibility with things I did from 2013 and 2016 – aside of changes of PHP over the years.

So, whilst upgrading Opensuse my Apache 2.4 and PHP installations get upgraded, too.

When one activates the php-Apache-module by “a2enmod php” the presently installed version of the php-Apache-module gets loaded – on Leap 15.2 this is provided by mod_php7.so corresponding to PHP 7.4.6. (You can verify this with the command “php -v”). Using “a2enmod php7” would give you no different result. (Opensuse, unfortunately, does not support the parallel installation of multiple, different PHP-versions. See the critics in the conclusion-section.)

All upgrades from Leap 42.2 to 42.3, 15.0 and eventually to Leap 15.1 went well; but this time the Apache service did not start. Trying “rcapache2 start” followed by “rcapache2 status” (or equivalent systemctl status apache2.service) gave me an error message:

….
[php7:crit] [pid 16004:tid 139738320386304] Apache is running a threaded MPM, but your PHP Module is not compiled to be threadsafe. You need to recompile PHP.
Jan 26 19:15:39 MySRV start_apache2[16004]: AH00013: Pre-configuration failed
Jan 26 19:15:39 MySrv systemd[1]: apache2.service: Main process exited, code=exited, status=1/FAILURE
….

I remembered vaguely that I had experimented with MPM and prefork settings on Leap 15.1, but recompile PHP? With what options? A strange error message!

Confusion when comparing to a seemingly similar Apache2 installation on a different system

I cross-checked with a similarly configured Apache2 installation (with PHP7, but different virtual domains) on one of my my Leap 15.2 workstations; let us call it “MyWS” and the upgraded server “MySRV”. Confusingly, Apache2 started on MyWS without any complaints.

You may suspect now that I may have used different Apache2 configuration directives in “/etc/sysconfig/apache2” or the files in folder “/etc/apache2/” and sub-folders – without being aware of it. But no: I checked and eliminated even small, irrelevant differences for tests. The problem on the server MySRV remained.

When I looked at the Apache2-related processes on the workstation MyWS with “ps aux | grep https” I saw that Apache2 had been started in a pre-forked way. As I was used to …. I got multiple forked processes like

....
wwwrun    7043  0.0  0.0 236180 10488 ?        S    10:41   0:00 /usr/sbin/httpd-prefork -DSYSCONFIG -DSSL -C PidFile /var/run/httpd.pid -C Include .... 
.....
wwwrun    7045  0.0  0.0 236180 10488 ?        S    10:41   0:00 /usr/sbin/httpd-prefork -DSYSCONFIG -DSSL -C PidFile /var/run/httpd.pid -C Include ... 
... 

The interesting thing here is the indicated program “/usr/sbin/httpd-prefork”. Note that it is NOT “usr/bin/httpd”; I shall come back to this point later. The process information was of course consistent with startup information provided by “rcapache2 status” or equivalently by “systemctl status apache2”.

It took me a while to find out what was going on and what I had to do. Whilst searching the Internet and experimenting with Apache I passed some interesting information about details of the Apache-installation on Opensuse which may be interesting for others, too.

MPM methods on Apache 2.4

The error message indicated that the whole problem had to do with MPMs – “Multi-Processing-Modules”, i.e. the ability of Apache to handle a multitude of requests concurrently. There are 3 MPM-methods and respective pre-compiled static modules or loadable DSO modules available with Apache 2.4. See e.g.:

digitalocean tutorials how-to-configure-apache-http-with-mpm-event-and-php-fpm-on-ubuntu-18-04-de,
http://httpd.apache.org /docs /current /en /mpm.html
and links therein.

I quote:

  • Pre-fork: A new process is created for each incoming connection reaching the server. Each process is isolated from the others, so no memory is shared between them, even if they are performing identical calls at some point in their execution. This is a safe way to run applications linked to libraries that do not support threading—typically older applications or libraries.
  • Worker: A parent process is responsible for launching a pool of child processes, some of which are listening for new incoming
    connections, and others are serving the requested content. Each process is threaded (a single thread can handle one connection) so one process can handle several requests concurrently. This method of treating connections encourages better resource utilization, while still maintaining stability. This is a result of the pool of available processes, which often has free available threads ready to immediately serve new connections.
  • Event: Based on worker, this MPM goes one step further by optimizing how the parent process schedules tasks to the child processes and the threads associated to those. A connection stays open for 5 seconds by default and closes if no new event happens; this is the keep-alive directive default value, which retains the thread associated to it. The Event MPM enables the process to manage threads so that some threads are free to handle new incoming connections while others are kept bound to the live connections. Allowing re-distribution of assigned tasks to threads will make for better resource utilization and performance.

In a doc at apache.org about performance scaling we furthermore read:

Httpd Configuration
The Apache 2.2 httpd is by default a pre-forking web server. When the server starts, the parent process spawns a number of child processes that do the actual work of servicing requests. But Apache httpd 2.0 introduced the concept of the Multi-Processing Module (MPM). Developers can write MPMs to suit the process- or threading- architecture of their specific operating system. Apache 2 comes with special MPMs for Windows, OS/2, Netware and BeOS. On unix-like platforms, the two most popular MPMs are Prefork and Worker. The Prefork MPM offers the same pre-forking process model that Apache 1.3 uses. The Worker MPM runs a smaller number of child processes, and spawns multiple request handling threads within each child process. In 2.4 MPMs are no longer hard-wired. They too can be exchanged via LoadModule. The default MPM in 2.4 is the event MPM.
….
The maximum number of workers, be they pre-forked child processes or threads within a process, is an indication of how many requests your server can manage concurrently. It is merely a rough estimate because the kernel can queue connection attempts for your web server. When your site becomes busy and the maximum number of workers is running, the machine doesn’t hit a hard limit beyond which clients will be denied access. However, once requests start backing up, system performance is likely to degrade.
Finally, if the httpd server in question is not executing any third-party code, via mod_php, mod_perl or similar, we recommend the use of mpm_event. This MPM is ideal for situations where httpd serves as a thin layer between clients and backend servers doing the real job, such as a proxy or cache.
… 
Selecting your MPM
The prime reason for selecting a threaded MPM is that threads consume fewer system resources than processes, and it takes less effort for the system to switch between threads. This is more true for some operating systems than for others. On systems like Solaris and AIX, manipulating processes is relatively expensive in terms of system resources. On these systems, running a threaded MPM makes sense.
On Linux, the threading implementation actually uses one process for each thread. Linux processes are relatively lightweight, but it means that a threaded MPM offers less of a performance advantage than in other environments.
 
Running a threaded MPM can cause stability problems in some situations For instance, should a child process of a preforked MPM crash, at most one client connection is affected. However, if a threaded child crashes, all the threads in that process disappear, which means all the clients currently being served by that process will see their connection aborted. Additionally, there may be so-called “thread-safety” issues, especially with third-party libraries. In threaded applications, threads may access the same variables indiscriminately, not knowing whether a variable may have been changed by another thread.
 
This has been a sore point within the PHP community. The PHP processor heavily relies on third-party libraries and cannot guarantee that all of these are thread-safe. The good news is that if you are running Apache on Linux, you can run PHP in the preforked MPM without fear of losing too much performance relative to the threaded option.

Ok, what do we learn from this? php_mod is not thread-safe. But regarding performance, using the “pre-fork” MPM method is not really a bad choice on a Linux system. Obviously, this method was used on the workstation. An interesting question, which I come back to later, is whether Opensuse offers dynamically loadable or static MPM modules via its repositories.

Information about the MPM used by your installation?

One (naively) expects that the Apache2 daemon binary should be callable by “/usr/sbin/httpd” on an Opensuse system. Actually, the command “which httpd” gives you “/usr/sbin/httpd” as an answer.

On the Internet you find hints that one can check the compiled and available MPMs by either “httpd -V” and/or “httpd -l“. On the upgraded server MySRV I got:

MySRV:~ # httpd -V
Server version: Apache/2.4.43 (Linux/SUSE)
Server built:   2020-11-17 11:05:32.000000000 +0000
Server's Module Magic Number: 20120211:92
Server loaded:  APR 1.6.3, APR-UTIL 1.6.1
Compiled using: APR 1.6.3, APR-UTIL 1.6.1
Architecture:   64-bit
Server MPM:     worker
  threaded:     yes (fixed thread count)
    forked:     yes (variable process count)
Server compiled with....
 -D APR_HAS_SENDFILE
 -D APR_HAS_MMAP
....

and

MySRV:~ # httpd -l
Compiled in modules:
  core.c
  mod_so.c
  http_core.c
  worker.c
  mod_unixd.c
  mod_systemd.c

These commands seemingly worked on my server even if the Apache2-service itself could not be started (e.g. by “systemctl start apache2”). The man-pages for “httpd” say that the “-l”-option lists modules which were compiled into the server, i.e. static modules.

What did I get on the workstation?

MyWS:~ # httpd -V
Server version: Apache/2.4.43 (Linux/SUSE)
Server built:   2021-01-27 09:11:10.000000000 +0000
Server's Module Magic Number: 20120211:92
Server loaded:  APR 1.6.3, APR-UTIL 1.6.1
Compiled using: APR 1.6.3, APR-UTIL 1.6.1
Architecture:   64-bit
Server MPM:     prefork
  threaded:     no
    forked:     yes (variable process count)

and

MyWS:~ # httpd -l
Compiled in modules:
  core.c
  mod_so.c
  http_core.c
  prefork.c
  mod_unixd.c
  mod_systemd.c

Regarding the actually running processes after Apache-startup on the workstation this result is not too surprising. However, if the listed modules are static then the “httpd” binary on the server must be another one than on the workstation. Why? And how did this come about ?

Another question is: How much can we rely on the information given by e.g. “httpd -l” in the sense that this information reflects facts about the running Apache2 service? We shall see that the answer depends …

Removing the PHP-module on both installations

To analyze further I first aligned the list of modules on both systems via the corresponding setting in “/etc/sysconfig/apache2” to:

APACHE_MODULES=”actions alias auth_basic authn_file authz_host authz_groupfile authz_core authz_user autoindex cgi dir env expires include log_config mime rewrite negotiation
setenvif ssl socache_shmcb userdir reqtimeout php7 authn_core version”

and then removed the “php7” module.

(Instead of editing the file “/etc/sysconfig/apache2” directly you could also use a2enmod or a2dismod, respectively, for changing the module configuration. But I recommend to keep some commented alternatives for in the file for fast changes between repeated tests of different constellations).

In addition I used the seemingly standard Opensuse setting of

APACHE_MPM=””

in “/etc/sysconfig/apache2” for both installations.

Well, as expected no significant changes occurred on the workstation. Apache2 started successfully again (without PHP support); multiple prefork-processes were launched. On the server, however, I now got a successful start, too; “rcapache2 status” showed me a bunch of related processes afterwards:

MySRV:~ # rcapache2 status
● apache2.service - The Apache Webserver
   Loaded: loaded (/usr/lib/systemd/system/apache2.service; enabled; vendor preset: >
   Active: active (running) since Sun 2021-01-31 11:45:26 CET; 23s ago
  Process: 8555 ExecStop=/usr/sbin/start_apache2 -DSYSTEMD -DFOREGROUND -k graceful->
 Main PID: 20667 (httpd-worker)
   Status: "Processing requests..."
    Tasks: 83
   CGroup: /system.slice/apache2.service
           ├─20667 /usr/sbin/httpd-worker -DSYSCONFIG -DSSL -C PidFile /var/run/httpd.pid ....
           ├─20675 /usr/sbin/httpd-worker -DSYSCONFIG -DSSL -C PidFile /var/run/httpd.pid ....
           ├─20676 /usr/sbin/httpd-worker -DSYSCONFIG -DSSL -C PidFile /var/run/httpd.pid ....
           ├─20677 /usr/sbin/httpd-worker -DSYSCONFIG -DSSL -C PidFile /var/run/httpd.pid ....
           └─20678 /usr/sbin/httpd-worker -DSYSCONFIG -DSSL -C PidFile /var/run/httpd.pid ....

Jan 31 11:45:26 MySRV systemd[1]: Starting The Apache Webserver...
Jan 31 11:45:26 MyServ start_apache2[20667]: AH00548: NameVirtualHost has no effect and>
Jan 31 11:45:26 MySRV systemd[1]: Started The Apache Webserver.

But you see the difference?
=> The programs started were of the type “/usr/bin/http-worker” !

Looking at the startup scripts for the Apache2 service

To dig a bit deeper I looked at the startup scripts used on Opensuse for the Apache2 service. One reason to do so was a comment in Opensuse’s “/etc/sysconfig/apache2” regarding the variable APACHE_MPM:

# MPM (multi-processing module) to use.
#
# Needed to determine with which MPM apache will run, as well as
# against which header files modules will be built. 
#
# If not set, the system will simply pick one of the installed MPMs.
#
# The implementation of the logic is in /usr/share/apache2/find_mpm,
# a script which can be used standalone as well if needed.

This indicates that some special Opensuse logic is used during startup. But the comment is also misleading: A script “find-mpm” can only be found in the folder “/usr/share/apache2/deprecated-scripts” on a Leap 15.1 and 15.2 systems – and the “script” itself is not used any more, as we shall see in a minute.

So, where do we find the startup-script for Apache? We do not have to look into the systemd configuration for the Apache2 service; the output of the command “rcapache2 status” already revealed it as
/usr/sbin/start_apache2“.

In it’s script-code we find a line for the execution of a helper script:

 . /usr/share/apache2/script-helpers

and later a call to a function “find_mpm”, followed by some logic for the name of a binary to be saved in the variable “apache_bin”:

#
# figure out correct apache2 binary (/usr/sbin/httpd-prefork,
# /usr/sbin/httpd-worker, etc.) and serverflags
#
find_mpm
if [ -n "$HTTPD_MPM" ]; then
r
    apache_bin="$HTTPD_SBIN_BASE-$HTTPD_MPM"
    if ! [ -x $apache_bin ]; then
        echo >&2 "$apache_bin-$APACHE_MPM is not a valid httpd binary."
        echo >&2 "Check your APACHE_MPM setting in /etc/sysconfig/apache2."
        exit 1
    fi
else
  echo >&2 "${warn}No Apache binary found. No MPM package installed? $norm"
  echo >&2 "Hint: install the apache2-prefork package, and try again."
  exit 1
fi

Obviously, we will find the function “find_mpm” in “/usr/share/apache2/script-helpers”. A look into this file shows that this really is the case. The outcome of this function influences the startup of Apache2; obviously, a different Apache variant is started for different (valid) values of the script-variable “HTTPD_MPM“:

apache_bin="$HTTPD_SBIN_BASE-$HTTPD_MPM"

($HTTPD_SBIN_BASE resolves to “/usr/sbin/httpd” – this is set at the beginning of “/usr/share/apache2/script-helpers”.)

Now, you can easily build your own test-script “mpm_test” including this function and print its output:

mpm_test-script:

#!/bin/bash

HTTPD_SBIN_BASE="/usr/sbin/httpd"
#
# loads sysconfig variables into environment
# return value in: APACHE_*
function load_sysconfig
{
  [ -n "$sysconfig_loaded" ] && return
  [ ! -f "$SYSCONFIG_FILE" ] && return

  . $SYSCONFIG_FILE

  export ${!APACHE_*} sysconfig_loaded=true
}
#
# finds prefered multiprocessing module
# return value in: HTTPD_MPM 
function find_mpm
{
  # load sysconfig variables if they weren't yet;
  # this has no effect when find_mpm is not called
  # from start_apache2
  load_sysconfig

  # try to read from sysconfig's APACHE_MPM
  HTTPD_MPM="$APACHE_MPM"
  # if empty, then choose one from installed
  if [ -z "$HTTPD_MPM" ]; then
      installed_mpms=""
      for i in $HTTPD_SBIN_BASE-*; do
          test -f $i || continue
          i=$(basename $i)
          i=${i#*-}
          installed_mpms="$installed_mpms $i"
      done
      # hardcoded preference here:
      for mpm in event worker prefork; do
        if [[ $installed_mpms =~ "$mpm" ]]; then
          HTTPD_MPM=$mpm
          break
        fi
      done
  fi

  # in case no 
  export HTTPD_MPM
}

# ***********************
# Action
# ***********************
find_mpm
echo $HTTPD_MPM

This script produced the following output on my server MySRV (with a running Apache2 without PHP support):

worker

On my workstation (and a running Apache2 with and without PHP support) I got:

prefork

Available and different Apache binaries

Analyzing the logic of our file, which is the same as in the beginning of the startup-script, we see that
“apache_bin” can become either “/usr/sbin/httpd-prefork” or “/usr/sbin/httpd-worker”. It depends on the value of the variable “APACHE_MPM” in “/etc/sysconfig/apache2”:

  • If the variable “APACHE_MPM” is not empty its value is used.
  • But if it is empty the result depends on the existence of one or multiple of the 3 files “httpd-event”, “httpd-worker” or “httpd-prefork” in the folder “/usr/sbin”.

You also see that the logic, which the guys from Opensuse implemented, chooses a MPM-specific variant of Apache2 according to the following fixed priority order:

event > worker > prefork    !

The “prefork” variant has the least priority – and is obviously only chosen if neither the file “httpd-worker” nor the file “httpd-event” exist – and if “APACHE_MPM” is not set otherwise.

A quick check showed that binaries “/usr/sbin/httpd-worker”, “/usr/sbin/httpd-event” did not exist on my workstation MyWS; there, only “/usr/sbin/httpd-prefork” was available. But on the server MySRV I actually had two variants available “/usr/sbin/httpd-prefork” and “/usr/sbin/httpd-worker”.

This fact together with the “find_mpm”-logic explained the different start-up results for Apache2 on my systems MySRV and MyWS !

But not the original cause and why Apache had worked on MySRV with Leap 15.1 before ….

Differences in the package installation?!

The existence of different binaries on MyWS and the server MySRV could only be due to a package difference, which must have been there already before the upgrade to Leap 15.2. A check in my logs and also a look into the package constellation with YaST revealed that on the Leap 15.1 version of the server-system MySRV I actually had installed the RPM “apache2-worker” at some point in time! For some experiments which I had totally forgotten. But not on the workstation. So, I had caused the trouble somehow by myself – without anticipating any consequences for a coming upgrade to Leap 15.2.

But wait a second – this does not explain why I did not run into trouble with Apache on Leap 15.1. Despite the installed package “apache_worker” and an empty “APACHE_MPM” in the sysconfig file….. Why did Apache2 start on Leap 15.1 – as my old log-files show? With “httpd-prefork” processes, by the way …

I had to suspect a difference in the “find_mpm” logic – if I really had left “Apache_MPM” empty on the server. The latter was the case as I saw in my backups. And indeed – on Leap 15.1 the logic of the function “find_mpm” in “/usr/share/apache2/script-helpers” was different:

The old “find_mpm” from Leap 15.1

function find_mpm
{
  [ -n "$mpm_found" ] && return

  # load sysconfig variables if they weren't yet
  load_sysconfig

  # try to read from sysconfig's APACHE_MPM
  HTTPD_MPM=$APACHE_MPM
  # if empty, then choose one from installed
  if [ -z "$HTTPD_MPM" ]; then
      # guess
      for i in $HTTPD_SBIN_BASE-*; do
          test -f $i || continue
          i=$(basename $i)
          i=${i#*-}
          installed_mpms=(${installed_mpms[*]} $i)
      done
      if [ -z "${installed_mpms[*]}" ]; then
          HTTPD_MPM=""
          return
      elif [ ${#installed_mpms[*]} = 1 ]; then
          HTTPD_MPM=${installed_mpms[*]}
      else
          case ${installed_mpms[*]} in
              *prefork*)      HTTPD_MPM=prefork;;
              *worker*)       HTTPD_MPM=worker;;
              *event*)        HTTPD_MPM=event;;
          esac
      fi
  fi

  export HTTPD_MPM mpm_found=true
}

A test shows that the “case”-part of this function gives us an answer “prefork” on a system with both the “prefork” and the “worker” variants of Apache installed!

This makes me wonder why, for heavens sake, the Opensuse guys did not insert a section into their Leap 15.2 Release Notes telling us something about this change!

The solution to my Apache/MPM/PHP problem

After all the information gathering above the very trivial solution to my start-up problem was as that I just had to set the variable “APACHE_MPM” in “/etc/sysconfig/apache2” to

APACHE_MPM=”prefork”

when I wanted a successful start of an Apache instance with PHP support.

What do the http-* files in “/usr/sbin/” look like for multiple MPM-installations?

What happens if we install multiple MPM-variants of Apache? From some parts of the above text you have already concluded that we really must speak of “variants” as the MPM-module is compiled in statically by the Opensuse guys.

On the server MySRV with the worker-MPM active (no PHP loaded) we can see this by using the command ”
apachectl -M”. See the man pages. The result is:

MySRV: ~# apachectl -M
Loaded Modules:
 core_module (static)
 so_module (static)
 http_module (static)
 mpm_worker_module (static)
 unixd_module (static)
 systemd_module (static)
 actions_module (shared)
 alias_module (shared)
 auth_basic_module (shared)
 authn_file_module (shared)
 authz_host_module (shared)
...
...
 

So, we expect multiple, relatively big variants of the “http”-binary under “/usr/sbin” whenever we install the packages “apache2-worker” and “apache2-event”. Which I did for fun and step-wise on the workstation MyWS. And indeed – after the installation of the RPM “apache2-worker” we see:

MyWS:~ # la /usr/sbin/ | grep httpd
lrwxrwxrwx  1 root root           22 Jan 31 16:27 httpd -> /usr/sbin/httpd-worker
-rwxr-xr-x  1 root root       616816 Jan 27 16:17 httpd-prefork
-rwxr-xr-x  1 root root       629200 Jan 27 16:17 httpd-worker
lrwxrwxrwx  1 root root           23 Jan 30 11:02 httpd2 -> /usr/sbin/httpd-prefork
lrwxrwxrwx  1 root root           13 Jan 27 16:17 httpd2-prefork -> httpd-prefork
lrwxrwxrwx  1 root root           12 Jan 27 16:17 httpd2-worker -> httpd-worker

and after an additional installation of the “event” version of Apache2 we get:

MyWS:~ # la /usr/sbin/ | grep httpd
lrwxrwxrwx  1 root root           21 Jan 31 16:29 httpd -> /usr/sbin/httpd-event
-rwxr-xr-x  1 root root       641552 Jan 27 16:17 httpd-event
-rwxr-xr-x  1 root root       616816 Jan 27 16:17 httpd-prefork
-rwxr-xr-x  1 root root       629200 Jan 27 16:17 httpd-worker
lrwxrwxrwx  1 root root           23 Jan 30 11:02 httpd2 -> /usr/sbin/httpd-prefork
lrwxrwxrwx  1 root root           11 Jan 27 16:17 httpd2-event -> httpd-event
lrwxrwxrwx  1 root root           13 Jan 27 16:17 httpd2-prefork -> httpd-prefork
lrwxrwxrwx  1 root root           12 Jan 27 16:17 httpd2-worker -> httpd-worker

At the end of both installations a post-installation script is running, which creates the binaries, the “httpd”-link and the various “httpd2”-links.

These files/links do not change when we restart Apache2 with a different MPM-setting in “/etc/sysconfig/apache2”.

Can we trust in “httpd -l”? What about “apachectrl -M”?

From the above information we must conclude that we cannot trust the information provided by “httpd -l” or “httpd -V” to give us reliable information for a running Apache2 service on Opensuse. It only shows us, what the most modern MPM among the available ones on the system is. Something similar is to say about “httpd2 -l”; it shows us the other end of the spectrum. Only, if we have just one MPM variant of Apache2 installed all information gets identical and consistent.

If you have a running Apache2 service, a really reliable information regarding the compiled in MPM module is given by “Apachectrl -M” or indirectly by “ps aux | grep httpd”, which reveals the type of the httpd-binary variant started.

Conclusion

Apache2 is provided on Opensuse with static compiled in MPM modules. If you want to choose between different MPMs you have to install related packages; they provide different variants of the httpd-binary which you find afterwards in the directory “/usr/sbin”.

The present Leap 15.2-specific startup-logic for the Apache-service is based an Opensuse-specific script-function “find_mpm” (witin the file “/usr/share/apache2/script-helpers”). According to its logic you MUST specify the MPM-variant via setting a variable in “/etc/sysconfig/apache2” if you want to or must break the fixed preference order “even > worker > prefork”.

Opensuse changed the logic of its MPM selection in the function “find-mpm” between Leap 15.1 and Leap 15.2. This change may, unfortunately, block a successful start of an upgraded Apache2-
service with PHP-support. (A bit of concise information on this point in the Release Notes would have been helpful.)

In case you use “mod_php” for PHP support by Apache2 and thus add “php7” to the list of Apache-modules to be loaded at start-up, you MUST set the “APACHE_MPM”-variable to “prefork” – at least if you follow the standard Opensuse concepts for the Apache service. The enforced choice is reasonable as PHP itself is not thread-safe.

Criticism, Opensuse Apache/PHP deficits and alternatives

As I criticize the handling of Apache2 on Opensuse I also want to add that I miss the flexibility other distributions offer regarding the parallel installation of multiple PHP-versions and a switch between them on the Apache-server. At least with a fast-cgi implementation.

I furthermore do not understand why Opensuse restricts the PHP version in their Update repository to just one and in the PHP-specific repository to the latest 2 versions. And even for these two PHP versions a developer is forced to install two virtual Opensuse machines with different Apache/PHP constellations, if he does not want to change the installed packages all the time.

The somewhat unflexible provision of Apache in combination with PHP and the MPM-logic may make you think about alternatives. There are various situations:

  • If you look beyond Apache I would recommend to setup Nginx.
  • If you look at a way to combine PHP with threaded MPM variants as “worker” or “event, try an installation based on FastCGI and PHP-FPM instead of using php_mod. But consider security aspects then.
  • If you look for support of multiple PHP-versions with Nginx and/or Apache install Debian or Ubuntu in a virtual machine (on a Opensuse host 🙂 ).

For those who like to experiment I have added some links in the last section of this post.

Next post

In the next post I shall cover the startup-problem of the “smartd”-service. This will be a much shorter story.

Links

Links related to the Apache MPM-topic
http://httpd.apache.org /docs /current /en /mpm.html
https://documentation.suse.com/ sles/15-SP2/ html/ SLES-all/cha-apache2.html
http://httpd.apache.org/ docs/ trunk/misc/ perf-scaling.html

Check available MPM-method on Apache
https://www.binarytides.com/check-which-mpm-multi-processing-module-apache-is-running/

Links regarding an alternative FastCGI installation with PHP-FPM
https://en.opensuse.org/ SDB:Apache_FastCGI_and_PHP-FPM_configuration
https://www.p4tchwork.de/apache-php-system-resource-usage/
https://blog.wappler.systems/opensuse-leap-15-apache-2-php-fpm-http2/

Flexible Apache and PHP installations on Debian derivatives
https://www.digitalocean.com/ community/ tutorials/ how-to-configure-apache-http-with-mpm-event-and-php-fpm-on-ubuntu-18-04-de
https://www.epidemiology.tech/apache-php7-4-wordpress-setup/
https://tecadmin.net/setup-apache-php-fpm-ubuntu-20-04/
https://www.interserver.net/ tips/kb/ change-php-version-apache-ubuntu/

Links regarding Nginx
https://techviewleo.com/install-php-on-opensuse-nginx-apache/
https://www.tecmint.com/install-lemp-nginx-php-mariadb-phpmyadmin-in-opensuse/
https://websiteforstudents.com/install-php-7-4-php-7-4-fpm-on-ubuntu-18-04-with-apache2-nginx/

HTTP 2 installation on Opensuse
https://en.opensuse.org/ SDB: Apache_HTTP_2_configuration

Set default target on Opensuse
https://forums.opensuse.org/ showthread.php/ 526792-How-to-change-default-runlevel

Plasma, Nvidia and Wayland
https://community.kde.org/ Plasma/ Wayland/ Nvidia

Deep Dreams of a CNN trained on MNIST data – III – catching dream patterns at smaller length scales

In the first two articles of this series on Convolutional Neural Networks [CNNs] and “Deep Dream” pictures

Deep Dreams of a CNN trained on MNIST data – II – some code for pattern carving
Deep Dreams of a CNN trained on MNIST data – I – a first approach based on one selected map of a convolutional layer

I introduced the concept of CNN-based “pattern carving” and a related algorithm; it consists of a small number of iterations over of a sequence of image manipulation steps (downscaling to a resolution the CNN can handle, detection and amplification of a map triggering pattern by a CNN, upscaling to the original size and normalization). The included pattern detection and amplification algorithm is nothing else than the OIP-algorithm, which I discussed in another article series on CNNs in this blog. The OIP-algorithm is able to create a pattern, which triggers a chosen CNN map optimally, out of a chaotic fluctuation background. The difference is that we apply such an OIP-algorithms now to a structured input image – the basis for a “dream”. The “carving algorithm” is just a simplifying variation of more advanced Deep Dream algorithms; it can be combined with simple CNNs trained on gray and low-resolution images.

In the last article I provided some code which creates dream like pictures with the help of carving by a CNN trained on MNIST data. Nice, but … In the original theory of “Deep Dreaming” people from Google and others applied their algorithms to a cascade of so called “octaves“: Octaves represent the structures within an image at different levels of resolution. Unfortunately, at first sight, such an approach seems to be beyond the capabilities of our MNIST CNN, because the latter works on a fixed and very coarse resolution scale of 28×28 pixels, only.

As we cannot work on different scales of resolution directly: Can we handle pattern detection and amplification on different length scales in a different way? Can we somehow extend the carving process to smaller length scales and a possible detection of smaller map-triggering patterns within the input image?
The answer in short is: Yes, we can.
We “just” have to replace working with “octaves” by looking at sub-segments of the images – and apply our “carving” algorithm with up and down-scaling to these sub-segments as well as to the full image during an iteration process. We shall test the effects of such an approach in this article for a single selected sub-area of the input image, before we apply it more thoroughly to different input images in further articles.

A first step towards Deep Dreams “dreaming” detail structures of a chosen image …

We again use the image of a bunch of roses, which we worked with in the last articles, as a starting point for the creation of a Deep Dream picture. The easiest way to define sub-areas of our (quadratic) input image is to divide it into 4 (quadratic) adjacent sub-segments or sub-images, each with half of the side length of the original image. We thus get two rows, each with two neighboring sub-areas of the size 280×280 px. Then we could apply the carving algorithm to one selected or to all of the sub-segments. But how do we combine such an approach with an overall treatment of the full image? The answer is that we have to work in parallel on both length-scales involved. We can do this within each cycle of down- and up-scaling according to the following scheme:

  • Step 1: Pick the input image IMG at the full working size (here 560 x 560 px) – and reshape it into a tensor suitable for
    Tensorflow 2 [TF2] functions.
  • Step 2: Down- and upscale the input image “IMG” (560×560 => 28×28 => 560×560) – with information loss
  • Step 3: Calculate the difference between the re-up-scaled image to the original image as a tensor for later detail corrections.
  • Step 4: Subdivide IMG into 4 quadrants (each with a size of of 280×280 px). Save the quadrants as sub-images S_IMG_1, …. S_IMG_4 – and create respective tensors.
  • Step 5: For all sub-images: Down- and up-scale the sub-images S_IMG_m (280×280 => 28×28 => 280×280) – with information loss.
  • Step 6: For all sub-images: Determine the differences between the re-upscaled sub-images and the original sub-images as tensors for later detail corrections.
  • Step 7: Loop A [around 4 to 6 iterations]
    • Step LA-1: Pick the down-scaled image IMG_d (of size 28×28 px) as the present input image for the CNN and the OIP-analysis
    • Step LA-2: Apply the OIP-algorithm for N epochs (20 < N < 40) to the downscaled image IMG_d (28x28)
    • Step LA-3: Upscale the resulting changed version of image IMG_d to the original IMG- size by bicubic interpolation => IMG_u
    • Step LA-4: Add the (constant) correction tensor for details to IMG_u.
    • (Step LA-5: Loop B over Sub-Images)
      • Step LB-1: Cut out the area corresponding to sub-image S_IMG_m from the changed full image IMG_u. Use it as Sub_IMG_m during the following steps.
      • Step LB-2: Downscale the present sub-image Sub_IMG_m to the size suitable for the CNN (here: 28×28) by bicubic interpolation => Sub_IMG_M_d.
      • Step LB-3: Apply the OIP-algorithm for N epochs (20 < N < 40) to the downscaled sub-image (28x28) Sub_IMG_m_d
      • Step LB-4: Upscale the changed Sub_IMG_m_d to the original size of Sub_Img_m by bicubic interpolation => Sub_IMG_m_u
      • Step LB-5: Add the (constant) correction tensor to the tensor for the upscaled sub-image Sub_IMG_m_u
      • Step LB-6: Replace the sub-image region of the upscaled full image IMG_u with the (corrected) sub-image Sub_IMG_m_u
      • Step LB-7: Downscale the new full image IMG_u again (to 28×28 => IMG_d) and use the resulting IMG_d for the next iteration of Loop A

The Loop B has been placed in brackets because we are going to apply the suggested technique only to a single one of the 4 sub-image quadrants in this article.

Code fragments for the sub-image preparation

The following code for a Jupyter cell prepares the quadrants of the original image IMG – here of a bunch of roses. The code is straight forward and easy to understand.

# ****************************
# Work with sub-Images  
# ************************

%matplotlib inline

import PIL
from PIL import Image, ImageOps

fig_size = plt.rcParams["figure.figsize"]
fig_size[0] = 16
fig_size[1] = 20
fig_3 = plt.figure(
3)
ax3_1_1 = fig_3.add_subplot(541)
ax3_1_2 = fig_3.add_subplot(542)
ax3_1_3 = fig_3.add_subplot(543)
ax3_1_4 = fig_3.add_subplot(544)
ax3_2_1 = fig_3.add_subplot(545)
ax3_2_2 = fig_3.add_subplot(546)
ax3_2_3 = fig_3.add_subplot(547)
ax3_2_4 = fig_3.add_subplot(548)
ax3_3_1 = fig_3.add_subplot(5,4,9)
ax3_3_2 = fig_3.add_subplot(5,4,10)
ax3_3_3 = fig_3.add_subplot(5,4,11)
ax3_3_4 = fig_3.add_subplot(5,4,12)
ax3_4_1 = fig_3.add_subplot(5,4,13)
ax3_4_2 = fig_3.add_subplot(5,4,14)
ax3_4_3 = fig_3.add_subplot(5,4,15)
ax3_4_4 = fig_3.add_subplot(5,4,16)
ax3_5_1 = fig_3.add_subplot(5,4,17)
ax3_5_2 = fig_3.add_subplot(5,4,18)
ax3_5_3 = fig_3.add_subplot(5,4,19)
ax3_5_4 = fig_3.add_subplot(5,4,20)

# size to work with 
# ******************
img_wk_size = 560 

# bring the orig img down to (560, 560) 
# ***************************************
imgvc = Image.open("rosen_orig_farbe.jpg")
imgvc_wk_size = imgvc.resize((img_wk_size,img_wk_size), resample=PIL.Image.BICUBIC)
# Change to np arrays 
ay_picc = np.array(imgvc_wk_size)
print("ay_picc.shape = ", ay_picc.shape)
print("r = ", ay_picc[0,0,0],    " g = ", ay_picc[0,0,1],    " b = " , ay_picc[0,0,2] )
print("r = ", ay_picc[200,0,0],  " g = ", ay_picc[200,0,1],  " b = " , ay_picc[200,0,2] )

# Turn color to gray 
#Red * 0.3 + Green * 0.59 + Blue * 0.11
#Red * 0.2126 + Green * 0.7152 + Blue * 0.0722
#Red * 0.299 + Green * 0.587 + Blue * 0.114

ay_picc_g = ( 0.299  * ay_picc[:,:,0] + 0.587  * ay_picc[:,:,1] + 0.114  * ay_picc[:,:,2] )  
ay_picc_g = ay_picc_g.astype('float32') 
t_picc_g  = ay_picc_g.reshape((1, img_wk_size, img_wk_size, 1))
t_picc_g  = tf.image.per_image_standardization(t_picc_g)

# downsize to (28,28)  
t_picc_g_28 = tf.image.resize(t_picc_g, [28,28], method="bicubic", antialias=True)
t_picc_g_28 = tf.image.per_image_standardization(t_picc_g_28)

# get the correction for the full image 
t_picc_g_wk_size = tf.image.resize(t_picc_g_28, [img_wk_size,img_wk_size], method="bicubic", antialias=True)
t_picc_g_wk_size = tf.image.per_image_standardization(t_picc_g_wk_size)

t_picc_g_wk_size_corr = t_picc_g - t_picc_g_wk_size
t_picc_g_wk_size_re   = t_picc_g_wk_size + t_picc_g_wk_size_corr


# Display wk_size orig images  
ax3_1_1.imshow(imgvc_wk_size)
ax3_1_2.imshow(t_picc_g[0,:,:,0], cmap=plt.cm.gray)
ax3_1_3.imshow(t_picc_g_28[0,:,:,0], cmap=plt.cm.gray)
ax3_1_4.imshow(t_picc_g_wk_size_re[0,:,:,0], cmap=plt.cm.gray)


# Split in 4 sub-images 
# ***********************

half_wk_size  = int(img_wk_size / 2)  

t_picc_g_1  = t_picc_g[:, 0:half_wk_size, 0:half_wk_size, :]
t_picc_g_2  = t_picc_g[:, 0:half_wk_size, half_wk_size:img_wk_size, :]
t_picc_g_3  = t_picc_g[:, half_wk_size:img_wk_size, 0:half_wk_size, :]
t_picc_g_4  = t_picc_g[:, half_wk_size:img_wk_size, half_wk_size:img_wk_size, :]

# Display wk_size orig images  
ax3_2_1.imshow(t_picc_g_1[0,:,:,0], cmap=plt.cm.gray)
ax3_2_2.imshow(t_picc_g_2[0,:,:,0], cmap=plt.cm.gray)
ax3_2_3.imshow(t_picc_g_3[0,:,:,0], cmap=plt.cm.gray)
ax3_2_4.imshow(t_picc_g_4[0,:,:,0], cmap=plt.cm.gray)


# Downscale sub-images 
t_picc_g_1_28 = tf.image.resize(t_picc_g_1, [28,28], method="bicubic", antialias=True)
t_picc_g_2_28 = tf.image.resize(t_picc_g_2, [28,28], method="bicubic", antialias=True)
t_picc_g_3_28 = tf.image.resize(t_picc_g_3, [28,28], method="bicubic", antialias=True)
t_picc_g_4_28 = tf.image.resize(t_picc_g_4, [28,28], method="bicubic", antialias=True)


# Display downscales sub-images  
ax3_3_1.imshow(t_picc_g_1_28[0,:,:,0], cmap=plt.cm.gray)
ax3_3_2.imshow(t_picc_g_2_28[0,:,:,0], cmap=plt.cm.gray)
ax3_3_3.imshow(t_picc_g_3_28[0,:,:,0], cmap=plt.cm.gray)
ax3_3_4.imshow(t_picc_g_4_28[0,:,:,0], cmap=plt.cm.gray)


# get correction values for upsizing 
t_picc_g_1_wk_half = tf.image.resize(t_picc_g_1_28, [half_wk_size,
half_wk_size], method="bicubic", antialias=True)
t_picc_g_1_wk_half = tf.image.per_image_standardization(t_picc_g_1_wk_half)
t_picc_g_2_wk_half = tf.image.resize(t_picc_g_2_28, [half_wk_size,half_wk_size], method="bicubic", antialias=True)
t_picc_g_2_wk_half = tf.image.per_image_standardization(t_picc_g_2_wk_half)
t_picc_g_3_wk_half = tf.image.resize(t_picc_g_3_28, [half_wk_size,half_wk_size], method="bicubic", antialias=True)
t_picc_g_3_wk_half = tf.image.per_image_standardization(t_picc_g_3_wk_half)
t_picc_g_4_wk_half = tf.image.resize(t_picc_g_4_28, [half_wk_size,half_wk_size], method="bicubic", antialias=True)
t_picc_g_4_wk_half = tf.image.per_image_standardization(t_picc_g_4_wk_half)
 
    
t_picc_g_1_corr =  t_picc_g_1 - t_picc_g_1_wk_half   
t_picc_g_2_corr =  t_picc_g_2 - t_picc_g_2_wk_half   
t_picc_g_3_corr =  t_picc_g_3 - t_picc_g_3_wk_half   
t_picc_g_4_corr =  t_picc_g_4 - t_picc_g_4_wk_half   

t_picc_g_1_re   = t_picc_g_1_wk_half + t_picc_g_1_corr 
t_picc_g_2_re   = t_picc_g_2_wk_half + t_picc_g_2_corr 
t_picc_g_3_re   = t_picc_g_3_wk_half + t_picc_g_3_corr 
t_picc_g_4_re   = t_picc_g_4_wk_half + t_picc_g_4_corr 

print(t_picc_g_1_re.shape)

# Display downscales sub-images  
ax3_4_1.imshow(t_picc_g_1_re[0,:,:,0], cmap=plt.cm.gray)
ax3_4_2.imshow(t_picc_g_2_re[0,:,:,0], cmap=plt.cm.gray)
ax3_4_3.imshow(t_picc_g_3_re[0,:,:,0], cmap=plt.cm.gray)
ax3_4_4.imshow(t_picc_g_4_re[0,:,:,0], cmap=plt.cm.gray)

ay_img_comp = np.zeros((img_wk_size,img_wk_size))
ay_t_comp   = ay_img_comp.reshape((1,img_wk_size, img_wk_size,1))
t_pic_comp  = ay_t_comp
t_pic_comp[0, 0:half_wk_size, 0:half_wk_size, 0]                       = t_picc_g_1_re[0, :, :, 0]
t_pic_comp[0, 0:half_wk_size, half_wk_size:img_wk_size, 0]             = t_picc_g_2_re[0, :, :, 0]
t_pic_comp[0, half_wk_size:img_wk_size, 0:half_wk_size, 0]             = t_picc_g_3_re[0, :, :, 0]
t_pic_comp[0, half_wk_size:img_wk_size, half_wk_size:img_wk_size, 0]   = t_picc_g_4_re[0, :, :, 0]

ax3_5_1.imshow(t_picc_g[0,:,:,0], cmap=plt.cm.gray)
ax3_5_2.imshow(t_pic_comp[0,:,:,0], cmap=plt.cm.gray)

 
Note the computation of the correction tensors; we prove their effectiveness below by displaying respective results for the full image, cut out and re-upscaled, corrected quadrants and replaced areas of the original image.

The results in the prepared image frames look like:

Good!

Some code for the carving algorithm – applied to the full image AND a single selected sub-image

Now, for a first test, we concentrate on the sub-image in the lower-right corner – and apply the above algorithm. A corresponding Jupyter cell code is given below:

# *************************************************************************
# OIP analysis on sub-image tiles (to be used after the previous cell)
# **************************************************************************
# Note: To be applied after previous cell !!!
# ******

#interactive plotting - will be used in a future version when we deal with "octaves" 
#%matplotlib notebook 
#plt.ion()

%matplotlib inline

# preparation of figure frames 
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
fig_size = plt.rcParams["figure.figsize"]
fig_size[0] = 12
fig_size[1] = 12
fig_4 = plt.figure(4)
ax4_1_1 = fig_4.add_subplot(441)
ax4_1_2 = fig_4.add_subplot(442)
ax4_1_3 = fig_4.add_subplot(443)
ax4_1_4 = fig_4.add_subplot(444)

ax4_2_1 = fig_4.add_subplot(445)
ax4_2_2 = fig_4.add_subplot(446)
ax4_2_3 = fig_4.add_subplot(447)
ax4_2_
4 = fig_4.add_subplot(448)

ax4_3_1 = fig_4.add_subplot(4,4,9)
ax4_3_2 = fig_4.add_subplot(4,4,10)
ax4_3_3 = fig_4.add_subplot(4,4,11)
ax4_3_4 = fig_4.add_subplot(4,4,12)

ax4_4_1 = fig_4.add_subplot(4,4,13)
ax4_4_2 = fig_4.add_subplot(4,4,14)
ax4_4_3 = fig_4.add_subplot(4,4,15)
ax4_4_4 = fig_4.add_subplot(4,4,16)

fig_size = plt.rcParams["figure.figsize"]
fig_size[0] = 12
fig_size[1] = 6
fig_5 = plt.figure(5)
axa5_1 = fig_5.add_subplot(241)
axa5_2 = fig_5.add_subplot(242)
axa5_3 = fig_5.add_subplot(243)
axa5_4 = fig_5.add_subplot(244)
axa5_5 = fig_5.add_subplot(245)
axa5_6 = fig_5.add_subplot(246)
axa5_7 = fig_5.add_subplot(247)
axa5_8 = fig_5.add_subplot(248)
li_axa5 = [axa5_1, axa5_2, axa5_3, axa5_4, axa5_5, axa5_6, axa5_7, axa5_8]

# Some general OIP run parameters 
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
n_epochs  = 30          # should be divisible by 5  
n_steps   = 6           # number of intermediate reports 
epsilon   = 0.01        # step size for gradient correction  
conv_criterion = 2.e-4  # criterion for a potential stop of optimization 

# image width parameters
# ~~~~~~~~~~~~~~~~~~~~~
img_wk_size   = 560    # must be identical to the last cell   
half_wk_size  = int(img_wk_size / 2) 

# min / max values of the input image 
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Note : to be used for contrast enhancement 
min1 = tf.reduce_min(t_picc_g)
max1 = tf.reduce_max(t_picc_g)
# parameters to deal with the spectral distribution of gray values
spect_dist = min(abs(min1), abs(max1))
spect_fact   = 0.85
height_fact  = 1.1

# Set the gray downscaled input image as a startng point for the iteration loop 
# ~~~~~~~~~~~~~----------~~~~--------------~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
MyOIP._initial_inp_img_data = t_picc_g_28

# ************
# Main Loop 
# ***********

num_main_iterate = 4
# map_index_main = -1       # take a cost function for a full layer 
map_index_main = 56         # map-index to be used for the OIP analysis of the full image  
map_index_sub  = 56         # map-index to be used for the OIP analysis of the sub-images   


for j in range(num_main_iterate):
   
    # ******************************************************
    # deal with the full (downscaled) input image 
    # ******************************************************
    
    # Apply OIP-algorithm to the whole downscaled image  
    # ~~~~~~~~~~~~~~~~~~~-----------~~~~~~~~~~~~~~~~~~~~~~~
    map_index_main_oip = map_index_main   # map-index we are interested in 

    # Perform the OIP analysis 
    MyOIP._derive_OIP(map_index = map_index_oip, 
                      n_epochs = n_epochs, n_steps = n_steps, 
                      epsilon = epsilon , conv_criterion = conv_criterion, 
                      b_stop_with_convergence=False,
                      b_print=False,
                      li_axa = li_axa5,
                      ax1_1 = ax4_1_1, ax1_2 = ax4_1_2)

    # display the modified downscaled image  (without and with contrast)
    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    t_oip_g_28  = MyOIP._inp_img_data
    ay_oip_g_28 = t_oip_g_28[0,:,:,0].numpy()
    ay_oip_g_28_cont = MyOIP._transform_tensor_to_img(T_img=t_oip_g_28[0,:,:,0], centre_move=0.33, fact=1.0)
    ax4_1_3.imshow(ay_oip_g_28_cont, cmap=plt.cm.gray)
    ax4_1_4.imshow(t_picc_g_28[0, :, :, 0], cmap=plt.cm.gray)
    
    # rescale to 560 and re-add details via the correction tensor  
    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    t_oip_g_wk_size        = tf.image.resize(t_oip_g_28, [img_wk_size,img_wk_size], 
                                             method="bicubic", antialias=True)
    t_oip_g_wk_size_re     = t_oip_g_wk_size + t_picc_g_wk_size_corr 
    
    # standardize to get an intermediate image 
    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
n    # t_oip_g_loop           = t_oip_g_wk_size_re.numpy()
    # t_oip_g_wk_size_re_std = tf.image.per_image_standardization(t_oip_g_wk_size_re)
    t_oip_g_wk_size_re_std = tf.image.per_image_standardization(t_oip_g_wk_size_re)
    t_oip_g_loop           = t_oip_g_wk_size_re_std.numpy()
    
    # contrast required as the reduction of irrelvant pixels has smoothed out the gray scale beneath 
    # ~~~~~~~~~~~~~~~~~
    # the added details at pattern areas = high level of whitened socket + relative small detail variation    
    t_oip_g_wk_size_re_std_plt = tf.clip_by_value(height_fact*t_oip_g_wk_size_re_std, -spect_dist*spect_fact, spect_dist*spect_fact)
    ax4_2_1.imshow(t_oip_g_wk_size[0,:,:,0], cmap=plt.cm.gray)
    ax4_2_2.imshow(t_oip_g_wk_size_re_std[0,:,:,0], cmap=plt.cm.gray)
    ax4_2_3.imshow(t_oip_g_wk_size_re_std_plt[0,:,:,0], cmap=plt.cm.gray)
    ax4_2_4.imshow(t_picc_g[0,:,:,0], cmap=plt.cm.gray)
   

    # *************************************************************
    # deal with a chosen sub-image - in this version: "t_oip_4_g"
    # **************************************************************
   
    # Cut out and downscale a sub-image region 
    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    t_oip_4_g     = t_oip_g_loop[:, half_wk_size:img_wk_size, half_wk_size:img_wk_size, :]
    t_oip_4_g_28  = tf.image.resize(t_oip_4_g, [28,28], method="bicubic", antialias=True)
    
    # use the downscaled sub-image as an input image to the OIP-algorithm 
    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    MyOIP._initial_inp_img_data = t_oip_4_g_28
    
    
    # Perform the OIP analysis on the sub-image
    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    map_index_sub_oip  = map_index_sub   # we may vary this in later versions 
    
    MyOIP._derive_OIP(map_index = map_index_sub_oip, 
                      n_epochs = n_epochs, n_steps = n_steps, 
                      epsilon = epsilon , conv_criterion = conv_criterion, 
                      b_stop_with_convergence=False,
                      b_print=False,
                      li_axa = li_axa5,
                      ax1_1 = ax4_3_1, ax1_2 = ax4_3_2)
    
    
    # display the modified sub-image without and with contrast  
    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    t_oip_4_g_28  = MyOIP._inp_img_data
    ay_oip_4_g_28 = t_oip_4_g_28[0,:,:,0].numpy()
    ay_oip_4_g_28_cont = MyOIP._transform_tensor_to_img(T_img=t_oip_4_g_28[0,:,:,0], centre_move=0.33, fact=1.0)
    ax4_3_3.imshow(ay_oip_4_g_28_cont, cmap=plt.cm.gray)
    ax4_3_4.imshow(t_oip_4_g[0, :, :, 0], cmap=plt.cm.gray)
    
    # upscaling of the present sub-image 
    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    t_pic_4_g_wk_half = tf.image.resize(t_oip_4_g_28, [half_wk_size,half_wk_size], method="bicubic", antialias=True)
    
    # add the detail correction to the manipulated upscaled sub-image
    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    t_pic_4_g_wk_half = t_pic_4_g_wk_half + t_picc_g_4_corr
    #t_pic_4_g_wk_half = tf.image.per_image_standardization(t_pic_4_g_wk_half)
    ax4_4_1.imshow(t_pic_4_g_wk_half[0, :, :, 0], cmap=plt.cm.gray)
    ax4_4_2.imshow(t_oip_4_g[0, :, :, 0], cmap=plt.cm.gray)
    
    # Overwrite the related region in the full image with the new sub-image  
    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    t_oip_g_loop[0, half_wk_size:img_wk_size, half_wk_size:img_wk_size, 0]  = t_pic_4_g_wk_half[0, :, :, 0]
    t_oip_g_loop_std = tf.image.per_image_standardization(t_oip_g_loop)
    
    #ax4_4_3.imshow(t_oip_g_loop[0, :, :, 0], cmap=plt.cm.gray)   
    ax4_4_3.imshow(t_oip_g_loop_std[0, :, :, 0], cmap=plt.cm.gray)   
    ay_oip_g_loop_std_plt = tf.clip_by_value(height_fact*t_oip_g_loop_std, -spect_dist*spect_fact, spect_dist*spect_fact)
    ax4_4_4.
imshow(ay_oip_g_loop_std_plt[0, :, :, 0], cmap=plt.cm.gray)
    
    # Downscale the resulting new full image and feed it into into the next iteration of the main loop
    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    MyOIP._initial_inp_img_data = tf.image.resize(t_oip_g_loop_std, [28,28], method="bicubic", antialias=True)
   

 

The result of these operations after 4 iterations is displayed in the following image:

We recognize again the worm-like shape resulting for map 56, which we have seen in the last article, already. (By the way: Our CNN’s map 56 originally is strongly triggered for the shapes of handwritten 9-digits and plays a significant role in classifying respective images correctly). But there is a new feature which appeared in the lower right part of the image – a kind of wheel or two closely neighbored 9-like shapes.

The following code compares the final result with the original image:

fig_size = plt.rcParams["figure.figsize"]
fig_size[0] = 12
fig_size[1] = 12
fig_7X = plt.figure(7)
ax1_7X_1 = fig_7X.add_subplot(221)
ax1_7X_2 = fig_7X.add_subplot(222)
ax1_7X_3 = fig_7X.add_subplot(223)
ax1_7X_4 = fig_7X.add_subplot(224)

ay_pic_full_cont = MyOIP._transform_tensor_to_img(T_img=t_oip_g_loop_std[0, :, :, 0], centre_move=0.46, fact=0.8)

ax1_7X_1.imshow(imgvc_wk_size)
ax1_7X_2.imshow(t_picc_g[0,:,:,0], cmap=plt.cm.gray)
ax1_7X_3.imshow(ay_oip_g_loop_std_plt[0, :, :, 0], cmap=plt.cm.gray)
ax1_7X_4.imshow(ay_pic_full_cont, cmap=plt.cm.gray)

Giving:

Those who do not like the relative strong contrast may do their own experiments with other parameters for contrast enhancement.

You see: Simplifying prejudices about the world, which get or got manifested in inflexible neural networks (or brains ?), may lead to strange, if not bad dreams. May also some politicians learn from this … hmmm, just joking, as this requires both an openness and capacity for learning … can’t await Jan, 21st …

Conclusion

Our “carving algorithm”, which corresponds to an amplification of traces of map-related patterns that a CNN detects in an input image, can also be used for sub-image-areas and their pixel information. We can therefore use even very limited CNNs trained on low resolution MNIST images to create a “Deep MNIST Dream” out of a chosen arbitrary input image. At least in principle.

Addendum 06.01.2021: There is, however, a pitfall – our CNN (for MNIST or other image samples) may have been trained and work on standardized image/pixel data, only. But even if the overall input image may have been standardized before operating on it in the sense of the algorithm discussed in this article, we should take into account that a sub-image area may contain pixel data which are far from a standardized distribution. So, our corrections for details may require more than just the calculation of a difference term. We may have to reverse intermediate standardization steps, too. We shall take care of this in forthcoming code changes.

In the next article we are going to extend our algorithm to multiple cascaded sub-areas of our image – and we vary the maps used for different sub-areas at the same time.

 

Manual requests to a https web-server with “openssl s_client”: do not forget the options “-ign_eof -crlf”

Sometimes you may need to analyze the behavior and responses of a web-server or a REST service to certain requests. And sometimes you are restricted to the command line of a Linux system (e.g. during penetration testing). Then you have to type and send HTTP commands in a direct manner. While this is trivial with telnet and HTTP-commands via n unencrypted connection on port 80, you must use a tool like openssl for HTTPS-servers using TLS tunnels.

A quick search on the Internet will show you that you should be able to use “openssl” on a Linux system in the following form:

openssl s_client -connect YOUR_TARGET_WEB_DOMAIN:443

Or – if you do not want to look at certificates and related CA chains in detail – with an additional option “-quiet”:

openssl s_client -quiet -connect YOUR_TARGET_WEB_DOMAIN:443

YOUR_TARGET_WEB_DOMAIN has to be replaced, of course, by a valid URI. For restricting the encryption to TLS V1.2 you would instead use:

openssl s_client -quiet -tls1_2 -connect YOUR_TARGET_WEB_DOMAIN:443

For some servers an additional option “-ign_eof” can be helpful: This hinders a connection to directly close when an “end of file” [EOF] may be reached (during a response). Meaning: The response will not be shown in some cases. The option “-quiet” triggers a “-ign_eof” behavior implicitly. But keep in mind that this option does not hinder any timeouts on the connection imposed by the server.

A problem with the interactive “openssl s_client” command-line on Linux systems

After entering the above commands at the command prompt of a Linux shell (e.g. bash) you will first see some information regarding the connection handshake and the establishment of the encryption tunnel. Then end up on a line where you can interactively enter HTTP commands. You expect to successfully enter commands in the following way:

We type “GET / HTTP/1.1.” (without the quotes)    =>    we press the “ENTER”-key    =>    we get a new line    =>    we type “Host: YOUR_TARGET_WEB_DOMAIN” (without the quotes)    =>    we press the “ENTER”-key twice.

You may try this sequence with “google.com”. This will work! You get, however, an information that the document has been moved to “www.google.com”. But at the new address everything is working, too.

So far so good! But let us try the given recipe with another domain: www.debian.org. Using the command line within “s_client” then leads to an error:

.....
    Start Time: 1609195616
    Timeout   : 7200 (sec)
    Verify return code: 0 (ok)
    Extended master secret: no
    Max Early Data: 0
---
read R BLOCK
GET / HTTP/1.1
HTTP/1.1 400 Bad Request
Date: Mon, 28 Dec 2020 22:47:04 GMT
Server: Apache
X-Content-Type-Options: nosniff
X-Frame-Options: sameorigin
Referrer-Policy: no-referrer
X-Xss-Protection: 1
Strict-Transport-Security: max-age=15552000
Content-Length: 291
Connection: close
Content-Type: text/html; charset=iso-8859-1

<!DOCTYPE HTML PUBLIC "-//IETF//DTD HTML 2.0//EN">
<html><head>
<title>400 Bad Request</title>
</head><body>
<h1>Bad Request</h1>
<p>Your browser sent a request that this server could not understand.<br />
</p>
<hr>
<address>Apache Server at www.debian.org Port 443</address>
</body></html>
closed

You do not even get a chance to enter the “Host: …” line!
The same is true for other servers – as the one supporting e.g. one of my own domains “anracon.de”.

A simple trick shows what the correct request format is and that it works …

Let us use a
simple trick with echo and a pipe on the Linux shell:

myself@mytux:~> (echo -ne "GET / HTTP/1.1\r\nHost: www.debian.org\r\n\r\n") | openssl s_client -tls1_2 -quiet  -connect www.debian.org:443

This leads to

depth=2 O = Digital Signature Trust Co., CN = DST Root CA X3
verify return:1
depth=1 C = US, O = Let's Encrypt, CN = Let's Encrypt Authority X3
verify return:1
depth=0 CN = www.debian.org
verify return:1
HTTP/1.1 200 OK
Date: Mon, 28 Dec 2020 23:09:01 GMT
Server: Apache
Content-Location: index.en.html
Vary: negotiate,accept-language,Accept-Encoding,cookie
TCN: choice
X-Content-Type-Options: nosniff
X-Frame-Options: sameorigin
Referrer-Policy: no-referrer
X-Xss-Protection: 1
Strict-Transport-Security: max-age=15552000
Upgrade: h2,h2c
Connection: Upgrade
Last-Modified: Sun, 27 Dec 2020 19:27:21 GMT
ETag: "36b1-5b777257b5a41"
Accept-Ranges: bytes
Content-Length: 14001
Cache-Control: max-age=86400
Expires: Tue, 29 Dec 2020 23:09:01 GMT
X-Clacks-Overhead: GNU Terry Pratchett
Content-Type: text/html
Content-Language: en

<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd">
<html lang="en">
<head>
  <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
  <title>Debian -- The Universal Operating System </title>
...
...
</div>
<!--/UdmComment-->
</div> <!-- end footer -->
</body>
</html>

 
After a timeout we get back our prompt.

So, we see that the server reacts properly for the end of line characters used, namely “\r\n” after “GET / HTTP/1.1” and after “Host: www.debian.org” plus after an empty line (leading to 2 “\r\n” at the end of the request).

But if we this change to

myself@mytux:~> (echo -ne "GET / HTTP/1.1\nHost: www.debian.org\r\n\r\n") | openssl s_client -tls1_2 -quiet  -connect www.debian.org:443

we run into a “HTTP/1.1 400 Bad Request” answer again
(Note the “\n” before “Host: …”)

Obviously, our “openssl s_client” interface works with LF-characters when pressing the “ENTER”-key on the command line. Well, the interface uses the typical Linux/Unix-style for EOLs …

By the way:
Some servers set very short timeouts for receiving all request lines – you may have difficulties with typing fast enough. Then the “trick” with an echo-command and a pipe is very useful to test the server none the less.

What is the reason for the unequal behavior of different servers?

The correct way to build HTTP-Requests is described e.g. in “www.tutorialspoint.com (http/http_requests.htm)“; I quote:

  • A Request-line
  • Zero or more header (General|Request|Entity) fields followed by CRLF
  • An empty line (i.e., a line with nothing preceding the CRLF) indicating the end of the header fields
  • Optionally a message-body

According to RFC7230 a HTTP-Request line has the following format:

request-line = method SP request-target SP HTTP-version CRLF

However, in section 3.5 the named RFC also says:

Although the line terminator for the start-line and header fields is the sequence CRLF, a recipient MAY recognize a single LF as a line terminator and ignore any preceding CR.

So, this explains the different behavior of some web-servers to the commands sent with “openssl s_client”.

What can we do with “openssl s_client” to enforce a CRLF as the end of lines when pressing the “ENTER”-key?

The ”
openssl s_client” has a lot of options. You find an overview at the following URI:
https://zoomadmin.com / HowToLinux / LinuxCommand / s_client
The required option for our purpose is “-crlf“; we thus arrive at:

openssl s_client <strong>-quiet -crlf</strong> -tls1_2 -connect YOUR_TARGET_WEB_DOMAIN:443

as the right way to bring RFC compliant web-servers to answer without error messages to manually sent HTTP requests within “openssl s_client”.

Prepare your terminal for long answers from the server

Some web-sites may present long web pages. The HTTP/HTML code sent to you as an answer may be pretty long. As most Linux terminals limit the scroll-able output length by default you should look out for settings that allow for long or infinite output within the terminal emulation of your choice. KDE’s “konsole” e.g. offers you an option to allow for scrolling output with unlimited length. A noteworthy side-effect is, however, that the contents is saved unencrypted in temporary files (for which you can define a location on your PC). So, be a bit careful when using such options.

Conclusion

Using CRLF is the RFC defined way to properly end HTTP request and header field lines. (… whether we from the Linux side may like this MS influenced definition or not …). To enforce an “openssl s_client” to interpret the signal from an “ENTER”-key as “CRLF” (instead of “LF”) we should use the option “-crlf” when opening “s_client”. The additional options “-ign_eof” or “-quiet” are useful to prevent a shutdown of the connection before the server’s answer is fully displayed.

Links

https://zoomadmin.com/HowToLinux/LinuxCommand/s_client
https://stackoverflow.com/questions/5757290/http-header-line-break-style
https://tools.ietf.org/html/rfc7230#section-3.5
https://www.tutorialspoint.com/http/http_requests.htm
https://www.tutorialspoint.com/http/http_quick_guide.htm

Deep Dreams of a CNN trained on MNIST data – II – some code for pattern carving

In the last article of this series

Deep Dreams of a CNN trained on MNIST data – I – a first approach based on one selected map of a convolutional layer

I presented some “Deep Dream” images of a CNN trained on MNIST data. Although one cannot produce spectacular images with a simple CNN based on gray scale, low resolution image data, we still got a glimpse of what the big guys do with more advanced CNNs and image data. But we have only started, yet ….

In the forthcoming articles we shall extend our abilities step by step beyond the present level. Our next goal is to work on different length scales within the image, i.e. we go down into sub-tiles of the original image, analyze the sub-structures there and include the results in the CNN’s dreams. But first we need to understand more precisely how we apply the image manipulation techniques that “carve” some additional dream-like optical elements, which correspond to map triggering features, into an arbitrary input image presented to the CNN.

Carving: A combination of low-resolution OIP-pattern creation and a reproduction of high-resolution details

I used the term “carved” intentionally. Artists or craftsmen who create little figures or faces out of wooden blocks by a steady process of carving often say that they only bring to the surface what already was there – namely in the basic fiber structures of the wood. Well, the original algorithm for the creation of pure OIP- or feature patterns systematically amplifies some minimal correlated pixel “structures” detected in an input image with random pixel values. The values of all other pixels outside the OIP pattern are reduced; sooner or later they form a darker and rather homogeneous background. The “detection” is based on (trained) filters which react to certain pixel constellations. We already encoded an adaption of a related method for CNN filter visualization, which Francois Chollet discusses in his book on “Deep learning with Python”, into a Python class in another article series of this blog (see the link section at the end of this article).

There is, however, a major difference regarding our present context: During the image manipulations for the visualization of Deep Dreams we do not want to remove all details of the original image presented to the DeepDream algorithm. Instead, we just modify it at places where the CNN detects traces of patterns which trigger a strong response of one or multiple maps of our CNN – and keep detail information alive elsewhere. I.e., we must counterbalance the tendency of our OIP algorithm to level out information outside the OIP structures.

How do we do this, when both our present OIP algorithm and the MNIST-trained CNN are limited to a resolution of (28×28) pixels, only?

Carving – amplification of certain OIP structures and intermittent replenishment of details of the original image

In this article we follow the simple strategy explained in the last article:

  • Preparation:
    • We choose a map.
    • We turn the colored original input image (about which the CNN shall dream) into a gray one.
    • We downscale the original image to a size of (28×28) by interpolation, upscale the result again by interpolating again (with loss) and calculate the difference to the original full resolution image (all interpolations done in a bicubic way). (The difference contains information on detail structures.)
  • Loop (4 times or so):
    • We apply the OIP-algorithm on the downscaled input image for a fixed amount of epochs with a small epsilon.
      The result is a slight emphasis of
      pixel structures resembling OIP patterns.
    • We upscale the result by bicubic interpolation to the original size.
    • We re-add the difference to the original image, i.e. we supply details again.
    • We downscale the result again.
  • Eventual image optimization: We apply some contrast treatment to the resulting image, whose pixel value range and distribution has gradually been changed in comparison to the original image.

This corresponds to a careful accentuation of an OIP-pattern detected by a certain map in the given structured pixel data of the input image (on coarse scales) – without loosing all details of the initial image. As said – it resembles a process of carving coarse patterns into the original image ..

Note: The above algorithm works on the whole image – a modified version working on a variety of smaller length scales is discussed in further articles.

What does this strategy look like in Python statements for Jupyter cells?

Start with loading required modules and involving your graphics card

We have to load python modules to support our CNN, the OIP algorithm and some plotting. In a second step we need to invoke the graphics card. This is done in two Jupyter cells (1 and 2). Their content was already described in the article

A simple CNN for the MNIST dataset – XI – Python code for filter visualization and OIP detection

of a parallel series in this blog. These initial steps are followed by an instantiation of a class “My_OIP” containing all required methods. The __init__() function of this class restores a Keras model of our basic CNN from a h5-file:

Jupyter cell 3

%matplotlib inline

# Load the CNN-model 
#  ~~~~~~~~~~~~~~~~~
imp.reload(myOIP)
try:
    with tf.device("/GPU:0"):
        MyOIP = myOIP.My_OIP(cnn_model_file = 'cnn_best.h5', layer_name = 'Conv2D_3')

except SystemExit:
    print("stopped")
    

Turn the original image into a gray one, downscale it and calculate difference-tensors for coarsely re-enlarged versions

In the next Jupyter cell we deal with the following initial image:

We choose a method (out of 3) to turn it into a gray scaled image. Afterwards we calculate the difference between an image tensor of the original (gray) image and a coarse, but softly rescaled image derived from a (28×28)-pixel version.

Jupyter cell 4

# Down-, gray- and Re-scale images, determine correction differences  
# *******************************************************************
import PIL
from PIL import Image, ImageOps

fig_size = plt.rcParams["figure.figsize"]
fig_size[0] = 15
fig_size[1] = 30
fig_1 = plt.figure(1)
ax1_1 = fig_1.add_subplot(531)
ax1_2 = fig_1.add_subplot(532)
ax1_3 = fig_1.add_subplot(533)
ax2_1 = fig_1.add_subplot(534)
ax2_2 = fig_1.add_subplot(535)
ax2_3 = fig_1.add_subplot(536)
ax3_1 = fig_1.add_subplot(537)
ax3_2 = fig_1.add_subplot(538)
ax3_3 = fig_1.add_subplot(539)
ax4_1 = fig_1.add_subplot(5,3,10)
ax4_2 = fig_1.add_subplot(5,3,11)
ax4_3 = fig_1.add_subplot(5,3,12)

# Parameters 
# ************
# size of the image to work with 
img_wk_size = 560 
# method to turn the original picture into a gray scaled one 
gray_scaling_
method = 2

# bring the orig img down to (560, 560) 
# ***************************************
imgvc    = Image.open("rosen_orig_farbe.jpg")
imgvc_g  = imgvc.convert("L")

# downsize with PIL 
# ******************
imgvc_wk_size    = imgvc.resize((img_wk_size,img_wk_size), resample=PIL.Image.BICUBIC)
imgvc_wk_size_g  = imgvc_wk_size.convert("L")
imgvc_28      = imgvc.resize(  (28,28), resample=PIL.Image.BICUBIC)
imgvc_g_28    = imgvc_g.resize((28,28), resample=PIL.Image.BICUBIC)

# Change to np array
ay_picc = np.array(imgvc_wk_size)

# Display orig and wk_size images  
# *********************************
ax1_1.imshow(imgvc)
ax1_2.imshow(imgvc_wk_size)
ax1_3.imshow(imgvc_g, cmap=plt.cm.gray)

# Apply 3 different methods to turn the image into a gray one 
# **************************************************************
#Red * 0.3 + Green * 0.59 + Blue * 0.11
#Red * 0.2126 + Green * 0.7152 + Blue * 0.0722
#Red * 0.299 + Green * 0.587 + Blue * 0.114

ay_picc_g1 = ( 0.3    * ay_picc[:,:,0] + 0.59   * ay_picc[:,:,1] + 0.11   * ay_picc[:,:,2] )  
ay_picc_g2 = ( 0.299  * ay_picc[:,:,0] + 0.587  * ay_picc[:,:,1] + 0.114  * ay_picc[:,:,2] )  
ay_picc_g3 = ( 0.2126 * ay_picc[:,:,0] + 0.7152 * ay_picc[:,:,1] + 0.0722 * ay_picc[:,:,2] )  

ay_picc_g1 = ay_picc_g1.astype('float32') 
ay_picc_g2 = ay_picc_g2.astype('float32') 
ay_picc_g3 = ay_picc_g3.astype('float32') 

# Prepare tensors
# *****************
t_picc_g1 = ay_picc_g1.reshape((1, img_wk_size, img_wk_size, 1))
t_picc_g2 = ay_picc_g2.reshape((1, img_wk_size, img_wk_size, 1))
t_picc_g3 = ay_picc_g3.reshape((1, img_wk_size, img_wk_size, 1))

t_picc_g1 = tf.image.per_image_standardization(t_picc_g1)
t_picc_g2 = tf.image.per_image_standardization(t_picc_g2)
t_picc_g3 = tf.image.per_image_standardization(t_picc_g3)

# Display gray-tensor variants    
# ****************************
ax2_1.imshow(t_picc_g1[0,:,:,0], cmap=plt.cm.gray)
ax2_2.imshow(t_picc_g2[0,:,:,0], cmap=plt.cm.gray)
ax2_3.imshow(t_picc_g3[0,:,:,0], cmap=plt.cm.gray)

# choose one gray img
# ******************
if   gray_scaling_method == 1:
    t_picc_g = t_picc_g1
elif gray_scaling_method == 2:
    t_picc_g = t_picc_g2
else:
    t_picc_g = t_picc_g3

# downsize to (28,28) and standardize
# *****************************************
t_picc_g_28_scd  = tf.image.resize(t_picc_g, [28,28], method="bicubic", antialias=True)
t_picc_g_28_std  = tf.image.per_image_standardization(t_picc_g_28_scd)
t_picc_g_28      = t_picc_g_28_std

# display 
ax3_1.imshow(imgvc_g_28, cmap=plt.cm.gray)
ax3_2.imshow(t_picc_g_28_scd[0,:,:,0], cmap=plt.cm.gray)
ax3_3.imshow(t_picc_g_28_std[0,:,:,0], cmap=plt.cm.gray)

# Upscale and get correction values  
# **********************************
t_picc_g_wk_size_scd     = tf.image.resize(t_picc_g_28, [img_wk_size,img_wk_size], method="bicubic", antialias=True)
t_picc_g_wk_size_scd_std = tf.image.per_image_standardization(t_picc_g_wk_size_scd)
t_picc_g_wk_size_corr    =  t_picc_g - t_picc_g_wk_size_scd_std   

# Correct and display 
# **********************
t_picc_g_wk_size_re   = t_picc_g_wk_size_scd_std + t_picc_g_wk_size_corr 
# Display
ax4_1.imshow(t_picc_g_wk_size_scd_std[0,:,:,0], cmap=plt.cm.gray)
ax4_2.imshow(t_picc_g[0,:,:,0],                 cmap=plt.cm.gray)
ax4_3.imshow(t_picc_g_wk_size_re[0,:,:,0],      cmap=plt.cm.gray)

 
The code is straightforward – a lot of visualization steps will later on just show tiny differences. For production you may shorten the code significantly.

In a first step we downscale the image to a work-size of (560×560) px with the help of PIL modules. (Note that 560 can be divided by 2, 4, 8, 16, 28). We then test out three methods to turn the image into a gray one. This step is required, because our MNIST-oriented CNN is trained for gray images, only. We
reshape the image data arrays to prepare for tensor-operations and standardize with the help of Tensorflow – and produce image tensors on the fly. Afterwards, we downscale to a size of (28×28) pixels; this is the size of the MNIST images for which the CNN has been trained. We upscale to the work-size again, compare it to the original detailed image and calculate respective differences of the tensor-components. These differences will later be used to re-add details to our manipulated images.

The following picture displays the resulting output:

Note: The last two images should not show any difference. The middle row and the displays how coarse the (28×28) resolution really is and how much information we have lost after upscaling to (560×560)px.

OIP-Carving and detail replenishment

We now follow the algorithm described above. This leads to a main loop with steps for downscaling, OIP-carving on the (28×28)-scale, upscaling and a supplementation of original details:

Jupyter cell 5

# **************************************************+
# OIP analysis (to be used after the previous cell)
# **************************************************+

# Prepare plotting 
# -------------------
#interactive plotting 
#%matplotlib notebook 
#plt.ion()

fig_size = plt.rcParams["figure.figsize"]
fig_size[0] = 15
fig_size[1] = 25
fig_2 = plt.figure(2)
ax2_1_1 = fig_2.add_subplot(531)
ax2_1_2 = fig_2.add_subplot(532)
ax2_1_3 = fig_2.add_subplot(533)
ax2_2_1 = fig_2.add_subplot(534)
ax2_2_2 = fig_2.add_subplot(535)
ax2_2_3 = fig_2.add_subplot(536)
ax2_3_1 = fig_2.add_subplot(537)
ax2_3_2 = fig_2.add_subplot(538)
ax2_3_3 = fig_2.add_subplot(539)
ax2_4_1 = fig_2.add_subplot(5,3,10)
ax2_4_2 = fig_2.add_subplot(5,3,11)
ax2_4_3 = fig_2.add_subplot(5,3,12)

fig_size = plt.rcParams["figure.figsize"]
fig_size[0] = 16
fig_size[1] = 8
fig_3 = plt.figure(3)
axa_1 = fig_3.add_subplot(241)
axa_2 = fig_3.add_subplot(242)
axa_3 = fig_3.add_subplot(243)
axa_4 = fig_3.add_subplot(244)
axa_5 = fig_3.add_subplot(245)
axa_6 = fig_3.add_subplot(246)
axa_7 = fig_3.add_subplot(247)
axa_8 = fig_3.add_subplot(248)
li_axa = [axa_1, axa_2, axa_3, axa_4, axa_5, axa_6, axa_7, axa_8]

fig_size = plt.rcParams["figure.figsize"]
fig_size[0] = 16
fig_size[1] = 4
fig_h = plt.figure(50)
ax_h_1 = fig_h.add_subplot(141)
ax_h_2 = fig_h.add_subplot(142)
ax_h_3 = fig_h.add_subplot(143)
ax_h_4 = fig_h.add_subplot(144)
li_axa = [axa_1, axa_2, axa_3, axa_4, axa_5, axa_6, axa_7, axa_8]

# list for img tensors 
li_t_imgs = []

# Parameters for the OIP run
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
map_index = 56         # map-index we are interested in 
# map_index = -1         # for this value the whole layer is taken; the loss function is averaged over the maps
n_epochs  = 20          # should be divisible by 5  
n_steps   = 6           # number of intermediate reports 
epsilon   = 0.01        # step size for gradient correction  
conv_criterion = 2.e-4  # criterion for a potential stop of optimization 

n_iter    = 5           # number of iterations for down-, upscaling and supplementation of details

# Initial image
# -----------------
MyOIP._initial_inp_img_data = t_picc_g_28

# Main iteration loop 
# ***********************
for i in range(n_iter): 

    # Perform the OIP analysis 
    MyOIP._derive_OIP(map_index = map_index, 
                      n_epochs = n_epochs, n_steps = n_steps, 
       
               epsilon = epsilon , conv_criterion = conv_criterion, 
                      b_stop_with_convergence=False,
                      b_print = False,
                      li_axa = li_axa,
                      ax1_1 = ax2_1_1, ax1_2 = ax2_1_2,
                      centre_move = 0.33, fact = 1.0)

    t_oip_c_g_28  = MyOIP._inp_img_data
    ay_oip_c_g_28 = t_oip_c_g_28[0,:,:,0].numpy()

    ay_oip_c_g_28_cont = MyOIP._transform_tensor_to_img(T_img=t_oip_c_g_28[0,:,:,0], centre_move=0.33, fact=1.0)
    ax2_1_3.imshow(ay_oip_c_g_28_cont, cmap=plt.cm.gray)

    # Rescale to 500
    t_oip_c_g_wk_size  = tf.image.resize(t_oip_c_g_28, [img_wk_size,img_wk_size], method="bicubic", antialias=True)
    #t_oip_c_g_500 = tf.image.per_image_standardization(t_oip_c_g_500)

    t_oip_c_g_wk_size_re     = t_oip_c_g_wk_size + t_picc_g_wk_size_corr 
    t_oip_c_g_wk_size_re_std = tf.image.per_image_standardization(t_oip_c_g_wk_size_re)

    # remember the img tensors of the iterations    
    li_t_imgs.append(t_oip_c_g_wk_size_re_std)
    
    # do not change the gray spectrum for iterations 
    MyOIP._initial_inp_img_data = tf.image.resize(t_oip_c_g_wk_size_re_std, [28,28], method="bicubic", antialias=True)
    
    # Intermediate plotting - if required and if plt.ion() is set 
    # ------------------------
    if b_plot_intermediate and j < n_iter-1: 
        ay_oip_c_g_wk_size = t_oip_c_g_wk_size[0,:,:,0].numpy()
        #display OIPs
        ax2_2_1.imshow(t_picc_g_28[0,:,:,0], cmap=plt.cm.gray)
        ax2_2_2.imshow(ay_oip_c_g_28,  cmap=plt.cm.gray)
        ax2_2_3.imshow(ay_oip_c_g_wk_size, cmap=plt.cm.gray)

        ax2_3_1.imshow(t_oip_c_g_wk_size[0,:,:,0], cmap=plt.cm.gray)
        ax2_3_2.imshow(t_oip_c_g_wk_size_re[0,:,:,0], cmap=plt.cm.gray)
        ax2_3_3.imshow(t_picc_g[0,:,:,0], cmap=plt.cm.gray)
    
        # modify the gray spectrum for plotting - cut off by clipping to limit the extend of the gray color map  
        min1 = tf.reduce_min(t_picc_g)
        min2 = tf.reduce_min(t_oip_c_g_wk_size_re)
        max1 = tf.reduce_max(t_picc_g)
        max2 = tf.reduce_max(t_oip_c_g_wk_size_re)
        # tf.print("min1 = ", min1, " :: min2 = ", min2, " :: max1 = ", max1, " :: max2 = ", max2 )

        #fac1 = min2/min1 * 0.95
        #fac2 = 1.0*fac1
        fac2 = 1.18
        fac3 = 0.92
        cut_lft = fac3*min1
        cut_rht = fac3*max1
        cut_rht = fac3*abs(min1)
        tf.print("cut_lft = ", cut_lft, "  cut_rht = ", cut_rht)

        t_oip_c_g_wk_size_re_std_plt = fac2 * t_oip_c_g_wk_size_re_std 
        t_oip_c_g_wk_size_re_std_plt = tf.clip_by_value(t_oip_c_g_wk_size_re_std_plt, cut_lft, cut_rht)
    
        ax2_4_1.imshow(t_oip_c_g_wk_size[0,:,:,0], cmap=plt.cm.gray)
        ax2_4_2.imshow(t_oip_c_g_wk_size_re_std[0,:,:,0], cmap=plt.cm.gray)
        ax2_4_3.imshow(t_oip_c_g_wk_size_re_std_plt[0,:,:,0], cmap=plt.cm.gray)


# Eventual plotting 
# ********************
ay_oip_c_g_wk_size = t_oip_c_g_wk_size[0,:,:,0].numpy()

#display OIPs
ax2_2_1.imshow(t_picc_g_28[0,:,:,0], cmap=plt.cm.gray)
ax2_2_2.imshow(ay_oip_c_g_28,  cmap=plt.cm.gray)
ax2_2_3.imshow(ay_oip_c_g_wk_size, cmap=plt.cm.gray)

ax2_3_1.imshow(t_oip_c_g_wk_size[0,:,:,0], cmap=plt.cm.gray)
ax2_3_2.imshow(t_oip_c_g_wk_size_re[0,:,:,0], cmap=plt.cm.gray)
ax2_3_3.imshow(t_picc_g[0,:,:,0], cmap=plt.cm.gray)

# modify the gray spectrum for plotting - cut off by clipping to limit the extend of the gray color map  
min1 = tf.reduce_min(t_picc_g)
min2 = tf.reduce_min(t_oip_c_g_wk_size_re)
max1 = tf.reduce_max(t_picc_g)
max2 = tf.reduce_max(t_oip_c_g_wk_size_re)
# tf.print("min1 = ", min1, " :: min2 = ", min2, " :: max1 = ", max1, " :: max2 = ", max2 )

#fac1 = min2/min1 * 0.95
#fac2 = 1.0*fac1
fac2 = 1.10
fac3 = 0.92
cut_lft = fac3*min1
cut_rht = fac3*
max1
#cut_rht = fac3*abs(min1)
tf.print("cut_lft = ", cut_lft, "  cut_rht = ", cut_rht)

t_oip_c_g_wk_size_re_std_plt = fac2 * t_oip_c_g_wk_size_re_std 
t_oip_c_g_wk_size_re_std_plt = tf.clip_by_value(t_oip_c_g_wk_size_re_std_plt, cut_lft, cut_rht)

ax2_4_1.imshow(t_oip_c_g_wk_size[0,:,:,0], cmap=plt.cm.gray)
ax2_4_2.imshow(t_oip_c_g_wk_size_re_std[0,:,:,0], cmap=plt.cm.gray)
ax2_4_3.imshow(t_oip_c_g_wk_size_re_std_plt[0,:,:,0], cmap=plt.cm.gray)
    
# Histogram analysis and plotting 
# *********************************    
# prepare the histograms of greay values for the saved imagetensors     
ay0 = np.sort( t_picc_g[0, :, :, 0].numpy().ravel() )
ay1 = np.sort( li_t_imgs[0][0,:,:,0].numpy().ravel() )
ay2 = np.sort( li_t_imgs[1][0,:,:,0].numpy().ravel() )
ay3 = np.sort( li_t_imgs[2][0,:,:,0].numpy().ravel() )
ay4 = np.sort( li_t_imgs[3][0,:,:,0].numpy().ravel() )
#print(ay1)

nx, binsx, patchesx = ax_h_1.hist(ay0, bins='auto')
nx, binsx, patchesx = ax_h_2.hist(ay2, bins='auto')
nx, binsx, patchesx = ax_h_3.hist(ay3, bins='auto')
nx, binsx, patchesx = ax_h_4.hist(ay4, bins='auto')

 

We first prepare some Matplotlib axes frames. Note that the references should be clearly distinguishable and unique across the cells of the Jupyter notebook.

Parameters
Then we define parameters for a DeepDream run. [The attentive reader has noticed that I now allow for a value of nmap_index = -1. This corresponds to working with a loss function calculated and averaged over all maps (nmap=-1) of a layer. The changed code of the class My_OIP is given at the end of the article. We come back to this option in a later post of this series. For the time being we just work with a chosen single map and its cost function.]

The most interesting parameters are “n_epochs” and “n_iter“. These parameters help to tune a careful balance between OIP-“carving” and “detail replenishment”. These are parameters you absolutely should experiment with!

Iterations
We first pick the (28×28) gray image prepared in the previous Jupyter cell as the initial image presented to the CNN’s input layer. We then start an iteration by invoking the “My_OIP” class discussed in the article named above. (The parameters “centre_move = 0.33” and “fact = 1.0” refer to an image manipulation, namely a contrast enhancement.)
After OIP-carving over the chosen limited number of epochs we apply a special treatment to the data for the OPI-image derived so far. In a side step we enhance the contrast contrast for displaying the intermediate results. We then rescale by bicubic interpolation from 28×28 to the working size 560×560, add the prepared correction data and standardize. We add the results to a list for a later histogram analysis.

The preparation of plotting (of an intermediate and eventual image) comprises a step to reduce the spectrum of pixel values. This step is required as the iterative image manipulation leads to an extension of the maximum and the minimum values, i.e. the flanks of the value spectrum, up to a factor of 1.7 – although we only find a few pixels in that extended range. If we applied a standard gray-scale color map to the whole wider range of values then we would get a strong reduction in contrast – for the sake of including some extreme pixel values. The remedy is to cut the spectrum off at the original maximum/minimum value.

The histogram analysis performed at the end in addition shows that the central spectrum has become a bit narrower – we can compensate this by some factors fac2 and fac3. Play around with them.

Results of the carving process

Results for map 57 are displayed in the following image:

The effect of the “carving” is clearly visible as an overlayed ghost like pattern apparition. Standardization of the final images makes no big difference – however, the discussed compensation for the widened spectrum of the eventual image is important.

The histogram analysis shows the flank widening and the narrowing of the central core of the pixel value spectrum:

The following image allow a comparison of the original gray image with the carved image for parameters [n_epochs=20, n_iter=5]:

Conclusion

Deep Dreams are based on special image analysis and manipulation techniques. We adapted some of the techniques discussed in literature for high resolution images and advanced CNNs in a special way and for a single deep layer map of a low resolution CNN. We called the required intermittent combination of down- and upscaling, enhancing OIP-patterns and re-adding details “carving“. This helped us to create a “dream” of a CNN – which was trained with a bunch low resolution MNIST images, only; the dream can be based on arbitrary quadratic input images (here: an image of a bunch of roses).

The technique presented in this article can be extended to include a pattern analysis on various sub-parts and length-scales of the original image. In the next blog post we shall have a look at some of the required programming steps.

My_OIP class – V0.65

I have modified the My_OIP code a bit since the last articles on CNNs. You find the new version below.

 
'''
Module to create OIPs as visualizations of filters of a simple CNN for the MNIST data set
@version: 0.65, 12.12.2020
@change: V0.5:  was based on version 0.4 which was originally created in Jupyter cells
@change: V0.6: General revision of class "my_OIP" and its methods
@change: V0.6: Changes to the documentation 
@change: V0.65: Added intermediate plotting of selected input fluctuation patterns 
@attention: General status: For experimental purposes only! 
@requires: A full CNN trained on MNIST data 
@requires: A Keras model of the CNN and weight data saved in a h5-file, e.g."cnn_MIST_best.h5". 
           This file must be placed in the main directory of the Jupyter notebooks.
@requires: A Jupyter environment - from where the class My_OIP is called and where plotting takes place 
@note: The description to the interface to the class via the __init__()-method may be incomplete
@note: The use of prefixes li_ and ay_ is not yet consistent. ay_ should indicate numpy arrays, li_ instead normal Python lists
@warning: This version has not been tested outside a Jupyter environment - plotting in GTK/Qt-environment may require substantial changes 
@status: Under major development with frequent changes 
@author: Dr. Ralph Mönchmeyer
@copyright: Simplified BSD License, 12.12.2020. Copyright (c) 2020, Dr. Ralph Moenchmeyer, Augsburg, Germany

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, this
   list of conditions and 
the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright notice,
   this list of conditions and the following disclaimer in the documentation
   and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
'''

# Modules to be imported - these libs must be imported in a Jupyter cell nevertheless 
# ~~~~~~~~~~~~~~~~~~~~~~~~
import time 
import sys 
import math
import os 
from os import path as path

import numpy as np
from numpy import save  # used to export intermediate data
from numpy import load

import scipy
from sklearn.preprocessing import StandardScaler
from itertools import product 

import tensorflow as tf
from tensorflow import keras as K
from tensorflow.python.keras import backend as B  # this is the only version compatible with TF 1 compat statements
#from tensorflow.keras import backend as B 
from tensorflow.keras import models
from tensorflow.keras import layers
from tensorflow.keras import regularizers
from tensorflow.keras import optimizers
from tensorflow.keras.optimizers import schedules
from tensorflow.keras.utils import to_categorical
from tensorflow.keras.datasets import mnist

from tensorflow.python.client import device_lib

import matplotlib as mpl
from matplotlib import pyplot as plt
from matplotlib.colors import ListedColormap
import matplotlib.patches as mpat 

from IPython.core.tests.simpleerr import sysexit

class My_OIP:
    '''
    @summary: This class allows for the creation and the display of OIP-patterns, 
              to which a selected map of a CNN-model and related filters react maximally
    
    @version: Version 0.6, 10.10.2020
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    @change: Revised methods 
    @requires: In the present version the class My_OIP requires: 
                * a CNN-model which works with standardized (!) input images, size 28x28,
                * a CNN-Modell which was trained on MNIST digit data,
                * exactly 4 length scales for random data fluctations are used to compose initial statistial image data 
                  (the length scales should roughly have a factor of 2 between them) 
                * Assumption : exactly 1 input image and not a batch of images is assumed in various methods 
    
    @note: Main Functions:     
    0) _init__()
    1) _load_cnn_model()             => load cnn-model
    2) _build_oip_model()            => build an oip-model to create OIP-images
    3) _setup_gradient_tape_and_iterate_function()        
                                    => Implements TF2 GradientTape to watch input data for eager gradient calculation
                                    => Creates a convenience function by the help of Keras to iterate and optimize the OIP-adjustments
    4) _oip_strat_0_optimization_loop():
                                     => Method implementing a simple strategy to create OIP-images, 
                                        based on superposition of random data on long range data (compared to 28 px) 
                                        The optimization uses "gradient ascent"
 to get an optimum output of the selected Conv map 
    6) _derive_OIP():                => Method used to start the creation of an OIP-image for a chosen map 
                                        - based on an input image with statistical random date 
    6) _derive_OIP_for_Prec_Img():   => Method used to start the creation of an OIP-image for a chosen map 
                                       - based on an input image with was derived from a PRECURSOR run, 
                                       which tests the reaction of the map to large scale fluctuations 
                                        
    7) _build_initial_img_data():    => Builds an input image based on random data for fluctuations on 4 length scales 
    8) _build_initial_img_from_prec():    
                                     => Reconstruct an input image based on saved random data for long range fluctuations 
    9) _prepare_precursor():         => Prepare a _precursor run by setting up TF2 GradientTape and the _iterate()-function 
    10) _precursor():                => Precursor run which systematically tests the reaction of a selected convolutional map 
                                        to long range fluctuations based on a (3x3)-grid upscaled to the real image size  
    11) _display_precursor_imgs():   => A method which plots up to 8 selected precursor images with fluctuations,
                                        which triggered a maximum map reaction 
    12) _transfrom_tensor_to_img():  => A method which allows to transform tensor data of a standardized (!) image to standard image data 
                                        with (gray)pixel valus in [0, 255]. Parameters allow for a contrast enhancement. 
    
    Usage hints 
    ~~~~~~~~~~~
    @note: If maps of a new convolutional layer are to be investigated then method _build_oip_model(layer_name) has to be rerun 
           with the layer's name as input parameter
    '''
    
    # Method to initialize an instantiation object 
    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    def __init__(self, cnn_model_file = 'cnn_MNIST_best.h5', 
                 layer_name = 'Conv2D_3', 
                 map_index = 0, 
                 img_dim = 28, 
                 b_build_oip_model = True  
                ): 
        '''
        @summary: Initialization of an object instance - read in a CNN model, build an OIP-Model 
        @note: Input: 
        ~~~~~~~~~~~~
        @param cnn_model_file:  Name of a file containing a fully trained CNN-model; 
                                the model can later be overwritten by self._load_cnn_model()
        @param layer_name:      We can define a layer name, which we are interested in,  already when starting;
                                the layer can later be overwritten by self._build_oip_model()
        @param map_index:       We can define a map, which we are interested in, already when starting;
                                A map-index is NOT required for building the OIP-model, but for the GradientTape-object 
        @param img_dim:         The dimension of the assumed quadratic images (2 for MNIST)

        @note: Major internal variables:
        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        @ivar _cnn_model:        A reference to the CNN model object
        @ivar _layer_name:       The name of convolutional layer 
                                 (can be overwritten by method _build_oip_model() ) 
        @ivar _map_index:        Index of the map in the chosen layer's output array 
                                 (can later be overwritten by other methods) 
        @ivar _r_cnn_inputs:     A reference to the input tensor of the CNN model 
                                 Could be a batch of images; but in this class only 1 image is assumed
        @ivar _layer_output:     Tensor with all maps of a certain layer    
        @
ivar _oip_submodel:     A new model connecting the input of the CNN-model with a certain map's (!) output
        @ivar _tape:             An instance of TF2's GradientTape-object 
                                    Watches input, output, loss of a model 
                                    and calculates gradients in TF2 eager mode 
        @ivar _r_oip_outputs:    Reference to the output of the new OIP-model = map-activation
        @ivar _r_oip_grads:      Reference to gradient tensors for the new OIP-model (output dependency on input image pixels)
        @ivar _r_oip_loss:       Reference to a loss defined on the OIP-output - i.e. on the activation values of the map's nodes;
                                 Normally chosen to be an average of the nodes' activations 
                                 The loss defines a hyperplane on the (28x28)-dim representation space of the input image pixel values  
        @ivar _val_oip_loss:     Loss value for a certain input image 
        @ivar _iterate           Reference toa Keras backend function which invokes the new OIP-model for a given image
                                 and calculates both loss and gradient values (in TF2 eager mode) 
                                 This is the function to be used in the optimization loop for OIPs
        
        @note: Internal Parameters controlling the optimization loop:
        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        @ivar _oip_strategy:        0, 1 - There are two strategies to evolve OIP patterns out of statistical data 
                                    - only the first strategy is supported in this version 
                                    Both strategies can be combined with a precursor calculation 
                                    0: Simple superposition of fluctuations at different length scales
                                    1: NOT YET SUPPORTED 
                                    Evolution over partially evolved images based on longer scale variations 
                                    enriched with fluctuations on shorter length scales 

        @ivar _ay_epochs:           A list of 4 optimization epochs to be used whilst 
                                    evolving the img data via strategy 1 and intermediate images 
        @ivar _n_epochs:            Number of optimization epochs to be used with strategy 0 
        @ivar _n_steps:             Defines at how many intermediate points we show images and report 
                                    during the optimization process 
        @ivar _epsilon:             Factor to control the amount of correction imposed by the gradient values of the oip-model 

        @note: Input image data of the OIP-model and references to it 
        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        @ivar _initial_precursor_img:  The initial image to start a precursor optimization with.
                                       Would normally be an image of only long range fluctuations. 
        @ivar _precursor_image:        The evolved image created and selected by the precursor loop 

        @ivar _initial_inp_img_data:  A tensor representing the data of the input image 
        @ivar _inp_img_data:          A tensor representig the 
        @ivar _img_dim:               We assume quadratic images to work with 
                                      with dimension _img_dim along each axis
                                      For the time being we only support MNIST images 
        
        @note: Internal parameters controlling the composition of random initial image data 
        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        @ivar _li_dim_steps:        A list of the intermediate dimensions for random data;
                                    these data are smoothly scaled to the image dimensions 
        @ivar _ay_facts:            A Numpy array of 4 factors to control the amount of 
                         
           contribution of the statistical variations 
                                    on the 4 length scales to the initial image
        
        @note: Internal variables to save data of a precursor run
        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
        @ivar _list_of_covs: list of long range fluctuation data for a (3x3)-grid covering the image area 
        @ivar _li_fluct_enrichments: [li_facts, li_dim_steps] data for enrichment with small fluctuations 
        
        
        @note: Internal variables for plotting
        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        @ivar _li_axa: A Python list of references to external (Jupyter-) axes-frames for plotting 
        
        '''    
        
        # Input data and variable initializations
        # ****************************************
        
        # the model 
        # ~~~~~~~~~~
        self._cnn_model_file = cnn_model_file
        self._cnn_model      = None 
        
        # the chosen layer of the CNN-model
        self._layer_name = layer_name
        # the index of the map in the layer array
        self._map_index  = map_index
        
        # References to objects and the OIP-submodel
        # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        self._r_cnn_inputs  = None # reference to input of the CNN_model, also used in the oip-model  
        self._layer_output  = None
        self._oip_submodel  = None
        
        # References to watched GradientTape objects 
        # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        self._tape          = None # TF2 GradientTape variable
        # some "references"
        self._r_oip_outputs = None # output of the oip-submodel to be watched 
        self._r_oip_grads   = None # gradients determined by GradientTape   
        self._r_oip_loss    = None # loss function
        # loss and gradient values (to be produced ba a backend function _iterate() )
        self._val_oip_grads = None
        self._val_oip_loss  = None
        
        # The Keras function to produce concrete outputs of the new OIP-model  
        self._iterate       = None
        
        # The strategy to produce an OIP pattern out of statistical input images 
        # ~~~~~~~~~~~~~~~~~~~~~~~~~--------~~~~~~
        # 0: Simple superposition of fluctuations at different length scales 
        # 1: Move over 4 interediate images - partially optimized 
        self._oip_strategy = 0
        
        # Parameters controlling the OIP-optimization process 
        # ~~~~~~~~~~~~~~~~~~~~~~~~~--------~~~~~~
        # number of epochs for optimization strategy 1
        self._ay_epochs    = np.array((20, 40, 80, 400), dtype=np.int32)
        len_epochs         = len(self._ay_epochs)
        
        # number of epochs for optimization strategy 0
        self._n_epochs     = self._ay_epochs[len_epochs-1]   
        self._n_steps      = 6   # divides the number of n_epochs into n_steps to produce intermediate outputs
        
        # size of corrections by gradients
        self._epsilon       = 0.01 # step-size for gradient correction
        
        # Input image-typess and references 
        # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        # precursor image
        self._initial_precursor_img = None
        self._precursor_img         = None # output of the _precursor-method 
        
        # The input image for the OIP-creation - a superposition of inial random fluctuations
        self._initial_inp_img_data  = None  # The initial data constructed 
        self._inp_img_data          = None  # The data used and varied for optimization 
        # image dimension
        self._img_dim               = img_dim   # = 28 => MNIST images for the time being 
        
        # Parameters controlling the setup of an initial image 
        # ~~~~~~~~~~~~~~~~~~~~~~~~~--------
~~~~~~~~~~~~~~~~~~~
        # The length scales for initial input fluctuations
        self._li_dim_steps = ( (3, 3), (7,7), (14,14), (28,28) ) # can later be overwritten 
        # Parameters for fluctuations  - used both in strategy 0 and strategy 1  
        self._ay_facts     = np.array( (0.5, 0.5, 0.5, 0.5), dtype=np.float32 )
        
        # Data of a _precursor()-run 
        # ~~~~~~~~~~~~~~~~~~~~~~~~~~~
        self._list_of_covs = None       # list of long range fluctuation data for a (3x3)-grid covering the image area 
        self._li_fluct_enrichments = None # = [li_facts, li_dim_steps] list with with 2 list of data enrichment for small fluctuations 
        # These data are required to reconstruct the input image to which a map reacted 
        
        # List of references to axis subplots
        # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        # These references may change from Jupyter cell to Jupyter cell and provided by the called methods
        self._li_axa = None # will be set by methods according to axes-frames in Jupyter cells 
        # axis frames for a single image in 2 versions (with contrast enhancement)
        self._ax1_1 = None
        self._ax1_2 = None 
        
        
        # Variables for Deep Dream Analysis 
        # ***********************************
        #     Here we need a stack of tiled sub-images of a large image, as we subdivide on different length scales 
        self._num_tiles = 1
        self._t_image_stack = None 
        
        # ********************************************************
        # Model setup - load the cnn-model and build the OIP-model
        # ************
        if path.isfile(self._cnn_model_file): 
            # We trigger the initial load of a model 
            self._load_cnn_model(file_of_cnn_model = self._cnn_model_file, b_print_cnn_model = True)
            # We trigger the build of a new sub-model based on the CNN model used for OIP search 
            self._build_oip_model(layer_name = self._layer_name, b_print_oip_model = True ) 
        else:
            print("<\nWarning: The standard file " + self._cnn_model_file + 
                  " for the cnn-model could not be found!\n " + 
                  " Please use method _load_cnn_model() to load a valid model")
            sys.exit()
        return
    
    
    #
    # Method to load a specific CNN model
    # **********************************
    def _load_cnn_model(self, file_of_cnn_model=None, b_print_cnn_model=True ):
        '''
        @summary: Method to load a CNN-model from a h5-file and create a reference to its input (image)
        @version: 0.2 of 28.09.2020
        @requires: filename must already have been saved in _cnn_model_file or been given as a parameter 
        @requires: file must be a h5-file 
        @change: minor changes - documentation 
        @note: A reference to the CNN's input is saved internally
        @warning: An input in form of an image - a MNIST-image - is implicitly assumed
        
        @note: Parameters
        -----------------
        @param file_of_cnn_model: Name of h5-file with the trained (!) CNN-model
        @param b_print_cnn_model: boolean - Print some information on the CNN-model 
        '''
        if file_of_cnn_model != None:
            self._cnn_model_file = file_of_cnn_model
        
        # Check existence of the file
        if not path.isfile(self._cnn_model_file): 
            print("\nWarning: The file " + file_of_cnn_model + 
                  " for the cnn-model could not be found!\n" + 
                  "Please change the parameter \"file_of_cnn_model\"" + 
                  " to load a valid model")
        
        # load the CNN model 
        self._cnn_model = models.load_model(self._cnn_model_file)
        
        # Inform about the model and its 
file
        # ~~~~~~~~~~~~~~~~~~~~~~~~~~~
        print("Used file to load a ´ model = ", self._cnn_model_file)
        # we print out the models structure
        if b_print_cnn_model:
            print("Structure of the loaded CNN-model:\n")
            self._cnn_model.summary()
        
        # handle/references to the models input => more precise the input image 
        # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        #    Note: As we only have one image instead of a batch 
        #    we pick only the first tensor element!
        #    The inputs will be needed for buildng the oip-model 
        self._r_cnn_inputs = self._cnn_model.inputs[0]  # !!! We have a btach with just ONE image 
        
        # print out the shape - it should be known from the original cnn-model
        if b_print_cnn_model:
            print("shape of cnn-model inputs = ", self._r_cnn_inputs.shape)
        
        return
    
    
    #
    # Method to construct a model to optimize input for OIP creation 
    # ***************************************
    def _build_oip_model(self, layer_name = 'Conv2D_3', b_print_oip_model=True ): 
        '''
        @summary: Method to build a new (sub-) model - the "OIP-model" - of the CNN-model by 
                  connectng the input of the CNN-model with one of its Conv-layers
        @version: 0.4 of 28.09.2020
        @change: Minor changes - documentation 
        @note: We need a Conv layer to build a working model for input image optimization
        We get the Conv layer by the layer's name 
        The new model connects the first input element of the CNN to the output maps of the named Conv layer CNN 
        We use Keras' models.Model() functionality 
        @note: The layer's name is crucial for all later investigations - if you want to change it this method has to be rerun 
        @requires: The original, trained CNN-model must be loaded and referenced by self._cnn_model 
        @warning: Only 1 input image and not a batch is implicitly assumed 
        
        @note: Parameters
        -----------------
        @param layer_name: Name of the convolutional layer of the CNN for whose maps we want to find OIP patterns
        @param b_print_oip_model: boolean - Print some information on the OIP-model 
        
        '''
        # free some RAM - hopefully 
        del self._oip_submodel
        
        # check for loaded CNN-model
        if self._cnn_model == None: 
            print("Error: CNN-model not yet defined.")
            sys.exit()
        
        # get layer name 
        self._layer_name = layer_name
        
        # We build a new model based on the model inputs and the output 
        self._layer_output = self._cnn_model.get_layer(self._layer_name).output
        # Note: We do not care at the moment about a complex composition of the input 
        # We trust in that we handle only one image - and not a batch
        
        # Create the sub-model via Keras' models.Model()
        model_name = "mod_oip__" + layer_name 
        self._oip_submodel = models.Model( [self._r_cnn_inputs], [self._layer_output], name = model_name)                                    

        # We print out the oip model structure
        if b_print_oip_model:
            print("Structure of the constructed OIP-sub-model:\n")
            self._oip_submodel.summary()
        return
    
    
    #
    # Method to set up GradientTape and an iteration function providing loss and gradient values  
    # *********************************************************************************************
    def _setup_gradient_tape_and_iterate_function(self, b_print = True, b_full_layer = False):
        '''
        @summary: Central method to watch input variables and resulting gradient changes 
        @version: 0.6 of 12.12.2020
n        @change: V0.6: We add the definition of a loss function for a whole layer 
        @note: For TF2 eager execution we need to watch input changes and trigger automatic gradient evaluation
        @note: The normalization of the gradient is strongly recommended; as we fix epsilon for correction steps 
               we thereby will get changes to the input data of an approximately constant order.
               This - together with standardization of the images (!) - will lead to convergence at the size of epsilon !
        @param b_full_layer: Boolean; define a loss function for a single map (map_index > -1) or the full layer (map_index = -1)
        '''   
        # Working with TF2 GradientTape
        self._tape = None

        # Watch out for input, output variables with respect to gradient changes
        # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        with tf.GradientTape() as self._tape: 
            # Input
            # ~~~~~
            self._tape.watch(self._r_cnn_inputs)
            # Output
            # ~~~~~~
            self._r_oip_outputs = self._oip_submodel(self._r_cnn_inputs)
            
            # only for testing
            # *******************
            #t_shape = self._r_oip_outputs[0, :, :, :].shape  # the 0 is important, else "Noen" => problems 
            #scaling_factor = tf.cast(tf.reduce_prod(t_shape), tf.float32) 
            #tf.print("scaling_factor: ", scaling_factor)
            #t_loss = tf.reduce_sum(tf.math.square(self._r_oip_outputs[0, :, :, :])) / scaling_factor
            
            # Define Loss 
            # ***********
            if (b_full_layer):
                # number of neurons over all layer maps 
                t_shape = self._r_oip_outputs[0, :, :, :].shape  # the 0 is important, else "Noen" => problems 
                scaling_factor = tf.cast(tf.reduce_prod(t_shape), tf.float32) # shape! eg. 3x3x128 for layer 3 
                tf.print("scaling_factor: ", scaling_factor)
                
                # The scaled total loss of the layer 
                self._r_oip_loss = tf.reduce_sum(tf.math.square(self._r_oip_outputs[0, :, :, :])) / scaling_factor
                
                # self._r_oip_loss = tf.reduce_mean(self._r_oip_outputs[0, :, :, self._map_index])
            else:
                self._r_oip_loss = tf.reduce_mean(self._r_oip_outputs[0, :, :, self._map_index])
            
            #self._loss = B.mean(oip_output[:, :, :, map_index])
            #self._loss = B.mean(oip_outputs[-1][:, :, map_index])
            #self._loss = tf.reduce_mean(oip_outputs[-1][ :, :, map_index])
            if b_print:
                print(self._r_oip_loss)
                print("shape of oip_loss = ", self._r_oip_loss.shape)
        
        # Gradient definition and normalization
        # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        self._r_oip_grads  = self._tape.gradient(self._r_oip_loss, self._r_cnn_inputs)
        print("shape of grads = ", self._r_oip_grads.shape)

        # normalization of the gradient - required for convergence 
        t_tiny = tf.constant(1.e-7, tf.float32)
        self._r_oip_grads /= (tf.sqrt(tf.reduce_mean(tf.square(self._r_oip_grads))) + t_tiny)
        #self._r_oip_grads /= (B.sqrt(B.mean(B.square(self._r_oip_grads))) + 1.e-7)
        #self._r_oip_grads = tf.image.per_image_standardization(self._r_oip_grads)

        # Define an abstract recallable Keras function as a convenience function for iterations 
        #     The _iterate()-function produces loss and gradient values for corrected img data 
        #     The first list of addresses points to the input data, the last to the output data 
        self._iterate = B.function( [self._r_cnn_inputs], [self._r_oip_loss, self._r_oip_grads] )
        
    
    #        
    # Method to optimize an emerging OIP out of statistical input 
image data 
    # (simple strategy - just optimize, no precursor, no intermediate pattern evolution, ..) 
    # ********************************
    def _oip_strat_0_optimization_loop(self, conv_criterion = 5.e-4, 
                                            b_stop_with_convergence = False,
                                            b_show_intermediate_images = True,
                                            b_print = True):
        '''
        @summary: Method to control the optimization loop for OIP reconstruction of an initial input image
                  with a statistical distribution of pixel values. 
        @version: 0.4 of 28.09.2020
        @changes: Minor changes - eliminated some unused variables
        @note:    The function also provides intermediate output in the form of printed data and images !.
        @requires: An input image tensor must already be available at _inp_img_data - created by _build_initial_img_data()
        @requires: Axis-objects for plotting, typically provided externally by the calling functions 
                  _derive_OIP() or _precursor()
        
        @note: Parameters:
        ----------------- 
        @param conv_criterion:  A small threshold number for (difference of loss-values / present loss value )
        @param b_stop_with_convergence: Booelan which decides whether we stop a loop if the conv-criterion is fulfilled
        @param b_show_intermediate_image: Boolean which decides whether we show up to 8 intermediate images 
        @param b_print: Boolean which decides whether we print intermediate loss values 
        
        @note: Intermediate information is provided at _n_steps intervals, 
               which are logarithmically distanced with respect to _n_epochs
               Reason: Most changes happen at the beginning 
        @note: This method produces some intermediate output during the optimization loop in form of images.
        It uses an external grid of plot frames and their axes-objects. The addresses of the 
        axis-objects must provided by an external list "li_axa[]" to self._li_axa[].  
        We need a seqence of _n_steps+2 axis-frames (or more) => length(_li_axa) >= _n_steps + 2). 
        
        @todo: Loop not optimized for TF 2 - but not so important here - a run takes less than a second 
        
        '''
        
        # Check that we already an input image tensor
        if ( (self._inp_img_data == None) or 
             (self._inp_img_data.shape[1] != self._img_dim) or 
             (self._inp_img_data.shape[2] != self._img_dim) ) :
            print("There is no initial input image or it does not fit dimension requirements (28,28)")
            sys.exit()

        # Print some information
        if b_print:
            print("*************\nStart of optimization loop\n*************")
            print("Strategy: Simple initial mixture of long and short range variations")
            print("Number of epochs = ", self._n_epochs)
            print("Epsilon =  ", self._epsilon)
            print("*************")

        # some initial value
        loss_old = 0.0
        
        # Preparation of intermediate reporting / img printing
        # --------------------------------------
        # Logarithmic spacing of steps (most things happen initially)
        n_el = math.floor(self._n_epochs / 2**(self._n_steps) ) 
        li_int = []
        if n_el != 0:
            for j in range(self._n_steps):
                li_int.append(n_el*2**j)
        else:  # linear spacing
            n_el = math.floor(self._n_epochs / (self._n_steps+1) ) 
            for j in range(self._n_steps+1):
                li_int.append(n_el*j)
        
        if b_print:
            print("li_int = ", li_int)
        
        # A counter for intermediate reporting  
        n_rep = 0
        
        # 
Convergence? - A list for steps meeting the convergence criterion
        # ~~~~~~~~~~~~
        li_conv = []
        
        
        # optimization loop 
        # *******************
        # counter for steps with zero loss and gradient values 
        n_zeros = 0
        
        for j in range(self._n_epochs):
            
            # Get output values of our Keras iteration function 
            # ~~~~~~~~~~~~~~~~~~~
            self._val_oip_loss, self._val_oip_grads = self._iterate([self._inp_img_data])
            
            # loss difference to last step - shuold steadily become smaller 
            loss_diff = self._val_oip_loss - loss_old 
            #if b_print:
            #    print("j: ", j, " :: loss_val = ", self._val_oip_loss, " :: loss_diff = ",  loss_diff)
            #    # print("loss_diff = ", loss_diff)
            loss_old = self._val_oip_loss
            
            if j > 10 and (loss_diff/(self._val_oip_loss + 1.-7)) < conv_criterion:
                li_conv.append(j)
                lenc = len(li_conv)
                # print("conv - step = ", j)
                # stop only after the criterion has been met in 4 successive steps
                if b_stop_with_convergence and lenc > 5 and li_conv[-4] == j-4:
                    return
            
            grads_val     = self._val_oip_grads
            #grads_val =   normalize_tensor(grads_val)
            
            # The gradients average value 
            avg_grads_val = (tf.reduce_mean(grads_val)).numpy()
            
            # Check if our map reacts at all
            # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            if self._val_oip_loss == 0.0 and avg_grads_val == 0.0 and b_print :
                print( "0-values, j= ", j, 
                       " loss = ", self._val_oip_loss, " avg_loss = ", avg_grads_val )
                n_zeros += 1
            
            if n_zeros > 10 and b_print: 
                print("More than 10 times zero values - Try a different initial random distribution of pixel values")
                return
            
            # gradient ascent => Correction of the input image data 
            # ~~~~~~~~~~~~~~~
            self._inp_img_data += self._val_oip_grads * self._epsilon
            
            # Standardize the corrected image - we won't get a convergence otherwise 
            # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            self._inp_img_data = tf.image.per_image_standardization(self._inp_img_data)
            
            # Some output at intermediate points 
            #     Note: We us logarithmic intervals because most changes 
            #           appear in the intial third of the loop's range  
            if (j == 0) or (j in li_int) or (j == self._n_epochs-1) :
                if b_print or (j == self._n_epochs-1):
                    # print some info 
                    print("\nstep " + str(j) + " finalized")
                    #print("Shape of grads = ", grads_val.shape)
                    print("present loss_val = ", self._val_oip_loss)
                    print("loss_diff = ", loss_diff)
                # show the intermediate image data 
                if b_show_intermediate_images: 
                    imgn = self._inp_img_data[0,:,:,0].numpy()
                    # print("Shape of intermediate img = ", imgn.shape)
                    self._li_axa[n_rep].imshow(imgn, cmap=plt.cm.get_cmap('viridis'))
                    # counter
                    n_rep += 1
        
        return
    
    
    #        
    # Standard UI-method to derive OIP from a given initial input image
    # ********************
    def _derive_OIP(self, map_index = 1, 
                          n_epochs = None, 
                          n_steps = 6, 
                          epsilon = 
0.01, 
                          conv_criterion = 5.e-4, 
                          li_axa = [], 
                          ax1_1 = None, ax1_2 = None, 
                          b_stop_with_convergence = False,
                          b_show_intermediate_images = True,
                          b_print = True,
                          centre_move = 0.33, fact = 1.0):
        '''
        @summary: Method to create an OIP-image for a given initial input image
                  This is the standard user interface for finding an OIP 
        @warning: This method should NOT be used for finding an initial precursor image 
                  Use _prepare_precursor() to define the map first and then _precursor() to evaluate initial images 
        @version: V0.5, 12.12.2020
        @change:  V0.4: Minor changes - added internal _li_axa for plotting, added documentation 
                  This method starts the process of producing an OIP of statistical input image data
        @change:  V0.5 : map_index can now have the value "-". Then a loss function for the whole layer (instead of a single map) is prepared
        @requires: A map index (or -1) should be provided to this method 
        @requires: An initial input image with statistical fluctuations of pixel values must have been created. 

        @warning:    This version only supports the most simple strategy - "strategy 0" 
        -------------    Optimize in one loop - starting from a superposition of fluctuations
                         No precursor, no intermediate evolutions

        @note: Parameters:
        -----------------
        @param map_index: We can and should chose a map here  (overwrites previous settings)
                          If map_index == -1 => We take a loss for the whole layer  
        @param n_epochs: Number of optimization steps  (overwrites previous settings) 
        @param n_steps:  Defines number of intervals (with length n_epochs/ n_steps) for reporting
                         standard value: 6 => 8 images - start image, end image + 6 intermediate 
                         This number also sets a requirement for providing (n_step + 2) external axis-frames 
                         to display intermediate images of the emerging OIP  
                         => see _oip_strat_0_optimization_loop()
        @param epsilon:  Size for corrections by gradient values
        @param conv_criterion: A small threshold number for convegenc (checks:  difference of loss-values / present loss value )
        @param b_stop_with_convergence: 
                         Booelan which decides whether we stop a loop if the conv-criterion is fulfilled
        @param _li_axa: A Python list of references to external (Jupyter-) axes-frames for plotting 
                 
        
        @note: Preparations for plotting: 
        We need n_step + 2 axis-frames which must be provided externally
        
        With Jupyter this can externally be done by statements like 

        # figure
        # -----------
        #sizing
        fig_size = plt.rcParams["figure.figsize"]
        fig_size[0] = 16
        fig_size[1] = 8
        fig_a = plt.figure()
        axa_1 = fig_a.add_subplot(241)
        axa_2 = fig_a.add_subplot(242)
        axa_3 = fig_a.add_subplot(243)
        axa_4 = fig_a.add_subplot(244)
        axa_5 = fig_a.add_subplot(245)
        axa_6 = fig_a.add_subplot(246)
        axa_7 = fig_a.add_subplot(247)
        axa_8 = fig_a.add_subplot(248)
        li_axa = [axa_1, axa_2, axa_3, axa_4, axa_5, axa_6, axa_7, axa_8]
        
        '''
        # Get input parameters
        self._map_index = map_index
        self._n_epochs  = n_epochs   
        self._n_steps   = n_steps
        self._epsilon   = epsilon
        
        # references to plot frames 
        self._li_axa = li_axa
        num_frames = len(li_
axa)
        if num_frames < n_steps+2:
            print("The number of available image frames (", num_frames, ") is smaller than required for intermediate output (", n_steps+2, ")")
            sys.exit()
        
        # 2 axes frames to display the final OIP image (with contrast enhancement) 
        self._ax1_1 = ax1_1
        self._ax1_2 = ax1_2
        
        # number of epochs for optimization strategy 0 
        # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        if n_epochs == None:
            len_epochs = len(self._ay_epochs)
            self._n_epochs   = self._ay_epochs[len_epochs-1]
        else: 
            self._n_epochs = n_epochs

        # Reset some variables  
        self._val_oip_grads = None
        self._val_oip_loss  = None 
        self._iterate       = None 

        # Setup the TF2 GradientTape watch and a Keras function for iterations 
        # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        # V0.5 => map_index == "-1"
        if self._map_index == -1: 
            self._setup_gradient_tape_and_iterate_function(b_print = b_print, b_full_layer=True)
            if b_print:
                print("GradientTape watch activated (map-based: map = ", self._map_index, ")")
        else: 
            self._setup_gradient_tape_and_iterate_function(b_print = b_print, b_full_layer=False)
            if b_print:
                print("GradientTape watch activated (layer-based: layer = ", self._layer_name, ")")
        
        '''
        # Gradient handling - so far we only deal with addresses 
        # ~~~~~~~~~~~~~~~~~~
        self._r_oip_grads  = self._tape.gradient(self._r_oip_loss, self._r_cnn_inputs)
        print("shape of grads = ", self._r_oip_grads.shape)
        
        # normalization of the gradient 
        self._r_oip_grads /= (B.sqrt(B.mean(B.square(self._r_oip_grads))) + 1.e-7)
        #self._r_oip_grads = tf.image.per_image_standardization(self._r_oip_grads)
        
        # define an abstract recallable Keras function 
        # producing loss and gradients for corrected img data 
        # the first list of addresses points to the input data, the last to the output data 
        self._iterate = B.function( [self._r_cnn_inputs], [self._r_oip_loss, self._r_oip_grads] )
        '''
            
        # get the initial image into a variable for optimization 
        self._inp_img_data = self._initial_inp_img_data

        # Start optimization loop for strategy 0 
        # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        if self._oip_strategy == 0: 
            self._oip_strat_0_optimization_loop( conv_criterion = conv_criterion, 
                                                b_stop_with_convergence = b_stop_with_convergence,  
                                                b_show_intermediate_images = b_show_intermediate_images,
                                                b_print = b_print
                                               )
        
        # Display the last OIP-image created at the end of the optimization loop
        # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        # standardized image 
        oip_img = self._inp_img_data[0,:,:,0].numpy()
        # transformed image 
        oip_img_t = self._transform_tensor_to_img(T_img=self._inp_img_data[0,:,:,0], 
                                                  centre_move = centre_move,
                                                  fact = fact)
        
        # display 
        ax1_1.imshow(oip_img, cmap=plt.cm.get_cmap('viridis'))
        ax1_2.imshow(oip_img_t, cmap=plt.cm.get_cmap('viridis'))
        
        return
    
    
    # 
    # Method to derive OIP from a given initial input image if map_index is already defined 
    # ********************
    def _derive_OIP_for_Prec_
Img(self, 
                          n_epochs = None, 
                          n_steps = 6, 
                          epsilon = 0.01, 
                          conv_criterion = 5.e-4, 
                          li_axa = [], 
                          ax1_1 = None, ax1_2 = None, 
                          b_stop_with_convergence = False,
                          b_show_intermediate_images = True,
                          b_print = True):
        '''
        @summary: Method to create an OIP-image for an already given map-index and a given initial input image
                  This is the core of OIP-detection, which starts the optimization loop  
        @warning: This method should NOT be used directly for finding an initial precursor image 
                  Use _prepare_precursor() to define the map first and then _precursor() to evaluate initial images 
        @version: V0.4, 28.09.2020
        @changes: Minor changes - added internal _li_axa for plotting, added documentation 
                  This method starts the process of producing an OIP of statistical input image data
        
        @note:    This method should only be called after _prepare_precursor(), _precursor(), _build_initial_img_prec() 
                  For a trial of different possible precursor images rerun _build_initial_img_prec() with a different index
        
        @requires: A map index should be provided to this method 
        @requires: An initial input image with statistical fluctuations of pixel values must have been created. 

        @warning:    This version only supports the most simple strategy - "strategy 0" 
        -------------    Optimize in one loop - starting from a superposition of fluctuations
                         no intermediate evolutions

        @note: Parameters:
        -----------------
        @param n_epochs: Number of optimization steps  (overwrites previous settings) 
        @param n_steps:  Defines number of intervals (with length n_epochs/ n_steps) for reporting
                         standard value: 6 => 8 images - start image, end image + 6 intermediate 
                         This number also sets a requirement for providing (n_step + 2) external axis-frames 
                         to display intermediate images of the emerging OIP  
                         => see _oip_strat_0_optimization_loop()
        @param epsilon:  Size for corrections by gradient values
        @param conv_criterion: A small threshold number for convegenc (checks:  difference of loss-values / present loss value )
        @param b_stop_with_convergence: 
                         Booelan which decides whether we stop a loop if the conv-criterion is fulfilled
        @param _li_axa: A Python list of references to external (Jupyter-) axes-frames for plotting 
                 
        
        @note: Preparations for plotting: 
        We need n_step + 2 axis-frames which must be provided externally
        
        With Jupyter this can externally be done by statements like 

        # figure
        # -----------
        #sizing
        fig_size = plt.rcParams["figure.figsize"]
        fig_size[0] = 16
        fig_size[1] = 8
        fig_a = plt.figure()
        axa_1 = fig_a.add_subplot(241)
        axa_2 = fig_a.add_subplot(242)
        axa_3 = fig_a.add_subplot(243)
        axa_4 = fig_a.add_subplot(244)
        axa_5 = fig_a.add_subplot(245)
        axa_6 = fig_a.add_subplot(246)
        axa_7 = fig_a.add_subplot(247)
        axa_8 = fig_a.add_subplot(248)
        li_axa = [axa_1, axa_2, axa_3, axa_4, axa_5, axa_6, axa_7, axa_8]
        
        '''
        # Get input parameters
        self._n_epochs  = n_epochs   
        self._n_steps   = n_steps
        self._epsilon   = epsilon
        
        # references to plot frames 
        self._li_axa = li_axa
        num_frames 
= len(li_axa)
        if num_frames < n_steps+2:
            print("The number of available image frames (", num_frames, ") is smaller than required for intermediate output (", n_steps+2, ")")
            sys.exit()
            
        # 2 axes frames to display the final OIP image (with contrast enhancement) 
        self._ax1_1 = ax1_1
        self._ax1_2 = ax1_2
        
        # number of epochs for optimization strategy 0 
        # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        if n_epochs == None:
            len_epochs = len(self._ay_epochs)
            self._n_epochs   = self._ay_epochs[len_epochs-1]
        else: 
            self._n_epochs = n_epochs
        
        # Note: No setup of GradientTape and _iterate(required) - this is done by _prepare_precursor 
            
        # get the initial image into a variable for optimization 
        self._inp_img_data = self._initial_inp_img_data

        # Start optimization loop for strategy 0 
        # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        if self._oip_strategy == 0: 
            self._oip_strat_0_optimization_loop( conv_criterion = conv_criterion, 
                                                b_stop_with_convergence = b_stop_with_convergence,  
                                                b_show_intermediate_images = b_show_intermediate_images,
                                                b_print = b_print
                                               )
        # Display the last OIP-image created at the end of the optimization loop
        # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        # standardized image 
        oip_img = self._inp_img_data[0,:,:,0].numpy()
        # transfored image 
        oip_img_t = self._transform_tensor_to_img(T_img=self._inp_img_data[0,:,:,0])
        
        # display 
        ax1_1.imshow(oip_img, cmap=plt.cm.get_cmap('viridis'))
        ax1_2.imshow(oip_img_t, cmap=plt.cm.get_cmap('viridis'))
        
        return
    
    
    # 
    # Method to build an initial image from a superposition of random data on different length scales 
    # ***********************************
    def _build_initial_img_data( self, 
                                 strategy = 0, 
                                 li_epochs    = [20, 50, 100, 400], 
                                 li_facts     = [0.5, 0.5, 0.5, 0.5],
                                 li_dim_steps = [ (3,3), (7,7), (14,14), (28,28) ], 
                                 b_smoothing = False,
                                 ax1_1 = None, ax1_2 = None):
        
        '''
        @summary: Standard method to build an initial image with random fluctuations on 4 length scales
        @version: V0.2 of 29.09.2020
        
        @note: This method should NOT be used for initial images based on a precursor images. 
               See _build_initial_img_prec(), instead.  
        
        @note: This method constructs an initial input image with a statistical distribution of pixel-values.
        We use 4 length scales to mix fluctuations with different "wave-length" by a simple approach: 
        We fill four square with a different numer of cells below the number of pixels 
        in each dimension of the real input image; e.g. (4x4), (7x7, (14x14), (28,28) <= (28,28). 
        We fill the cells with random numbers in [-1.0, 1.]. We smootly scale the resulting pattern 
        up to (28,28) (or what ever the input image dimensions are) by bicubic interpolations 
        and eventually add up all values. As a final step we standardize the pixel value distribution.          
        
        @warning: This version works with 4 length scales, only. 
        @warning: In the present version th eparameters "strategy " and li_epochs have no effect 
        
        @note: Parameters:
        -----------------
r
        @param strategy:  The strategy, how to build an image (overwrites previous settings) - presently only 0 is supported 
        @param li_epochs: A list of epoch numbers which will be used in strategy 1 - not yet supported 
        @param li_facts:  A list of factors which the control the relative strength of the 4 fluctuation patterns 
        @param li_dim_steps: A list of square dimensions for setting the length scale of the fluctuations  
        @param b_smoothing: Parameter which builds a control image   
        @param ax1_1: matplotlib axis-frame to display the built image 
        @param ax1_2: matplotlib axis-frame to display a second version of the built image 
        
        '''
        
        # Get input parameters 
        # ~~~~~~~~~~~~~~~~~~
        self._oip_strategy = strategy               # no effect in this version 
        self._ay_epochs    = np.array(li_epochs)    # no effect in this version 
        
        # factors by which to mix the random number fluctuations of the different length scales 
        self._ay_facts     = np.array(li_facts)
        # dimensions of the squares which simulate fluctuations 
        self._li_dim_steps = li_dim_steps
        
        # A Numpy array for the eventual superposition of random data 
        fluct_data = None
        
        
        # Strategy 0: Simple superposition of random patterns at 4 different wave-length
        # ~~~~~~~~~~
        if self._oip_strategy == 0:
            
            dim_1_1 = self._li_dim_steps[0][0] 
            dim_1_2 = self._li_dim_steps[0][1] 
            dim_2_1 = self._li_dim_steps[1][0] 
            dim_2_2 = self._li_dim_steps[1][1] 
            dim_3_1 = self._li_dim_steps[2][0] 
            dim_3_2 = self._li_dim_steps[2][1] 
            dim_4_1 = self._li_dim_steps[3][0] 
            dim_4_2 = self._li_dim_steps[3][1] 
            
            fact1 = self._ay_facts[0]
            fact2 = self._ay_facts[1]
            fact3 = self._ay_facts[2]
            fact4 = self._ay_facts[3]
            
            # print some parameter information
            # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            print("\nInitial image composition - strategy 0:\n Superposition of 4 different wavelength patterns")
            print("Parameters:\n", 
                 fact1, " => (" + str(dim_1_1) +", " + str(dim_1_2) + ") :: ", 
                 fact2, " => (" + str(dim_2_1) +", " + str(dim_2_2) + ") :: ", 
                 fact3, " => (" + str(dim_3_1) +", " + str(dim_3_2) + ") :: ", 
                 fact4, " => (" + str(dim_4_1) +", " + str(dim_4_2) + ")" 
                 )
            
            # fluctuations
            fluct1 =  2.0 * ( np.random.random((1, dim_1_1, dim_1_2, 1)) - 0.5 ) 
            fluct2 =  2.0 * ( np.random.random((1, dim_2_1, dim_2_2, 1)) - 0.5 ) 
            fluct3 =  2.0 * ( np.random.random((1, dim_3_1, dim_3_2, 1)) - 0.5 ) 
            fluct4 =  2.0 * ( np.random.random((1, dim_4_1, dim_4_2, 1)) - 0.5 ) 
            
            # Scaling with bicubic interpolation to the required image size
            fluct1_scale = tf.image.resize(fluct1, [28,28], method="bicubic", antialias=True)
            fluct2_scale = tf.image.resize(fluct2, [28,28], method="bicubic", antialias=True)
            fluct3_scale = tf.image.resize(fluct3, [28,28], method="bicubic", antialias=True)
            fluct4_scale = fluct4
            
            # superposition
            fluct_data = fact1*fluct1_scale + fact2*fluct2_scale + fact3*fluct3_scale + fact4*fluct4_scale
            
        
        # get the standardized plus smoothed and unsmoothed image 
        # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        #  TF2 provides a function performing standardization of image data function
        fluct_data_unsmoothed = tf.image.per_image_
standardization(fluct_data) 
        fluct_data_smoothed   = tf.image.per_image_standardization(
                                    tf.image.resize( fluct_data, [28,28], method="bicubic", antialias=True) )
        
        if b_smoothing: 
            self._initial_inp_img_data = fluct_data_smoothed
        else:
            self._initial_inp_img_data = fluct_data_unsmoothed
        
        # Display of both variants => there should be no difference 
        # ~~~~~~~~~~~~~~~~~~~~~~~~
        img_init_unsmoothed = fluct_data_unsmoothed[0,:,:,0].numpy()
        img_init_smoothed   = fluct_data_smoothed[0,:,:,0].numpy()
        ax1_1.imshow(img_init_unsmoothed, cmap=plt.cm.get_cmap('viridis'))
        ax1_2.imshow(img_init_smoothed, cmap=plt.cm.get_cmap('viridis'))
        
        print("Initial images plotted")
        
        return self._initial_inp_img_data    


    # 
    # Method to build an initial image from a superposition of a PRECURSOR image with random data on different length scales 
    # ***********************************
    def _build_initial_img_from_prec( self, 
                                 prec_index = -1,
                                 strategy = 0, 
                                 li_epochs    = (20, 50, 100, 400), 
                                 li_facts     = (1.0, 0.0, 0.0, 0.0, 0.0),
                                 li_dim_steps = ( (3,3), (7,7), (14,14), (28,28) ), 
                                 b_smoothing = False,
                                 b_print = True, 
                                 b_display = False, 
                                 ax1_1 = None, ax1_2 = None):
        
        '''
        @summary: Method to build an initial image based on a Precursor image => Input for _derive_OIP_for_Prec_Img()
        @version: V0.3, 03.10.2020
        @changes: V0.2: Minor- only documentation and comparison of index to length of _li_of_flucts[] 
        @changes: V0.3: Added Booleans to control the output and display of images 
        @changes: V0.4: Extended the reconstruction part / extended documentation
        
        @note: This method differs from _build_initial_img_data() as follows:
                * It uses a Precursor image as the fundamental image 
                * The data of the Precursor Image will be reconstructed from a (3x3) fluctuation pattern and enrichments
                * This method adds even further fluctuations if requested 
        @note: This method should be called manually from a Jupyter cell 
        @note: This method saves the reconstructed input image into self._initial_inp_img_data
        @note: This method should be followed by a call of self._derive_OIP_for_Prec_Img()
        
        @requires: Large scale fluctuation data saved in self._li_of_flucts[]
        @requires: Additional enrichment information in self._li_of_fluct_enrichments[]
        
        @param prec_index: This is an index ([0, 7[) of a large scale fluctuation pattern which was saved in self._li_of_flucts[] 
                           during the run of the method "_precursor()". The image tensor is reconstructed from the fluctuation pattern. 
        @warning: We support a maximum of 8 selected fluctuation patterns for which a map reacts 
        @warning: However, less precursor patterns may be found - so you should watch for the output of _precursor() before you run this method
        
        @param li_facts:  A list of factors which the control the relative strength of the precursor image vs. 
                          4 additional fluctuation patterns 
        @warning: Normally, it makes no sense to set li_facts[1] > 0 - because this will destroy the original large scale pattern 
        
        @note: For other input parameters see _build_initial_img_data()
        '''
        
        # Get input parameters 
        self._oip_strategy 
= strategy
        self._ay_facts     = np.array(li_facts)
        self._ay_epochs    = np.array(li_epochs)
        self._li_dim_steps = li_dim_steps
        
        fluct_data = None
        
        # Reconstruct an precursor image from a saved large scale fluctuation pattern (result of _precursor() ) 
        # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        length_li_cov = len(self._li_of_flucts)
        if prec_index > -1: 
            if length_li_cov < prec_index+1:
                print("index too large for number of detected patterns (", length_li_cov, ")")
                sys.exit()
            cov_p = self._li_of_flucts[prec_index][1][1]
            
            fluct0_scale_p = tf.image.resize(cov_p, [28,28], method="bicubic", antialias=True)
            
            # Add fluctuation enrichments - if saved [ len(self._li_fluct_enrichments) > 0 ]
            if len(self._li_fluct_enrichments) > 0:
                # Scaling enrichment flucts with bicubic interpolation to the required image size
                fluct1_p = self._li_fluct_enrichments[2][0] 
                fluct2_p = self._li_fluct_enrichments[2][1] 
                fluct3_p = self._li_fluct_enrichments[2][2] 
                fluct4_p = self._li_fluct_enrichments[2][3] 
                
                fact0_p = self._li_fluct_enrichments[0][0]
                fact1_p = self._li_fluct_enrichments[0][1]
                fact2_p = self._li_fluct_enrichments[0][2]
                fact3_p = self._li_fluct_enrichments[0][3]
                fact4_p = self._li_fluct_enrichments[0][4]
                
                fluct1_scale_p = tf.image.resize(fluct1_p, [28,28], method="bicubic", antialias=True)
                fluct2_scale_p = tf.image.resize(fluct2_p, [28,28], method="bicubic", antialias=True)
                fluct3_scale_p = tf.image.resize(fluct3_p, [28,28], method="bicubic", antialias=True)
                fluct4_scale_p = fluct4_p
                
                fluct_scale_p = fact0_p*fluct0_scale_p \
                 + fact1_p*fluct1_scale_p + fact2_p*fluct2_scale_p \
                 + fact3_p*fluct3_scale_p + fact4_p*fluct4_scale_p
                 
            else: 
                fluct_scale_p = fluct0_scale_p
            
            # get the img-data 
            fluct_data_p  = tf.image.per_image_standardization(fluct_scale_p)     
            fluct_data_p  = tf.where(fluct_data_p > 5.e-6, fluct_data_p, tf.zeros_like(fluct_data_p))
            self._initial_inp_img_data = fluct_data_p
            self._inp_img_data         = fluct_data_p
            
        
        # Strategy 0: Simple superposition of the precursor image AND additional patterns at 4 different wave-length
        # ~~~~~~~~~~
        if self._oip_strategy == 0:
            
            dim_1_1 = self._li_dim_steps[0][0] 
            dim_1_2 = self._li_dim_steps[0][1] 
            dim_2_1 = self._li_dim_steps[1][0] 
            dim_2_2 = self._li_dim_steps[1][1] 
            dim_3_1 = self._li_dim_steps[2][0] 
            dim_3_2 = self._li_dim_steps[2][1] 
            dim_4_1 = self._li_dim_steps[3][0] 
            dim_4_2 = self._li_dim_steps[3][1] 
            
            fact0 = self._ay_facts[0]
            fact1 = self._ay_facts[1]
            fact2 = self._ay_facts[2]
            fact3 = self._ay_facts[3]
            fact4 = self._ay_facts[4]
            
            # print some parameter information
            # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            if b_print:
                print("\nInitial image composition - strategy 0:\n Superposition of 4 different wavelength patterns")
                print("Parameters:\n", 
                     fact0, " => precursor image \n", 
                     fact1, " => (" + str(dim_1_1) +", " + str(dim_1_2) + ") :: ", 
                  
   fact2, " => (" + str(dim_2_1) +", " + str(dim_2_2) + ") :: ", 
                     fact3, " => (" + str(dim_3_1) +", " + str(dim_3_2) + ") :: ", 
                     fact4, " => (" + str(dim_4_1) +", " + str(dim_4_2) + ")" 
                     )
            
            # fluctuations
            fluct1 =  2.0 * ( np.random.random((1, dim_1_1, dim_1_2, 1)) - 0.5 ) 
            fluct2 =  2.0 * ( np.random.random((1, dim_2_1, dim_2_2, 1)) - 0.5 ) 
            fluct3 =  2.0 * ( np.random.random((1, dim_3_1, dim_3_2, 1)) - 0.5 ) 
            fluct4 =  2.0 * ( np.random.random((1, dim_4_1, dim_4_2, 1)) - 0.5 ) 
            
            # Scaling with bicubic interpolation to the required image size
            fluct1_scale = tf.image.resize(fluct1, [28,28], method="bicubic", antialias=True)
            fluct2_scale = tf.image.resize(fluct2, [28,28], method="bicubic", antialias=True)
            fluct3_scale = tf.image.resize(fluct3, [28,28], method="bicubic", antialias=True)
            fluct4_scale = fluct4
            
            # superposition with the already calculated image 
            fluct_data = fact0 * self._initial_inp_img_data  \
                         + fact1*fluct1_scale + fact2*fluct2_scale \
                         + fact3*fluct3_scale + fact4*fluct4_scale
            
        
        # get the standardized plus smoothed and unsmoothed image 
        # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        #    TF2 provides a function performing standardization of image data function
        fluct_data_unsmoothed = fluct_data 
        fluct_data_smoothed   = tf.image.per_image_standardization(
                                    tf.image.resize( fluct_data, [28,28], 
                                                     method="bicubic", antialias=True) )
        
        if b_smoothing: 
            self._initial_inp_img_data = fluct_data_smoothed
        else:
            self._initial_inp_img_data = fluct_data_unsmoothed
        
        # There should be no difference 
        img_init_unsmoothed = fluct_data_unsmoothed[0,:,:,0].numpy()
        img_init_smoothed   = fluct_data_smoothed[0,:,:,0].numpy()
        
        if b_display:
            ax1_1.imshow(img_init_unsmoothed, cmap=plt.cm.get_cmap('viridis'))
            ax1_2.imshow(img_init_smoothed, cmap=plt.cm.get_cmap('viridis'))
            print("Initial images plotted")
        
        return self._initial_inp_img_data    
    
    #
    # Method to prepare a Precursor run which checks a variety of large scale fluctuations for optimum actvation  
    # ***********************************
    def _prepare_precursor(self, map_index = 120, b_print = False):
        '''
        @summary: A method to prepare a Precursor run by setting up GradientTape and the _ierate() function for an optimization loop
        @version: 0.2, 30.09.2020
        @changes: Minor - adaption to _setup_gradient_tape_and_iterate_function() instead of the obsolete _setup_gradient_tape()
        @requires: A loaded CNN-Model and an already built OIP-model 
        @requires: A valid map-index as input parameter 
        
        @note: This method sets up the GradientTape and _iterate only once 
               - it will not be done again during the investigation of thousands of different input images (with large scale fluctuations) 
               which we investigate during the _precursor()-run. 
        
        @param map_index: Index selecting a map for the CNN layer defined previously by _build_oip_model()
        @param b_print: Boolean - decides about intemediate output 
        
        '''
        
        # Get input parameters 
        # ~~~~~~~~~~~~~~~~~~~~~~
        self._map_index = map_index
        # Rest some variables  
        self._val_oip_grads = None
        self._val_oip_loss  = None 
 
       self._iterate       = None 
        
        # Setup the TF2 GradientTape watch
        # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        self._setup_gradient_tape_and_iterate_function(b_print = b_print)
        if b_print:
            print("GradientTape watch activated and _iterate()-function defined")
    
    
    # Method to prepare a Precursor run which checks a variety of large scale fluctuations for optimum actvation
    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    def _precursor(self, 
                   li_pre_val=[0.2, 0.5, 0.8], 
                   num_epochs=10, 
                   loss_limit = 0.5, 
                   b_print = True, 
                   b_with_enriched_fluct = False, 
                   li_facts     = (1.0, 0.0, 0.0, 0.0, 0.1),
                   li_dim_steps = ( (3,3), (7,7), (14,14), (28,28) ), 
                   b_check = True, 
                   fig_a = None, 
                   li_axa = None,
                   # the next parameters are new in V0.6
                   b_show_test_input_images = False,
                   fig_test_img = None,
                   ax_b = None, 
                   interval_test_img = 1):

        '''
        @summary: Method to investigate thousands of input images with large scale fluctuations for the reaction of a specified map (filters)
                  and a given number of epochs in pattern creation
        @version: 0.6, 12.12.2020
        @changes: V0.5: Minor - documentation, skipped some superfluous statements 
        @changes: V0.6: Added intermediate plot output of constructed input images (a fig_a-reference must have been provided!) 
        @note: We select the 8 most dominant images - or less, if there are fewer input images which trigger the map 
        @requires: Previous run of _prepare_cursor() with a definition of the map-index
        @note: We vary only 3 given pixel values on (3x3) grids (19683 possibilities)
        @note: The optimization loop is completely done within this method - due to performance reasons
        
        @param li_pre_val: A list of three scaled pixel values between ]0, 1[ which shall be combined in (3x3)-fluctuation patterns
        
        @param num_epochs: The number of epochs used in the optimization loop for pattern creation
        @note: It is worthwile to experiment with the number of epochs - the (3x3)-pattern selection may change !!!
        
        @param loss_limit: Threshold of loss for which we register a fluctuation image as relevant 
        
        @param b_print: Boolean - controls printout for intermediate results - useful to see map response 
        @param b_check: Check the response of map for the first saved pattern - check the image reconstruction at the same time
        
        @note: Parameters to enrich the (3x3)-large scale fluctuation with a small scale pattern 
        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        @param b_with_enriched_fluct: Boolean - controls whether we enrich the long-range pattern with other additional patterns 
        @param li_facts:  A list of factors which the control the relative strength of the precursor image vs. 
                          4 additional fluctuation patterns 
        @param li_dim_steps: A list of square dimensions for setting the length scale of the fluctuations  
        
        @note: Parameters for plotting  
        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        @param fig_a: reference to the 8 axes-frame 
        @param li_axa: A Python list of references to external (Jupyter-) axes-frames for plotting 
        
        Plots for intermediate test images 
        Note that for inernediate image plotting the cell must have a command line with metacommand 
        %matplotlib notebook
        plt.ion()
        @param b_show_test_input_images = Boolean; show intermediate 
fluctuation patterns 
        @param fig_test_img: reference to external image plot frame. 
        @param ax_b: Rference to axis-frame for the intermediate display of test-images  
        
        @note: Preparations for plotting: 
        ---------------------------------
        fig_b : a simple figure for an image
        fig_a: We need 8 axis-frames which must be provided externally
        With Jupyter this can externally be done by statements like 

        # figure
        # -----------
        fig_b_pre, ax_b_pre = plt.subplots(1, figsize=(3,3))
        fig_b_pre.canvas.draw()
        
        fig_a_pre = plt.figure(2, figsize=(10,5))
        axa_1 = fig_a_pre.add_subplot(241)
        axa_2 = fig_a_pre.add_subplot(242)
        axa_3 = fig_a_pre.add_subplot(243)
        axa_4 = fig_a_pre.add_subplot(244)
        axa_5 = fig_a_pre.add_subplot(245)
        axa_6 = fig_a_pre.add_subplot(246)
        axa_7 = fig_a_pre.add_subplot(247)
        axa_8 = fig_a_pre.add_subplot(248)
        li_axa = [axa_1, axa_2, axa_3, axa_4, aab_5, axa_6, axa_7, axa_8]
        fig_a_pre.canvas.draw()
        
        
        '''
        # Internal parameter for number of selected input patterns 
        num_selected = 8 
        # check if length of li_axa is sufficient
        if li_axa == None or len(li_axa) < num_selected:
            print("Error: The length of the provided list with axes-frames for plotting must be at least ", num_selected )
            sys.exit()
        
        # get required exernal params 
        # ~~~~~~~~~~~~~~~~~~~~~~~~
        self._n_steps   = 2           # only a dummy 
        self._epsilon   = 0.01        # only a dummy  
        
        # number of epochs for optimization strategy 0 
        if num_epochs == None:
            len_epochs = len(self._ay_epochs)
            self._n_epochs   = self._ay_epochs[len_epochs-1]
        else: 
            self._n_epochs = num_epochs  
        
        # Create a fixed random flutution pattern which we later can overlay to the long range fluctuation patterns 
        # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        self._li_dim_steps = li_dim_steps # dimensions for fluctuations 
        self._ay_facts     = np.array(li_facts)
        
        dim_1_1 = self._li_dim_steps[0][0] 
        dim_1_2 = self._li_dim_steps[0][1] 
        dim_2_1 = self._li_dim_steps[1][0] 
        dim_2_2 = self._li_dim_steps[1][1] 
        dim_3_1 = self._li_dim_steps[2][0] 
        dim_3_2 = self._li_dim_steps[2][1] 
        dim_4_1 = self._li_dim_steps[3][0] 
        dim_4_2 = self._li_dim_steps[3][1] 
        
        fact0 = self._ay_facts[0]
        fact1 = self._ay_facts[1]
        fact2 = self._ay_facts[2]
        fact3 = self._ay_facts[3]
        fact4 = self._ay_facts[4]
        
        # Create fluctuation patterns for enrichment 
        fluct1 =  2.0 * ( np.random.random((1, dim_1_1, dim_1_2, 1)) - 0.5 ) 
        fluct2 =  2.0 * ( np.random.random((1, dim_2_1, dim_2_2, 1)) - 0.5 ) 
        fluct3 =  2.0 * ( np.random.random((1, dim_3_1, dim_3_2, 1)) - 0.5 ) 
        fluct4 =  2.0 * ( np.random.random((1, dim_4_1, dim_4_2, 1)) - 0.5 ) 
        
        li_fluct = [fluct1, fluct2, fluct3, fluct4]
        
        # Scaling with bicubic interpolation to the required image size
        fluct1_scale = tf.image.resize(fluct1, [28,28], method="bicubic", antialias=True)
        fluct2_scale = tf.image.resize(fluct2, [28,28], method="bicubic", antialias=True)
        fluct3_scale = tf.image.resize(fluct3, [28,28], method="bicubic", antialias=True)
        fluct4_scale = fluct4

        
        # Create cartesian product of combinatoric possibilities for a (3x3)-grid of long range fluctuations
        # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        cp = list(product(li_pre_val, repeat=9))
      
  num = len(cp) 
        print ("We test ", num, " possibilities for a (3x3) fluctuations ")
        
        # Prepare lists to save parameter data for the fluctuation pattern 
        # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        # Intermediate list to save relevant long scale fluctuations 
        d_cov = {}
        # list to save parameters for an enrichments of the large scale pattern with small fluctuations
        if b_with_enriched_fluct:
            self._li_fluct_enrichments = [li_facts, li_dim_steps, li_fluct]
        else:
            self._li_fluct_enrichments = []
            
        # Loop to check for relevant fluctuations => Loop over all combinations
        # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        for i_cp in range(num): 
            
            # create the value distribution 
            cov = np.array(cp[i_cp])
            cov = cov.reshape(1,3,3,1) 

            # create basic image to investigate
            # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            fluct_scale = tf.image.resize(cov, [28,28], method="bicubic", antialias=True) 
            
            # enrich with additional smallscale fluctuations 
            if b_with_enriched_fluct: 
                # superposition with the already calculated image 
                fluct_scale = fact0 * fluct_scale \
                             + fact1*fluct1_scale + fact2*fluct2_scale \
                             + fact3*fluct3_scale + fact4*fluct4_scale
            
            #standardization of the image data 
            fluct_data  = tf.image.per_image_standardization(fluct_scale)     
            # eliminatng very small values - prove to be helpful in many cases 
            fluct_data = tf.where(fluct_data > 5.e-6, fluct_data, tf.zeros_like(fluct_data))
            
            # save image data 
            self._initial_inp_img_data = fluct_data
            self._inp_img_data         = fluct_data
            
            # Display the present input image 
            if b_show_test_input_images and i_cp%interval_test_img == 0: 
                img_b = fluct_data[0, :,:, 0].numpy()
                ax_b.imshow(img_b, cmap=plt.cm.get_cmap('viridis'))
                fig_test_img.canvas.draw()
            
            # optimization loop 
            # ~~~~~~~~~~~~~~~~~
            for j in range(self._n_epochs):
                # Get output values of our Keras iteration function 
                self._val_oip_loss, self._val_oip_grads = self._iterate([self._inp_img_data])
                # gradient ascent => Correction of the input image data 
                self._inp_img_data += self._val_oip_grads * self._epsilon
                # Standardize the corrected image - we won't get a convergence otherwise 
                self._inp_img_data = tf.image.per_image_standardization(self._inp_img_data)
                
            # Check if we have a loss value > 0.x and save the (3x3)-data
            # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            if self._val_oip_loss > loss_limit:
                d_cov[i_cp] = [self._val_oip_loss, cov] 
                if b_print:
                    tf.print("i = ", i_cp, " loss = ", self._val_oip_loss)

        # We restrict the number to those 8 with highest loss 
        # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        if len(d_cov) > 0: 
            self._li_of_flucts = sorted(d_cov.items() , reverse = True,  key=lambda x: x[1][0])
            # print("num of relevant covs = ", len(self._li_of_flucts), len(d_cov))   
            print("\nnum of relevant covs = ", len(self._li_of_flucts))   
            self._li_of_flucts = self._li_of_flucts[:num_selected].copy()
            #save( 'li_of_flucts.npy', np.array(self._li_of_flucts, dtype=np.float32) )
   
         save( 'li_of_flucts.npy', self._li_of_flucts)
            # save the enrichment-setting 
            save('li_of_cov_enrichments.npy',self._li_fluct_enrichments)
            
            
            # Check of map really reacted - and check the reconstruction mechanism 
            # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            if b_check: 
                print("check of map reaction to first selected image")
                cov_t = self._li_of_flucts[0][1][1]
                #print("cov_t-shape = ", cov_t.shape)
                #cov_del = cov_t - cov
                #print("cov_del =\n", cov_del)
                
                fluct0_scale_t = tf.image.resize(cov_t, [28,28], method="bicubic", antialias=True) 
                        
                # Add fluctuation enrichments - if saved
                if b_with_enriched_fluct:
                    
                    # Scaling enrichment flucts with bicubic interpolation to the required image size
                    fluct1_t = self._li_fluct_enrichments[2][0] 
                    fluct2_t = self._li_fluct_enrichments[2][1] 
                    fluct3_t = self._li_fluct_enrichments[2][2] 
                    fluct4_t = self._li_fluct_enrichments[2][3] 
                    
                    fact0_t = self._li_fluct_enrichments[0][0]
                    fact1_t = self._li_fluct_enrichments[0][1]
                    fact2_t = self._li_fluct_enrichments[0][2]
                    fact3_t = self._li_fluct_enrichments[0][3]
                    fact4_t = self._li_fluct_enrichments[0][4]
                    
                    fluct1_scale_t = tf.image.resize(fluct1_t, [28,28], method="bicubic", antialias=True)
                    fluct2_scale_t = tf.image.resize(fluct2_t, [28,28], method="bicubic", antialias=True)
                    fluct3_scale_t = tf.image.resize(fluct3_t, [28,28], method="bicubic", antialias=True)
                    fluct4_scale_t = fluct4_t
                    
                    fluct_scale_t = fact0_t*fluct0_scale_t \
                                 + fact1_t*fluct1_scale_t + fact2_t*fluct2_scale_t \
                                 + fact3_t*fluct3_scale_t + fact4_t*fluct4_scale_t
                    
                else: 
                    fluct_scale_t = fluct0_scale_t
                
                #standardization
                fluct_data_t  = tf.image.per_image_standardization(fluct_scale_t)     
                fluct_data_t  = tf.where(fluct_data_t > 5.e-6, fluct_data_t, tf.zeros_like(fluct_data_t))
                self._initial_inp_img_data = fluct_data_t
                self._inp_img_data         = fluct_data_t
                self._precursor_img        = fluct_data_t
                
                # optimization loop 
                for j in range(self._n_epochs):
                    self._val_oip_loss, self._val_oip_grads = self._iterate([self._inp_img_data])
                    self._inp_img_data += self._val_oip_grads * self._epsilon
                    self._inp_img_data = tf.image.per_image_standardization(self._inp_img_data)
                print("loss for 1st selected img = ", self._val_oip_loss )

            # show the imgs 
            # ~~~~~~~~~~~~~~
            self._display_precursor_imgs(li_axa = li_axa, fig_a = fig_a)
            
        # list contains no patterns 
        else:
            print("No image found !")
    
    
    # Method to display initial fluctuation images identified as objects for OIP creation 
    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    def _display_precursor_imgs(self, li_axa = None, fig_a = None):
        '''
        @summary: Method to display up to 8 selected precursor images 
        @version: 0.2, 02.10.2020
        @change: Only some documentation 
        @note: We first 
reconstruct the image from saved data of the large scale fluctuations 
        @note: We then display the images in externally delivered axes-frames of matplotlib
        @requires: A filled set of valid fluctuation patterns in self._li_of_flucts[][][] - determined by a run of _precursor()
        @requires: Information on fluctuation enrichments in self._li_fluct_enrichments[][] - determined by a _precursor() run
        @requires: A set of axes-frames for plotting - preferably defined in a Jupyter cell calling thi smethod 
        @note: Parameters for plotting  
        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        @param _li_axa: A Python list of references to external (Jupyter-) axes-frames for plotting 
         '''
        # length of _li_of_flucts[] vs. length of li_axa
        len_cov = len(self._li_of_flucts)
        if li_axa == None or len(li_axa) < len_cov:
            print("Error: The length of the provided list with axes-frames for plotting must be at least ", len_cov )
            sys.exit()
        
        # Loop to reconstruct and display the found precursor images 
        # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        for j in range(len(self._li_of_flucts)):
            print(j, "loss = ", self._li_of_flucts[j][1][0])
            
            # reconstruct the image from the data of the precursor run 
            # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            cov = self._li_of_flucts[j][1][1]
            fluct0_scale_t  = tf.image.resize(cov, [28,28], method="bicubic", antialias=True)        
            
            # add enrichments if defined 
            if len(self._li_fluct_enrichments) > 0: 
                # Scaling enrichment flucts with bicubic interpolation to the required image size
                fluct1_t = self._li_fluct_enrichments[2][0] 
                fluct2_t = self._li_fluct_enrichments[2][1] 
                fluct3_t = self._li_fluct_enrichments[2][2] 
                fluct4_t = self._li_fluct_enrichments[2][3] 
                
                fact0_t = self._li_fluct_enrichments[0][0]
                fact1_t = self._li_fluct_enrichments[0][1]
                fact2_t = self._li_fluct_enrichments[0][2]
                fact3_t = self._li_fluct_enrichments[0][3]
                fact4_t = self._li_fluct_enrichments[0][4]
                
                fluct1_scale_t = tf.image.resize(fluct1_t, [28,28], method="bicubic", antialias=True)
                fluct2_scale_t = tf.image.resize(fluct2_t, [28,28], method="bicubic", antialias=True)
                fluct3_scale_t = tf.image.resize(fluct3_t, [28,28], method="bicubic", antialias=True)
                fluct4_scale_t = fluct4_t
                
                fluct_scale_t = fact0_t*fluct0_scale_t \
                             + fact1_t*fluct1_scale_t + fact2_t*fluct2_scale_t \
                             + fact3_t*fluct3_scale_t + fact4_t*fluct4_scale_t
            else:
                fluct_scale_t = fluct0_scale_t                 
            
            fluct_datx  = tf.image.per_image_standardization(fluct_scale_t)     
            fluct_dat  = tf.where(fluct_datx > 5.e-6, fluct_datx, tf.zeros_like(fluct_datx))
            img = fluct_dat[0, :,:, 0].numpy()
            li_axa[j].imshow(img, cmap=plt.cm.get_cmap('viridis'))
        
        fig_a.canvas.draw()
        return
    
    
    # Method to transform an img tensor into a standard image - used for contrast enhancemant  
    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    def _transform_tensor_to_img(self, T_img = None, centre_move = 0.33, fact = 1.0):
        '''
        @summary: Method to transform (standardized) tensor img data into standard img array data 
                  and to apply contrast enhancement(> V0.2)
        @change: V0.2: Additional 
parameers to control contrast enhancement. 
                 We scale the the value distribution, move its center and apply some clipping  
        @note: Clipping is used to remove pixel values outside [0, 255] 
        @version: 0.2, 12.10.2020
        
        @requires: A defined TF2/Keras backend 
        
        @param T_img: The TF2 or keras tensor for image data 
        @param centre_move: A decimal to move the centre of the pixel value distribution  
        
        
        '''
        ay_x = T_img.numpy()  # floating point array  
        
        maxi_o    = np.max(T_img)
        avg_o     = np.mean(T_img)
        mini_o    = np.min(T_img)
        std_dev_o = np.std(T_img)
        print("\nInfos on pixel value distribution during contrast enhancement: ") 
        print("\max_orig = ", maxi_o, " :: avg_orig = ", avg_o, " :: min_orig: ", mini_o) 
        print("std_dev_orig = ", std_dev_o)

        # the following operation should have no effect on standardized images
        ay_x -= ay_x.mean()
        ay_x /= (ay_x.std() + B.epsilon())
        
        maxi    = np.max(ay_x)
        avg     = np.mean(ay_x)
        mini    = np.min(ay_x)
        std_dev = np.std(ay_x)
        print("max_ay = ", maxi, " :: avg_ay = ", avg, " :: min_ay: ", mini) 
        print("std_dev_ay = ", std_dev)
        
        div = fact * 0.5 * ( abs(maxi_o) + abs(mini_o) )
        print("div = ", div)
        ay_x /= div          # scaling  
        ay_x += centre_move  # moving the data center - 0.5 would move of the distribution into [0,1]
        #                    # smaller values emphasize light/dark contratss 
        maxi = np.max(ay_x)
        avg  = np.mean(ay_x)
        mini = np.min(ay_x)
        std_dev = np.std(ay_x)
        print("max_fin = ", maxi, " :: avg_fin = ", avg, " :: min_fin: ", mini) 
        print("std_dev_fin = ", std_dev)
        
        ay_x = np.clip(ay_x, 0, 1)
        
        ay_x *= 255
        ay_x_img = np.clip(ay_x, 0, 255).astype('uint8')
        
        maxi = np.max(ay_x_img)
        avg  = np.mean(ay_x_img)
        mini = np.min(ay_x_img)
        std_dev = np.std(ay_x_img)
        print("max_img = ", maxi, " :: avg_img = ", avg, " :: min_img: ", mini) 
        print("std_dev_img = ", std_dev, "\n")
        
        
        return ay_x_img


 

Links

A simple CNN for the MNIST dataset – XI – Python code for filter visualization and OIP detection